Identification of optimal parameter combinations for the emergence of bistability

Bistability underlies cellular memory and maintains alternative differentiation states. Bistability can emerge only if its parameter range is either physically realizable or can be enlarged to become realizable. We derived a general rule and showed that the bistable range of a reaction parameter is maximized by a pair of other parameters in any gene regulatory network provided they satisfy a general condition. The resulting analytical expressions revealed whether or not such reaction pairs are present in prototypical positive feedback loops. They are absent from the feedback loop enclosed by protein dimers but present in both the toggle-switch and the feedback circuit inhibited by sequestration. Sequestration can generate bistability even at narrow feedback expression range at which cooperative binding fails to do so, provided inhibition is set to an optimal value. These results help to design bistable circuits and cellular reprogramming and reveal whether bistability is possible in gene networks in the range of realistic parameter values.


Introduction
Bistability is the co-existence of two stable equilibria. It is a hallmark of self-organization in dynamical systems, and has a wide range of occurrence and applications [1][2][3].
Bistability creates two distinct cell states or phenotypes in genetically identical cells. It can also store memory of past stimuli by hysteresis and promote cell-cycle oscillations [4,5]. These effects underlie cellular differentiation, adaptation to varying environments by bet hedging, and may have practical implications for cellular reprogramming and antibiotic resistance [6][7][8][9][10].
The presence of positive feedback is a necessary condition of bistability [11], and it is frequently encountered in networks that control the above processes [12,13]. The existence of positive feedback is however not a sufficient condition for bistability; the feedback components must also participate in reactions that generate ultrasensitive responses. More precisely formulated, the emergence of bistability requires the following condition: if the feedback loop is opened, the output response f of the resulting openloop must be a sigmoidal function of the input w with a logarithmic sensitivity w S f larger than one; in other words, w ( ) f has to be an ultrasensitive response [14,15]; f ln ln f . This sufficient condition for bistability is satisfied even in the simplest one-gene feedback loop provided the self-activating transcriptional factor (TF) participates in an ultrasensitive reaction. This is the case when the TF binds cooperatively to the promoter of its own gene or the TF forms dimers to be active.
Bistability can function robustly if its parameter range is broad enough. The logarithmic sensitivity can provide clues about the bistable range of a parameter, as evidenced by the positive correlation between cooperativity and the size of the bistable domain in the parameter space [16][17][18]. w S f is proportional to the degree of cooperativity (Hill-coefficient). With moderately cooperative binding of the TF to DNA, bistability is restricted to a narrow range of TF affinity; this range broadens as cooperativity increases [17].
This intuitive relationship may be blurred for other types of ultrasensitive reactions, including protein homo-dimerization, sequestration and multistep phosphorylation [19][20][21]. For example, reducing the production rate of a dimeric TF increases sensitivity but also lowers the concentration of the active TF quadratically, and hence the feedback intensity.
We hypothesized that the trade-off between such effects may maximize the bistable range of a parameter. We sought to identify general conditions for this maximization in a multi-dimensional parameter space. The derived general relations were then tested on prototypical feedback loops incorporating three different ultrasensitive reactions: homo-dimerization, inhibition by sequestration and cooperative binding.
Identifying the maximal bistable range of a parameter permits the robust design of bistable networks and provides answer to the question as to whether bistability is possible in endogenous networks when only the topology and few experimental data are available, which is often the case.

