Stability analysis of an autocatalytic protein model

Life is maintained by complex interactions of genes and proteins.¹ The first step toward understanding such a system is to understand its steady-state behavior. Depending on the network structure, the gene network may exhibit multistability, that is, there may be more than one steady state.^2–16 Multistability forms the basis for explaining the response of the cell to the environment, and it is also important for understanding stem cell differentiation into different tissues.^17,18 The existence of a positive feedback loop in the regulatory network is necessary for multistability.^12–15 

The genetic regulatory system can be modeled in terms of the probability distribution of molecule numbers, which follows a master equation.^{2–5,7–9,19–31} At the limit where the number of molecules is large and the fluctuation of molecule numbers becomes negligible, the system can be described in terms of the molecule density, which follows a system of first-order differential equations.^{2,10–16,32,32–34} This set of differential equations describes the time-dependent production and decay of various RNA and protein molecules, as well as binding of transcription factors to DNA that stimulate or inhibit gene transcription. In this case, a general method for detecting multistability for a large class of systems with positive feedback loops has been proposed.¹⁰ However, it is of interest to perform an explicit analysis of such a network, when possible, to obtain more detailed information on the behavior of the system.

The simplest biological system with a positive feedback loop is a self-regulatory genetic circuit. In a self-regulatory genetic circuit, a protein acts as a transcription factor for its own gene, thus positively regulating its own production. According to central dogma, there should be at least three components to this self-regulatory circuit: DNA, RNA, and the protein itself. Even for the simplest system consisting of these three components, obtaining an exact solution to the chemical master equation is a nontrivial task; thus, various approximations have been employed. The intermediate step of RNA production is often omitted such that protein production from the gene is approximated as a one-step process.^2,3 A fast-equilibrium DNA state, in which the binding and unbinding of the protein to the DNA is assumed to be more rapid than all other processes, also is frequently employed as an approximation.^2,11,35–37 In this case, the DNA component can be eliminated from the equations, which is a special case of a quasi-steady-state approximation in which fast degrees of freedom are removed from the equations.^19–31 

Here, we retain the assumption that the system can be described by a set of first-order differential equations, but we explicitly include both the RNA and DNA components in the system of equations. In fact, the dissociation rate of transcription factors from DNA is on the order of 10⁻³ ∼ 10 s⁻¹,^38–42 whereas mRNA transcriptions are initiated at rates of approximately 10⁻³ ∼ 1 per second^43,44 and proteins are translated at rates on the order of 10 amino acids per second from one ribosome.^45–48 If we assume that the number of ribosomes attached to an mRNA molecule is on the order of 1 ∼ 10 at any given time,⁴⁷ and if we also assume a protein size of 10² ∼ 10³ amino acids, then we can estimate the overall translation rate for one mRNA molecule as 10⁻² ∼ 1 molecules per second. These hand-waving arguments show that it is likely that there are examples in which the dissociation rate of a transcription factor from DNA is not much larger than the transcription or translation rate. In addition, because the transcription and translation rates are of a similar order, it is more accurate to describe protein production as a two-step rather than a one-step process.

However, we still neglect the time scales for the oligomerization of the protein molecules, transport of mRNA from the gene to a ribosome, and diffusion of the protein back to the gene.^1,2,8,9,16 In fact, the oligomerization of protein molecules are characterized by second-order rate constants in the range of 0.5 × 10⁶ ∼ 5 × 10⁶ M⁻¹s⁻¹,⁴⁹ and consequently the oligomerization rate is much higher than the frequency of transcriptional and translational events if the protein concentration is much larger than 2 × 10⁻⁹ ∼ 2 × 10⁻⁶ M. The diffusions of molecules in the nucleoplasm and cytoplasm exhibit almost the same characteristics, and the diffusion coefficients are measured to be 0.01 ∼ 1 μm²s⁻¹ and 1 ∼ 100 μm²s⁻¹, for mRNA and proteins, respectively.^50–57 Therefore, although the diffusion of mRNA is the rate-limiting step, the corresponding rate is much higher than the frequencies of transcription and translation events if the distance between the gene and region of translation is significantly less than 0.1 ∼ 10 μm. This condition is easily realized for prokaryotes, in which the nucleoid and ribosomes are colocalized,⁵⁷ as well as for an eukaryotic cell if the endoplasmic reticulum is located close enough to the nucleus.

