TRANSITION BETWEEN MONOSTABILITY AND BISTABILITY OF A GENETIC TOGGLE SWITCH IN ESCHERICHIA COLI

. In this paper, we investigate a genetic toggle switch in Escherichia Coli , which models an artiﬁcial double-negative feedback loop with two mu- tually repressors. This model is a planar diﬀerential system with three parameters, one of which is an integer power n ≥ 1, in the case that repressors 1 and 2 multimerize with n and 1 subunits respectively and its equilibria are decided by a polynomial of degree n +1. Since one hardly solves such a polyno- mial equation, a known result on bistability was given by omitting some small terms under the assumption that the promoters are strong and the expression ratio between the ON state and the OFF state is large. In this paper, determining distribution of zeros qualitatively for the polynomial of high degree, we analytically discuss on the system without the assumption and completely give qualitative properties for all equilibria, which corrects the known result of bistability. Furthermore, we prove that there may occur in the system a codimension 2 bifurcation, called cusp bifurcation, which is a collision of two saddle-node bifurcations and manifests the transition between bistability and monostability. We exhibit the global dynamics of repressors in various cases by analyzing equilibria at inﬁnity and proving nonexistence of closed orbits.

1. Introduction. Biological cells are extremely sophisticated systems which can execute a lot of complex behaviors, for example, response with an appropriate change dynamically as monitoring their environment. This attracts great attention to investigate those complex behaviors, one of which is genetic circuit. Explanation to its components of genetic circuit with DNA and RNA can be found from [6,10,14,19,21,26]. Along with the advance of genomics and genetic engineering, a subject known as synthetic biology ( [16,30]) is developed to construct synthetic gene circuits which were not found in nature. One of synthetic gene circuits is toggle switch, constructed in Escherichia coli by Gardner, Cantor and Collins ( [10]) in 2000. Toggle switch is composed of two repressors and two constitutive promoters, where each promoter is inhibited by the repressor transcribed by the opposing promoter, as shown in Figure 1. In this switch, it is expected to flip between two stable states (high state and low state ( [10,22])), which refers to as a doublenegative feedback loop ( [7,8,22]). It forms a synthetic and addressable cellular unit for storing "memory" and involves many fields such as biotechnology, biocomputing and gene therapy. In [10] the change of two repressors with respect to time in toggle switch was modeled by the system of two differential equations where u and v present the concentrations of repressor 1 and repressor 2 respectively, a and b are the effective rates of synthesis of repressor 1 and repressor 2 respectively, and β and γ are the cooperativities of repression of promoter 2 and promoter 1 respectively. In the biological sense, u, v, a, b, β and γ are all positive.  The main attention of [10] was paid to bistability ( [1,7,12,20,23]) of system (1), a state having two stable equilibria, which can store memory as shown in [18,27,32]. Although it is difficult to solve zeros from the right hand sides of system (1), the authors of [10] plotted numerically the horizontal isocline H : v = b/(1 + u γ ) and the vertical isocline V : u = a/(1 + v β ) as in Figures 3 and 4 for β > 1 and γ > 1, showing that system (1) has either exactly one equilibrium or three equilibria. As surveyed in [31], the authors of [10] assumed to work on strong promoters (i.e., parameters a and b are both large) with large expression ratio between the ON state and the OFF state. In such a case, the coordinates of equilibrium A 2 : (u 2 , v 2 ) satisfy that u 2 v 2 and therefore were approximated as u 2 ≈ a and v 2 ≈ b/a γ and the coordinates of A 0 : (u 0 , v 0 ) were similarly approximated as u 0 ≈ a/b β and v 0 ≈ b, which enables them to obtain the condition to ensure a vanished determinant D of the Jacobian matrix of system (1) at A 2 (or A 0 ). Further, they omitted the relatively small term log(βγ)/(βγ) from the equality log a = log b γ + log(βγ) βγ (or log b = log a β + log(βγ) βγ ), an equivalent form of the above condition, and reduced the condition approximately to the following Thus, (2) determines those lines in the (log b, log a)-coordinate as shown in Figure 2, which partition the first quadrant into regions for D > 0 or D < 0 and therefore give parameter conditions for bistability and monostability. In 2007, Xiu ( [33]) investigated the system ) ηγ − v, a modification of (1) mentioned in [15,25], with a tuple p := (a, b, β, γ, η, K) of six parameters and an input [IP T G]. Considering p to be a definite choice with a small random change, i.e., p = p (I +δY ), where p := (156.25, 15.6, 2.5, 1, 2.0015, 2.9618 ×10 −5 ), δ = 0.1, I is the unit matrix and Y is a 6×6 matrix composed by 6 random row vectors which distributed in the six-dimensional cube (−1, 1) 6 uniformly, they numerically computed for the solution curve of the normalized v, which is shown in error bars with respect to varying [IP T G]-input, and exhibited a good agreement between numerical error bars and experimental error bars.
As observed above, one hardly simulate the existence of exact two equilibria numerically, the transitional case between one equilibrium and three ones, although the continuity of the vector field makes sure of the existence of two equilibria. Actually, the transitional case as well as the corresponding values of parameters is important to the change of stability and geometric phase structures, which leads to an investigation to bifurcations of equilibria, showing how equilibria arise and how the stability of equilibria changes. Such an investigation needs a more complete discussion to parameters a, b, β and γ than the approximative analysis with coordinates u 0 , v 0 , u 2 and v 2 and condition (2). As shown above, equilibria of system (1) are determined by the equation of fractional power (u − a)(u γ + 1) β + b β u = 0.
One hardly obtains a coordinate for any one of those equilibria even if β and γ are restricted to be integers m ≥ 1 and n ≥ 1 respectively, not mentioning the computation of eigenvalues at those equilibria.
In this paper, for convenience in computation, we study equation (1) with integers β = 1 and γ = n, i.e., the system where a > 0 and b > 0. As shown in the beginning of next section, finding equilibria is reduced to finding real zeros of a polynomial of degree n + 1 with symbolic coefficients, but one hardly solve such a polynomial. Our strategy is to abandon the routine of qualitative analysis with exact coordinates and eigenvalues but use inequalities to determine signs of some quantities. In section 2, making monotonic interval partition in the domain and applying elimination for complicated terms, we determine the distribution of real zeros for the polynomial of high degree, which gives conditions of parameters a, b and n for system (3) to have exactly one, two and three equilibria separately. In section 3 we give the qualitative properties for those equilibria, applying the same techniques as above to handle systems of polynomials of high degree given by the determinant, trace, and discriminant of the Jacobian matrix, in which coordinates of equilibria are not solved yet. In section 4, we further prove that there may occur in the system a codimension 2 bifurcation, called cusp bifurcation ([17, pp.301-307]), which is a collision of two saddle-node bifurcations and manifests the transition between bistability and monostability. In section 5, summarizing the above results, analyzing equilibria at infinity and proving nonexistence of closed orbits, we exhibit the global dynamics of repressors in various cases. We finish the paper by returning to the problem of bistability in section 6, in which we can conclude Theorems 2.1 and 5.1 to give a condition for bistability even without the assumption of [31] that the promoters are strong and the expression ratio between the ON state and the OFF state is large. In the case of [31], our theorems show a smaller size of the bistable region, which corrects the boundaryΥ ± of bistability, as shown in Figure 15. Our correction is demonstrated by numerical simulation (Figures 16 and 17) with parameters chosen betweenΥ + and Υ + as well as betweenΥ − and Υ − .
2. Distribution of equilibria. As a routine, equilibria of system (3) are determined by the right-hand sides of (3) being equal to zero, i.e., the system However, it is difficult to find real roots of the equation of degree n + 1 which is deduced by eliminating v in (4). Actually, we need to find zeros of Φ on the interval (0, a) because in the first quadrant the first equation of (4) implies The simplest one is that n = 1, where the polynomial Φ is quadratic, i.e., Φ(u) = u 2 − (a − b − 1)u − a = 0, and has two different real zeros Clearly, only u + lies on the interval (0, a). In this case system (3) has exactly one equilibrium A : (u 0 , v 0 ) in the first quadrant, where u 0 = u + and v 0 = a/u + − 1. We claim that A is a stable node. In fact, computing the Jacobian matrix of system (3) at A, we get the trace and the determinant Difficulties mainly come from the case that n ≥ 2, where we hardly solve zeros from equation (5) of degree n + 1. Those difficulties make the usual routine of equilibrium analysis unavailable. We abandon pursuing the determination of signs for the determinant, trace and discriminant at accurate coordinates of equilibria, but estimate those quantities near equilibria. For convenience, we need the following notations: Theorem 2.1. For n ≥ 2, the following assertions hold in the first quadrant: Proof. As known at the beginning of this section, the equilibria of system (3) are determined by positive zeros of Φ. Obviously, Φ has at least one zero in the interval (0, a) because of the continuity of Φ and the inequalities but it is hard to find all zeros for an irreducible polynomial of degree n+1 in general.
In the case that either 0 < a < a * (b) or a = a * (b), we have Φ (u) ≥ 0 for all u ∈ (0, a) and Φ has at most one zero, as shown in Figures 5 and 6. It follows that Φ increases strictly on (0, a), implying the uniqueness of zero of Φ in (0, a) if it exists. Additionally, the zero of Φ exists because of (8). Letũ * denote the zero, which determines a unique equilibrium of system (3). Although the degree n + 1 of the polynomial Φ prevents from solving the value ofũ * , it is still possible to determine which of sub-intervals (0,ũ 2 ] and [ũ 2 , a) the zeroũ * locates in. This can be done by checking the sign of Φ(ũ 2 ) because of (8). Note that Φ(ũ 2 ) = aF(a, b), where 1/n as denoted in (7) if b > b 0 . Moreover, Φ(ũ 2 ) > 0 (or < 0) if a < a 0 (b) (or a > a 0 (b)). Furthermore, one can prove that a * (b) = a 0 (b) (or <, or >) if and only if b = b * (or >, or <), where b * is given in (7). It implies the critical case that Φ(ũ 2 ) = Φ (ũ 2 ) = 0 (13) if and only if a = a * (b) and b = b * . For convenience, use the notation Thus, we discuss in three cases: In case (i), from (12) we see that Φ(ũ 2 ) < 0. It implies thatũ * lies in the interval (ũ 2 , a), giving the results of the theorem in case (i) except a > a * (b).
We only need to return those situations (S1)-(S5) on signs of Φ(u 0i ) (i = 1, 2) to relations of parameters a, b and n, which gives conditions for different numbers of equilibria. For this purpose, we compute those values of a, b and n in the critical cases (S2) and (S4), i.e., respectively. These are the same requirement as in (13). We hardly do the same routine as in [9,13,28], solving u 0i as a function of a, b and n from one equation of (22), because the two equations are of degrees n + 1 and n separately. Thus we cannot substitute u 0i as a function of a, b and n in the other equation to give conditions of a, b and n under which the system has exactly one, two or three equilibria. However, the above discussion of case (iii) suggests an effective elimination, i.e., using u n−1 0i = (b + 1)/(an − (n + 1)u 0i ) (obtained from the second equation of (22)), to replace the terms of u n−1 0i in the first equation of (22), which gives where Ξ is the quadratic function defined just below (16).
Knowing the discriminant ∆ of Ξ in (17), we discuss in the three cases: (C1) In case (C1), Ξ has no real zeros, implying that Ξ(u) > 0 for all u ∈ (0, a). It follows from (23) that Φ(u 0i ) < 0 for i = 1 and 2, which is exactly the situation (S1). This shows that if 0 < b < b * and a > a * (b) then Φ has a unique simple zero u * , which satisfies by the continuity of Φ thatũ 2 < u * < a, the same results as given in (i) and (ii) of the theorem. Since the case that b ≤ b * and 0 < a ≤ a * (b) is completely discussed in the paragraphs of cases (i) and (ii) below (13), we conclude that results (i) and (ii) of the theorem are proved except that b = b * and a > a * (b), which however is covered in the following case (C2).
We leave the proof of the lemma to the end of this section. Since the assumption a > a * (b) was required in the "other case" of case (iii), as shown at the beginning of the paragraph of (20), and the assumption b > b * was required in (C3) of the "other case", applying Lemma 2.2, we obtain the following conclusions: If a > a + (b), the zeros u 01 and u 02 of Φ distribute as in case (C3-1-1). By (24), we have Ξ(u 0i ) > 0, i = 1, 2, implying from (23) that Φ(u 0i ) < 0, i = 1, 2, i.e., the situation (S1) listed below (21). It implies by the statement just below the list of (S1)-(S5) that Φ has exactly one zero, denoted by u * , which lies in (ζ + , a) by the continuity of Φ. Therefore, in this case system (1) has a unique equilibrium, as stated in (iii-1) of the theorem.
As a consequence, we have completed our discussion in the 'other' case, i.e., a > a * (b), as mentioned in the paragraph of (20), and completed the proof of the theorem.
We end this section with the proof of Lemma 2.2.
3. Properties of equilibria. In section 3 we gave conditions of parameters for exact numbers of equilibria but did not obtain coordinates of those equilibria, which prevents us from determining qualitative properties of those equilibria by computing determinant and trace of the Jacobian matrix at each of them in the usual routine.
In what follows, we only use the interval distributions of those equilibria obtained in last section to give their qualitative properties.
Theorem 3.1. For the system (3) with n ≥ 2, equilibria obtained in Theorem 2.1 have the following properties: In case (i) of Theorem 2.1, i.e., 0 < b ≤ b 0 and a > 0, A * is a stable node. In case (ii-1) of Theorem 2.1, i.e., b 0 < b ≤ b * and either 0 < a < a 0 (b) or a > a 0 (b), A * is a stable node. In case (ii-2), i.e., b 0 < b ≤ b * and a = a 0 (b), A * is a stable node if b 0 < b < b * and A * is a degenerate stable node if b = b * . In case (iii-1) of Theorem 2.1, i.e., b > b * and either 0 < a < a − (b) or a > a + (b), A * is a stable node. In case (iii-2), i.e., b > b * and either a = a − (b) or a = a + (b), A * is a stable node and A 0 is a saddle-node. In case (iii-3), i.e., b > b * and a − (b) < a < a + (b), A 1 and A 3 are both stable nodes and A 2 is a saddle.
Proof. Let E : (ũ,ṽ) present a general equilibrium of system (3), whereṽ = b/(ũ n + 1) by (4). Although the coordinateũ is unknown yet, we can formally compute the Jacobian matrix of system (44) at E, i.e., where Ξ is the quadratic function defined just below (16). Clearly, since 0 <ũ < a as required in (6) and 0 <ũ < a/(b + 1) by (15). In order to determine the sign of D(ũ), we cannot compute its exact value as usual because we did not obtain precise coordinates of E in section 2. However, the monotonic interval whereũ locates gives enough information for the sign of D(ũ).
The above discussion indicates two degenerate situations: (ii-2D) b = b * and a = a * , and (iii-2) b > b * and a = a ± (b), which are included in subcases (ii-2) and (iii-2) respectively and in which D(ũ) = 0. In both situations exact one eigenvalue vanishes since T (ũ) < 0 by (32). In what follows, we give qualitative properties for the two degenerate equilibria.
In situation (iii-2), b > b * and a = a ± (b). Since the discussion for a = a − (b) is similar, we concentrate to properties of A 0 for a = a + (b). By Theorem 2.1 (iii-2), the coordinates of A 0 are given by (bn + b − 2n − )a + (b)/(2n(b + 1)), (bn − b + )/(2n) , where defined in (7) depends on (b, n). The eigenvalues are 0 and −bn − b + , which is negative. Expanding system (44) at A 0 , we get where a + presents a + (b) for simplicity. Using an invertible linear transformation , we normalize the linear part of (37) and obtain where we compute < 0 because b > b * = 4n/(n − 1) 2 as given in (7). It implies that A 0 is a saddle-node. This proves the results of Theorem 3.1 in case (iii-2), and therefore the proof of the theorem is completed.

