First passage times and minimum actions for a stochastic minimal bistable system

Chinese Journal of Physics 59 (2019) 220–230

Contents lists available at ScienceDirect

Chinese Journal of Physics journal homepage: www.elsevier.com/locate/cjph

First passage times and minimum actions for a stochastic minimal bistable system

T

⁎,a

Hongwei Yina,b, Xiaoqing Wen a b

School of Sciences, Nanchang University, Nanchang 330031, PR China Numerical Simulation and High-Performance Computing Laboratory, Nanchang University, Nanchang 330031, PR China

A R T IC LE I N F O

ABS TRA CT

Keywords: Bistability Chemical master equation Mean first passage time Minimum actions

Bistability is an ubiquitous phenomenon in biological systems, and always plays important roles in cell division, differentiation, cancer onset, apoptosis and so on. However, stochastic fluctuations in bistable systems are still hard to understand. To address this issue, we propose a chemical master equation model for a minimal bistable system, which underlies generally bistable systems. For this master equation model, we mainly focus on the mean first passage times (MFPTs) by respectively using Gillespie algorithm and an approximation method of the large deviation theory, and does on minimum actions along optimal transition paths from OFF to ON states by the large deviation theory. Further, we find that for this stochastic system the MFPTs have different change tendencies compared to the corresponding minimum actions. Our results of this minimal stochastic model can also well understand more general bistable systems.

1. Introduction Bistability is extremely important to understand basic phenomena of cellular function, such as decision-making processes in cell cycle progression, cell differentiation, apoptosis, switch in gene express, cancer onset and prion diseases [1–6]. A common viewpoint is that bistable switches are usually produced by positive feedback loops in signal transduction networks in [7]. Bistability with outstanding biological importance always attracts the attention of theoretical biologists, who addresses various bistable models in real systems. For example, a sufficiently strong external signal can switch on a self-amplifying process, leading to expression of the corresponding target genes. Even though there is no any persistent signal, this process can retain this state because of the hysteresis. However, comlexity of these models due to stronger non-liner terms urges researchers to seek some analytical or approximated methods for identifying bistability, see [8–10]. To identify the bistable system, the simple approach is the identification of necessary structural conditions [11]. Some different minimal bistable systems are gradually uncovered, for instance, cell polarity [12], G protein signalling [13] and mass-preserving chemical reaction network [14]. The minimal bistable systems are underlying to really bistable systems, because some real bistable systems are constructed through adding some feedback loops, such as positive and negative feedback loops and feedforward loops etc, to the minimal bistable systems. Therefore, research on the minimal bistable systems is important to understand more general bistable systems. In order to distinguish minimal bistable systems, Wilhelm in [11] proposed some criteria: (1) minimal numbers of reactants; (2) minimal numbers of reactions, and (3) minimal number of terms in the ODEs. According to these criteria, Wilhelm proposed a minimal universal bistable system, whose reaction process is below,

⁎

Corresponding author. E-mail addresses: [email protected] (H. Yin), [email protected] (X. Wen).

https://doi.org/10.1016/j.cjph.2019.02.009 Received 13 September 2018; Received in revised form 28 December 2018; Accepted 14 February 2019 Available online 19 February 2019 0577-9073/ © 2019 Published by Elsevier B.V. on behalf of The Physical Society of the Republic of China (Taiwan).

Chinese Journal of Physics 59 (2019) 220–230

H. Yin and X. Wen k1

k2

k3

k4

S + Y → 2X , 2X → X + Y , X + Y → Y + P , X → P ,

(1.1)