In the case of an eukaryote, extra time (approximately 180 ms) is required for mRNA to be transported through a nucleic pore. The corresponding rate is approximately 5 s⁻¹, which is not significantly higher than the frequency of transcription and translation events.⁵⁸ This extra time for mRNA transport should be treated as part of the transcriptional or translational time delay. We also neglect the time delays caused by the transcription and translation processes.^1,2,8,9,16 During transcription, mRNA is elongated at a rate of 10 ∼ 10² nucleotides per second,^48,59–68 corresponding to 10⁻³ ∼ 10⁻¹ molecules per second for a protein size on the order of 10² ∼ 10³ amino acids, even if we consider the transcription of the coding region only. As mentioned above, the protein is synthesized at a rate on the order of 10 amino acids per second from one ribosome,^45–48 corresponding to 10⁻² ∼ 10⁻¹ molecules per second for a protein size on the order of 10² ∼ 10³ amino acids, which suggests that the transcriptional and translational time delay may not be negligible in real situations. We will discuss the implications of such time delays later in this paper.

In the next sections, we introduce our model and analyze the stability of the fixed points. We then compare our results with previous results based on a fast-equilibration approximation,^11,35–37 and discuss similarities and differences.

The simplest system with a positive feedback loop is a self-regulatory genetic circuit in which a protein X binds to its gene to act as a transcription factor for its own species. Let us denote the inactive and active genes as D and D^∗, respectively. We also assume that m molecules of proteins must bind to the inactive gene to activate the gene and initiate transcription. The active gene D^∗ produces the RNA molecule R, which in turn produces molecules of protein X via translation. Both the RNA and protein degrade over time. These processes can be described by the following reactions:^{2,3,11,35–37}

The rate equations of these processes are given as:

The condition that the DNA concentration is a constant,

was imposed to eliminate [D] from the first and third equations. If one considers an average behavior of a large population of cells, then [D] and [D^∗] are proportional to the average number of genes in the off and on states, which also are proportional to the probabilities that a gene is in the off and on state.

Eq. (2) can be simplified by rescaling time as

and defining the rescaled variables as

and the parameters

where μ, ϵ, and τ are chosen to satisfy

The resulting set of equations

is the object of analysis in this work.

The fixed points of the system (8) are obtained by solving $x ̇ = y ̇ = z ̇ =0$⁠. Note that x = y = z = 0 is a fixed point of Eq. (8). The condition for the existence of other fixed points, for m = 1 and m = 2, have been published previously.^35,36 Such conditions for general values of m ≥ 1 are provided in  Appendix A. As shown in  Appendix A, if

then there are two additional fixed points for m > 1 whose x values are given by the roots of the equations

For m = 1, the condition Eq. (9) should be interpreted as

for which there is one additional fixed point whose x value is given by the root of Eq. (10). The values of y and z are then obtained by solving $y ̇ = z ̇ =0$⁠. The stability of each fixed point is analyzed by examining the signs of the three eigenvalues of the linearized set of equations, which are obtained by expanding Eq. (8) around the given fixed point.

Note that one can first solve $z ̇ =0$ in Eq. (8) to obtain the reduced set of equations

which shares the same set of fixed points as the original set in (Eq. (8)). The two-variable system described by Eq. (12) has been investigated in detail in the literature.^35–37 Such a dimensional reduction is based on the assumption that $z ̇ =0$ is always satisfied, which corresponds to the case in which the DNA equilibrates over a time scale much shorter than that of RNA or protein. In fact, according to Tikhonov’s theorem, the solutions of Eq. (8) converge to the solutions of the dimensionally reduced equations in Eq. (12) in the limit γ → ∞ as long as z, which is given by the last equation in the set of Eq. (12), is a stable fixed point of

where x is treated as a constant and the initial value z(0) lies in the domain of attraction of Eq. (13) with x = x(0).^69,70 In fact, because Eq. (13) is linear, z^∗ = x^m/(1 + x^m) is a unique fixed point. By expanding Eq. (13) around z^∗, with z = z^∗ + δz, we obtain

