## Analytic stability analysis of three-component self-regulatory genetic circuit

## Julian Lee ∗

Department of Bioinformatics and Life Science, Soongsil University, Seoul, Korea (Dated: August 12, 2014)

## Abstract

A self-regulatory genetic circuit, where a protein acts as a positive regulator of its own production, is known to be a simplest form of biological network with a positive feedback loop. Although at least three components, DNA, RNA, and the protein, are required to form such a circuit, the stability analysis of fixed points of the self-regulatory circuit has been performed only after reducing the system into a two-component system consisting of RNA and protein only, assuming a fast equilibration of the DNA component. Here, the stability of fixed points of the three-component positive feedback loop is analyzed by obtaining eigenvalues of full three dimensional Hessian matrix. In addition to rigorously identifying the stable fixed points and the saddle points, detailed information can be obtained, such as the number of positive eigenvalues near a saddle point. In particular, complex eigenvalues is shown to exist for sufficiently slow binding and unbinding of the auto-regulatory transcription factor to DNA, leading to oscillatory convergence to the steady states, a novel feature unseen in the two-dimensional analysis.

PACS numbers: 87.10.Ed, 87.16.ad, 87.16.Yc, 87.18.Cf

∗ jul@ssu.ac.kr

## I. INTRODUCTION

Life is maintained by complex interaction of genes and proteins. The genes are turned on and off according to environmental change, so the gene network can be considered as a Turing machine with input and output. The genetic regulatory system is often modelled as a system of first-order differential equations that 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 the transcription of genes[1-12]. The first step toward understanding such a system is to understand its steady state behavior. Depending on the network structure, the gene network exhibits multistability, meaning there can be more than one steady state of the protein concentrations[4-6, 9-13]. Multistability not only forms the basis for explaining the response of the cell to the environment, but can also be used to model epigenetics, where information other than genetic one are passed on to the next generation of a cell[6]. It has been shown that the existence of a positive feedback loop in the regulatory network is a necessary condition for a multistability[9-12]. Although general method for detecting multistability for a large class of systems with positive feedback loop has been proposed[5], it is of interest to perform an explicit analysis of a network if one can, where one can obtain more detailed information on the behavior of the system near the fixed point, such as the numbers of real and complex eigenvalues, and those with positive real parts in the case of a saddle point.

The simplest biological system with a positive feedback loop is a self-regulatory genetic circuit where a protein acts as a transcription factor for its own gene, thus positively regulating its own production. According to the central dogma, there should be at least three components, DNA, RNA, and the protein. However, fixed-point analysis for such a simple self-regulatory genetic circuit has been performed only for two-dimensional system[19, 20, 22] or one-dimensional system[6] which are obtained after dimensional reductions of the original set of equations, and it has been concluded that if the binding of the transcription factor to the DNA is cooperative, in the sense that more than two protein molecules are required to be bound in order to stimulate the production of their species, then there can be more than one possible steady-state concentrations of these proteins.

Although fixed-points obtained from a dimensionally reduced system of equations are quite valid, the stability analysis of the fixed point in the dimensionally reduced system

of equations is not justified: they are obtained under the assumption that some of the components, usually taken to be the DNA, are equilibrated to their steady-state values. The assumption of the fast-equilibration itself has to be tested, by showing that the eliminated variables are stable with respect to their steady-state values so that such a limit can indeed be taken.

In this work, I consider the three-component auto-regulatory system consisting of DNA, RNA, and obtain analytic form for the eigenvalues of the Hessian at the fixed points, using the cubic root formula. Not only I show that this system can indeed exhibit bistability, but I can also explicitly investigate whether and how many of the eigenvalues have positive real values at the saddle point. I also find that the eigenvalues at the stable fixed point can have imaginary parts, leading to oscillatory convergence of the concentrations to the steady state, a novel feature that cannot be seen in the two-dimensional analysis.

## II. THE CHEMICAL EQUATIONS

The simplest system with a positive feedback loop is a self-regulatory genetic circuit where 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. Also assume that m molecules of proteins have to bind to the inactive gene in order to make it active and turn on the transcription. The active gene D ∗ produces the RNA molecule R , which in turn produces molecules of the protein X via translation mechanism. Both RNA and the protein degrade with time. This processes can be described by the following reactions:

$$D + m X \overset { k _ { 1 } } { \Longrightarrow } D ^ { * } \\ D ^ { * } \overset { a } { \rightarrow } D ^ { * } + R, \quad R \overset { b } { \rightarrow } \varnothing \\ R \overset { c } { \rightarrow } X + R, \quad X \overset { d } { \rightarrow } \varnothing.$$

The rate equation for this process is

$$[ \begin{matrix} 1 \end{matrix}$$

$$\ell ^ { \quad } \\ [ \dot { X } ] & = - d [ X ] + c [ R ] + m k _ { 2 } [ D ^ { * } ] - m k _ { 1 } [ D ] [ X ] ^ { m } \\ & = - d [ X ] + c [ R ] + m ( k _ { 1 } [ X ] ^ { m } + k _ { 2 } ) [ D ^ { * } ] - m k _ { 1 } T [ X ] ^ { m } \\ [ \dot { R } ] & = - b [ R ] + a [ D ^ { * } ] \\ [ \dot { D } ^ { * } ] & = k _ { 1 } [ D ] [ X ] ^ { m } - k _ { 2 } [ D ^ { * } ] \\ & = - ( k _ { 1 } [ X ] ^ { m } + k _ { 2 } ) [ D ^ { * } ] + k _ { 1 } T [ X ] ^ { m } \\ e the condition that the concentration of DNA is a constant.$$

where the condition that the concentration of DNA is a constant,

$$[ D ] + [ D ^ { * } ] = T$$

was used to eliminate [ D ] from the first and the third equation. Here, the concentration [ X ] denotes the number density of the species X divided by some fixed reference number density. Rigorously, it is the expectation value of the concentration. The distinction is often neglected when the fluctuation is small, but it is clear that the stochastic effect for DNA cannot be neglected, since typical gene copy number encoding a given protein is of order one per cell. Assuming there is only one copy of the gene, [ D ] and [ D ∗ ], the expected concentrations of active and inactive genes, are also equal to the probabilities that a gene is in the off or on state, respectively, divided by the volume and the reference number density. As shown in the appendix A, the condition for the validity of Eq.(2) is

$$\langle N _ { D } N _ { X } ^ { m } \rangle = \langle N _ { D } \rangle \langle N _ { X } \rangle ^ { m }$$

where N D and N X denote the numbers of the protein X and the active genes, respectively, and the brackets denote the expectation values. This condition will hold true if the fluctuation of the gene and that of the proteins are independent. In particular, it will hold if the number of protein X is large enough so that the relative fluctuation of its number is negligible. However, although this condition may hold at a certain point of time, the fluctuations of N X may grow in time and make the equation Eq.(2) invalid for later times. The full stochastic master equation or Fokker-Planck equation will have to analyzed to find the time-scale that the Eq.(4) is valid [13-18], perhaps by numerical method, which is out of scope of the current manuscript. Here, I will simply assume the validity of Eq.(4) and analyze the consequence of the equation (2).

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

$$t \rightarrow \tau t,$$

and defining the rescaled variables,

$$x = [ X ] / \epsilon, y = [ R ] / \mu, z = [ D ^ { * } ] / T$$

and the parameters

$$\alpha = b \tau, \beta = d \tau, \gamma = \tau k _ { 2 }, \nu = m \gamma T / \epsilon$$

where µ , /epsilon1 , and τ are chosen to satisfy

$$k _ { 1 } \epsilon ^ { m } / k _ { 2 } = c \mu \tau \epsilon ^ { - 1 } = T \tau a \mu ^ { - 1 } = 1,$$

to get the rescaled set of equations

$$\dot { x } & = - \beta x + y + \nu ( z + z x ^ { m } - x ^ { m } ) \\ \dot { y } & = - \alpha y + z \\ \dot { z } & = - \gamma ( 1 + x ^ { m } ) z + \gamma x ^ { m }$$

$$- \prime$$

which is the object of analysis in this work.

Note that the in the first equation of Eq.(9), the first, second and the third terms describe the change of the protein number due to its degradation, production, and its binding and unbinding to the gene. Since only m molecules of the protein bind to the gene, in most cases of interest where the total number of proteins is much larger than m , the binding and unbinding process will have negligible contribution to the change of the total protein number. This can also be seen from the fact that ν is proportional to T , the concentration of gene which is usually much smaller than [ R ] and [ X ]. Therefore ν is expected to be negligible in most practical applications. However, for the sake of generality, I will keep ν in the rest of the article. Also note that the ν dependent term becomes non-negligible near origin when only a few molecules of proteins are present, but in this limit the stochastic fluctuation of [ X ] cannot be neglected and the deterministic equation (2) will no more be valid anyway(See also the appendix A).

## III. THE DIMENSIONAL REDUCTION AND THE TWO-DIMENSIONAL SYSTEM

The fixed points of the system (9) is obtained by solving ˙ x = ˙ = ˙. y z One first solves ˙ z = 0 to obtain

$$z = \frac { x ^ { m } } { 1 + x ^ { m } }.$$

Substituting into Eq.(9), we get

$$\dot { x } & = - \beta x + y \\ \dot { y } & = - \alpha y + \frac { x ^ { m } } { 1 + x ^ { m } } \\ \cdots \cdots$$

The two-variable system described by Eq.(11) is the one investigated in the literature in detail[19, 20, 22] . 1 It is clear that although such a dimensional reduction is allowed for obtaining fixed points, one has to analyze the original system Eq.(9) in order to correctly analyze the stability of the fixed points. Eq(11) is based on the assumption that ˙ z = 0 is always satisfied. This corresponds to the case when DNA equilibrates on a time scale much faster than that of the RNA or the protein. However, it has to be checked that the fixed point is stable with respect to the perturbation orthogonal to the nullcline ˙ z = 0, which is the main topic of interest in this work.

In order to get the fixed points, we first solve the equation for the nullcline ˙ y = 0, to get

$$y = \frac { x ^ { m } } { \alpha ( 1 + x ^ { m } ) }$$

and substitute into the nullcline equation for x , ˙ x = 0, which after some manipulation becomes

$$x \left ( \alpha \beta x ^ { m } - x ^ { m - 1 } + \alpha \beta \right ) & = 0 & ( 1 3 ) \\ \text{e fixed points are obtained by solving Eq.(13), and those of $y$ and $z$ by}$$

The values of x for the fixed points are obtained by solving Eq.(13), and those of y and z by Eq.(12) and Eq.(10). Since x = 0 a solution of (13), it is always a fixed point. Depending

1 Although a three-dimensional system has been considered in ref.[22], whose fixed points coincide with those of Eq(11), it corresponds to the case where the first equation in Eq(11) is replaced by two linear equations, a rather unnatural situation where protein x produces another protein that acts as a transcriptor for x . The resulting equation is totally different from Eq(9) where only one species of protein is present and only the assumption of instantaneous equilibration of DNA is dropped.

on the values of α and β , there can be a positive real fixed point ˜ that makes the second x factor of Eq.(13) zero:

$$\alpha \beta \tilde { x } ^ { m } - \tilde { x } ^ { m - 1 } + \alpha \beta = 0$$

The stability of the two-dimensional system is analyzed by linearizing the Eq.(11) around a fixed point x ∗ :

The details of the procedure are reviewed in appendix B. Here I just summarize the results: i) m = 1

