Realization of tristability in a multiplicatively coupled dual-loop genetic network

Multistability is a crucial recurring theme in cell signaling. Multistability is attributed to the presence of positive feedback loops, but the general condition and essential mechanism for realizing multistability remain unclear. Here, we build a generic circuit model comprising two transcription factors and a microRNA, representing a kind of core architecture in gene regulatory networks. The circuit can be decomposed into two positive feedback loops (PFLs) or one PFL and one negative feedback loop (NFL), which are multiplicatively coupled. Bifurcation analyses of the model reveal that the circuit can achieve tristability through four kinds of bifurcation scenarios when parameter values are varied in a wide range. We formulate the general requirement for tristability in terms of logarithmic gain of the circuit. The parameter ranges for tristability and possible transition routes among steady states are determined by the combination of gain features of individual feedback loops. Coupling two PFLs with bistability or one NFL with a bistable PFL is most likely to generate tristability, but the underlying mechanisms are largely different. We also interpret published results and make testable predictions. This work sheds new light on interlinking feedback loops to realize tristability. The proposed theoretical framework can be of wide applicability.

Here, we present a detailed derivation of the dimensionless differential equations. The model comprises two transcription factors (TFs) and one microRNA (miRNA), constituting two interlinked feedback loops.

Transcriptional dynamics
The transcription process is characterized by a simple two-state model-the promoter of a gene can be bound by n TFs or unbound. We assume that the transitions between unbound (U ) and bound (B) states of the promoter are both a one-step process: According to the law of mass action, the rates of forward and backward reactions are where [TF] is the concentration of TFs, and k B and k U denote the rate constants associated with binding and dissociating, respectively. The evolution of the probability for the promoter bound by TFs, P B , is governed by the following chemical master equation: where P U is the probability of the promoter unbound by TFs, i.e., P B + P U = 1. The stationary distributions for the unbound and bound states of the promoter obey Two different states of the promoter allow for different rates of transcription conducted by RNA polymerase. k 0 and k 1 denote the transcription rate when the promoter is unbound or bound, respectively. Given the binding and dissociating rates of TFs are much faster than the elongation rate, the average transcription ratek obeys with K = n √ kU kB describing the binding affinity of TF to the promoter. k 1 > k 0 when TF activates transcription, whereas k 1 < k 0 when TF inhibits transcription. If transcription is initialized only after TF binds to the promoter, k 0 should be zero, leading tok If transcription is completely prevented after TF binds to the promoter, k 1 should be zero, resulting in Of note, the Hill coefficient is not necessarily a positive integer. The above formulas were widely used to describe the transcription dynamics [1, 2].

Dynamics of microRNA-mediated translational repression
miRNA is a class of short non-coding RNAs with about 22 nucleotides. It can function as a probe to detect the target mRNA by fully or partially binding to its specific sequences in 3' UTR. Meanwhile, as part of miRNAinduced silencing complex (miRISC), miRNA can lead the miRISC to deal with its target mRNA by degrading it or preventing it from being translated. No matter how miRISC processes mRNA, the subsequent translation of mRNA will be prevented, which is equivalent to a decrease in the effective amount of mRNA. Thus, we consider the process of regulation by miRNA as the following reaction: miRNA actually functions as an enzyme during its lifetime. Instead of describing this reaction by the Michaelis-Menten kinetics, we take it as a second-order reaction for simplicity. Thus, the reaction rate is k Mm [M ][m].

Equations for the model
Combining the above two processes, we write down the equations governing the dynamics of the circuit: Given that the lifetime of mRNA is much shorter than those of proteins [3] and miRNA [4], the amount of mRNA can rapidly reach an equilibrium. Thus, the right-hand sides of Eqs. (8) and (10) are set to zero, leading to Consequently, we obtain the following equations where k X = km 1 βX γm 1 and k Y = km 2 βY γm 2 denote the effective production rates of X and Y, respectively, and k Z = kZm 1 γm 1 represents the effective strength of inhibition of Y mRNA translation. The dimensionless time, variables and parameters are defined as follows: γZ . The equations in dimensionless forms are: Similarly, the dimensionless equations for the TF-TF and TF-M-TF loops are separately and