Because the coefficient of δz on the right-hand side of the equation is negative, we observe that z^∗ is a stable fixed point. Therefore, z^∗ is the global stable fixed point of Eq. (13) for any given x. Hence, Tikhonov’s theorem can be invoked, and the full set of equations in Eq. (8) reduces to that in Eq. (12) in the limit of γ → ∞. The consequence of Tikhonov’s theorem applied to the system in Eq. (8) is shown in Fig. 1, where the numerical solution of Eq. (8) for m = 1, α = β = 1/4, and ν = 0 is displayed in the x − z plane for various values of γ and compared with the limiting solution at γ → ∞ given by Tikhonov’s theorem.

If we analyze the flow near the fixed point, then we can conclude that two of the three eigenvalues of the linearized equations around a fixed point (x^∗, y^∗, z^∗) will approach those obtained from Eq. (12):^35–37 

The remaining eigenvalue will approach that obtained from Eq. (13), which governs the fast dynamics of z:

where we substituted x = x^∗ in Eq. (13), because we are considering an infinitesimal perturbation around x^∗.

The analysis of the reduced system Eq. (12) is reviewed in  Appendix B. For m = 1, when there is an additional positive fixed point, this positive fixed point is a stable fixed point and x = 0 is an unstable fixed point. Otherwise, x = 0 is a stable fixed point. However, for m > 1, x = 0 is always a stable fixed point, and if there are two additional positive fixed points x_± with x₋ < x₊, then x₋ is a saddle point and x₊ is a stable fixed point, leading to bistability.

We have observed that the original set of equations given by Eq. (8) is reduced to the set in Eq. (12) in the limit of γ → ∞; the stability of this set of equations has been analyzed previously.^35–37 However, to examine the stability of the fixed points of Eq. (8) for small values of γ, one has to perform the exact analysis of the linearized equations. Such equations are obtained by writing x = x_∗ + δx, y = y_∗ + δy, and z = z_∗ + δz:

The stability of this set of equations is determined by examining the real parts of the eigenvalues of the matrix A₃, which are obtained by solving the characteristic equation:

where

The case of m > 1 and x_∗ = y_∗ = z_∗ = 0 is the easiest to analyze. In this case, $β ̄ =β$⁠, $γ ̄ =γ$⁠, $ν ̄ =ν$⁠, and ζ = 0, and the solutions to the characteristic equation Eq. (18) are

all of which are negative. Therefore, for m > 1, the origin is always a stable fixed point. The other cases remain to be analyzed.

In general, a fixed point is stable if the real parts of the three eigenvalues of Eq. (18) are negative, whereas a fixed point is a saddle point if some of the eigenvalues have positive real parts. Denoting the three eigenvalues as λ_i(i = 1, 2, 3), we observe that

The first equation in the set of Eq. (21) implies that at least one eigenvalue has a negative real part, which is also a necessary condition for a fixed point to be stable.

If $∏ i λ i >0$⁠, the realnesses and signs of the eigenvalues are either (−, − , + ) or (C, C^∗, + ), where + and − denote positive and negative real roots, respectively, and C, C∗ denotes a conjugate pair of complex roots. In addition, because at least one eigenvalue must have a negative real part due to the first equation in the set of Eq. (21), the real parts of the complex conjugate roots are negative.

However, if $∏ i λ i <0$⁠, the realnesss and signs are either (−, − , − ) or (−, C, C^∗). Note that (−, + , + ) is excluded due to the second equation in the set of Eq (21). If λ₃ < 0 and λ₁, λ₂ > 0, then

where the last inequality follows from the first equation in the set of Eq. (21), which is

Given that

we obtain from Eq. (22)

which is an obvious contradiction. Therefore, the signs of the eigenvalues cannot be (−, + , + ). The only remaining task is to determine the sign of the real parts of the complex roots for the case $∏ i λ i <0$⁠, which is achieved by utilizing the analytic form for the roots of the cubic equation in Eq. (18).⁷¹ The details are provided in  Appendix C, where it is shown that even if there is a conjugate pair of complex roots, $∏ i λ i <0$ implies that their real parts are all negative.

In summary, we obtain a saddle point for $∏ i λ i >0$ in which the three roots of the Hessian are of the form (−p, − q, r) or (−p + iq, − p − iq, r) with p, q, r > 0. In contrast, we obtain a stable fixed point for $∏ i λ i <0$ in which the roots are of the form (−p, − q, − r) or (−p + iq, − p − iq, − r) with p, q, r > 0. $∏ i λ i =0$ corresponds to the marginal case in which the roots are of the form (−p, − q, 0) or (−p + iq, − p − iq, 0) with p, q > 0. We did not obtain pure imaginary eigenvalues in our analysis; therefore, there is no Hopf bifurcation in this model.