$$\begin{pmatrix} \dots \dots \\ \delta \dot { x } \\ \delta \dot { y } \end{pmatrix} = \mathbf A _ { 2 } \begin{pmatrix} \delta x \\ \delta y \end{pmatrix} = \begin{pmatrix} - \beta & 1 \\ \frac { m x _ { * } ^ { m - 1 } } { ( 1 + \dot { x } _ { * } ^ { m } ) ^ { 2 } } & - \alpha \end{pmatrix} \begin{pmatrix} \delta x \\ \delta y \end{pmatrix} \\ \text{ils of the procedure are reviewed in appendix B. Here I just summarize the results:}$$

x = 0 is a unique stable fixed point when αβ &lt; 1. For αβ &gt; 1, x = 0 becomes an unstable fixed point, and there is an additional nonzero stable fixed point which is a root of Eq.(14). αβ = 1 is the marginal case where these two fixed points merge.

## ii) m&gt; 1

x = 0 is always a stable fixed point. There are additional nonzero fixed points x ± with x -&lt; x + that correspond to the roots of Eq.(14) for where x -and x + are a saddle point and a stable fixed point, respectively. Therefore, x = 0 is the unique stable fixed point for ( m -1 mαβ ) m &lt; m -1 whereas we get bistability for ( m -1 mαβ ) m &gt; m -1. Again, ( m -1 mαβ ) m = m -1 corresponds to the marginal case where x + and x -merge into one.

$$\cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \\ \left ( \frac { m - 1 } { m \alpha \beta } \right ) ^ { m } > m - 1, \\ \text{saddle- noint. and a stable fixed noint. respectively} \quad \text{Therefore}$$

In order to compare with those obtained from the three-dimensional analysis, we write the eigenvalues of the two-dimensional Hessian A 2 in Eq.(15):