where reaction rates ki, i = 1, …, 4 are constant, and S and P respectively denote substrates and wasters. X and Y are reactants. Wilhelm in [11] analyzed the bistability of the ODE model of (1.1). It is well known that stochastic fluctuations or noise, in micro systems, are common phenomena and unavoidable, stemming from either thermal effects or environment variability. For bistable systems, the fluctuations are viewed as a kind of energy for inducing switches between bistable states [15–17]. These noises result from the fluctuations of external and internal systems. For the external noise, the systems can be described by the Langevin equations. This kind of systems can be easily analyzed by Ito’s integral [18]. When the number of the reaction pieces is lower, the internal noise plays a leading role to the systems, and can be given by the form of the chemical master equation [19]. However, this kind of equations has almost no the analytical solution. Previous studies for the minimal reaction process (1.1) focused more on investigation of its corresponding ODE model, but less on stochastic switches of this system due to internal noise. On the other hand, the mean first passage time (MFPT), which refers to the mean of a random time when a random variable first passes a threshold value, is a key index of quantifying stochastic switches. It has been often used to examine dynamical behaviors of stochastic multistable systems and explain decision-making of cells [16,20,21]. Wang et al. in [22] investigated the mean extinction time for a metapopulation system induced by non-Gaussian multiplicative noise, and by the fast descend method and the path integral approach. Effect of time delay on a bistability system is examined by using delay Fokker-Planck equation and projection operator method, see [23]. Guo et al.in [24,25] showed the MFPTs of Gaussian noise. Besides, they studied information entropy and stochastic resonance induced by external noise in [26,27]. Besides, transition paths for stochastic bistable systems are another interesting topic. When the events with very little likelihood occur, a rough estimate for the probability that the trajectory lies in a given region are obtained by the large deviation theory (LDT), which is the field of mathematics concerned with the analysis of rare events and statistical outliers. Ever since the LDT proposed in [28], the LDT has been widely used to intuitive illustration of stochastically dynamic systems, such as population models [29–33], gene regulatory networks [34–36] and epidemic models [37–40]. This theoretical method is a methodology to understand the metastability of stochastic systems, and to identify the optimal transition paths between metastable states and compute the transition rates, etc. In this paper, we study a chemical master equation model of the minimal bistable system (1.1), which can better describe intrinsic fluctuations resulting from the lower numbers of reaction species. For this model, we first investigate the MFPTs of the transition from OFF to ON states for different reaction parameters by applying respectively the Gillespie algorithm and the large deviation principle. Further, for these parameters we research optimal transition paths from OFF to ON states and minimum actions along these pathes through the geometric minimum action method (gMAM) [41]. We show relationships between the MFPTs and minimum actions with respect to the different reaction parameters, respectively. 2. The chemical master equation model According to the reaction processes (1.1), we will propose a model in the chemical master equation. First, denote the state of the system n = (n1, n2)T with n1 and n2 representing the numbers of the species X and Y, respectively. Interactions among these reaction species change the variable n in steps of size rj in the jth reaction. For the four reactions in (1.1), the state change vectors rj respectively are r1 = (2, −1)T , r2 = (−1, 1)T , r3 = (−1, 0)T and r4 = (−1, 0)T . Once the jth reaction fires, the state of the system n would be updated to n + rj . The reaction propensity functions are determined by each term and the system size V, where V −1 can characterizes the magnitude of intrinsic fluctuations [42]. The four propensity functions, respectively, are

W1 (n) = k1 n2, W2 (n) = k2 n1 (n1 − 1)/ V , W3 (n) = k3 n1 n2/ V , W4 (n) = k 4 n1. We choose these forms for the propensities because the stochastic process ni (t ), i = 1, 2 will tend to the deterministic process when the volume size V tends to infinity in this scaling. Let P(n, t) represent the probability of the system at the state n and on the time t. The probability master equation is governed by

∂P (n , t ) = ∂t

4

∑ [Wj (n − rj) P (n − rj, t ) − Wj (n) P (n, t )]. (2.1)

j=1

Define a vector x = (x1, x2), x1 = n1/ V and x2 = n2/ V , and let wj (x) = Wj (V x) = Wj (n) . Then, we can rewrite (2.1) as

∂P (x, t ) = ∂t

4

∑ [wj (x − rj/V ) P (x − rj/V , t ) − wj (x) P (x, t )]. (2.2)

j=1

This form of the master equation suggests that 1/V ≪ 1 when V ≫ 1, and hence a Taylor expansion around x of wi and P is legitimate at least far from the boundaries. This is known as the Kramers–Moyal expansion [43]. Truncating to the second order in V, one can obtain a Fokker–Planck equation 2

∂P (x, t ) ∂ 1 = −∑ [Ai (x) P (x, t )] + ∂t ∂x i 2V i=1

2,2

∑ i, k = 1

∂2 [Bi, k (x) P (x, t )], ∂x i ∂xk

where 221

(2.3a)

Chinese Journal of Physics 59 (2019) 220–230

H. Yin and X. Wen

Ai (x) =