4.
Bifurcations. In last section a transition is exhibited between a monostable state and a bistable state. In this section we show how those transitions arise. Theorem 3.1 indicates that system (44) has exactly two degenerate situations, i.e., (ii-2D) b = b * and a = a * , and (iii-2) b > b * and a = a ± (b), as given in the paragraph just above (34). The relation between (ii-2D) and (iii-2) can be seen from the limit lim b→b * a ± (b) = a * , which was shown in the proof of Theorem 2.1.
Near those situations, transition of states may occur and will be displayed by the following discussion of bifurcations.
Theorem 4.1. For n ≥ 2, a stable node and an unstable saddle arise from a saddlenode bifurcation of system (3) at A 0 as (a, b) varies across the curve , as shown in Figure 14. The curves Υ ± meet at point C : (a * , b * ) and divide the first quadrant of the (a, b)-plane into the region B := {(a, b) ∈ R 2 : b > b * and a − (b) < a < a + (b)}, called the bistability region, and the outside one M, called the monostability region. When (a, b) varies from (a * , b * ) to Υ ± or B, a codimension 2 bifurcation called the cusp bifurcation happens, producing an unstable saddle-node as (a, b) ∈ Υ + (or Υ − ) and a stable node and an unstable saddle as (a, b) ∈ B.
Theorem 5.1. System (3) has no closed orbits. Moreover, in case a = a ± (b) and b > b * , system (3) has an orbit linking equilibria A 0 and A * ; in case a − (b) < a < a + (b) and b > b * , system (3) has an orbit linking A 1 and A 2 and an orbit linking A 2 and A 3 ; in other cases, system (3) is globally asymptotically stable.
In order to prove the theorem, we need to investigate equilibria at infinity and closed orbits. For convenience, we use the time-scaling t = τ (v + 1)(u n + 1) to simplify system (3) as du dτ = a(u n + 1) − u(v + 1)(u n + 1), a polynomial system orbitally equivalent to system (3).
Lemma 5.2. For any positive integer n, system (44) has two equilibria I u and I v at infinity in the first quadrant, where I u lies on the u-axis and I v lies on the v-axis. I u is a saddle-node with a one-dimensional unstable manifold tangent to the u-axis at infinity. I v is a saddle-node with a one-dimensional unstable manifold tangent to the v-axis at infinity if n = 1. For n ≥ 2, I v is of fully null degeneracy and there is a unique orbit leaving from I v , which goes in the direction of the positive v-axis. No orbits connect to I u and I v in other directions.
Proof. Applying the Poincaré transformation u = 1/z and v = x/z to project all points at infinity but the two on the v-axis to the x-axis of the new (x, z)-coordinate, we change system (44) into the form dx d = −axz + bxz n + bz n+1 − axz n+1 := X(x, z), where d = dτ /z n , so that we only need to consider equilibria of (45) on the x-axis. From the equations X(x, 0) = 0 and Z(x, 0) = 0, we obtain a unique zero x = 0, showing that system (45) has a unique equilibrium on the x-axis. This equilibrium, denoted by O 1 , locates at the origin of (x, z)-plane, and therefore corresponds to an equilibrium of system (44) at infinity on the positive u-axis, denoted by I u . In order to determine if there is an equilibrium at infinity on the v-axis, denoted by I v , we apply another Poincaré transformation u = y/z and v = 1/z to change system (44) into the form dy d = ay n z − byz n + az n+1 − byz n+1 := Y (y, z), where d = dτ /z n , so that we only need to consider if system (46) has an equilibrium at the origin O 2 of the (y, z)-plane because O 2 corresponds to the point at infinity on the positive v-axis. One can check that this is true because Y (0, 0) = Z 2 (0, 0) = 0.  (45) and (47) Concerning I u , we note that system (45) at O 1 has eigenvalues 0 and 1. Using the invertible linear transformation ν = −x, ω = x + z, we rewrite (45) as Since Ψ(ν, ω) has a common factor ω, we easily see that ω = 0 is a center manifold. Restricted to this center manifold, system (47) is reduced to the equation dν/d = −aν 2 − aν n+2 , showing that O * 1 is a saddle-node with a saddle sector on the side of the positive ν-axis. Further, we claim that the branch of the unstable manifold on the side of the positive ω, denoted by W u + , lies in the second quadrant of (ν, ω)plane, as shown in Figure 13(a). In fact, we know that W u + is tangent to the ω-axis at O * 1 and, on the other hand, dν/d = Φ(0, ω) = −bω n+1 < 0 on the positive ω-axis, implying that the positive half ω-axis lies in the saddle sector. This proves the claim. Correspondingly, by the linear transformation as mentioned just before (47), on the (x, z)-plane system (45) has a center manifold z = −x and an unstable manifold tangent to the z-axis at O 1 . By the claim, the branch of the unstable manifold which corresponds to W u + lies in the first quadrant of (x, z)-plane (see Figure 13(b)). It follows that system (45) has at most one orbit near O 1 in the first quadrant. Actually, system (45) has exactly one orbit near O 1 in the first quadrant, which leaves from O 1 , because dx/d = bz n+1 > 0 by (45), restricted to the positive half z-axis. This gives the results of Lemma 5.2 on I u .
Concerning I v , we discuss equilibrium O 2 of system (46) in the two cases: n = 1 and n ≥ 2. The case n = 1 is simple, where system (46) has eigenvalues 0 and 1 and can be simplified by the linear transformation y = −ν 1 , z = ν 1 + ω 1 as Similarly to I u , we can reduce the system to the center manifold ω 1 = 0 and see that the origin O 2 is a saddle-node, which gives the result of Lemma 5.2 on I v for n = 1. For n ≥ 2, O 2 is of fully null degeneracy, i.e., the Jacobian at O 2 is a zero matrix. Using the Briot-Bouquet transformation (y, z) → (y, η 1 ), where y = y and z = η 1 y, to desingularize O 2 , we change system (46) into the form The common factor y in the first equation of (48) shows that the η 1 -axis is an orbit. Hence, by the uniqueness, no orbits intersect the η 1 -axis. This shows by the geometric sense of the Briot-Bouquet transformation that in the first quadrant of the (y, z)-plane system (48) has no orbits connecting to O 2 in a direction other than the positive z-axis. In order to determine if there is an orbit connecting to O 2 in the direction of the positive z-axis, we use another Briot-Bouquet transformation (z, y) → (z, η 2 ), where z = z and y = η 2 z, which rewrites system (46) as One can check that the origin O * 2 of the (z, η 2 )-plane is a saddle of system (49) with eigenvalues 1 and −1. Note that dη 2 /dt = az > 0 if we restrict (49) to the positive z-axis. This shows that a branch of the unstable manifold of the saddle lies in the first quadrant of the (z, η 2 )-plane and therefore, by the geometric sense of the Briot-Bouquet transformation, in the first quadrant of the (y, z)-plane system (46) has exactly one orbit near O 2 , which leaves from O 2 in the direction of the z-axis. This implies the results on I v for n ≥ 2. Hence, the proof of Lemma 5.2 is competed. Now, we are ready to prove the main theorem of this section.
In case a = a ± (b) and b > b * , the system has exactly two equilibria A 0 and A * by Theorem 2.1. Moreover, A * is a stable node and A 0 is a saddle-node by Theorem 3.1. Since Q(u, 0) = b/(u n + 1) > 0 on the positive u-axis and P (0, v) = a/(v + 1) on the positive v-axis by (4), all orbits initiating either the u-axis or the v-axis enter the first quadrant. Additionally, we see from (45), the system reduced by the Poincaré transformation in the proof of Lemma 5.2, that Z 1 (x, 0) = x > 0 on the positive half x-axis, which implies by Lemma 5.2 that all orbits in the first quadrant are eventually bounded. Thus, the branch of a center manifold at A 0 which lies in the saddle sector links with A * .
In case a − (b) < a < a + (b) and b > b * , by Theorems 2.1 and 3.1 the system has exactly three equilibria A 1 , A 2 and A 3 , the middle one A 2 of which is a saddle and the other two A 1 and A 3 are both stable nodes. It follows that the two branches of the unstable manifold at A 2 link with A 1 and A 3 separately for the same reason as in the last paragraph.
In remaining cases, system (3) has exactly one equilibrium A * , which is a stable node, showing that system (3) is globally asymptotically stable. Thus, the proof of Theorem 5.1 is completed.
Summarizing Theorems 2.1-5.1 and some assertions given in their proofs, we see that system (3) has exactly one equilibrium as (a, b) lies in M ∪ C, exactly two equilibria as (a, b) lies on Υ ± , exactly three equilibria as (a, b) lies in B, as shown in Figure 14. Conclusions. Clearly, saddle-node bifurcation occurs on the curves Υ + and Υ − and cusp bifurcation occurs at the point C in Figure 14. These bifurcations show that the genetic toggle switch changes from monostability to bistability as (a, b) enters B (the angular sector between Υ + and Υ − ) from M (outside B).
For an intuitive observation, we plot the broader curves Υ ± for n = 5 in Figure 15. They actually correspond in the (log a, log b)-coordinate to the lines log a = β log b and log a = log b/γ respectively, which were given in [10] as shown in Figure 2 and can be plotted precisely in the case n = 5 as the finer curvesΥ + : a = b and Υ − : a = b 1/5 in Figure 15. As indicated in [10] and our Introduction, the curveΥ + (resp.Υ − ) being the critical conditions of transition between monostability and bistability were obtained approximatively in the special case that the degenerate equilibrium E : (u, v) lies closely near the v-axis, i.e., u v (resp. the u-axis, i.e., v u) under the assumption that both a and b are large. In our paper we obtain the curves Υ ± accurately without their assumption and prove that they are both saddle-node bifurcation curves. Therefore, our results correct the boundaryΥ of bistable region. Choosing (b, a) = (3, 2.8), which lies betweenΥ + and Υ + , and (b, a) = (3, 1.8), which lies betweenΥ − and Υ − , we make numerical simulations with the software PPLANE, showing that the system is monostable rather than bistable in Figures 16 and 17. This demonstrates that the region between Υ + and Υ − in Figure 15 is true but was not found in [10].
In comparison with the work [2] by the iGEM Duke Team in 2013, which plots intersection between the horizontal isocline and the vertical isocline numerically to shows that system (1) has exactly three equilibria for the choice (a, b, β, γ) = (4, 4, 3, 3) and a unique equilibrium for the choice (a, b, β, γ) = (3, 3, 1, 2), our results do not specify any choice but generally give conditions for exact number of equilibria analytically.
The genetic toggle switch is a synthetically constructed circuit in Escherichia coli to control gene expression. Although some dynamical properties such as multistability and oscillations have been found in the natural genetic circuit, researchers want to build genetic circuits artificially for desired properties. One of desired properties in a gene-regulatory network is the switching behavior between two equilibria. Our paper presents global dynamics of system (3) and gives a condition for bistability. In the bistable state (seen in Figure 15) , the stable manifold of the saddle A 2 divides the first quadrant into two basins of attraction. On the other hand, the unstable manifold of the saddle A 2 , being heteroclinic connections with A 1 and A 3 , makes a separatrix of the two basins so that the upper (lower) half-basin above (below) the separatrix preserves evolution inside the half-basin. The two basins make a mechanism for such a switching and the bifurcation values given in our Theorem 4.1 provide accurate parameter choices which lead to bistability and regulate the capacity to achieve switching behavior.
In this paper we only consider system S(1, n) for convenient computation. Because of the symmetry of (1), one can discuss system S(m, 1) similarly. Our discussion on system (3) exhibits a routine for the general system S(m, n) or system (1) but the computation will be more complicated.