Cusp points delimit the bistable range of two parameters
First, we examined how bifurcation points delimit the bistable range in the parameter space. For this purpose, two key parameters were varied in a simple positive feedback loop with cooperative binding. The system is described by the following equation at equilibrium: The Hill number n is larger than one if the binding is cooperative. g is the protein degradation rate constant. K is the equilibrium dissociation constant of the TF-DNA binding reaction, and b P is the basal production rate. The ratio of the basal ( ) b P to the maximal ( ) V m production rate,b = b V , P m / is a dimensionless parameter and can be used to compare different feedback loops. The inverse of this parameter is the feedback expression (dynamic) range. Figure 1 illustrates the equilibrium manifold of system (1). For low values of b, the manifold is S-shaped and has two turning points, which represent fold (saddle-node) bifurcation points. The two fold bifurcation points determine the bistability range of K in which system (1) has three distinct equilibria. As b increases, the two fold points move toward each other and collide at a cusp bifurcation point.
Since the intersection of the fold curves is tangential [22], the cusp point defines the maximum of parameter ranges in a two-parameter space in which bistability can occur (figure 2(a), right).
The condition for bistability is defined in terms of logarithmic sensitivity of the open-loop response. Therefore, it is important to see how bifurcations arise in this context (see also appendix A). The intersection determines the number of equilibria and their values. If w ( ) f intersects with the equivalence line tangentially, the system undergoes a fold (saddle-node) bifurcation (figure 2(a), table A1). When b is increased and K is adjusted so that the equivalence line intersects w ( ) f at its inflection point, then these values of b and K define the position of the cusp bifurcation point (figure 2(a)).

Absence of extrema on the fold-curves
The cusp points represent the extremal values of the bistable domains since the fold curves are monotone with respect to the parameters, i.e. they do not turn back (figures 1 and 2). If the fold-bifurcation curves were not monotone, they might delimit the extrema of bistable parameter domains. To explore this possibility, we defined a general condition for the existence of non-monotone fold curves.
The locus of fold points satisfies C (1) and C (2) in table A1. The fold curve has an extremum in a 2-dimensional parameter space a a { } , i k if the following expression equals zero, plane (dark purple) are indicated. Between the fold bifurcation curves, the system has three equilibria and is bistable.
where C stands for C (1) and C (2) . This is the case if Thus, an extremum in the fold curve can arise if the open-loop response is non-monotone with respect to a reaction parameter. Most well-characterized reactions in gene networks display monotone responses to variations in the relevant parameters. An exception may be the bellshaped response, which is the product of Hillfunctions for activation and repression. The bellshaped response can arise in promoters subject to transcriptional interference or feedforward regulation, and is characterized by a low noise behavior [23,24]. It has a non-monotone dependence on the state variable (figure 2(b), left). Similarly, ¶ ¶ = f K 0 / for some parameter values because the exchange of w and K preserves the response function [23]. However, ¶ ¶ = f K 0 / and the fold condition w cannot be met simultaneously (see appendix B).
Consequently, there is no extremum in the fold curves. In fact, the overall bistability range is very similar to that of the simple cooperative feedback despite the extremum in the bell-shaped response function (figure 2, right).
Thus, the cusp point is generally expected to delimit the range of two parameters, in which bistability is permitted.

Derivation of general conditions for the existence of extrema in the locus of cusp points
In order to detect the parameters that can maximize the permitted bistable range of another parameter, we set up the following equations, which establish the necessary conditions for the existence of a cusp bifurcation point, upon eliminating three degrees of freedom, The Jacobi determinant in the numerator of equation (4) is given as (see also appendix C): Assuming the genericity condition w ¶ ¶ ¹ f 0 3 3 is satisfied, a maximum exists if the second factor in (5) is zero. This is equivalent to: The choice of parameters for α j and α k is made easier by identifying a specific parameter that satisfies w a w a = ( ) ( ( ) ) f f g ; ; 1 .
C C α C is covariant with the state variable; therefore, we term it covariant parameter, and g is an arbitrary smooth function with non-zero derivative at a .
C It can be shown that =  Next, we applied condition (8) to analyze feedback prototypes in different parameter spaces. We always include the parameter pair representing the feedback expression range and binding affinity of the TF to DNA. This facilitates the consistent comparison of the different feedback loops.