2 ⎧− k2 x1 + (k2/ V − k3 x2 − k 4 ) x1 + 2k1 x2 , i = 1, kx ⎨ k2 x12 − 2 1 − k1 x2 , i = 2, V ⎩

(2.3b)

and 2 ⎧ (k2 x1 + (k3 x2 + k 4 ) x1 + 4 k1 x2) V − k2 x1 , ⎪ V ⎪ 2 − 2 k Vx + k x Vk x − 2 1 1 2 2 1 ⎪ , ⎪ 2 V Bi, k (x) = 2 ⎨ −Vk2 x1 − 2 k1 Vx2 + k2 x1 , ⎪ 2V ⎪ ⎪ k1 Vx2 + k2 x1 (Vx1 − 1) , ⎪ V ⎩

i = k = 1, i = 1, k = 2, i = 2, k = 1, i = k = 2.

(2.3c)

As the volume V extends to infinity, we recover the mean-field equation given by the Liouville equation

∂P (x, t ) ∂ [(−k2 x12 + 2k1 x2 − k3 x1 x2 − k 4 x1) P (x1, x2 , t )] = − ∂t ∂x1 ∂ [(k2 x12 − k1 x2) P (x1, x2 , t )], − ∂x2

(2.4)

whose solution, for the initial condition P (x, 0) = δ (x − x 0), has the simple form, where δ is the Dirac delta function. This form of P implies that x(t), the relative component densities, is the solution to the mean-field system:

dx1 = 2k1 x2 − k2 x12 − k3 x1 x2 − k 4 x1, dt

(2.5a)

dx2 = k2 x12 − k1 x2 . dt

(2.5b)

3. Results 3.1. Analysis of the mean-field system The mean-field system (2.5) always has the equilibrium point (0,0) denoted by E0. Let k1c = 4k3 k 4/ k2 . By some calculation, one can know that if k1 < k1c then the system (2.5) only has one equilibrium point E0; if k1 = k1c , the system (2.5) has one positive equilibrium ⎛ 2k k ⎞ point Ec = ⎜ 4 , 4 ⎟, but for k1 > k1c , the system (2.5) has another two positive equilibrium points. Thus, k1c is a critical parameter. To ⎝ k2 k3 ⎠ identify the bifurcation type of the parameter k1c , we first calculate Jacobin matrix of the model (2.5) at Ec

⎡−6k 4 6k3 k 4 ⎤ ⎢ k2 ⎥ Ac = D f(Ec , k1c ) = ⎢ ⎥. 4 ⎢ 4k 4 − k3 k 4 ⎥ k2 ⎥ ⎢ ⎦ ⎣ Obviously, Ac has a simple eigenvalue λ = 0 with the eigenvector V = (k3, k2)T . Its transpose matrix AcT also has the zero eigenvalue with the corresponding eigenvector W = (2, 3)T . We use f = (f1 , f2 )T to denote the vector of the reaction field of (2.5), where f1 (f2) is the right hand side of Eq. (2.5a) (Eq. (2.5b)). One can calculate and obtain

f k1 (Ec , k1) = (2k 4/ k3, −k 4/ k3)T , and

−2k2 −k3 −k3 0 ⎤ . Bc = D 2 f(Ec , k1c ) ⎡ ⎢ 2k2 0 0 0⎥ ⎣ ⎦

(3.1)

Further, after some complex calculation, one has

WT f k1 (Ec , k1) = k 4/ k3, WT [Bc (V , V)] = −2k32.

(3.2)

According to Sotomayor theorem in [44], we can know that (Ec , k1c ) is a saddle-node bifurcation point. Here, we take k1 = 8, k2 = k3 = 1 and k 4 = 1.5, the second and third steady states are E1 = (2, 0.5) and E2 = (6, 4.5), respectively. Fig. 1A shows the bistability of the system (2.5) induced by the parameter k1 because of the saddle-node bifurcation. Similarly, we can get the saddlenode bifurcations induced by the other parameters k2, k3 and k4, see Fig. 1B–D. 222

Chinese Journal of Physics 59 (2019) 220–230

H. Yin and X. Wen

A

B

C

D

Fig. 1. The saddle-node bifurcations of the model (2.5). The blue lines denote the stable steady states, but the yellow lines do unstable steady states. In (A), k2 = k3 = 1 and k 4 = 1.5. In (B), k1 = k3 = 1 and k 4 = 1.5. In (C), k1 = k2 = 1 and k 4 = 1.5. In (D), k1 = k2 = k3 = 1.

3.2. Optimal transition paths Next, we will discuss optimal transition paths and minimum actions for the model (2.1) from OFF state (the low numbers of X reaction species) to ON one (the large numbers of X reaction species). The LDT gives a rough estimate for the probability that the trajectory Xε(t), t ∈ [0, T], T < ∞, of the random dynamical system lies in a small neighborhood around a given path ψ ∈ C(0, T), where C(0, T) denotes the space of all continuous functions mapping from [0, T] into n in [45]. This theory asserts that, for δ and ε sufficiently small,

 ⎧ sup X ε (t ) − ψ (t ) ≤ δ ⎫ = exp( −ε −1ST (ψ)), ⎬ ⎨ ⎭ ⎩0 ≤ t ≤ T

(3.3)

where sup is short for supremum, which means the least upper bound of a subset.  denotes the probability conditioned on and assume ψ (0) = x . The action can be written as

ST (ψ) =

T ⎧ ∫0 L (ψ, ψ˙ ) dt , if ψ ∈ C (0, T ) is absolutely continuous and the integral converges, ⎨ + ∞, otherwise, ⎩

X ε (0)

=x

(3.4)

where the Lagrangian L(x, y) is given by

L (x , y ) = sup (〈y, θ〉 − H (x , θ)).

(3.5)

θ ∈ n

223

Chinese Journal of Physics 59 (2019) 220–230

H. Yin and X. Wen

Fig. 2. Minimum actions along optimal transition paths of transitions from OFF to ON. (A) k2 = k3 = 1 and k 4 = 1.5. (B) k1 = k3 = 1 and k 4 = 1.5. (C) k1 = k2 = 1 and k 4 = 1.5. (D) k1 = k2 = k3 = 1.

Here ⟨ · , · ⟩ denotes the Euclidean scalar product in n and H(x, θ) is a Hamiltonian, x is the generalized coordinate, y is the time derivative of x and θ is generalized momentum. Here, our aim is to construct exponential estimates for the probability densities to observe paths from OFF to ON states. We look for the solution of (2.1) by proposing a large deviation form of the stationary distribution, i.e., the stationary distribution for (2.2) is 4

Pst (x)=

∑ [wj (x − rj/V ) Pst (x − rj/V ) − wj (x) Pst (x)]. (3.6)

j=1

Pst (x) = C exp(−VS (x))

(3.7)

where C is a normalized factor and S(x1, x2, t) is a quantity known as the action [30,46–49]. Putting (3.7) into (2.2). After some ∂S algebra, one arrives at a Hamilton–Jacobic equation ∂t = −H . At leading order of V, the Hamilton–Jacobi equation has the form H (x, ∂ x1 S, ∂ x2 S ) = 0, where H, known as the effective Hamiltonian, is given as

H = k1 x2 (e2 p x1 − p x 2 − 1) + k2 x 2 (e−p x1 + p x 2 − 1) + k3 x1 x2 (e−p x1 − 1) + k 4 x1 (e−p x1 − 1),

(3.8)

where the conjugate momentum px1 = ∂ x1 S (x) and px2 = ∂ x2 S (x) . The solutions to H (x, p) = 0 with p = (px1 , px2 ), are the zero224

Chinese Journal of Physics 59 (2019) 220–230

H. Yin and X. Wen

A

B

C

D

E

F

G

H

I

Fig. 3. The optimal transition path and momentums. A–C: k1 = 8.52 and other parameters are taken as in Fig. 2A. D–F: k3 = 0.0684 and other parameters are taken as in Fig. 2C. G–I: k 4 = 0.025 and other parameters are taken as in Fig. 2D.

energy curves of the system. At least one solution is the optimal path where the action S is minimum. This solution corresponds to the path that maximizes the transition probability. The associated Hamiltons equations are found by the relations

⎧ dx1 = ∂H , dpx1 = − ∂H , ⎪ ∂px1 dt ∂x1 ⎪ dt dpx2 ⎨ dx2 ∂H ∂H = =− , . ⎪ ⎪ dt ∂px2 dt ∂x2 ⎩

(3.9)

whose detailed forms are as follows

225

Chinese Journal of Physics 59 (2019) 220–230

H. Yin and X. Wen

A

B

C

D

Fig. 4. The MFPTs from OFF state to ON one for the different parameters. (A) k2 = k3 = 1 and k 4 = 1.5. (B) k1 = k3 = 1 and k 4 = 1.5. (C) k1 = k2 = 1 and k 4 = 1.5. (D) k1 = k2 = k3 = 1.

dx1 = −k2 x12 e−p x1 + p x 2 − (e−p x1 x2 k3 + e−p x1 k 4) x1 + 2 k1 x2 e2 p x1 − p x 2 , dt dpx1 = −2 k2 x1 (e−p x1 + p x 2 − 1) − k3 x2 (e−p x1 − 1) − k 4 (e−p x1 − 1), dt dx2 = −k1 x2 e2 p x1 − p x 2 + k2 x12e−p x1 + p x 2 , dt dpx2 = −k1 (e2 p x1 − p x 2 − 1) − k3 x1 (e−p x1 − 1). dt Hamiltons equations (3.9) can describe the systems dynamics in the four-dimensional phase space. The x1 − x2 dynamics along the deterministic line of px1 = 0 and px2 = 0, are described by

dx1 ∂H |p = 0 , = dt ∂px1 x1

dx2 ∂H |p = 0 , = dt ∂px2 x 2

which is the rescaled mean-field rate equations associated with the mean-field equations (2.5). We numerically calculate the minimum actions, Smin, along the optimal transition paths for the different reaction parameters, referring to Fig. 2, by using the geometric minimum action method in [41,50,51]. In Fig. 2A, we take the parameters k2 = k3 = 1 and k 4 = 1.5, and find that Smin gets a globally minimum value at k1 = 8.52 . For this case, we further plot the optimal path from OFF to ON 226

Chinese Journal of Physics 59 (2019) 220–230

H. Yin and X. Wen

A

B

C

D

Fig. 5. The MFPTs approximated by (3.10). (A)–(D) The diagrams of τ with respect to the parameters k1, k2, k3 and k4. The other parameters are respectively taken as in Fig. 4.

in Fig. 3A. In addition, Fig. 3B and C respectively show the curves for the momentum px1 as a function of x1 and for the momentum px2 as a function of x2, both of them being not zero. Fig. 2B shows that the function Smin increases as k2 with k1 = k3 = 1 and k 4 = 1.5, implying that no globally minimum values for k2 except for its boundary point. However, in Fig. 2C with k1 = k2 = 1 and k 4 = 1.5, Smin as a function of k3 reaches a globally minimum value at k3 = 0.0684 . For the parameters in this case, we respectively plot the optimal transition path from OFF to ON in Fig. 3D, the momentum px1 as a function of x1 in Fig. 3E, and the momentum px2 as a function of x2 in Fig. 3F, where px1 and px2 are not zero. In Fig. 2D with k1 = k2 = k3 = 1, we also find that there exists a globally minimum value for Smin at k 4 = 0.025. Similarly, for this case, we fix the parameters and plot Fig. 3G–I. Fig. 3G displays the optimal transition path from OFF to ON, Fig. 3H does a curve of the momentum px1 as a function of x1, and Fig. 3I does a curve of the momentum px2 as a function of x2. Likewise, the momentums px1 in Fig. 3H and px2 in Fig. 3I are not zero. From these numerical results and the meaning of the parameters k1, k3 and k4, we can know that the production and transformation of the piece X can form some optimal paths from OFF to ON states. 3.3. Mean first passage times In this subsection, we will study the MFPTs of Eq. (2.1) in the bistable parameter regimes, from statistical views, by performing 227

Chinese Journal of Physics 59 (2019) 220–230

H. Yin and X. Wen

A

B

C

D

Fig. 6. The residuals of the functions τ for the parameters ki, i = 1, 2, 3, 4 . The other parameters are taken as in Fig. 4.

numerical simulation for Eq. (2.1) by the Gillespie algorithm [19]. Since the bistable regimes of (2.1) strongly depend on the parameters ki, i = 1, …, 4 (see Fig. 1), the MFPTs are respectively simulated under the different parameters, see Fig. 4. Fig. 4A with k2 = k3 = 1 and k 4 = 1.5 is a curve of the MFPT against k1, and shows that the MFPT decrease with the parameter k1 except for the bifurcation point k1c . This relationship is due to the role of k1 as the production rate for X in the first reaction. Larger k1 implies more production of X, resulting that the MFPT of X is shorter. But, this law does not hold for k2, see Fig. 4B with k1 = k3 = 1 and k 4 = 1.5, in which there exists a maximum value of the MFPT at k2 = 9.714 . In addition, the MFPT increases as k3, referring to Fig. 4C with k1 = k2 = 1 and k 4 = 1.5, but it decreases as k4, referring to Fig. 4D with k1 = k2 = k3 = 1. These results in Fig. 4 imply that the MFPTs tightly depend on the reaction parameters in the system (2.1). In the above discussion, we mainly applied statistical methods for numerical results of (2.1) to address the MFPTs of Eq. (2.1). Next, we will use an approximation method of the LDP to further examine the MFPTs of Eq. (2.1) and uncover the relationship between the LDP and the MFPTs. Given a quasi-potential Φ, the MFPT τ to escape from the fixed point at the OFF state takes the general Arrhenius form [52]

τ∼

2π exp{V [Φ(x 0 ) − Φ(x OFF )]}, |Φ″ (x 0)|Φ″ (x OFF )

(3.10)

where Φ(xOFF) denotes the quasi-potential at the OFF state, which is defined by − ln[P (x OFF )], Φ(x0) defined by − ln[P (x 0)] on an unstable fixed point of (2.5). Φ(x 0) − Φ(x OFF ) is the value for the action along the optimal path from xOFF to x0. Here, we perform 228

Chinese Journal of Physics 59 (2019) 220–230

H. Yin and X. Wen

Gillespie algorithm for (2.1), computer statistical probability of (2.1) at each state, and calculate Φ(x0) and Φ(xOFF). We demonstrate τ according to (3.10) in Fig. 5. This figure shows that the changing tendencies of τ are similar to the corresponding MFPTs in Fig. 4 for the parameters ki, i = 1, 2, 3, 4, respectively. To examine the errors of the approximated functions τ, in Fig. 6 we respectively plot residuals of τ for the different parameters, compared to the accurate MFPTs of numerical simulations. Fig. 6 shows that Eq. (3.10) can approximate well to the accurate MFPTs from the numerical simulation for (2.2). 4. Conclusion In this paper, we have proposed the master equation model for the minimum bistable system, which is fundamental for some real bistable systems, such as biochemical reaction systems, gene regulatory networks and so on, and analyzed the saddle-node bifurcations for the corresponding mean-field model. Besides, the MFPTs and optimal transition paths have been respectively examined by applying the Gillespie algorithm and the LDT, which are classical and valid to address stochastic dynamics. Via performing theoretical analysis and numerical simulations, we have found that for the reactions in which the MFPT can arrive at a maximum value, the corresponding minimum action on the transition path from OFF to ON is not globally optimal, comparing Fig. 4B with Fig. 2B, and that for the reactions in which the minimal actions have the globally optimal values, but the MFPTs have no maximum values, see Figs. 4 and Fig. 2. Therefore, the scalars of the MFPT and minimum action present opposite characteristics for the stochastic bistable system (1.1). These results is interesting and significant. However, this result is limited in the minimal bistable system in this paper. Whether or not it will still be sure for other bistable systems is uncertain, and will be answered in our future works. Acknowledgments We are grateful to the anonymous referee for his or her valuable comments, which lead to a substantial improvement in the contents of this paper. This work was supported by grant nos 61563033 (H.Y) and 11563005 (X.W)from the National Natural Science Foundation of China. References [1] F. Ortega, J.L.G. s, F. Mas, B.N. Kholodenko, M. Cascante, Bistability from double phosphorylation in signal transduction, FEBS J. 273 (2006) 3915–3926. [2] X. Chen, Y.-M. Kang, Y.-X. Fu, Switches in a genetic regulatory system under multiplicative non-gaussian noise, J. Theor. Biol. 435 (2017) 134–144. [3] S.N. Semenov, L.J. Kraft, A. Ainla, M. Zhao, M. Baghbanzadeh, V.E. Campbell, K. Kang, J.M. Fox, G.M. Whitesides, Autocatalytic, bistable, oscillatory networks of biologically relevant organic reactions, Nature 537 (2016) 656–660. [4] C. Gupta, J.E.M.L. Opez, W. Ott, K.I.C.S.J.I.A.C.E.C. Fi, M.R. Bennett, Transcriptional delay stabilizes bistable gene networks, Phys. Rev. Lett. 111 (2013) 058104. [5] S. Rati, E. Roberts, Gradient sensing by a bistable regulatory motif enhances signal amplification but decreases accuracy in individual cells, Phys. Biol. 13 (2016) 036003. [6] X.-P. Zhang, Z. Cheng, F. Liu, W. Wang, Linking fast and slow positive feedback loops creates an optimal bistable switch in cell signaling, Phys. Rev. E 76 (2007) 031924. [7] U. Alon, An Introduction to Systems Biology: Design Principles of Biological Circuits, CRC Press, 2006. [8] M. Assaf, E. Roberts, Z. Luthey-Schulten, Determining the stability of genetic switches: explicitly accounting for mrna noise, Phys. Rev. Lett. 106 (2011) 24810224. [9] N. Friedman, L. Cai, X.S. Xie, Linking stochastic dynamics to population distribution: an analytical framework of gene expression, Phys. Rev. Lett. 97 (2006). [10] Y.T. Lin, C.R. Doering, Gene expression dynamics with stochastic bursts: construction and exact results for a coarse-grained model, Phys. Rev. E 93 (2016) 022409. [11] T. Wilhelm, The smallest chemical reaction system with bistability, BMC Syst. Biol. 3 (2009) 90. [12] Y. Mori, A. Jilkine, L. Edelsteinkeshet, Wave-pinning and cell polarity from a bistable reaction-diffusion system, Biophys. J. 94 (2008). 3684–97 [13] D. Csercsik, K.M. Hangos, G.M. Nagy, A simple reaction kinetic model of rapid (g protein dependent) and slow (β-arrestin dependent) transmission, J. Theor. Biol. 255 (2008). 119–28 [14] A. Shiu, The Smallest Multistationary Mass-Preserving Chemical Reaction Network, Springer Berlin Heidelberg, 2010. [15] G.H. Peter, Yen, Phenotypic switching of populations of cells in a stochastic environment, J. Stat. Mech. 2018 (2018) 023501. [16] Y. Xu, J. Feng, J. Li, H. Zhang, Lvy noise induced switch in the gene transcriptional regulatory system, Chaos 23 (2013) 083101. [17] B. Nowakowski, A.L.K. ski, Stochastic transitions between attractors in a tristable thermochemical system: competition between stable states, React. Kinet. Mech. Catal. (2017) 1–11. [18] D.J. Higham, An algorithmic introduction to numerical simulation of stochastic differential equations, SIAM Rev. 43 (2001) 525–546. [19] D. Higham, Modeling and simulating chemical reactions, SIAM Rev. 50 (2008) 347–368. [20] S. Mayank, A. Surendran, Estimation of mean first passage time for bursty gene expression, Phys. Biol. 13 (2016) 036004. [21] Y.T. Lin, T. Galla, Bursting noise in gene expression dynamics: linking microscopic and mesoscopic models, J. R. Soc. Interface 13 (2016). [22] K.-K. Wang, Y.J. Wang, S.-H. Li, J.-C. Wu, Mean extinction time and stability for a metapopulation system subjected to correlated gaussian and non-gaussian noises, Chin. J. Phys. 54 (2016) 205–215. [23] T. Yang, R. Liu, C. Zeng, W. Duan, Effects of time delay and non-gaussian noise on the dynamics of a perceptual bistability, Chin. J. Phys. 55 (2017) 275–288. [24] Y.-F. Guo, B. Xi, Y.-J. Shen, J.-G. Tan, Mean first-passage time of second-order and under-damped asymmetric bistable model, Appl. Math. Model. 40 (2016) 9445–9453. [25] Y.-F. Guo, B. Xi, F. Wei, J.-G. Tan, The mean first-passage time in simplified Fitzhugh Nagumo neural model driven by correlated non-gaussian noise and gaussian noise, Mod. Phys. Lett. B 32 (2018) 1850339. [26] Y. Guo, Y. Shen, J. Tan, Stochastic resonance in a piecewise nonlinear model driven by multiplicative non-gaussian noise and additive white noise, Commun. Nonlinear Sci. Numer. Simul. 38 (2016) 257–266. [27] Y.-F. Guo, Y.-J. Shen, J.-G. Tan, Effects of colored noise and external periodic force on the time derivative of information entropy for a damped harmonic oscillator, Appl. Math. Comput. 252 (2015) 20–26. [28] M.I. Dykman, E. Mori, J. Ross, P.M. Hunt, Large fluctuations and optimal paths in chemical kinetics, J. Chem. Phys. 100 (1994) 5735–5750. [29] B. Meerson, P.V. Sasorov, Wkb theory of epidemic fade-out in stochastic populations, Phys. Rev. E 80 (2009) 041130. [30] M. Assaf, B. Meerson, Extinction of metastable stochastic populations, Phys. Rev. E 81 (2010) 021116.

229

Chinese Journal of Physics 59 (2019) 220–230

H. Yin and X. Wen

[31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42] [43] [44] [45] [46] [47] [48] [49] [50] [51] [52]

M. Khasin, B. Meerson, E. Khain, L.M. Sander, Minimizing the population extinction risk by migration, Phys. Rev. Lett. 109 (2012) 138104. B. Meerson, N.R. Smith, Extinction of oscillating populations, Phys. Rev. E 93 (2016) 032109. A. Michael, B. Meerson, Wkb theory of large deviations in stochastic populations, J. Phys. A. 50 (2017) 263001. J.M. Horowitz, R.V. Kulkarni, Stochastic gene expression conditioned on large deviations, 2017. S. Chaudhury, Modeling the effect of transcriptional noise on switching in gene networks in a genetic bistable switch, J. Biol. Phys. 41 (2015) 235–246. J. Newby, Bistable switching asymptotics for the self regulating gene, J. Phys. A. 48 (2015) 185001. E. Forgoston, S. Bianco, L. Shaw, I. Schwartz, Maximal sensitive dependence and the optimal path to? Epidemic extinction, Bull. Math. Biol. 73 (2011) 495–514. D.S. nchez Taltavull, A. Vieiro, T.A. n, Stochastic modelling of the eradication of the HIV-1 infection by stimulation of latently infected cells in patients under highly active anti-retroviral therapy, J. Math. Biol. 73 (2016) 919–946. C. Hanshuang, F. Huang, Epidemic extinction in a generalized susceptible-infected-susceptible model, J. Stat. Mech. 2017 (2017) 013204. G.T. Nieddu, L. Billings, J.H. Kaufman, E. Forgoston, S. Bianco, Extinction pathways and outbreak vulnerability in a stochastic Ebola model, J. R. Soc. Interface 14 (2017). M. Heymann, E. Vanden-Eijnden, The geometric minimum action method: a least action principle on the space of curves, Commun. Pure Appl. Math. 61 (2008) 1052–1117. N.G. van Kampen, Stochastic Processes in Physics and Chemistry, North holland, 2007. C.W. Gardiner, Stochastic Methods : A Handbook for the Natural and Social Sciences, Springer, New York, 2009. L. Perko, Differential Equations and Dynamical Systems, Springer, Berlin, 2001. M.I. Freidlin, A.D. Wentzell, Random perturbations of dynamical systems, Springer, New York, 1984. D. Kessler, N. Shnerb, Extinction rates for fluctuation-induced metastabilities: a real-space wkb approach, J. Stat. Phys. 127 (2007) 861–886. B.E. Shay, Michael, Rare events in stochastic populations under bursty reproduction, J. Stat. Mech. 2016 (2016) 113501. M. Bauver, E. Forgoston, L. Billings, Computing the optimal path in stochastic dynamical systems, Chaos 26 (2016) 083101. M. Parker, A. Kamenev, B. Meerson, Noise-induced stabilization in population dynamics, Phys. Rev. Lett. 107 (2011) 180603. C. Lv, X. Li, F. Li, T. Li, Energy landscape reveals that the budding yeast cell cycle is a robust and adaptive multi-stage process, PLoS Comput. Biol. 11 (2015) e1004156. C. Lv, L. Xiaoguang, Constructing the energy landscape for genetic switching system driven by intrinsic noise, PLoS One 9 (2014) e88167. P.C. Bressloff, Stochastic Processes in Cell Biology, Springer International Publishing, 2014.

230