Now we analyze the stability of each fixed point by determining the sign of $∏ i λ i$⁠. For m = 1 and x_∗ = 0, we have

Therefore, when αβ > 1, such that the origin is the unique fixed point, $∏ i λ i <0$ such that the origin is a stable fixed point. In contrast, when αβ < 1, such that there is an additional nonzero fixed point, $∏ i λ i >0$ and the origin becomes a saddle point.

For additional non-zero fixed points,

where Eq. (10) was used in the transformation from the first equation to the second equation. For m = 1, the expression in Eq. (10) is negative, such that the additional positive fixed point is stable whenever it exists. For m > 1, it remains to determine the sign of

As shown in  Appendix D, when there are two additional positive fixed points x_± with x₋ < x₊, then g(x₋) > 0 and g(x₊) < 0. That is, x₋ is a saddle point and x₊ is a stable fixed point.

Thus, we rigorously identified the stability of the fixed points in the three-dimensional system of Eq. (8), which is summarized below:

i) m = 1

(x, y, z) = (0, 0, 0) is a unique stable fixed point when αβ < 1. For αβ > 1, (x, y, z) = (0, 0, 0) becomes an unstable fixed point, and there is an additional nonzero stable fixed point whose x-value is obtained as a root of Eq. (10). αβ = 1 is the marginal case in which these two fixed points merge. Given that there is a unique stable fixed point regardless of the parameters, there is no bistability in the system for m = 1. In fact, extending x, y, and z to the nonphysical negative region, we can observe that (x, y, z) = (0, 0, 0) is a saddle point for αβ > 1, where the stable manifold exists in the nonphysical region. At αβ = 1, the saddle point merges with the stable fixed point, and the additional saddle reappears in the nonphysical region for αβ < 1, whose x-value is obtained as a root of Eq. (10). Therefore, we observe transcritical bifurcation for m = 1 if we extend (x, y, z) to negative values.

ii) m > 1

(x, y, z) = (0, 0, 0) is always a stable fixed point. There are additional nonzero fixed points whose x values are x_± with x₋ < x₊ that correspond to the roots of Eq. (10) for

where x₋ and x₊ are a saddle point and a stable fixed point, respectively. Therefore, (x, y, z) = (0, 0, 0) is the unique stable fixed point for αβ > (m − 1)^1−1/m/m, whereas we observe bistability for αβ < (m − 1)^1−1/m/m. Again, αβ = (m − 1)^1−1/m/m corresponds to the marginal case in which the stable fixed point x₊ and the saddle point x₋ merge into one fixed point. These fixed points disappear for αβ > (m − 1)^1−1/m/m. Therefore, for m > 1 we observe a saddle-node bifurcation.

The stability of the fixed points obtained in the previous section is consistent with that obtained from the analysis of the reduced set of equations in Eq. (12), as reviewed in  Appendix B, which is not surprising for large values of γ. As mentioned earlier, we expect that two of the three eigenvalues obtained by linearizing the set of equations in Eq. (8) around a given fixed point will approach $λ ± ( 2 d )$ in Eq. (15) and that the remaining eigenvalue approaches $λ 0 ( fast )$ in Eq. (16) as γ → ∞, due to Tikhonov’s theorem. Because $λ 0 ( fast ) =−γ ( 1 + x ∗ m )$ is negative, the stability condition will be the same as those obtained from the reduced set of equations if γ is large enough. In fact, we can explicitly check from the explicit expressions for the eigenvalues (See Eq. (C2) in  Appendix C) that

for large values of γ. The nontrivial result is that the stability condition remains the same even for small values of γ. In fact, the condition of bistability does not involve the parameter γ. However, the detailed behavior of the fixed points does depend on γ. Although the limiting values $λ ± ( 2 d )$ and $λ 0 ( fast )$ are all real, there exists the possibility that some of the eigenvalues are complex for a finite value of γ. As explained in the previous section, complex eigenvalues appear if R² − 4Q³ > 0. Indeed, there exists a region of small γ where complex eigenvalues appear, which leads to an oscillatory convergence to a fixed point that is not observed in the reduced set of equations Eq. (12).