Detection of reaction pairs that maximize bistability in the one-gene loop with sequestration
First, we examined a positive feedback circuit inhibited by sequestration. Sequestration is a common but often hidden source of non-linearity that can significantly alter system behavior, which underscores the importance of this prototype of bistable feedback loops [19,20,25]. In this feedback circuit, the TF (A) binds to the promoter and activates the production of its own mRNA (M); A can also be sequestered by an inhibitor (I) into a complex (C), which cannot bind to the promoter ( figure 3(a)). The promoter response is described by a Hill-function. The system is described by:  To further simplify the system before opening the loop, it was reduced to two variables and nondimensionalized: The SI units are given in brackets. Upon opening, the following equation is obtained: This opening corresponds to 'breaking' A into two parts. In the resulting open-loop system, two different components arise from A. The first part serves as the input w and the sole reaction it imitates is the binding to the promoter. The second part, A' is retained in all the other reactions and becomes the output of the response function ( figure 3(a)).
The resulting open-loop response and the bifurcation diagram (figure 3(b)) are similar to that of the simple positive feedback loop (figure 2(a)) when the feedback expression range, b -, A 1 and the binding affinity of the TF to the DNA, k, are varied. On the other hand, the loop with the sequestration has two other parameters and they may have non-trivial impact on system behavior. When b I , the inhibitor production rate, is increased, the open-loop response is lowered (figure 3(c)). At the same time, the maximal value of w S p displays a non-monotone response, pointing to the possibility that there is an optimal value of b I that maximizes the bistable range of b A .
For the above opening (12), the covariant parameter is k. Thus / Applying the conditions for the double maxima, C (1) , C (2) and equation (8), to the response function constructed from equation (12) reveals that all three possible pairs of parameter combinations of k have joint maxima.
has the most permissive conditions: In order to focus on sequestration, no cooperativity was included in the analysis, = n 1. Thus, sequestration was the sole source of ultrasensitivity.
The four reduced parameters contain multiple reaction parameters that often belong to consecutive or related reactions and therefore, they describe a specific property, exemplified by the feedback expression range, b -. A 1 When a reduced parameter contains unrelated reaction parameters, as well, it is possible to focus on a dominant or unique reaction parameter by varying it, while keeping the others fixed. In this way, c reflects the maximal production rate of the transcriptional activator; k reflects the TF-DNA binding affinity and b I the inhibitor production rate. For The optimal inhibitor production rate that maximizes the bistable range of b A is b = 27.18 I when c = 35 (figure 4), revealing that the maximal bistable feedback expression range is attained when the inhibitor production rate nearly reaches the level of the maximal activator production rate. To visualize the maximal bistable range, we plotted the locus of cusp points as b I was varied ( figure 4). It is interesting to explore how realistic values of b I can promote bistability (appendix E). We explored two limiting scenarios with respect to the feedback expression range, b -. If b -A 1 is large, then bistability is possible with very low inhibitor production rate or affinity. For example, for b = -10 , 14). Since b I contains parameters for the rate of production of the inhibitor and for its binding to the activator, the following combinations were considered.
If transcription rate is low (table E1), then l = 0.28 min −1 (appendix F). Consequently, the inhibitor affinity, has to be less than 350 nM for bistability to emerge ( figure 5(a)). On the other hand, when transcription rate is high even a very weak binding ( = ) K 35 mM inh is sufficient to support bistability. This points to the possibility that a highly expressed inhibitor protein can generate bistability even if the binding is weak, potentially representing non-specific protein-protein interaction, which may be overlooked in standard experimental conditions. Is bistability possible if the feedback expression range is narrow, e.g. b = -4?  5(b)). Assuming weak and strong transcription rates, 5.8 and 590 nM inhibitor affinities are needed, respectively, to generate bistability. This indicates that bistability can emerge in the range of typical protein binding constants even with a narrow, fourfold, dynamic range of feedback expression.