LINEAR STABILITY ANALYSIS AND BIFURCATION ANALYSIS
In this section, we analyze the linear stability of the systems of individual and interlinked feedback loops at steady states. We also give the necessary conditions for saddle-node (SN) bifurcation and Hopf bifurcation.

Linear stability analysis
. η is chosen as the control parameter of the system. According to Eqs. (18)-(20), the steady-state equations for the coupled system are: The steady-state values of x, y and z for a given η are denoted by x s , y s and z s . Note that these equations can be reduced to the following single equation: where the symbol • stands for the composition of functions. Obviously, the solution to this equation is the steady-state value of x, i.e., x s . Moreover, given any x s ∈ [0, +∞), the corresponding η, y s and z s can be calculated from Eqs.
(29), (27) and (28), respectively. Thus, we can use x s ∈ [0, +∞) to represent any steady state for different η. In the following, the subscript of x s is omitted and x stands for the steady-state value of X.
To determine the linear stability of the system at steady states, we investigate the roots of the characteristic equation of the Jacobian matrix associated with the ordinary differential equations evaluated at the steady states. For the coupled system, the characteristic equation for Eqs. (18)-(20) is where Note that η is replaced with the function of x according to Eq. (29), i.e., Thus, all the above coefficients only depend on x. If all roots of Eq. (30) have a negative real part, the steady state is stable.
If any root has a positive real part, the steady state is unstable. Similarly, the characteristic equation for the TF-TF loop is with B 1 = 1 θ y + 1, The characteristic equation for the TF-M-TF loop is with

Bifurcation analysis
For SN bifurcation to occur, the characteristic equation must have a root λ = 0. That is, the constant term of the characteristic equation should be zero. For the coupled system, the necessary condition is For the TF-TF and TF-M-TF loops, the necessary conditions are separately and For Hopf bifurcation to occur, the characteristic equation must have a pair of pure imaginary roots λ = ±iω, where ω is a positive real number. For the coupled system, the necessary condition is For the TF-M-TF loop, the necessary condition is Because the TF-TF loop is always a PFL, Hopf bifurcation is impossible. From Eqs. (36) and (37), the occurrence of Hopf bifurcation depends on θ y and θ z . By properly adjusting θ y and θ z , the condition for Hopf bifurcation can be never satisfied for any x.

LOGARITHMIC GAIN
In this section, we introduce the logarithmic gain of the systems of individual and coupled feedback loops based on several previous studies [5][6][7][8]. A single feedback loop can sometimes be regarded as self-regulation of a specific component via other components. Consider a feedback loop composed of N components C i (i = 1, 2, · · · , N ), each of which can only regulate its adjacent downstream component. The dimensionless concentration of C i is denoted by c i (i = 1, 2, · · · , N ) and the strictly monotonic function describing the regulation of C i+1 by C i is denoted by F i with i = 1, 2, · · · , N − 1 (that of C 1 by C N is denoted by F N ). We take this feedback loop as the self-regulation of C 1 via the other components. Assume the dynamics of this feedback loop are governed by the following ordinary differential equations: where η is the control parameter and 1 θi is the degradation rate of C i normalized to that of C 1 . At steady states, all the derivatives equal zero, thus yielding the following N algebraic equations: By eliminating c i (i = 2, · · · , N − 1), we have the equation for c 1 : Of note, the function T is the composition of functions, each describing one of regulatory interactions in the feedback loop, and the way of composition reflects the propagation of these interactions. Thus, the function T indicates how C 1 regulates itself via the feedback loop. In this sense, we term T a transfer function. Owing to the strict monotonicity of F i , T is a strictly monotonically increasing function for PFLs, or a strictly monotonically decreasing function for NFLs. We then define the logarithmic gain (G, called gain for short) as the derivative of the transfer function T with respect to c 1 in log-log coordinates: Because of the strict monotonicity of T , G is always positive for PFLs or negative for NFLs. To formulate the necessary condition for SN bifurcation, we calculate the determinant of the system's Jacobian matrix: where c i (i = 1, 2, · · · , N ) satisfy Eq. (39). For SN bifurcation to occur, the determinant must be zero, thus yielding According to the chain rule, this equation can be rewritten as This is just the necessary condition for SN bifurcation in terms of the gain. Based on the above definitions, the transfer function and gain of the TF-TF loop are and respectively. For the TF-M-TF loop, the transfer function and gain are and Similarly, we can calculate the transfer function and gain of coupled feedback loops as follows: and Notably, i.e., the gain of coupled feedback loops is just the sum of those of individual feedback loops. The necessary condition for Hopf bifurcation in the coupled system in terms of gains is which also depends on the timescales. For the TF-M-TF loop, the necessary condition is