For example, if m = 1 with α = β = 1/8, γ = 1/2⁹, and ν = 0, then there exists a nonzero stable fixed point at x_∗ = 63, y_∗ = 63/8, and z_∗ = 63/64, with

and the corresponding eigenvectors

By taking the linear combination with complex coefficients such that the values of x, y, and z are real, we can obtain the asymptotic behavior near the fixed points:

with constants A, B, and ϕ that depend on the initial condition. The expression Eq. (33) does indeed exhibit oscillatory convergence toward the fixed point. As another example, the time-series of (x, y, z) near the positive stable fixed point for m = 2, α = β = 0.70, γ = 0.28, and ν = 0 is shown in Fig. 2, which exhibits oscillatory convergence.

The boundary between the region of oscillatory convergence and non-oscillatory convergence around the stable positive fixed point is shown in Fig. 3 for m = 2, ν = 0, and α = β. It was not possible to scan all possible regions of the parameter space, but the occurrence of oscillatory convergence at small values of γ seems to be a generic feature. A simple intuitive argument as to why such an oscillatory convergence occurs is unavailable at the moment. It has been shown that, in the single-cell limit, damped oscillations become sustained stochastic oscillations whose observability depends on the strength of the macroscopic damping.⁷² 

In this work, we analyzed the stability of the fixed points of the simplest auto-regulatory system, consisting of DNA, RNA, and the protein product that positively regulates its own production upon binding to the gene. After linearizing the equation around the fixed points, the eigenvalues of the three-dimensional Hessian were obtained using the root formula of the cubic equation. The result rigorously confirms the well-known results of the two-dimensional analysis as follows: when one molecule of the transcription factor is sufficient to activate the gene, there is no bistability, whereas for the case in which binding of more than one molecule of the transcription factor is required for activation, there is bistability for a certain range of parameters.

We found that the condition for bistability remains identical to that obtained from the dimensionally reduced system, even when the DNA binding and unbinding are slow, where Tikhonov’s theorem cannot be invoked. We found, however, that the detailed behavior of the system near the fixed points is dependent on the rates of DNA binding and unbinding. In particular, oscillatory convergence may occur when DNA binding and unbinding are slow, depending on the values of the other parameters, which is a behavior that is not observed in the limit of fast equilibration. However, the biological significance of such an oscillation remains unclear at the present time.

The current model retains some of the approximations that have been used previously, such as the instant oligomerization of the protein molecules,^1,2,8,9,16 neglect of diffusion times of various molecules and transcriptional and translational time delays,^1,2,8,9,16 and neglect of stochastic fluctuations.^{2,10–16,32–34} However, the current model also omits an oft-used approximation, the fast equilibration of protein and DNA. Therefore, the current model can be considered a relatively accurate description of the autocatalytic protein system in a parameter regime that cannot be covered by previous approaches. As discussed in the introduction, there are ranges of biological parameters over which most of the assumptions employed in the current model are valid. The transcriptional and translational time delays require a more detailed discussion here. The overall rates of transcription and translation, denoted as a and c, respectively, in Eq. (2), determine the characteristic time scale of the current model. The rate a (c) is the frequency of initiation of transcription (translation) events, which depends on the concentration of RNA polymerase (ribosome) molecules. These concentrations were implicitly assumed to be constants in the current model. The time scale t_bind ≡ a⁻¹ (≡c⁻¹) corresponds to the average time between the binding of the RNA polymerase (ribosome) molecule to the active DNA (mRNA) molecule. Let us denote the time spent on transcription (translation) as t_delay. The time delay can be neglected only if t_bind > t_delay. By denoting the number of RNA polymerase (ribosome) molecules bound to an active DNA (mRNA) averaged over time as n, we can observe that the condition for the neglect of time delay can be written as n ≪ 1. As discussed in the introduction, this condition is usually not satisfied in real situations. For example, during active translation, many ribosomes are almost always bound to mRNA,⁴⁷ implying n > 1. Similar arguments can be made for transcription. In particular, eukaryotic genes are large, due to the presence of introns; thus, the time spent on transcription is longer than would be estimated from the coding region only and can easily reach tens of minutes. In addition, the time spent for splicing adds an extra 1 ∼ 10 minutes to the transcriptional time delay.^73–75 These estimates suggest that the transcriptional and translational time delays, which have been neglected in the current work and most previous work, may play an important role.

Let us conclude with a brief discussion on how the results may be affected by the inclusion of time delays. To incorporate time delays, the equations Eq. (8) should be modified to delay-differential equations as follows:

where τ₁ and τ₂ denote the transcriptional and translational time delays, respectively. It is obvious that the time delay does not play a role in stationary solutions due to the lack of time dependence. Therefore, the result regarding the number of fixed points for a given parameter set does not change. However, the result regarding the stability of the fixed points may change if we include a time delay, because we consider the dynamics of the system after introducing a small perturbation around a fixed point. Note that only saddle-node and transcritical bifurcations were found for the model described by the set of equations in Eq. (8). Although oscillatory convergence exists for some region of the parameter set, the Hopf bifurcation is ruled out because pure imaginary eigenvalues do not appear around any fixed point. However, the Hopf bifurcation is found quite often in systems of delay-differential equations.^76–79 Therefore, we may expect that some of the bifurcations found in the current model may be changed to the Hopf bifurcation once Eq. (8) is changed to Eq. (34), the detailed analysis of which is left for future studies.

This study was supported by Grant No. HI12C0675 of the Korea Healthcare Technology R&D Project through the Korea Health Industry Development Institute (KHIDI), funded by the Ministry of Health and Welfare, Republic of Korea.

As explained in the main text, x = 0 is always the fixed point of the set of equations in Eq. (8). Additional non-negative fixed points are obtained by solving Eq. (10), which we rewrite here as

The analysis for m = 1 and m = 2, where the analytic solution to Eq. (A1) can be obtained easily, has been reported previously.^35–37 The new content in this appendix is the analysis for general values of m ≥ 1.

When m = 1, it is obvious that there is one positive root to the equation (A1) if and only if

Note that αβ = 1 is the marginal case with root at x = 0.

To investigate the existence of real roots of Eq. (A1) for m > 1, we first locate the extremum of f(x) for x ≥ 0 by taking its derivative and setting it to zero:

which yields one extremum at

and the other at x = 0 for m > 2. Given that f^(j)(0) = 0 for j < m − 1 and f^(m−1)(0) = − (m − 1)!/αβ < 0, where f⁽ⁿ⁾ denotes the n-th derivative of f(x), we can observe that x = 0 for m > 2 is a local maximum. In contrast, given that

we can observe that x₁ is the global minimum of f(x). Given that f(0) = 1 > 0 and $lim x → ∞ f ( x ) =∞>0$⁠, Eq. (A1) will have two distinct positive roots if and only if the value of f(x) is negative at the global minimum x₁:

which can be rewritten as

In contrast, there will be no real roots if the inequality is inverted. When the inequality is replaced by the equality, we obtain the marginal case in which the two roots coincide.

Let us consider the stability of the two-dimensional system Eq. (12) around the fixed points. Most of this analysis can be found in the previous literature,^35–37 but an explicit form of the eigenvalues has been provided to compare with those obtained from the three-dimensional system in the limit of ultra-fast binding and unbinding of the protein to the gene. In addition, the analysis is performed for general integer values of m instead of just m = 1 and m = 2. By considering a small deviation of x and y around a fixed point (x_∗, y_∗) of Eq. (12), with x = x_∗ + δx and y = y_∗ + δy, we obtain a linearized equation

Solving the characteristic equation for the matrix A,

we obtain the eigenvalue

Given that the expression inside the square root is non-negative, we can observe that the eigenvalues are real.

From the first line of Eq. (B3), we can observe that if

then the absolute value of the square root term is smaller than the first term, which is negative; therefore, the signs of both eigenvalues are negative, and we obtain a stable fixed point. However, if

then the signs of the eigenvalues λ_± are positive and negative, respectively, and we obtain a saddle point. It is straightforward to observe that the equality corresponds to the marginal case in which one of the eigenvalues is zero. This general result can be used to analyze the stability of specific fixed points.

As discussed in  Appendix A, x_∗ = 0 is a unique fixed point in this case. The condition αβ > 1 is the same as that in Eq. (B4) with x_∗ = 0; thus, this fixed point is stable.

It is clear that x_∗ = 0 is now a saddle point, because Eq. (B5) now holds for x_∗ = 0. There is an additional fixed point that is a solution to Eq. (A1):

We can observe that

Therefore, the condition in Eq. (B4) is satisfied, and we obtain a stable fixed point.

To summarize, for m = 1, when the RNA and protein is degraded fast enough such that αβ > 1, x_∗ = 0 is the unique stable fixed point of Eq. (12). Eventually, the concentrations of both species go to zero. In contrast, when the degradation is slow enough such that αβ < 1, then x_∗ = 0 is a saddle point, and an additional stable nonzero fixed point $x ∗ = 1 α β −1$ appears. We see that there is no bistability for m = 1.

In this case, x_∗ = 0 is always a stable fixed point because

As discussed in  Appendix A, when Eq. (A7) holds, we have additional nonzero fixed points that are the two distinct real roots x₋ < x₊ of Eq. (A1). After multiplying $αβ ( 1 + x ∗ m ) 2$ to both sides of Eq. (B4) and Eq. (B5), we can observe that either a stable fixed point or saddle point can be obtained, depending on whether the sign of the expression

is positive or negative. Note that Eq. (A1) was used in the transformation from the first expression to the second expression. Therefore, we only have to examine the sign of

Let us define

We can observe that this function is zero when

Note that the positive fixed points x_± are obtained as the roots of the equation

That is,

By evaluating f(x) at x = x₀, we can obtain

where the last inequality follows from the condition in Eq. (A7). Therefore, we can observe that

Given that g(x) is a monotonically decreasing function and g(x₀) = 0, it immediately follows that g(x₋) > 0 and g(x₊) < 0. Therefore, x₋ is a saddle point and x₊ is a stable fixed point. To summarize, when the RNA and protein are degraded fast enough that

then x = 0 is a unique sable fixed point. When the degradation is slow enough such that the inequality is inverted, then we have additional nonzero fixed points x₋ < x₊, where x₋ is a saddle point and x₊ is a stable fixed point. Given that x = 0 and x = x₊ are both stable fixed points, we obtain bistability.

In this appendix, we prove that even if two roots of the characteristic equation in Eq. (18) are complex, the real parts of the complex roots are negative as long as $∏ i λ i <0$ and the other conditions given in Eq. (21) hold true. We use the fact that the three roots of a cubic equation

are given by⁷¹ 

where

When the coefficients p₀, p₁, and p₂ are all real, we obtain either three real roots or one real root and a conjugate pair of complex roots, depending on the sign of R² − 4Q³. When R² − 4Q³ ≤ 0, its square root is either imaginary or zero, and consequently the second and third terms of each line of Eq. (C2) form a conjugate pair such that the three roots are real. In contrast, if R² − 4Q³ > 0, then the second and third terms of λ_± are linear combinations of distinct real quantities with complex coefficients e^±2πi/3; thus, λ_± form a conjugate pair of complex roots, which is the case of interest here.

For the characteristic equation in Eq. (18), we obtain

Consequently, $R 2 − 4 Q 3 <|R|$ when R² − 4Q³ > 0, such that the sign of $R ± R 2 − 4 Q 3 2 3$ is the same as that of R. After some tedious algebra, it can also be shown that

Given that the coefficient of R on the right-hand side of the equation is positive and the remaining term is negative, we can observe that if R² − 4Q³ > 0 such that there are complex roots, then R > 0. Given that the real part of the complex roots is:

which is the sum of three negative numbers, the real part of the complex roots is negative. Therefore, we conclude that the fixed point is stable under the condition in Eq. (21) and $∏ i λ i <0$⁠, even if some of the eigenvalues are complex.

As explained in the main text, the stability of a positive fixed point x^∗ for m > 1 is determined by the sign of

evaluated at x = x^∗. Note that g(x) vanishes at

Also note that the positive fixed points x_± are obtained as the roots of the equation

That is,

By evaluating f(x) at x = x₀, we can obtain

We can confirm that the condition for the existence of two distinct positive fixed points x_±, in Eq. (29), can be written as

Therefore, whenever two distinct positive fixed points of Eq. (8) exist, then Eqs. (D4) and (D6) hold. Because f(0) > 0 and $lim x → ∞ f ( x ) >0$⁠, we can observe that x₋ < x₀ < x₊. Given that g(x) is a monotonically decreasing function of x, and g(x₀) = 0, we can observe that g(x₋) > 0 and g(x₊) < 0. Therefore, x₋ is a saddle point and x₊ is a stable fixed point.