$$\begin{array} { c } \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \ dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \\ \lambda _ { \pm } = \frac { 1 } { 2 } \left [ - \alpha - \beta \pm \sqrt { \left ( \alpha - \beta \right ) ^ { 2 } + 4 \frac { m x _ { * } ^ { m - 1 } } { ( 1 + x _ { * } ^ { m } ) ^ { 2 } } \right ] & & ( 1 7 ) \\ s \text{ the fixed point around which we linearized the equation.  Note that since $x_{*}\geq$} \end{array}$$

where x ∗ is the fixed point around which we linearized the equation. Note that since x ∗ ≥ 0, the quantity inside the square root is non-negative, and therefore the eigenvalues are real. In case of stable fixed point where both of the eigenvalues are negative, this leads to asymptotically exponential convergence:

$$x ( t ) & \to x _ { * } + u _ { + x } e ^ { - \lambda _ { + } t } + u _ { - x } e ^ { - \lambda _ { - } t } \\ y ( t ) & \to y _ { * } + u _ { + y } e ^ { - \lambda _ { + } t } + u _ { - y } e ^ { - \lambda _ { - } t }$$

as t →∞ , where u ± are the eigenvectors corresponding to λ ± , with u ± x and u ± y denoting their x and y components. However, as will be shown later, the analysis of the full threedimensional system Eq.(9) reveals that oscillatory convergence to a stable fixed point is also possible when the conversion between the active and inactive gene is slow enough.

## IV. THE STABILITY ANALYSIS OF THE THREE-DIMENSIONAL SYSTEM

As noted earlier, the stability analysis of the two-dimensional system given by Eq.(11) is rather incomplete. Considering a small deviation of x , y , and z around a fixed point ( x , y ∗ ∗ , z ∗ ) of Eq (9), with x = x ∗ + δx y , = y ∗ + δy , and z = z ∗ + δz , we get a linearized equation

$$\begin{pmatrix} \delta \dot { x } \\ \delta \dot { y } \\ \delta \dot { z } \end{pmatrix} = \mathbf A _ { 3 } \begin{pmatrix} \delta x \\ \delta y \\ \delta z \end{pmatrix} = \begin{pmatrix} - \beta - \nu m x _ { * } ^ { m - 1 } ( 1 - z _ { * } ) & 1 & \nu ( 1 + x _ { * } ^ { m } ) \\ 0 & - \alpha & 1 \\ \gamma m x _ { * } ^ { m - 1 } ( 1 - z _ { * } ) & 0 & - \gamma ( 1 + x _ { * } ^ { m } ) \end{pmatrix} \begin{pmatrix} \delta x \\ \delta y \\ \delta z \end{pmatrix} \\ \text{Solving the characteristic equation for the matrix } \mathbf A _ { 3 },$$

Solving the characteristic equation for the matrix A 3 ,

$$| \lambda - \mathbf A _ { 3 } | & = ( \lambda + \alpha ) ( \lambda + \beta + \nu m x _ { * } ^ { m - 1 } ( 1 - z _ { * } ) ) ( \lambda + \gamma ( 1 + x _ { * } ^ { m } ) ) \\ & \quad - \gamma m x _ { * } ^ { m - 1 } ( 1 - z _ { * } ) - ( 1 + x _ { * } ^ { m } ) ( \lambda + \alpha ) \gamma \nu m x _ { * } ^ { m - 1 } ( 1 - z _ { * } ) \\ & = \lambda ^ { 3 } + ( \alpha + \bar { \beta } + \bar { \gamma } ) \lambda ^ { 2 } + ( \alpha \bar { \beta } + \bar { \beta } \bar { \gamma } + \bar { \gamma } \alpha - \zeta \bar { \nu } ) \lambda + ( \alpha \bar { \beta } \bar { \gamma } - \zeta ( 1 + \alpha \bar { \nu } ) ) = 0, ( 2 0 ) \\.$$

where

$$\bar { \beta } & \equiv \beta + m x _ { * } ^ { m - 1 } ( 1 - z _ { * } ) \nu = \beta + \frac { \nu m x _ { * } ^ { m - 1 } } { x _ { * } ^ { m } + 1 } \\ \bar { \gamma } & \equiv \gamma ( 1 + x _ { * } ^ { m } ) \\ \zeta & \equiv \gamma m x _ { * } ^ { m - 1 } ( 1 - z _ { * } ) = \gamma \frac { m x _ { * } ^ { m - 1 } } { 1 + x _ { * } ^ { m } } \\ \bar { \nu } & \equiv \nu ( 1 + x _ { * } ^ { m } ) \\ \intertext { > 1 \text{ and } } r \ = \nu \ = \varepsilon \ = \theta \text{ is the easiest to analyze } \ \text{$\bar{B} = \beta$} \ \bar { \varepsilon } = \varepsilon$$

The case of m&gt; 1 and x ∗ = y ∗ = z ∗ = 0 is the easiest to analyze. In this case, ¯ = β β , ¯ = γ γ , ¯ = ν ν , ζ = 0, and the solutions to the characteristic equation Eq.(20) are easily shown to be

$$\lambda = - \alpha, - \beta, - \gamma,$$

all of which are negative. Therefore, for m &gt; 1, the origin is always a stable fixed point, as was the case for the two-dimensional system. It remains to analyze the other cases.

In general, the three roots of a cubic equation

$$\lambda ^ { 3 } + p _ { 2 } \lambda ^ { 2 } + p _ { 1 } \lambda + p _ { 0 } = 0$$

are given by[21]

$$\lambda _ { 0 } & = - \frac { p _ { 2 } } { 3 } + \frac { 1 } { 3 } \sqrt { \frac { R + \sqrt { R ^ { 2 } - 4 Q ^ { 3 } } } { 2 } } + \frac { 1 } { 3 } \sqrt { \frac { R - \sqrt { R ^ { 2 } - 4 Q ^ { 3 } } } { 2 } } \\ \lambda _ { \pm } & = - \frac { p _ { 2 } } { 3 } + \frac { e ^ { \mp 2 \pi i / 3 } } { 3 } \sqrt { \frac { R + \sqrt { R ^ { 2 } - 4 Q ^ { 3 } } } { 2 } } + \frac { e ^ { \pm 2 \pi i / 3 } } { 3 } \sqrt { \frac { R - \sqrt { R ^ { 2 } - 4 Q ^ { 3 } } } { 2 } } \\. &.$$

where

$$Q & \equiv p _ { 2 } ^ { 2 } - 3 p _ { 1 } \\ R & = - 2 p _ { 2 } ^ { 3 } + 9 p _ { 1 } p _ { 2 } - 2 7 p _ { 0 }.$$

When the coefficients p 0 , p 1 , and p 2 are all real, we either get three real roots, or one real root and a conjugate pair of complex roots, depending on the sign of

$$R ^ { 2 } - 4 Q ^ { 3 } = - 2 7 \Delta$$

where

$$\Delta \equiv 1 8 p _ { 0 } p _ { 1 } p _ { 2 } - 4 p _ { 2 } ^ { 3 } p _ { 0 } + p _ { 2 } ^ { 2 } p _ { 1 } ^ { 2 } - 4 p _ { 1 } ^ { 3 } - 2 7 p _ { 0 } ^ { 2 }.$$

When ∆ &gt; 0, √ -27∆ = √ R 2 -4 Q 3 is imaginary, and consequently the second and the third terms of each line of Eq.(24) form conjugate pair, so all three roots are real. On the other hand, if ∆ ≤ 0, λ 0 is real, but the second and the third terms of λ ± are linear combinations of distinct real quantities with complex coefficients e ± 2 πi/ 3 , so λ ± form a conjugate pair of complex roots.

For the characteristic equation Eq.(20) we get

$$R & = - 2 ( \alpha + \bar { \beta } + \bar { \gamma } ) ^ { 3 } + 9 ( \alpha + \bar { \beta } + \bar { \gamma } ) ( \alpha \bar { \beta } + \bar { \beta } \bar { \gamma } + \bar { \gamma } \alpha - \zeta \bar { \nu } ) \\ & \quad - 2 7 \left ( \alpha \bar { \beta } \bar { \gamma } - \zeta ( 1 + \alpha \bar { \nu } ) \right ) \\ & = ( \alpha + \bar { \beta } - 2 \bar { \gamma } ) ( \bar { \beta } + \bar { \gamma } - 2 \alpha ) ( \bar { \gamma } + \alpha - 2 \bar { \beta } ) + 2 7 \zeta - 9 ( \bar { \beta } + \bar { \gamma } - 2 \alpha ) \zeta \bar { \nu }$$

and

$$\text{and} \\ \Delta & = - \frac { R ^ { 2 } - 4 Q$$

which will be used for further analysis below.

## A. Limit of instantaneous equilibration of the gene

The limit of instantaneous equilibration of the gene is obtained by taking the limit of γ →∞ . The leading order expansions of R and ∆ are

$$R & = - 2 \gamma ^ { 3 } + 3 ( \alpha + \bar { \$$

We see that since the coefficient of γ 4 in the expansion for ∆ is positive, ∆ &gt; 0 and consequently all the eigenvalues are real when γ is large enough. In this case, we also have

$$\begin{array} {$$

and consequently

$$\text{and consequently} \\ \left ( \frac { 1 } { 2 } ( R \pm$$

Therefore, from Eq.(24) we have

$$\lambda _ { 0 } = - \frac { \alpha + \bar { \beta } + \$$

which is real and negative, and goes to -∞ as γ → ∞ . This shows that the trajectory of the system zaps onto the nullcline ˙ z = 0 as γ →∞ , justifying the two-dimensional system Eq.(11) in such a limit. We also have

$$e \text{ system zaps onto the nullcline } \tilde { z } = 0 \text{ as } \gamma \rightarrow \infty, \text{justifying the two-dimensional system} \\ \mu ( 1 1 ) \text{ in such a limit.  We also have} \\ & \quad \text{Re} \left [ e ^ { 2 \pi i / 3 } \left ( \frac { R \pm \sqrt { - 2 7 \Delta } } { 2 } \right ) ^ { 1 / 3 } \right ] \\ = & \text{Re} \left [ \left ( \frac { - 1 + i \sqrt { 3 } } { 2 } \right ) \left ( - \gamma ( 1 + x _ { ^ { * } } ^ { m } ) + \frac { 1 } { 2 } ( \alpha + \beta ) - \frac { m \nu x _ { ^ { * } } ^ { m - 1 } } { ( 1 + x _ { ^ { * } } ^ { m } ) ^ { 2 } } \\ & \quad \pm i \frac { \sqrt { 3 } } { 2 } \sqrt { ( \alpha - \beta ) ^ { 2 } + \frac { 4 x _ { ^ { * } } ^ { m - 1 } } { ( 1 + x _ { ^ { * } } ^ { m } ) ^ { 2 } } } + O ( \gamma ^ { - 1 } ) \right ) \right ] \\ = & \quad \gamma ( 1 + x _ { ^ { * } } ^ { m } ) - \frac { \alpha + \beta } { 4 } + \frac { m \nu x _ { ^ { * } } ^ { m - 1 } } { 2 ( 1 + x _ { ^ { * } } ^ { m } ) } \\ & \quad \mp \frac { 3 } { 4 } \sqrt { ( \alpha - \beta ) ^ { 2 } + 4 \frac { x _ { ^ { * } } ^ { m - 1 } } { ( 1 + x _ { ^ { * } } ^ { m } ) ^ { 2 } } } + O ( \gamma ^ { - 1 } ), \\ \text{ading to} \\ \lambda. - \frac { \alpha + \beta + \tilde { \gamma } } { \rightharpoonup } \frac { 2 } { R _ { a } } \left [ \left ( - 1 + i \sqrt { 3 } } \right \| \left ( \frac { 1 } { 1 } \right | R + \cdot \sqrt { - 2 7 \Delta } \right \| ^ { 1 / 3 } \right ]$$

leading to which agree with Eq.(17).

$$\text{leading to} \\ \lambda _ { \pm } & = - \frac { \alpha + \bar { \beta } + \bar { \gamma } } { 3 } + \frac { 2 } { 3 } \, \text{Re} \left [ \left ( \frac { - 1 + i \sqrt { 3 } } { 2 } \right ) \left ( \frac { 1 } { 2 } ( R \mp \sqrt { - 2 7 \Delta } ) \right ) ^ { 1 / 3 } \right ] \\ & = \frac { 1 } { 2 } \left [ - \alpha - \beta \pm \sqrt { ( \alpha - \beta ) ^ { 2 } + \frac { 4 m x _ { * } ^ { m - 1 } } { ( 1 + x _ { * } ^ { m } ) ^ { 2 } } } \right ] + O ( \gamma ^ { - 1 } ) \\ \text{which agree with Eq.(17).}$$

## B. General stability analysis of the fixed points

A fixed point is stable if the real parts of three eigenvalues are all negative, whereas it is a saddle point if some of the them have positive real parts. By utilizing the root formula for a cubic equation, one can not only determine the stability of the fixed points, but also obtain more detailed information such as the imaginary parts and the signs of the eigenvalues.

Denoting the three eigenvalues as λ i ( i = 1 2 3), we note that , ,

$$\text{Denoting the three eigenvalues as } \lambda _ { i } ( i = 1, 2, 3 ), \, \text{we note that} \\ \sum _ { i } \lambda _ { i } & = - \alpha - \beta - \gamma ( 1 + x _ { * } ^ { m } ) - \nu \frac { m x _ { * } ^ { m - 1 } } { 1 + x _ { * } ^ { m } } < 0 \\ \sum _ { i < j } \lambda _ { i } \lambda _ { j } & = \alpha \beta + \beta \gamma ( 1 + x _ { * } ^ { m } ) + \gamma ( 1 + x _ { * } ^ { m } ) \alpha + \alpha \frac { \nu m x _ { * } ^ { m - 1 } } { 1 + x _ { * } ^ { m } } > 0 \\ \prod _ { i } \lambda _ { i } & = \gamma \left ( \frac { m x _ { * } ^ { m - 1 } } { 1 + x _ { * } ^ { m } } - \alpha \beta ( 1 + x _ { * } ^ { m } ) \right ). \\ \text{e first line implies that at least one eigenvalue has a negative real part, which is also}$$

The first line implies that at least one eigenvalue has a negative real part, which is also a necessary condition for a fixed point being stable. When ∏ i λ i &gt; 0, the realities and signs of the eigenvalues are either ( - -, , +) or ( C, C , ∗ +), where + and -denote positive and negative real roots respectively, and C, C ∗ denote a conjugate pair of complex roots. Also, since at least one eigenvalue must have negative real part due to the first equation of Eqs.(36), the real parts of the complex conjugate roots are negative.

On the other hand, if ∏ i λ &lt; i 0, the realities and signs are either ( - - -, , ) or ( -, C, C ∗ ). Note that ( -, + +) is excluded due to the second line of Eq (36). , If λ 3 &lt; 0 and λ , λ 1 2 &gt; 0, we have

$$\lambda _ { 1 } \lambda _ { 2 } > | \lambda _ { 3 } | ( \lambda _ { 1 } + \lambda _ { 2 } ) > ( \lambda _ { 1 } + \lambda _ { 2 } ) ^ { 2 }$$

where the last inequality follows from the first line of Eq.(36), which is

$$| \lambda _ { 3 } | > \lambda _ { 1 } + \lambda _ { 2 }.$$

Since we get from Eq.(37)

$$( \lambda _ { 1 } - \lambda _ { 2 } ) ^ { 2 } = ( \lambda _ { 1 } + \lambda _ { 2 } ) ^ { 2 } - 4 \lambda _ { 1 } \lambda _ { 2 } \geq 0,$$

$$\lambda _ { 1 } \lambda _ { 2 } > 4 \lambda _ { 1 } \lambda _ { 2 }$$

which is an obvious contradiction. Therefore, the signs of the eigenvalues cannot be ( -, + +). , The signs of the real parts of the complex roots remain to be determined. Note that two of the roots are complex if and only if -27∆ = R 2 -4 Q &gt; 3 0. Also note that

$$Q & = \alpha ^ { 2 } + \bar { \beta } ^ { 2 } + \bar { \gamma } ^ { 2 } - \alpha \bar { \beta } - \bar { \beta } \bar { \gamma } - \bar { \gamma } \alpha + 3 \bar { \nu } \zeta \\ & = \frac { 1 } { 2 } \left [ ( \alpha - \bar { \beta } ) ^ { 2 } + ( \bar { \beta } - \bar { \gamma } ) ^ { 2 } + ( \bar { \gamma } - \alpha ) ^ { 2 } \right ] + 3 \bar { \nu } \zeta > 0, \quad \text{(41)} \\ \text{issequently } \sqrt { R ^ { 2 } - 4 Q ^ { 3 } } < R, \, \text{so that the sign of } [ \frac { R \pm \sqrt { R ^ { 2 } - 4 Q ^ { 3 } } } { 2 } ] ^ { 1 / 3 } \text{ is the same as that}$$

and consequently √ R 2 -4 Q &lt; R 3 , so that the sign of [ R ± √ R 2 -4 Q 3 2 ] 1 3 / is the same as that of R . From Eqs.(28) and (29), we see that

$$\Delta = - 2 \zeta R + 2 7 \zeta ^ { 2 } + \left ( 4 \zeta \bar { \nu } + \left ( \bar { \beta } - \bar { \gamma } \right ) ^ { 2 } \right ) \left ( \zeta \bar { \nu } - \left ( \alpha - \bar { \beta } \right ) \left ( \alpha - \bar { \gamma } \right ) \right ) ^ { 2 }. \quad \ \ ( 4 2 ) \\ \intertext { since the coefficient of $R$ is negative and the remaining term is positive, we see that if $\Delta < 0$}$$

Since the coefficient of R is negative and the remaining term is positive, we see that if ∆ &lt; 0 so that there are complex roots, then R &gt; 0. Since the real part of the complex roots are

$$\cos \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \min \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \max \text{$\quad $1$} \\ \text{which are sum of three negative numbers, the real parts of the complex roots are negative.}$$

which are sum of three negative numbers, the real parts of the complex roots are negative.

In summary, we get a saddle point for ∏ i λ i &gt; 0 where the three roots of the Hessian are of the form ( -p, -q, r ) or ( -p + iq, - -p iq, r ) ( p, q, r &gt; 0). On the other hand we get a stable fixed point for ∏ i λ i &lt; 0 where the roots are of the form ( -p, -q, -r ) or ( -p + iq, - -p iq, -r ) ( p, q, r &gt; 0). ∏ i λ i = 0 corresponds to the marginal case where the roots are of the form ( -p, -q, 0) or ( -p + iq, - -p iq, 0) ( p, q &gt; 0).

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

$$\prod _ { i } \lambda _ { i } & = \gamma ( 1 - \alpha \beta ) \\ \circ \text{ that the origin is the unique fixed point. } \Pi \, \lambda _ { i } < 0. \text{ so it is a}$$

Therefore, when αβ &gt; 1 so that the origin is the unique fixed point, ∏ i λ i &lt; 0, so it is a stable fixed point, whereas when αβ &lt; 1 so that there is an additional nonzero fixed point, ∏ i λ &gt; i 0 and the origin becomes a saddle point. For additional nonzero fixed points, we have

$$\text{nonzero fixed points, we have} \\ \prod _ { i } \lambda _ { i } & = \gamma ( \frac { m x _ { * } ^ { m - 1 } } { 1 + x _ { * } ^ { m } } - \alpha \beta ( 1 + x _ { * } ^ { m } ) ) \\ & = \gamma ( \frac { m x _ { * } ^ { m - 1 } } { 1 + x _ { * } ^ { m } } - x _ { * } ^ { m - 1 } ) \\ & = \frac { \gamma x _ { * } ^ { m - 1 } } { 1 + x _ { * } ^ { m } } ( m - 1 - x _ { * } ^ { m } ) \\ & & 1 3$$

where Eq.(14) was used when going from the first to the second line. For m = 1, it is negative, so the additional nonzero fixed point is stable whenever it exists. For m &gt; 1, it remains to determine the sign of

$$m - 1 - x _ { * } ^ { m }.$$

Note that the function vanishes at

Substituting x 0 into we get

$$f ( x _ { 0 } ) = m - \frac { 1 } { \alpha \beta } ( m - 1 ) ^ { ( m - 1 ) / m }.$$

Whenever the nonzero fixed points exist, the condition for which is Eq.(16), f ( x 0 ) &lt; 0. Since f (0) &gt; 0 and lim x →∞ f ( x ) &gt; 0, we see that x -&lt; x 0 &lt; x + . Since g x ( ) is a monotonically decreasing function of x , we see that g x ( -) &gt; 0 and g x ( + ) &lt; 0. Therefore, 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 equations Eq.(9). For m = 1, when the degradation rates are large compared to production rates, αβ &gt; 1, then x ∗ = y ∗ = z ∗ = 0 is the unique stable fixed point. When αβ &lt; 1, the origin becomes an unstable fixed point, and the unique stable fixed point appears at nonzero values of the concentrations, at ( x , y ∗ ∗ , z ∗ ) = (1 /αβ -1 1 , /α -β, 1 -αβ ). Since there is a unique stable fixed point regardless of the parameters, we have no bistability for m = 1. On the other hand, for m&gt; 1, x ∗ = y ∗ = z ∗ = 0 is a stable fixed point regardless of the parameter values. For (( m -1) /mαβ ) m &lt; m -1, this is the unique stable fixed point. For (( m -1) /mαβ ) m &gt; m -1 we get an additional stable fixed point x + and a saddle point x -which are the roots to the equation (14), leading to bistability.

$$g ( x ) \equiv m - 1 - x ^ { m }$$

$$x = x _ { 0 } \equiv ( m - 1 ) ^ { 1 / m }.$$

$$f ( x ) \equiv x ^ { m } - \frac { x ^ { m - 1 } } { \alpha \beta } + 1,$$

## C. Complex eigenvalues and the oscillatory behavior

It is to be noted that, in contrast to the two-dimensional system, the eigenvalues of the full three dimensional systems are allowed to have complex values, in the region of the parameter space where ∆ &lt; 0. To be specific, let us consider the Hessian around a stable fixed point. The eigenvector u 0 of A 3 for real eigenvalue λ 0 = -r ( r &gt; 0) is also a real vector, but when ∆ &lt; 0 so that λ ± for a conjugate pair, the corresponding eigenvectors u ± also form a conjugate pair. Denoting λ ± = - ± p iq ( p &gt; 0), we construct real vectors

$$u _ { 1 } & = u _ { + } + u _ { - } \\ u _ { 2 } & = \frac { 1 } { i } ( u _ { + } - u _ { - } )$$

so that

$$\mathbf A _ { 3 } u _ { 1 } = - p u _ { 1 } - q u _ { 2 }, \quad \mathbf A _ { 3 } u _ { 2 } = - p u _ { 2 } + q u _ { 1 }.$$

The parameter p is the rate of approach to the fixed point along the plane spanned by u 1 and u 2 , whereas q is the rate of oscillation between the vectors u 1 and u 2 .

In the limit of t →∞ , the behavior of the variables x , y , and z are thus

$$\left ( \begin{array} { c } x ( t ) \\ y ( t ) \\ z ( t ) \end{array} \right ) = \left ( \begin{array} { c } x _ { * } + A u _ { 0 x } e ^ { - r t } + B e ^ { - p t } \left ( u _ { 1 x } \cos ( q t + \phi ) - u _ { 2 x } \sin ( q t + \phi ) \right ) \\ y _ { * } + A u _ { 0 y } e ^ { - r t } + B e ^ { - p t } \left ( u _ { 1 y } \cos ( q t + \phi ) - u _ { 2 y } \sin ( q t + \phi ) \right ) \\ z _ { * } + A u _ { 0 z } e ^ { - r t } + B e ^ { - p t } \left ( u _ { 1 z } \cos ( q t + \phi ) - u _ { 2 z } \sin ( q t + \phi ) \right ) \end{array} \right ) \\ \text{where } u _ { i j } \text{ denotes the } j - t h \text{ component of the eigenvector } u _ { i }, \text{ and the coefficients } A, B, \text{ and}$$

where u ij denotes the j -th component of the eigenvector u i , and the coefficients A B , , and the phase φ depend on the initial condition.

As an example, let us consider m = 1 with α = β = 1 8, / γ = 1 2 , / 9 ν = 0. We have a nonzero stable fixed point at x ∗ = 1 /αβ -1 = 63, y ∗ = x / α ∗ ( (1 + x ∗ )) = 63 8, / z ∗ = x / ∗ (1 + x ∗ ) = 63 / 64. Since ¯ = 0, ν ¯ = 1 β / 8, ¯ = 1 γ / 8, ζ = 1 2 / 15 , we have

$$\Phi ^ { * } \quad \Phi ^ { * } \, \cdots \$$

leading to one real root and two complex roots for the characteristic equation Eq.(20):

$$\begin{array} {$$

with the corresponding eigenvectors

$$\mathbf u _ { 0 } ^ { T } = ( 3 2, 1, 1 / 3$$

leading to

$$\mathfrak { u } _ { 1 } ^ { T } = ( 1 2 8, -$$

and the asymptotic behavior

$$& \text{and the asymptotic behavior} \\ & \left ( \begin{array}$$

This behavior can be understood using nullclines, just as in the case of two-dimensional system[19]. Note that in the case of the three-dimensional system, a nullcline is a twodimensional surface. Assuming ν = 0 for simplicity, we see from eq.(9) that the nullclines for x y , , and z are the surfaces y = βx , z = αy , and z = x m x m +1 , respectively. Linearizing the equation around the fixed points, we get the planes:

$$\delta y & = \beta \delta x \\ \delta z & = \alpha \delta y \\ \delta z & = \frac { m x _ { ^ { * } } ^ { m - 1 } } { ( x _ { ^ { * } } ^ { m } + 1 ) ^ { 2 } } \delta x. \\ 5 8 ) \text{ with } A & = \phi = 0 \text{ and } B = 1 \text{ is shown in Fig.1. The projection}$$

Part of the trajectory Eq.(58) with A = φ = 0 and B = 1 is shown in Fig.1. The projection to x y -, z y -, and x z -planes are shown in Figure 1 (a),(b), and (c), which may be considered as front, left, and bottom view. The fixed point is denoted by a plus symbol at the origin. The projection of x -nullcline, y -nullcline, and the z -nullcline are indicated in Figures(a),(b), and (c) as dashed lines. The intersections of the trajectory with x -nullcline, y -nullcline, and the z -nullcline are indicated with the open square, open circle, and open triangle respectively, in the figures where the projection of the corresponding nullclines cannot be drawn as lines.

Initially the concentrations of the protein X is decreasing, that of RNA is increasing, and the probability of the gene being in the active state is increasing. As can be shown in

FIG. 1. Projection of the asymptotic trajectory to (a) x y -plane (front view), (b) z y -plane (left view), and (c) x z -plane (bottom view). The projections of nullclines are shown in dashed lines whenever possible. See the text for more details. 17

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0502v2-figure-b4664be0-45f4-46dd-8782-513d3f96db48.jpg)

Figure 1 (c), eventually the trajectory crosses the z -nullcline, indicating that the protein level has decreased to a level where probability of the gene being in the active state begins to decrease. The intersection point is denoted by triangles in Figure 1 (a) and (b). The crossing of the x and y -nullcline is invisible with this resolution, and the protein level x keeps decreasing and the RNA level y keeps increasing until the trajectory goes very close to the fixed point. Figures (d),(e), and (f) are the magnification of (a),(b) and (c) near the fixed point. As can be seen from Fig.1 (d), the trajectory indeed crosses the x -nullcline, indicating that the RNA has been increased to a level so that it can produce more proteins than being degraded, so the protein level x begins to increase. The intersection is denoted as open squares in figures (e) and (f). The crossing of y -nullcline is invisible in this resolution, and the concentrations of both protein and RNA increase, and the probability of the gene being in the active state decreases, until the system goes closer to the fixed point. Figures (g), (h), and (i) are the magnification of (d),(e) and (f) near the fixed point. As can be seen from Fig.1 (h), the trajectory indeed crosses the y -nullcline, indicating that the probability of gene being in the active state has decreased so much that the production of RNA is less than that being degraded, and consequently the RNA concentration begins to decrease. Similar cyclic behavior will repeat in smaller scales as the trajectory converges to the fixed point.

Similar oscillatory behavior exists near a saddle point when the Hessian has complex roots, although the concentrations will eventually diverge away from the fixed point and the linear approximation will break down.

## V. CONCLUSION

In this work, I analyzed the stability of the fixed points of the simplest auto-regulatory system, consisting of DNA, RNA, and the protein product that positively regulate its production upon binding to the gene. After linearizing the equation around 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 from the two-dimensional analysis: When Hill coefficient is one there is only one stable fixed point, that with zero concentrations when the decay rates of the RNA and the protein product are large compared to their production rate, and that with nonvanishing concentration if the production rates

are large compared to the decay rates, zero concentrations being an unstable fixed point in this case. On the other hand, for Hill coefficient larger than one, the vanishing concentration is always a stable fixed point, and an additional nonzero stable fixed point appears when the production rates are large compared to the decay rates.

Although the general conclusion about the stability of the fixed points are the same as the two-dimensional analysis, we could obtain more detailed information such as the number of positive eigenvalues near a saddle point, and the number of complex eigenvalues near an arbitrary fixed point. In particular, we could see that the approach to the fixed point in the three-dimensional system can exhibit oscillatory behavior for certain values of the parameters, a property that cannot be observed in the two-dimensional system. In fact, the two-dimensional system can be justified only in the limit of fast equilibration of gene, in which case the eigenvalues are all real. On the other hand, the novel oscillatory behavior appears when binding and unbinding of protein to the DNA is slow.

Although the analysis is much more general than the two-dimensional analysis, the deterministic equation (2) is valid only when the concentration of protein is large enough so that its relative fluctuation can be neglected. Eventually the fluctuations will grow with time and will no more be negligible. Also, the stochastic effect cannot be neglected near the origin where the number of proteins and RNA are expected to be small. Therefore, the full master equation or Fokker-Planck equation will have to be considered for more accurate analysis[13-18]. Moreover, the current formalism neglects the time delay between the transcription and translation, and that between the translation and the positive regulation. The current formalism is expected to work well for the regulation process in prokaryotes where all the processes happen in more or less the same region in space, and will be less accurate for Eukaryote where the process of transcription and translation happen at distant locations. The effect of time delay due to finite diffusion coefficients will have to be considered for a more accurate description of such systems[17, 18].

## VI. ACKNOWLEDGEMENTS

This work was supported by the National Research Foundation of Korea, funded by the Ministry of Education, Science, and Technology (NRF-2012M3A9D1054705).

## Appendix A: Deviation of the rate equations

Let us denote the numbers of proteins, RNA molecules, and active genes as N X , N Y , and N Z respectively. We assume there is only one copy of the gene, so N Z = 0 or N Z = 1. Denoting the probability for N X = x , N Y = y , and N Z = z as P x, y, z ( ), the master equations corresponding to the process Eq.(1) is

$$\dot { P } ( x, y, z ) & = \tilde { k } _ { 1 } P ( x + m, y, 1 - z ) ( x + m ) ^ { m } z - \tilde { k } _ { 1 } P ( x, y, z ) x ^ { m } ( 1 - z ) \\ & \quad - k _ { 2 } P ( x, y, z ) z + k _ { 2 } P ( x - m, y, 1 - z ) ( 1 - z ) + a P ( x, y - 1, z ) z - \ a P ( x, y, z ) z \\ & \quad + c P ( x - 1, y, z ) y - c P ( x, y, z ) y + b P ( x, y + 1, z ) ( y + 1 ) - b P ( x, y, z ) y \\ & \quad + d P ( x + 1, y, z ) ( x + 1 ) - d P ( x, y, z ) x$$

where ˜ k 1 ≡ k / c V 1 ( 0 ) m is used as the rate constant for binding of the protein X to the gene, where V is the volume of the region where the process is taking place and c 0 is a standard reference concentration. Usage of ˜ k 1 here ensures that we get binding rate constant k 1 when the equations are expressed in terms of [ X ], [ Y ], and [ Z ]. The equations describing the time evolution of the expectation values of N X , N Y , and N Z can be obtained by multiplying Eq.(A1) by x , y , and z and taking summation:

$$\langle \dot { N } _ { \zeta } \rangle = \sum _ { x, y, z } z \dot { P } ( x, y, z ) = & \sum _ { x, y, z } [ \tilde { k } _ { 1 } P ( x + m, y, 1 - z ) ( x + m ) ^ { m } z ^ { 2 } - \tilde { k } _ { 1 } P ( x, y, z ) x ^ { m } z ( 1 - z ) \\ & - k _ { 2 } P ( x, y, z ) z ^ { 2 } + k _ { 2 } P ( x - m, y, 1 - z ) ( 1 - z ) z ] \\ = & \tilde { k } _ { 1 } \sum _ { x, y, z } P ( x, y, z ) x ^ { m } ( 1 - z ) - k _ { 2 } \sum _ { x, y, z } P ( x, y, z ) z \\ & - \tilde { k } _ { 1 } \sum _ { x < m } P ( x, y, z ) x ^ { m } ( 1 - z ) ^ { 2 } \\ = & \tilde { k } _ { 1 } \langle N _ { X } ^ { m } ( 1 - N _ { z } ) \rangle - k _ { 2 } \langle N _ { Z } \rangle \\ & - \tilde { k } _ { 1 } \sum _ { x < m } P ( x, y, z ) x ^ { m } ( 1 - z )$$

$$\langle \dot { N } _ { Y } \rangle = \sum _ { x, y, z } y \dot { P } ( x, y, z ) & = a \sum _ { x, y, z } [ P ( x, y, z ) ( y + 1 ) z - P ( x, y, z ) y z ] \\ & \quad + b \sum _ { x, y, z } [ P ( x, y, z ) y ( y - 1 ) - P ( x, y, z ) y ^ { 2 } ] \\ & \quad - \tilde { k } _ { 1 } \sum _ { x < m } P ( x, y, 1 - z ) x ^ { m } y z ^ { 2 } \\ & = a \sum _ { x, y, z } P ( x, y, z ) z - b \sum _ { x, y, z } P ( x, y, z ) y \\ & \quad - \tilde { k } _ { 1 } \sum _ { x < m } P ( x, y, z ) x ^ { m } y ( 1 - z ) ^ { 2 } \\ & = a \langle N _ { Z } \rangle - b \langle N _ { Y } \rangle \\ & \quad - \tilde { k } _ { 1 } \sum _ { x < m } P ( x, y, z ) x ^ { m } y ( 1 - z ) \\ \\ \langle \dot { N } _ { X } \rangle = \sum _ { x } \dot { P } ( x, y, z ) = \sum [ \tilde { k } _ { 1 } P ( x + m, y, 1 - z ) ( x + m ) ^ { m } x z - \tilde { k } _ { 1 } P ( x, y, z ) x ^ { m + 1 } ( 1 - z )$$

$$\overline { x < m } \\ \langle \dot { N } _ { X } \rangle = \sum _ { x, y, z } x \dot { P } ( x, y, z ) & = \sum _ { x, y, z } [ \tilde { k } _ { 1 } P ( x + m, y, 1 - z ) ( x + m ) ^ { m } x z - \tilde { k } _ { 1 } P ( x, y, z ) x ^ { m + 1 } ( 1 - z ) \\ & \quad - k _ { 2 } P ( x, y, z ) x z + k _ { 2 } P ( x - m, y, 1 - z ) x ( 1 - z ) \\ & \quad + c P ( x - 1, y, z ) x y - c P ( x, y, z ) x y \\ & \quad + d P ( x + 1, y, z ) ( x + 1 ) x - d P ( x, y, z ) x ^ { 2 } ] \\ & = - m \tilde { k } _ { 1 } \sum _ { x, y, z } P ( x, y, z ) x ^ { m } z + m k _ { 2 } \sum _ { x, y, z } P ( x, y, z ) z \\ & \quad + c \sum _ { x, y, z } P ( x, y, z ) y - d \sum _ { x, y, z } P ( x, y, z ) x \\ & \quad - \sum _ { x < m } \tilde { k } _ { 1 } P ( x, y, 1 - z ) x ^ { m } ( x - m ) z \\ & = - m \tilde { k } _ { 1 } \langle N ^ { m } _ { X } N _ { Z } \rangle + m k _ { 2 } \langle N _ { Z } \rangle + c \langle N _ { Y } \rangle - d \langle N _ { X } \rangle \\ & \quad - \sum _ { x < m } \tilde { k } _ { 1 } P ( x, y, z ) x ^ { m } ( x - m ) ( 1 - z ) \\ \text{where we used the fact $z^{2}=z$ and $z(1-z)=0$ since $z=0$ or $z=1$.  When $\langleN_{X}$) is so large} \\ \text{that its fluctuation is small relative to the average value $\langleN_{X}$, we may make approximations}$$

where we used the fact z 2 = z and z (1 -z ) = 0 since z = 0 or z = 1. When 〈 N X 〉 is so large that its fluctuation is small relative to the average value 〈 N X 〉 , we may make approximations

$$\langle N _ { X } ^ { m } N _ { Z } \rangle & \simeq \langle N _ { X } \rangle ^ { m } \langle N _ { Z } \rangle, \\ \langle N _ { X } ^ { m } ( 1 - N _ { Z } ) \rangle & = \langle N _ { X } ^ { m } \rangle - \langle N _ { X } ^ { m } N _ { Z } \rangle \simeq \langle N _ { X } ^ { m } \rangle - \langle N _ { X } \rangle ^ { m } \langle N _ { Z } \rangle = \langle N _ { X } \rangle ^ { m } \langle 1 - N _ { Z } \rangle. \quad ( A 5 )$$

After dividing by c V 0 , last terms in Eqs.(A2), (A3), and (A4) becomes negligible in the limit of V →∞ . Defining [ A ] ≡ 〈 N A 〉 / c V ( 0 ), the rate equation (2) is now derived.

## Appendix B: Fixed points

In this section, I first review how the fixed points of the system (11), or equivalently those of Eq.(9), are determined. In order to get the fixed points, we first solve the equation for the nullcline ˙ y = 0, to get

$$y = \frac { x ^ { m } } { \alpha ( 1 + x ^ { m } ) }$$

and substitute into the nullcline equation for x , ˙ x = 0, which after some manipulation becomes

$$x \left ( \alpha \beta x ^ { m } - x ^ { m - 1 } + \alpha \beta \right ) & = 0 & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & && & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & &  & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & &
 & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & - & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & ; & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & \\ & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & &. & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & } & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & + & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & _ & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & &, & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & = & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & \ & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & b & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & ^ & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & | & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & ( & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & &$$

We see that x = 0 is always a fixed point, and there can be additional fixed points which are the roots of the equation

$$f ( x ) \equiv x ^ { m } - \frac { 1 } { \alpha \beta } x ^ { m - 1 } + 1 = 0$$

for certain values of α and β .

When m = 1, there is one nonzero real root to the equation (B3) if and only if

$$\alpha \beta < 1$$

For m &gt; 1, there can be up to two real roots. In order to investigate the existence of such roots, we first locate the extremum of f ( x ), by taking its derivative and setting it to zero:

$$x ^ { m - 2 } \left ( m x - \frac { m - 1 } { \alpha \beta } \right ) = 0 \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \ \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \cdots \, \quad \quad \quad$$

yielding one extremum at

$$x _ { 1 } = \frac { m - 1 } { m \alpha \beta }$$

and the other one at x = 0 if m &gt; 2. Since f ( j ) (0) = 0 for j &lt; m -1 and f ( m -1) (0) = -( m -1)! /αβ &lt; 0, where f ( n ) denotes the n -th derivative of f ( x ), we see that x = 0 for m&gt; 2 is a local maximum. On the other hand, since

$$\text{cal maximum. On the other hand, since} \\ f ^ { ( 2 ) } ( x _ { 1 } ) & = m ( m - 1 ) \left ( \frac { m - 1 } { m \alpha \beta } \right ) ^ { m - 2 } \\ & \quad - \frac { ( m - 1 ) ( m - 2 ) } { \alpha \beta } \left ( \frac { m - 1 } { m \alpha \beta } \right ) ^ { m - 3 } \\ & = m \left ( \frac { m - 1 } { m \alpha \beta } \right ) ^ { m - 2 } > 0, \\ & \quad 2 2$$

we see that x 1 is the global minimum of f ( x ). Since f (0) = 1 &gt; 0 and lim x →∞ f ( x ) = ∞ &gt; 0, Eq. (B3) will have two distinct real roots if and only if the value of f ( x ) negative at the global minimum x 1 :

$$f ( x _ { 1 } ) & = \left ( \frac { m - 1 } { m \alpha \beta } \right ) ^ { m } - \frac { 1 } { \alpha \beta } \left ( \frac { m - 1 } { m \alpha \beta } \right ) ^ { m - 1 } + 1 \\ & = - \left ( \frac { 1 } { m - 1 } \right ) \left ( \frac { m - 1 } { m \alpha \beta } \right ) ^ { m } + 1 \\ & \quad < 0$$

which can be rewritten as

$$\left ( \frac { m - 1 } { m \alpha \beta } \right ) ^ { m } > m - 1.$$

## Appendix C: The stability analysis of the two-dimensional system

Let us consider the stability of the two-dimensional system Eq.(11) around the fixed points. Most of the material can be found in previous literature[19, 20, 22], but I write down an explicit form of eigenvalues so that they can be compared with those obtained from the three-dimensional system in the limit of ultra-fast binding and unbinding of the protein to the gene. Considering a small deviation of x and y around a fixed point ( x , y ∗ ∗ ) of Eq (11), with x = x ∗ + δx y , = y ∗ + δy , we get a linearized equation

$$\begin{pmatrix} \delta x \\ \delta y \end{pmatrix} = \mathbf A \begin{pmatrix} \delta x \\ \delta y \end{pmatrix} = \begin{pmatrix} - \beta & 1 \\ \frac { m x _ { ^ { m - 1 } } ^ { m - 1 } } { ( 1 + x _ { ^ { * } } ^ { m } ) ^ { 2 } } & - \alpha \end{pmatrix} \begin{pmatrix} \delta x \\ \delta y \end{pmatrix} \\ \text{te characteristic equation for the matrix } A,$$

Solving the characteristic equation for the matrix A ,

$$| \lambda - \mathbf A | & = ( \lambda + \alpha ) ( \lambda + \beta ) - \frac { m x _ { * } ^ { m - 1 } } { ( 1 + x _ { * } ^ { m } ) ^ { 2 } } \\ & = \lambda ^ { 2 } + ( \alpha + \beta ) \lambda + \alpha \beta - \frac { m x _ { * } ^ { m - 1 } } { ( 1 + x _ { * } ^ { m } ) ^ { 2 } } = 0, \\ e \text{eigenvalue}$$

we get the eigenvalue

$$\text{$tne eigenvalue} \\ \lambda _ { \pm } & = \frac { 1 } { 2 } \left [ - \alpha - \beta \pm \sqrt { ( \alpha + \beta ) ^ { 2 } - 4 \alpha \beta + 4 \frac { m x _ { * } ^ { m - 1 } } { ( 1 + x _ { * } ^ { m } ) ^ { 2 } } } \right ] \\ & = \frac { 1 } { 2 } \left [ - \alpha - \beta \pm \sqrt { ( \alpha - \beta ) ^ { 2 } + 4 \frac { m x _ { * } ^ { m - 1 } } { ( 1 + x _ { * } ^ { m } ) ^ { 2 } } } \right ]. \\ & 2 3$$

Since the expression inside the square root is nonnegative, we see that the eigenvalues are real.

From the first line of Eq.(C3), we see that if

$$\frac { m x _ { * } ^ { m - 1 } } { ( 1 + x _ { * } ^ { m } ) ^ { 2 } } < \alpha \beta, & & ( \text{C4} ) \\ \text{more about form is email for the fact form which is not written}$$

the absolute value of the square root term is smaller than the first term which is negative, so signs of both eigenvalues are negative and we get a stable fixed point. On the other hand, if

$$\frac { m x _ { * } ^ { m - 1 } } { ( 1 + x _ { * } ^ { m } ) ^ { 2 } } > \alpha \beta, & & ( \mathbb { C } 5 ) \\ \intertext { o d o r m i n d } \text{br the same root term } \alpha \text{ both entities and}$$

the sign of the eigenvalues are determined by the square root term, so both positive and negative signs appear, and we get a saddle point. It is straightforward see that the equality corresponds to the marginal case where one of the eigenvalues becomes zero. This general result can be used for analyzing the stability of the specific fixed points.

- 1. m = 1
- a. αβ &gt; 1

As discussed in the appendix B, x ∗ = 0 is a unique fixed point in this case. The condition αβ &gt; 1 is the same as (C4) with x ∗ = 0, so this is a stable fixed point.

$$b. \ \alpha \beta < 1$$

It is clear that x ∗ = 0 is now a saddle point, since now (C5) holds with x ∗ = 0. There is an additional fixed point which is a solution to the equation (B3):

$$x _ { * } = \frac { 1 } { \alpha \beta } - 1.$$

We see that

$$\frac { m x _ { * } ^ { m - 1 } } { ( 1 + x _ { * } ^ { m } ) ^ { 2 } } = \frac { 1 } { ( 1 + x _ { * } ) ^ { 2 } } = ( \alpha \beta ) ^ { 2 } < \alpha \beta. \text{and} \, \text{and} \, \text{and} \, \text{and} \, \text{and} \, \text{and} \, \text{and} \, \text{and} \, \text{and} \, \text{and} \, \text{and} \, \text{and} \, \text{and} \, \text{and} \, \text{and} \, \text{and} \, \text{and} \, \text{and} \, \text{and} \, \text{and} \, \text{and} \, \text{and} \, \text{and} \, \text{and} \, \text{and} \, \text{and} \, \text{and} \, \text{and} \, \text{and} \, \text{or} \, \text{or} \, \text{or} \, \text{or} \, \text{or} \, \text{or} \, \text{or} \, \text{or} \, \text{or} \, \text{or} \, \text{or} \, \text{or} \, \text{or} \, \text{or} \, \text{or} \, \text{or} \, \text{or} \, \text{or} \, \text{or} \, \text{or} \, \text{or} \, \text{or} \, \text{or} \, \text{or} \, \text{or} \, \text{or} \, \text{or} \, \text{or} \, \text{or} \,\text{or} \, \text{or} \, \text{or} \, \text{or} \, \text{or} \, \text{or} \, \text{or} \, \text{or} \, \text{or} \, \text{or} \, \text{or} \, \text{or} \, \text{or} \, \text{or} \, \text{or} \, \text{or} \, \text{or} \, \text{or} \, \text{or} \, \text{or} \, \text{or} \, \text{or} \, \text{or} \, \text{or} \, \text{or} \, \text{or} \, \text{or} \, \text{or} \, \text{or},$$

Therefore, the condition (C4) is satisfied and we get a stable fixed point.

To summarize, for m = 1, when the RNA and protein is degraded fast enough so that αβ &gt; 1, x ∗ = 0 is the unique stable fixed point of (11). Eventually the concentrations

of both species go to zero. On the other hand, when the degradation is slow enough so that αβ &lt; 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.

## 2. m&gt; 1

In this case, x ∗ = 0 is always a stable fixed point since

$$\frac { m x _ { * } ^ { m - 1 } } { ( 1 + x _ { * } ^ { m } ) ^ { 2 } } = 0 < \alpha \beta.$$

As discussed in the appendix B, when Eq.(B9) holds, we have additional nonzero fixed points which are the two distinct real roots x -&lt; x + of the equation (B3). After multiplying αβ (1 + x m ∗ ) 2 to both sides of Eq.(C4) and Eq.(C5), we see that we get stable fixed point or saddle point depending on whether the sign of the expression

$$\alpha \beta m x _ { \pm } ^ { m - 1 } - ( \alpha \beta ) ^ { 2 } ( 1 + x _ { \pm } ^ { m } ) ^ { 2 } = x _ { \pm } ^ { m - 1 } ( m \alpha \beta - x _ { \pm } ^ { m - 1 } )$$

is negative or not, where Eq.(B3) was used in going from the second to the third expression.

Therefore, we only have to examine the sign of

$$m \alpha \beta - x _ { \pm } ^ { m - 1 }.$$

Let us define

$$g ( x ) \equiv m \alpha \beta - x ^ { m - 1 }.$$

we see that this function is zero when

$$x = x _ { 0 } \equiv ( m \alpha \beta ) ^ { 1 / ( m - 1 ) }.$$

On the other hand, we have

$$f ( x _ { 0 } ) = ( m \alpha \beta ) ^ { m / ( m - 1 ) } - ( m - 1 ) < 0$$

where the last inequality follows from the condition (B9). Therefore, we see that

$$x _ { - } < x _ { 0 } < x _ { + }.$$

Since g x ( ) is a monotonically decreasing function, it immediately follows that g x ( -) &gt; 0 and g x ( + ) &lt; 0. Therefore, x -is a saddle point and x + is a stable fixed point. To summarize, when RNA and the protein are degraded fast enough so that

$$\left ( \frac { m - 1 } { m \alpha \beta } \right ) ^ { m } < m - 1,$$

x = 0 is a unique sable fixed point. When the degradation is slow enough so that the inequality is inverted, then we have additional nonzero fixed points x -&lt; x + , where x -is a saddle point and x + is a sable fixed point. Since x = 0 and x = x + are both stable fixed points, we get bistability.

- [1] H. de Jong, J. Comput. Biol. 9 , 67 (2002).
- [2] P. Smolen, D. A. Baxter, and J. H. Byrne, Bull. Math. Biol. 62 , 247 (2000).
- [3] S. Pigolotti, S. Krishina, and M. H. Jensen, Proc. Natl. Acad. Sci. U.S.A. 104 , 6533 (2007).
- [4] N. Radde, Bioinformatics 26 , 2874 (2010).
- [5] D. Angeli, F. E. Ferrel, and E. D. Sontag, Proc. Natl. Acad. Sci. U.S.A. 101 , 1822 (2004).
- [6] A. D. Keller, J. Theor. Biol. 170 , 175 (1994).
- [7] T. olde Scheper, D. Klinkenberg, C. Pennartz, and J. van Pelt, J. Neurosci. 19 , 40 (1999).
- [8] P. Smolen, D. A. Baxter, and J. H. Byrne, J. Neurosci. 21 , 6644 (2001).
- [9] E. H. Snoussi, J. Biol. Sys. 6 , 3 (1998).
- [10] E. Plahtye, T. Mestl, and W. S. Omholt, J. Biol. Sys. 3 , 409 (1995).
- [11] J. L. Gouz´, J. Biol. Sys. e 6 , 11 (1998).
- [12] O. Cinquin and J. Demongeot, J. Theor. Biol. 216 , 229 (2002).
- [13] H. Ge and H. Qian, Phys. Rev. Lett. 103 , 148103 (2009).
- [14] H. Ge and H. Qian, Chaos 22 , 023140 (2012).
- [15] H. Ge and H. Qian, J. R. Soc. Interface 8 , 107 (2011).
- [16] A. Lipshtat, A. Loinger, N. Q. Balaban, and O. Biham, Phys. Rev. Lett. 96 , 188101 (2006).
- [17] W. Wu and J. Wang, J. Chem. Phys. 139 , 121920 (2013).
- [18] W. Wu and J. Wang, J. Phys. Chem. B 117 , 12908 (2013).
- [19] S. H. Strogatz, Nonlinear Dynamics and Chaos (Westview Press, Perseus Book Publishing, Cambridge, MA, 1994).

- [20] J. S. Griffith, Mathematical Neurobiology (Academic Press, New York, 1971).
- [21] G. V. Milovanovi´, D. S. Mitrinovi´, Th. M. Rassias, c c Topics in Polynomials: extremal Problems, Inequalities, Zeros (World Scientific, Singapore, 1994).
- [22] J. S. Griffith, J. Theoret. Biol. 20 , 209 (1968).
Julian Lee
Analytic stability analysis of three-component self-regulatory genetic circuit
A self-regulatory genetic circuit, where a protein acts as a positive
regulator of its own production, is known to be a simplest form of biological
network with a positive feedback loop. Although at least three components, DNA,
RNA, and the protein, are required to form such a circuit, the stability
analysis of fixed points of the self-regulatory circuit has been performed only
after reducing the system into to a two-component system consisting of RNA and
protein only, assuming a fast equilibration of the DNA component. Here, the
stability of fixed points of the three-component positive feedback loop is
analyzed by obtaining eigenvalues of full three dimensional Hessian matrix. In
addition to rigorously identifying the stable fixed points and the saddle
points, detailed information can be obtained, such as the number of positive
eigenvalues near a saddle point. In particular, complex eigenvalues is shown to
exist for sufficiently slow binding and unbinding of the auto-regulatory
transcription factor to DNA, leading to oscillatory convergence to the steady
states, a novel feature unseen in the two-dimensional analysis.