Estimation of the width of the bistability region for a single PFL
In this section, we estimate the width of a bistability region (W B ) for a single PFL as a function of the maximum gain (G max > 1) and the width of the effective concentration range (W ). The G max is defined as the global maximum of G. The W is defined as the logarithmic difference between w 2 and w 1 , i.e., W = ln(w 2 /w 1 ) (w 2 > w 1 ) with w 1 and w 2 denoting the x values at half points (i.e., G(w 1 ) = G(w 2 ) = G max /2). The W B is defined as the logarithmic difference between η 1 and η 2 , i.e., W B = ln η1 η2 , where η 1 and η 2 are the upper and lower thresholds for the bistability region, respectively. The x values at η 1 and η 2 are denoted by b 1 and b 2 , respectively. According to Eq. (43), G(b 1 ) = G(b 2 ) = 1 (Fig. 1). To explicitly present the relationship between η and x, we define the transfer function in units of η,T : Thus, The width of the bistability region can be expressed as G can be expressed in terms ofT : which is equivalent to the expression of G in terms of T (x). Thus, That is, W B only depends on the part of G(x) 1.
To estimate the integral and (ln b 2 − ln b 1 ), we use an isosceles triangle to replace the bell-shaped part of G(x). This triangle has a height of G max and a base length of 2W (Fig. 1). Then, the integral is approximately It is easy to see Substituting Eqs. (58) and (59) into Eq. (57), we have the following estimate:

Estimation of the widths of two separate bistability regions for the coupled system
Similar to Eq. (57), the widths of two separate bistability regions, W BC1 = ln(η 1 /η 2 ) and W BC2 = ln(η 3 /η 4 ), can be determined by two bell-shaped parts of the G C (x) curve (more exactly, the parts of G C (x) 1), i.e., and where c 11 , c 12 , c 21 and c 22 are the x values at and c 11 < c 12 < c 21 < c 22 . Because of the 'bandpass' feature and large R, the two bell-shaped parts of the G C (x) curve can be approximately represented by G 1 (x) and G 2 (x), respectively, leading to and Notably, the right-hand sides of Eqs. and

Estimation of the width of the intermediate branch of stable steady states for coupled PFLs
Here, we estimate the width of the intermediate branch of stable steady states, W I , in the case of x 2m > x 1m . Similarly, where c 12 and c 21 (c 12 < c 21 ) are the x values at the lower (η 2 ) and upper (η 3 ) thresholds for the intermediate branch, respectively. That is, G C (c 12 ) = G C (c 21 ) = 1. For a large R, two bell-shaped parts of the G C curve can maintain the features of the G 1 and G 2 curves. Thus, these parts can be approximated by G 1 and G 2 , and we have c 12 ≈ b 12 and c 21 ≈ b 21 , with G 1 (b 12 ) = G 2 (b 21 ) = 1. Using the same approximation method mentioned above, we replace the bell-shaped parts of G 1 and G 2 with two isosceles triangles, respectively. These two triangles have heights of G 1max and G 2max and base lengths of 2W 1 and 2W 2 , respectively. As shown in Fig. 2, the integral in Eq. (61) approximately equals the negative total area of two small gray triangles, i.e., ).
As for (ln c 12 − ln c 21 ), it can be estimated by (ln b 12 − ln b 21 ). Furthermore, and Thus, W I is approximately

Estimation of the thresholds of R for tristability
The widths of two separate bistability regions, W BC1 and W BC2 , can be estimated at: and respectively. If W I > W BC1 + W BC2 , two separate bistability regions appear in two ranges of η; if W I < W BC1 + W BC2 , a tristability region is admitted. Thus, when W I = W BC1 + W BC2 , the approximate upper threshold of R, R 2,appr is given by If the local minimum of G C (x) between two local maxima is below 1.0, there are four SN bifurcation points, and tristability is possible; if the local minimum is elevated above 1.0, only two SN bifurcation points are admitted, and tristability is forbidden. Thus, the local minimum determines whether tristability is possible. As shown in Fig. 3, two isosceles triangles have an overlap when R is relatively small. P is the intersection of two sides of two triangles, and |P O| is the distance of P to the x-axis. We set |QO| = 2|P O| as the estimate of the local minimum of G C (x). Since we have When |QO| = 1, we get the approximate lower threshold of |R| for tristability: ). (80)

THE NUMERICAL METHOD FOR CHANGING ONLY R
In this section, we develop a numerical method in which we only change x 2m while keeping G 1max , x 1m , W 1 , G 2max , and W 2 unchanged. In this way can we merely alter R = ln(x 2m /x 1m ). Without altering the feature of the TF-TF loop, we only adjust two interactions: g(z) = 1 1+δz , and u(y) = n 3 , which contain four parameters. Adjusting n 3 , α and K in u(y) simultaneously can change x 2m and ensure that G 2max and K 2 are fixed.
First, we design two functions F G2max (n 3 , α, K) and F W2 (n 3 , α, K) to numerically calculate G 2max and W 2 from a set of (n 3 , α, K), respectively. Second, we numerically solve α and K for different n 3 from the following two equations: are the fixed values for G 2max and W 2 , respectively. These sets of (n 3 , α, K) can lead to identical G 2max and K 2 but different x 2m , thus generating distinct R.

THE CONDITIONS FOR DIFFERENT TYPES OF BIFURCATION DIAGRAMS AND DERIVATION OF THE THRESHOLDS OF R FOR TRISTABILITY IN SIMPLIFIED PIIN SYSTEMS
Here we construct several reduced P II N systems with piecewise linear gain functions which are a simplification of the nonlinear bell-shaped gain curves of the original systems. For PFL, the gain curve is an isosceles triangle with a height of G 1max and base length of 2W 1 , while the inverted bell-shaped gain curve for NFL is an isosceles triangle with a height of |G min | and base length of 2W 2 . Given the PFL should have a bifurcation diagram of type II, G 1max > 1. Next, we derive the conditions for different types of bifurcation diagrams with those piecewise linear gain curves.

Case 1:
x 2m x G B 0 1 x 1m x 2m , the G C (x) curve has only two intersections with G(x) = 1 for any R (Fig. 4). Thus, the coupled system always admits the bifurcation diagram of type II.
x 1m x 2m If W1 W2 > r 1 , the G C curve intersects with G(x) = 1 four times within a limited range of R (Fig. 5C). The lower threshold of |R| for four intersections, R ′ 3 , equals G1max+G2min−1 G1max W 1 when the middle minimum of G C just equals 1.0 (Fig. 5B), while the upper threshold, R ′ 4 , is G1max−1 G1max W 1 − W 2 when one of the local maxima of G C equals 1.0 (Fig.  5D). For |R| < R ′ 3 or |R| > R ′ 4 , only two intersections are allowed and the coupled system can admit bistability rather than tristability.

Case 2:
x 2m , the G C curve can have at most two intersections with G(x) = 1, and thus bistability is allowed but tristability is forbidden (Fig. 6). When |R| approaches 0, the maximum of G C becomes less than 1 due to G 1max + G 2min < 1, and bistability also disappears (Fig. 6A). The threshold of R for bistability is R W 2 when the maximum of G C equals 1 (Fig. 6B). If |R| < R ′ 5 , only the bifurcation diagram of type I is allowable; if |R| > R ′ 5 , type II is admitted and bistability is possible (Fig. 6C).
G1max−1 , two local maxima in the G C curve are less than 1.0 for R = 0 (Fig. 7A). When R is non-zero, one maximum rises but the other drops. Thus, the G C curve still has at most two intersections with G(x) = 1. However, the threshold of R for bistability is R G1max W 1 when the larger of the two maxima of G C equals 1 (Fig. 7B). If |R| < R ′ 6 , only the bifurcation diagrams of type I can be produced; if |R| > R ′ 6 , type II is admitted and bistability is allowable (Fig. 7C).
If W1 W2 > r 2 , the G C curve can have four intersections with G(x) = 1 when |R| approaches 0, indicating the possibility of tristability (Fig. 8A). However, when |R| is much larger, there are only two intersections, and thus only bistability x 2m x G C 0 1 FIG. 7: Illustration for the simplified G1(x), G2(x) and GC(x) curves with G1max + G2min < 1 and r1 < W 1 x 2m x G B 0 1 x 1m x 2m x G C 0 1 FIG. 8: Illustration for the simplified G1(x), G2(x) and GC(x) curves with G1max + G2min < 1 and W 1 W 2 > r2. Only x2m is varied for different panels. arises (Fig. 8C). The threshold of R, used to discriminate between these two cases, is R G1max W 1 − W 2 when the smaller maximum of G C equals 1 (Fig. 8B). If |R| > R ′ 7 , only the bifurcation diagrams of type II are allowed; otherwise, either type III or type IV-VII is admitted.
It is noted that r 3 > r 2 always holds for G 1max + G 2min < 1. As mentioned above, the relationship between W1 W2 and r 3 determines whether the bifurcation diagram is of type III or IV-VII. Thus, if W1 W2 > r 3 , type IV-VII is allowed and tristability is possible; if r 2 < W1 W2 < r 3 , only type III is admitted.

ANALYSIS OF TWO OTHER TYPES OF MULTIPLICATIVELY COUPLED FEEDBACK LOOPS
Although the results in the main text are based on a particular type of coupled feedback loops, the analysis method is also applicable to at least two types of coupled feedback loops given the loops are multiplicatively coupled.

A mutual inhibitory loop combined with a self-regulation loop
FIG. 9: Schematic of a mutual inhibitory loop combined with a self-regulation loop. The system can be decomposed into two separate loops. The arrow-and bar-headed lines denote activation and inhibition, respectively.
Such a motif often contains two components denoted by X and Y . Their concentrations are x and y, respectively.
The kinetic equations for this motif are as follows: where f (x) describes the autoregulation of X, while g(y) and h(x) characterize the mutual regulation between X and Y .
For two individual loops in the system, the transfer functions are and Thus, the corresponding logarithmic gains are and respectively. For the coupled feedback loops, the transfer function and gain are separately and Denote the steady-state value by x. The characteristic equation for the coupled loops is with By setting λ = 0, we obtain the necessary condition for SN bifurcation:

Two mutual inhibitory loops sharing one component
This type of motif usually contains three components. Denote them by X, Y and Z, and their concentrations by x, y and z, respectively. The kinetic equations for this motif are as follows: (96) where f (y) and h(x) describe the mutual regulation between X and Y , while g(z) and u(x) characterize that between X and Z.
For two individual loops in the system, the transfer functions are and Thus, the corresponding logarithmic gains are and respectively. For the coupled loops, the transfer function and gain are separately and Denote by x the steady-state value. The characteristic equation for the coupled feedback loops is where By setting λ = 0, we obtain the necessary condition for SN bifurcation:   Table S1 The sampling ranges of the parameters n 1 n 2 n 3 δ K α κ PP 1-10 1-10 1-10 0.1-100 0.01-100 1-100 0.1-100 PN 1-10 1-10 1-10 0.1-100 0.01-100 0.01-1 0.1-100 Parameter sets are sampled uniformly on a logarithmic scale in a 7-dimensional parameter space using the Latin hypercube sampling method. Parameters θ z and θ y are kept at 1.0. and P II denote the individual PFLs with bifurcation diagrams of type I and II, respectively. P I P I , P I P II , P II P I and P II P II denote different kinds of dual-PFL systems, while P I N and P II N denote two kinds of coupled PFL and NFL systems. The first and second items denote the types of TF-TF loop and TF-M-TF loop, respectively. -0.56 9.67 4.50 0.0517 10 0.7 *For each system, the values of n 3 , α and K are changed such that only x 2m is altered while the other characteristic quantities, i.e., G 1max , x 1m , W 1 , G 2max (G 2min ), and W 2 are fixed. The parameters besides those shown in the table are always set as follows: n 1 =2, n 2 =2, δ=10, and θ y =1.