Analysis of the one-gene loop with protein homo-dimerization
Next, we examined if maxima exist in the locus of cusp points for a one-gene model with protein homodimerization ( figure 6). The TF binds to the promoter only as dimer C. To focus on the effect of homodimerization, the dimer binds to the promoter non-cooperatively. The system was non-dimensionalized and then opened (appendix G). The main reduced parameters are dominated by the following reactions: c by the maximal rate of activator production, k q by the binding strength in the homodimerization, and k d by the TF-DNA affinity. b -1 is the feedback expression range. The covariant parameter is k d .
No positive real (physical) solution was found for any of the possible parameter pairs k a { } , d k that satisfies the conditions for the cusp points and equation (8). Thus, the cusp points delimiting the bistability domains lack extrema for positive parameter values. This can be illustrated by varying the dissociation rate of the dimer into the monomer (figure 6). At strong dimerization (k = 10, the bistability domain is small. This is because the TF concentration has to be below K dim to have sufficiently high logarithmic sensitivity in the open-loop response. Since the concentration of the TF is around 2000 nM in the higher stable equilibrium, the ultrasensitivity necessary for bistability is possible only if the basal expression is around 100 times less than the maximal expression, resulting in concentrations of around 20 nM ( figure 6). Indeed, bistability appears only when b < -10 . 2 The TF-DNA binding affinity is low (higher than 200 nM) in this bistable domain.
A reduction of the dimerization expands the bistable range of b (figure 6). Above » K 2000 nM dim (i.e. belowk = ) 0.1 , q the expansion of the bistable range levels off. Indeed, if k q is set to zero by eliminating the dimer decay, b approaches 1/8, which is the same value as the positive feedback loop with a cooperative promoter, = n 2 (appendix G.2). Thus, in this range of low dimerization strengths (high K dim , low κ q ), a further weakening of dimerization does not increase ultrasensitivity and the lower dimer concentration is compensated by the stronger binding of the TF to DNA. This explains the absence of maximum in the locus of cusp points. Having no optimal value for the dimerization binding, the feedback loop can generate a robust bistability in the face of variations in K dim . In order to assess how the affinities affect the bistability range we explored the relationship between k 1 and k 2 at the cusp points. There is an inverse relation between k 1 and k 2 in the bistable domains. At k = 5, 1 the bistable domain is narrow ( figure 7(a)). As k 1 is decreased, k 2 increases and the bistable range of b 1 expands. At the extremum in b , 1 the two affinities have similar values ( figure 7(a)).
To assess the above relation more quantitatively, symmetric conditions were enforced in the reaction scheme. The Hill coefficients were set equal. Secondly, β 1 and b 2 were required to be equal at the extremum of the cusp points (47). With these two constraints, the two affinities were indeed equal at the maximal point (48). This confirms that the broadest bistable range in β 1 and b 2 is attained when the affinities are equal. It also hints to why the values of k 1 and k 2 are similar even when the two basal production rates are different.
It is interesting to observe that the bistability range of β 1 declines very rapidly when k 1 increases beyond the optimal value ( figure 7(a)). This behavior is even more pronounced when the Hill coefficient assumes larger values. In this case, a new extremum arises, which delimits k 1 at low affinities. Thus, the extremum is surrounded by very asymmetric range of k 1 values in which bistable domains arise ( figure 7(b)).
The extrema in this feedback loop can be compared to those with other forms of binding reactions. In both homo-and heterodimerization two proteins interact. Thus, cooperative binding to two sites in the promoter permits a consistent comparison. If the cooperativity of the binding to the two sites is large then the Hill coefficient, = n 2 [27]. For = n 2, the extremal values in the normalized basal expression is 1/8 (49), which is the same as the limiting value in homodimerization but less than 1/3, the value found for the loop with sequestration with realistic parameters (13).

Optimal parameter values for the broadest bistable parameter range in the two-gene loop with mutual repression
Positive feedback can be also realized by two mutually inhibiting repressors, also known as the toggle-switch [16].
The locus of cusp points reveals a behavior similar to that of the two-gene loop with activators (appendix I, figure (8)). The main difference is that in this case the affinities of the two transcriptional factors are similar in all the bistable domains. Weak binding affinity of one repressor entails weak binding affinity of the other repressor in the bistable domains.

Discussion
3.1. Maximizing sensitivity with respect to two parameters maximizes the bistable range of a third parameter The logarithmic sensitivity provides a practical tool to predict and analyze bistability. Sensitivity can be easily obtained from the open-loop response fitted to experimental data.
Maximizing sensitivity was shown to specify the broadest dynamic range of a single parameter in a simple positive feedback loop with cooperative binding [17], similar to figure (1). When such an approach was applied for more complex feedback loops, a high statistical correlation was seen between the maximal sensitivity and bistability range [28]. It has remained unclear how to formulate a condition so that the correlation becomes full (i.e. one) and what analytical relations can be obtained for that.
In this work, we derived a general rule that stipulates that the largest bistable range of a parameter, i.e. extremum of a bistable domain, (equation (8)), is attained when the sensitivity of the response function has a maximum against two other parameters. We then verified whether such parameter pairs were present in specific network topologies, by obtaining analytical relations (e.g. equation (13)). This allows drawing general conclusions on specific network topologies since the entire parameter space was captured. These relations revealed that the prototypical feedback enclosed by protein dimers has no such  reaction pairs. Conversely, loops with sequestration and two-gene loops with cooperative binding do have such reaction pairs despite the distinct nature of the sequestration and cooperative binding to DNA. Also, loops having similar reactions, such as homo-and heterodimerization, (sequestration) fall in distinct categories with respect to the existence of maxima in bistable regions.
This suggests that not the chemical-kinetic nature of a specific reaction but the interaction of all reactions determines whether or not such parameter pairs are elicited. Most complex networks that display bistability are derivatives of a few prototypical network topologies and reactions [28][29][30][31][32]. It will be interesting to explore how the number of maximizing reaction parameter pairs varies as the complexity of a network increases.

Maximizing bistability and the design of bistable networks
In order to design a bistable gene network, an ultrasensitive reaction must be included. Furthermore, the system parameters have to assume values within the bistable parameter range. If the accessible parameter range does not overlap with the bistable range one can alter the reaction parameter physically or biologically to broaden the accessible range ( figure (3)). However, there are often limits to how much the accessible range can be widened. For example, the feedback expression range is often delimited by the lowest possible value of basal transcription rate, also known as leakage ( figure 1(b)) [26]. In some promoters, the basal expression has counter-intuitive determinants. For example, if the number of activator binding sites is reduced in the promoter the leakage increases despite the fact that the maximal expression is reduced [26]. In this case it is very difficult to broaden the range of the feedback expression.
If it is not possible to enlarge the accessible range of a reaction parameter in the cells, one alternative is to expand the bistable range. Our results provide support for this possibility: the bistable range of a parameter can be expanded by tuning other parameters to enable the emergence of bistability. Even if the range of particular parameter in a system is severely constrained, many others may be broadly tuned.
The conditions we derived in this work identify the pairs of parameters that maximize the bistable range of a parameter under consideration (8). For example, knowing the affinity of the inhibitor to the TF and the maximal transcription rate of the positive feedback loop permits the calculation of the optimal inhibitor production rate. Tuning this parameter to its optimal value makes it possible for bistability to occur even if the feedback expression range is very narrow ( figure 5). When the parameter is tuned away from this optimal value the bistable feedback expression range shrinks rapidly.
In contrast, the loop with homodimerization does not have optimal values, and the bistable feedback expression range remains relatively constant and close to the limiting value as soon as sufficiently weak dimerization strength is reached ( figure 6). This indicates that bistable feedback expression range of the loop with homodimerization is not highly sensitive to variations in the dimerization strength, reflecting the absence of double maxima, which typically makes the bistable ranges sensitive to parameter variations.
In the two-gene loops with cooperative binding, the two TFs must have similar or equal affinities to maximize the bistable feedback expression range. For the mutual repression loop, the affinities are similar also in other bistable domains, where the affinities are not optimal ( figure 8). Interestingly, the opposite relation is observed in loop with mutual activation. The further away the value of an affinity is shifted from the optimal value, the more different the affinities of the two TF have to be to generate bistability (figure 7).

Detection of bistability in endogenous networks
Endogenous bistable networks are subject to constraints similar to that encountered in the design of the networks. Further constraints arise during evolution when a reaction controls different biological functions or a single function in different environments with conflicting demands on the reaction parameters [33][34][35].
Our results may shed light on the prevalence of the different types of feedback loops. Interestingly, the mutually activating and repressing loops have identical maximal bistable feedback expression ranges (figures 7 and 8). That may explain why both network topologies have been identified in the model organism budding yeast [29].
The conditions for the extremal points also help to assess whether bistability is possible in endogenous networks when few data are available, which is typically the case. Most of the time, the network topology is identified and gene expression data can be easily measured, which can also define feedback expression range. For this reason, we compared the extremal value of the feedback expression range in prototypical feedback loops. Interestingly, heterodimerization (sequestration) outperforms both cooperativity and homodimerization since only sequestration permits the occurrence of bistability if the feedback expression range is less than eight.
Moreover, our results suggest that unexplained bistability may arise in positive feedback loops with no obvious ultrasensitive reaction because even very weak, and hence unidentified, interaction with inhibitors can generate bistability by sequestration.
Determining the limits of bistability in networks is particularly important in networks that control cellular differentiation [36]. is converted to the following one-dimensional expression: The open-loop response function w a ( ) f ; is then given by replacing x in a ( ) f x; with w: It is clear that the equilibrium condition C (1) for the open-loop re-closes the loop.
The bifurcation conditions for the multidimensional closed-loop system [22] are also given in table A1.

Appendix B. Fold points for the feedback loop with the bell-shaped response
The positive feedback incorporating a promoter with the bell-shaped response is described at equilibrium by Thus, the open-loop response function is given by The bell-shaped response is the product of Hillfunctions for an activator and a repressor. It can be seen that the exchange of K and w does not change the equation and the response function has maxima with respect to both K and w.
The fold condition and ¶ ¶ = f K 0 / have to be satisfied simultaneously for the fold curve to display a maximum:   Table A1. Conditions for fold (C (1) and C (2) ) and cusp (C (1) , C (2) and C (3) ) bifurcations. The brackets denote the inner product.
Multi-dimensional a One-dimensional with w a q is the null eigenvector of the Jacobian matrix, ⋅ = J q 0 so that it is normed to one, < > = q q , 1 .H is the Hessian matrix, and  q is the adjoint null eigenvector of the Jacobian Appendix D. w S f as a function of the covariant parameter If there exists a parameter α c for which the response function can be re-written as at the cusp point. The proof can be seen by applying constraints C (1) , C (2) and C (3) to equation (22): Then the sensitivity of the α c -derivative of the response function is

Appendix E. Biologically relevant ranges of parameter values
Specific model parameters were representative of reactions in yeast and they fall into a realistic range of values, as detailed below. The typical yeast mRNA half-lives range between 1 and 10 min (decay rate constant, μ=0.07-0.7 min −1 ) [26,37]. The average protein half-life is 120 min but the half-life of some regulators is shorter by an order of magnitude [38], which results in typical decay rate constants of μ=0.006-0.06 min −1 .
The translation rate p t ranges from 6 to 20 min −1 [39].
Most transcriptional factors and their regulators have mean concentrations between 200 and 2000 molecules/cell [40]. For nuclear proteins, this corresponds to 100-1000 nM, assuming a nuclear volume of 2 μm 3 [37]. Based on these data, the typical RNA production (transcription rate, p R ) and the lumped protein production rates (translation rate, p P ) can be calculated (see table E1).
The relation between p R and p P is given by The equilibrium dissociation constant of proteins ranges between 0.1 and 1000 nM [20]. The typical association rate of proteins is around k a = 1.66×10 6 M −1 s −1 =0.1 nM −1 min −1 [41].  with the following parameters: Parameters at the extremal points can be expressed analytically: