Early signatures of regime shifts in gene expression dynamics

Complex dynamical systems, ranging from ecosystems and the climate to gene regulatory networks and financial markets, are known to exhibit abrupt shifts from one dynamical regime to a contrasting one at critical parameter values [1–4]. Examples of sudden regime shifts include the collapse of vegetation when semi-arid conditions prevail, the transition from a clear lake to a turbid one, sudden changes in fish/wildlife populations [4, 5], distinct changes in the climate and patterns of oceanic circulation [1, 6], financial markets undergoing global crashes [1], spontaneous systemic failures such as asthma attacks [7] or epileptic seizures [8] etc. In a gene regulatory network, a sudden transition may occur from one stable gene expression state to another at a critical parameter value [9–12]. The induction of the lac operon in E. coli results in the synthesis of the protein β-galactosidase required for breaking up sugar molecules and releasing energy to the cell. There is now experimental evidence [9] that an abrupt transition from the uninduced (β-galactosidase level low) state to the induced (β-galactosidase level high) state of the lac operon occurs at a critical value of an inducer concentration.

The regime/state shifts in the examples mentioned above are mostly a consequence of the fold-bifurcation (or the fold-catastrophe), well characterized in dynamical systems theory [1, 2, 6, 13]. Figure 1 illustrates a specific type of fold-bifurcation based on bistability and hysteresis, which provides a physical understanding of the features associated with sudden regime shifts. The plot represents the steady states of a dynamical system versus a specific parameter. The state of the dynamical system at time t is defined by the magnitudes of the relevant variables at t. In the steady state, the net rate of change in the magnitude of a variable is zero so that there is no further time evolution. The solid lines in figure 1 denote stable steady states separated by a branch (dotted line) of unstable steady states. The steady state curve is folded backwards giving rise to bistability, i.e., the existence of two stable steady states in the shaded region. The parameter values u₁ and u₂ represent the bifurcation points at which the abrupt regime shifts from one stable steady state to another occur. On crossing the bifurcation point, the system loses bistability and becomes monostable. The transition from one branch of stable steady states to the other is not reversible but exhibits hysteresis. This implies that if a transition takes place at a bifurcation point, say the upper one, the reverse transition from the upper to the lower branch occurs at the lower bifurcation point and not at the upper bifurcation point itself. Bistability owes its origin to the presence of one or more positive feedback loops governing the system dynamics with the added condition that the dynamics be sufficiently nonlinear [10–12]. Each stable steady state has its own basin of attraction in the state space, i.e., the space of all states. In the case of more than one stable steady state, the location of the initial state of the system in a basin of attraction determines the steady state attained by the system in the course of time.

Zoom In Zoom Out Reset image size

The stability of a steady state indicates that the system returns to the state on being weakly perturbed from it. This is shown by arrow directions in figure 1. If the perturbation is sufficiently strong so that a transition takes place from one basin of attraction to the other, then a switch occurs between the stable steady states even before the bifurcation point is reached. The closer the system is to the bifurcation point, the lower the magnitude of the perturbation needed to bring about the switch. In the example of figure 1, the branch of unstable steady states constitutes the border between the two basins of attraction. The distance of a stable steady state from the corresponding unstable steady state is a measure of the resilience (robustness against perturbation) of the system. The resilience is gradually weakened as the system approaches a bifurcation point. The dynamics of natural systems, in general, have a stochastic component due to the presence of random external influences and the inherently probabilistic nature of the processes involved in the dynamics. Consider the dynamics of gene expression in a gene regulatory network. Gene expression consists of two major steps: transcription and translation during which messenger RNA (mRNA) and protein molecules respectively are synthesized. The biochemical events (e.g., the binding of a RNA polymerase at the promoter region of the DNA to initiate transcription) constituting gene expression are inherently stochastic in character resulting in fluctuations (noise) around mean mRNA and protein levels [14, 15]. Extrinsic influences such as variability in the number of regulatory molecules also contribute to the noise. Instead of a single steady state protein level, as in the case of deterministic time evolution, one now has a steady state probability distribution (SSPD) in the protein levels reflecting the stochastic nature of the time evolution. In the presence of low/moderate amounts of noise, the physical picture underlying sudden regime shifts still remains valid. As in the case of applied perturbations, fluctuations of sufficiently strong magnitude can bring about regime shifts before the bifurcation point is crossed. A number of early signatures of regime shifts have been proposed so far [1, 6, 16, 17] in the scenario of the fold-bifurcation. These are the critical slowing down (CSD) and its associated effects, namely, rising variances and the lag-1 autocorrelation function as the critical transition point is approached, increased skewness in the SSPD and the presence of flickering transitions between the stable steady states. The utility of such signals in the cases of ecological and financial systems [1] and complex deteriorating diseases [18] is significant, especially in developing appropriate risk-aversion/management strategies. Also, in the absence of physical models capturing the essential dynamics, one can obtain quantitative measures of the early signatures by analyzing time-series data.

In this paper, we compute quantities providing early signatures of abrupt regime shifts in gene expression dynamics. We specifically focus on the phenomenon of competence development in the soil bacteria B. subtilis involving binary gene expression, i.e., the existence of two expression states. We consider two alternative models of competence development based on bistable and excitable dynamics. In section 2, the two models are described and their different dynamics highlighted. In section 3, we define the quantitative measures of the early signatures of sudden regime shifts. Section 4 contains the results of our computations on these early signatures. Section 5 contains concluding remarks.

Microorganisms such as bacteria have to cope with a number of stresses during their lifetime. The bacteria adopt a number of strategies to enhance their chances of survival under changed circumstances [12, 19]. One of these is the generation of phenotypic heterogeneity (no genetic change) in the bacterial population so that a fraction of the population, even if small, may survive under stressful conditions. In the B. subtilis population, only a part of the population, in which the level of a key transcription factor ComK is high, develops competence [12, 20]. The rest of the population is in the so-called vegetative state. The ComK protein acts as a master regulator activating the transcription of several genes including those essential for the uptake of foreign DNA from the environment. The incorporation of the new DNA into the bacterial genome confers favorable traits on the bacterial subpopulation with a high ComK level, enabling the subpopulation to adapt to stress. The ComK activity resulting in the development of competence is in turn controlled by a host of other proteins. The core module of the complex regulatory network consists of an autoregulatory positive feedback loop in which the ComK proteins promote their own production. The positive feedback gives rise to binary gene expression in the cell population, i.e., two distinct subpopulations with low and high ComK levels, respectively. Two independent experiments [21, 22] provide confirmation that an autoregulatory positive feedback loop of ComK production is by itself sufficient to establish binary gene expression in a bacterial population. Considering deterministic time evolution, the steady state ComK level versus an appropriate gene expression parameter exhibits a hysteresis curve, resulting from the fold-bifurcation, similar to the one shown in figure 1 [23]. In this scenario, if the cells in a population are prepared to be in the same initial state, all the cells evolve to the same final state giving rise to a homogeneous cell population. The generation of heterogeneity in the form of two distinct subpopulations is brought about by fluctuation-driven transitions between the low and high ComK expression states, the fluctuations being associated with the ComK levels. This gives rise to the experimentally observed [12, 20] bimodal distribution in the ComK levels, i.e., a distribution with two prominent peaks. There is now experimental evidence [24] that a lower fraction of the B. subtilis population develops competence with reduced noise in the low ComK level. An alternative physical mechanism underlying competence development has been proposed by Süel et al [25] in terms of excitability in the dynamics of the genetic circuit. The excitable core module consists of both positive and negative feedback loops which bring about a transient activation to the high ComK state. In this case, there is only one stable steady state (low ComK level) and two unstable steady states, the lower of which, in terms of the ComK level, sets a threshold for switching [25]. Fluctuations in the low ComK level activate the switch to expression states in the neighborhood of the state with a high ComK level which constitutes an unstable steady state. The transient activation is followed by an ultimate return to the stable low ComK state. At any instant of time, the population divides into two subpopulations with the low and high ComK levels respectively, signifying a different origin for the bimodal distribution. Quantitative time lapse fluorescence microscopy provides experimental evidence [25] of the probabilistic and transient differentiation into competence.

We first consider model 1 in which the protein product of a single gene autoactivates its own production via a positive feedback loop. The genetic circuit is shown in figure 2(a) and, as mentioned earlier, constitutes a core module of the complex network resulting in competence development. The proteins synthesized by the comK gene form dimers. The dimer molecules bind at the promoter region of the DNA and activate gene expression. The gene also synthesizes proteins at a basal level. A kinetic scheme of the model is as follows [23]:

$$\begin{eqnarray} &&P_{2}+G \;{\mathop{\rightleftharpoons}^{k_{1}}_{k_{2}}}\;GP_{2} \; {\mathop {\rightleftharpoons}^{k_{a}}_{k_{d}}}\; {G^{*}}\quad\nonumber \\ &&G \; \mathop{\rightarrow}^{J_{0}} P\quad\nonumber \\ &&G^{*} \mathop{\rightarrow}^{J_{1}} P\quad\nonumber \\ &&P+P \mathop{\rightleftharpoons}_{K_{e}} P_{2}\qquad\nonumber \\ &&P \mathop{\rightarrow}^{k_{p}} \phi . \qquad \end{eqnarray} \tag{ 1 }$$

The gene can be in two possible states: inactive (G) and active (G*). Proteins are synthesized with a rate constant J₁(J₀) in the state G*(G) with J₀ ≪ J₁. The synthesized proteins dimerize with K_e as the equilibrium dissociation constant. The protein dimer P₂ binds to the gene in the state G and forms the complex GP₂, with k₁ and k₂ being the rate constants for binding and unbinding. The complex GP₂ in turn is activated to the state G*, with the rate constants k_a and k_d being the activation and deactivation rate constants, respectively. The synthesized proteins are degraded with a rate constant, k_p, ϕ denoting the degradation product. As shown by Karmakar and Bose [23], the kinetic scheme displayed in equation (1) can be mapped onto a simpler scheme

$$\begin{eqnarray} &&G \; {\mathop{\rightleftharpoons}^{k^{\prime }_{a}}_{k^{\prime }_{d}}} \;G^{*} \nonumber \\ &&G \; {\mathop{\rightarrow}^{J_{0}}} \;P \nonumber \\ &&G^{*} \;{\mathop{\rightarrow}^{J_{1}}} \;P \nonumber \\ &&P+P \;{\mathop{\rightleftharpoons}_{K_{e}}}\; P_{2} \nonumber \\ &&P\; {\mathop{\rightarrow}^{k_{p}}} \; \phi . \end{eqnarray} \tag{ 2 }$$

The effective activation and deactivation rate constants k'_a(x) and k'_d are

$$\begin{equation} k_{a}^{\prime }(x)=k_{a}\frac{(x/k_{s})^{2}}{1+(x/k_{s})^{2}},\;k_{d}^{\prime }=k_{d}. \end{equation} \tag{ 3 }$$

In equation (3), x denotes the protein concentration and $k_{s}=\sqrt{\frac{k_{2}}{k_{1}}K_{e}}$. In the simplified scheme, the rate of change of the protein concentrations is given by

$$\begin{equation} \frac{{\rm d}x}{{\rm d}t}=\frac{J_{1}k_{a}^{\prime }}{k_{a}^{\prime }+k_{d}}+\frac{J_{0}k_{d}}{k_{a}^{\prime }+k_{d}}-k_{p}x. \end{equation} \tag{ 4 }$$

The steady state condition $\frac{{\rm d}x}{{\rm d}t}=0$ yields three solutions in a specific parameter regime corresponding to two stable steady states separated by an unstable steady state. The model dynamics exhibit the fold-bifurcation of the type shown in figure 1. The rate constants J₀, J₁, and k_a serve as the bifurcation parameters. Figure 2(b) shows the genetic circuit (model 2), proposed in [25], as governing the dynamics of competence development via excitability. The circuit contains the autoregulatory positive feedback loop of figure 2(a). In addition, there is a negative feedback loop in which the ComK protein represses the expression of the comS gene, whereas the ComS protein inhibits the degradation of ComK by the MecA–ClpP–ClpC complex. The repression of comS by ComK is, however, not well established under wild-type expression conditions [20]. The dynamics of the model are described in terms of the differential equations [25]

$$\begin{equation} \frac{{\rm d}K}{{\rm d}t}=a_{k}+\frac{b_{k}\: K^{n}}{k_{0}^{n}+K^{n}}-\frac{K}{1+K+S} \end{equation} \tag{ 5 }$$

$$\begin{equation} \frac{{\rm d}S}{{\rm d}t}=\frac{b_{s}}{1+(K/k_{1})^{p}}-\frac{S}{1+K+S} . \end{equation} \tag{ 6 }$$

The variables K and S denote the concentrations of the ComK and ComS proteins, respectively, a_k and b_k represent the basal and fully activated rates of ComK synthesis and k₀ is the ComK concentration needed for 50% activation. The Hill coefficients n and p are indicative of the cooperativities of ComK autoactivation and ComS repression, respectively. The expression rate of ComS has maximal value b_s and is half-maximal at K = k₁. The nonlinear degradation terms are a consequence of the MecA complex-mediated degradation with a competitive inhibition of the degradation by ComS.

Zoom In Zoom Out Reset image size

The deterministic dynamics of model 1 are governed by the rate equation shown in equation (4). A simple stochastic version of the model has been studied in [23]. In this model, the only stochasticity considered is associated with the random transitions of the gene between the inactive (G) and active (G*) states. The rest of the processes undergo deterministic time evolution. The simple model yields bimodal protein distributions in the steady state in certain parameter regions. In this paper, the stochastic dynamics of the model are investigated using the formulations based on the Langevin and Fokker–Planck (FP) equations. The steady state analysis of a bistable system in these formalisms is described in [26, 28–31]. The one-variable Langevin equation (LE) including both multiplicative and additive noise terms is given by

$$\begin{equation} \frac{{\rm d}x}{{\rm d}t}=f(x)+g_{1}(x)\varepsilon (t)+\Gamma (t), \end{equation} \tag{ 7 }$$

where ε(t) (multiplicative noise) and Γ(t) (additive noise) are Gaussian white noises with mean zero and correlations:

$$\begin{eqnarray} &&\fl \langle \varepsilon (t)\varepsilon (t^{\prime })\rangle = 2d_{1}\delta (t-t^{\prime })\nonumber \\ &&\fl \langle \Gamma (t)\Gamma (t^{\prime })\rangle = 2d_{2}\delta (t-t^{\prime })\nonumber \\ &&\fl \langle \varepsilon (t)\Gamma (t^{\prime })\rangle = \langle \Gamma (t)\varepsilon (t^{\prime })\rangle =2\lambda _{d}\sqrt{d_{1}d_{2}}\delta (t-t^{\prime }). \end{eqnarray} \tag{ 8 }$$

In equation (8), d₁ and d₂ denote the strengths of the noises ε(t) and Γ(t) respectively and λ_d is the degree of correlation between them. The first term in equation (7) describes the deterministic dynamics. In the case of model 1, f(x) is given by the right-hand side expression in equation (4). The FP equation is a rate equation for the probability distribution P(x, t), obtained from the LE as [28, 36]

$$\begin{equation} \frac{\partial P(x,t)}{\partial t}= -\frac{\partial }{\partial x}[A(x)P(x,t)]+\frac{\partial ^{2}}{\partial x^{2}}[B(x)P(x,t)], \end{equation} \tag{ 9 }$$

where

$$\begin{equation} A(x)=f(x)+d_{1}g_{1}(x)g_{1}^{\prime }(x)+\lambda _{d}\sqrt{d_{1}d_{2}}g_{1}^{\prime }(x) \end{equation} \tag{ 10 }$$

and

$$\begin{equation} B(x)=d_{1}[g_{1}(x)]^{2}+2\lambda _{d}\sqrt{d_{1}d_{2}}g_{1}(x)+d_{2} . \end{equation} \tag{ 11 }$$

The SSPD is given by [28–30]

$$\begin{eqnarray} &&\fl P_{\rm st}(x) = \frac{N}{B(x)}\exp \left[\int ^{x}\frac{A(x)}{B(x)}\,{\rm d}x\right]\nonumber \\ &&{\hphantom{\fl P_{\rm st}(x) }}=\frac{N}{\lbrace d_{1}[g_{1}(x)]^{2}+2\lambda _{d}\sqrt{d_{1}d_{2}}g_{1}(x)+d_{2}\rbrace ^{\frac{1}{2}}}\nonumber \\ &&\times \exp \left[\int ^{x}\frac{f(x^{\prime })\,{\rm d}x^{\prime }}{d_{1}[g_{1}(x^{\prime })]^{2}+2\lambda _{d}\sqrt{d_{1}d_{2}}g_{1}(x^{\prime })+d_{2}}\right],\nonumber \\ &&\; \end{eqnarray} \tag{ 12 }$$

where N is the normalization constant. Equation (12) can be rewritten in the form

$$\begin{equation} P_{\rm st}(x)=N{\rm e}^{-\phi _{F}(x)} \end{equation} \tag{ 13 }$$

with

$$\begin{eqnarray} &&\fl \phi _{F}(x)= \frac{1}{2}\ln \Big[d_{1}[g_{1}(x)]^{2}+2\lambda _{d}\sqrt{d_{1}d_{2}}g_{1}(x)+d_{2}\Big]\nonumber \\ &&-\int ^{x}\frac{f(y)\,{\rm d}y}{d_{1}[g_{1}(y)]^{2}+2\lambda _{d}\sqrt{d_{1}d_{2}}g_{1}(y)+d_{2}}, \end{eqnarray} \tag{ 14 }$$

where ϕ_F(x) defines the 'stochastic potential' of the dynamics. Once the SSPD is determined, quantities such as the variance and skewness, to be defined below, which provide early signatures of regime shifts can be determined.

We next define the quantitative measures of the early signatures of regime shifts. The first signature is that of the CSD, which is a distinguishing feature of the dynamics close to a bifurcation point [16, 17, 32]. This involves a progressively larger relaxation time, as the bifurcation point is approached, to an attractor of the dynamics (say, a stable steady state) after being weakly perturbed from it. The physical origin of this phenomenon can be understood in terms of the stability landscape U(x). The rate equation governing the dynamics of model 1 (equation (4)) can be written as $\frac{{\rm d}x}{{\rm d}t}=f(x)=-\frac{\partial U(x)}{\partial x}$. Figure 3(a) shows the phase diagram of the model in the k_a–k_d plane with a region of bistability separating two regions of monostability. The other parameter values are J₁ = 1, J₀ = 0.01 and k_p = 0.03 in appropriate units. Figure 3(b) shows the three stability landscapes at three successively decreasing values of k_a as one approaches the bifurcation point in the direction of the arrow from within the region of bistability. At point 1 (k_a = 1.5), the states corresponding to the minima of U(x) represent the stable steady states. The location of the 'ball' represents the state of the dynamical system. As one progresses from point 1 to point 2 (k_a = 1.35), the steady state with a high value of x becomes less stable. The associated basin of attraction becomes flatter with reduced size so that it takes a longer time for a perturbed state (ball shifted from the minimum position) to relax back to the stable steady state. The relaxation time diverges at the bifurcation point where the stable steady state is on the verge of losing stability (point 3, k_a = 1.344). A weak perturbation pushes the ball to the minimum with a low value of x, i.e., the system does not relax back to the high x state. One can define a return time T_R which provides a measure of the time taken by the dynamical system to regain a stable steady state after being weakly perturbed from it. Let x_st be the stable steady state value of x and η(t) = x(t) − x_st be the small deviation from the steady state value under weak perturbation. The deterministic rate equation is given by $\frac{{\rm d}x}{{\rm d}t}=f(x).$ On Taylor expanding f(x) around x = x_st and keeping terms up to the order of η(t), one obtains

$$\begin{equation} \frac{{\rm d}\eta }{{\rm d}t}=\lambda \eta (t),\;\lambda =\left.\frac{{\rm d}f(x)}{{\rm d}x}\right|_{x=x_{\rm st}} . \end{equation} \tag{ 15 }$$

The solution of equation (15) is

$$\begin{equation} \eta (t)=\eta (0)\,{\rm e}^{\lambda t}, \end{equation} \tag{ 16 }$$

where η(0) is the initial value of η(t) at time t = 0. The sign of the parameter λ indicates the stability of the steady state. If λ is <0 (>0), then the steady state is stable (unstable). We designate λ as the stability parameter. It is a well-known result from dynamical systems theory that at a bifurcation point, the stability parameter λ, associated with the steady state losing stability, becomes zero [13, 16, 17, 32]. In the case of model 1, one can check that λ = 0 at the two bifurcation points. The return time $T_{R}=\frac{1}{\left|\lambda \right|}$ thus diverges at a bifurcation point. If the dynamical system is described by more than one variable, then the stability of a steady state is determined by the eigenvalues of the Jacobian matrix governing the dynamics of the perturbed system [13]. Let λ_max be the real part of the dominant eigenvalue of the Jacobian matrix (for stability, the real parts of all the λs should be negative). The dominant eigenvalue is the one with the largest real part. The return time T_R is given by $T_{R}=\frac{1}{\left|\lambda _{{\rm max}}\right|}$. Examples of the experimental observations of the CSD include the transition from the G2 growth phase to the mitotic phase of the eukaryotic cell division cycle [31], a direct observation of the CSD in a laboratory population of budding yeast before population collapse occurs at a critical experimental condition [33] and the demonstration of the CSD in a population of cyanobacteria [34].

Zoom In Zoom Out Reset image size

In the presence of noise, the variance σ² is determined from the fluctuation–dissipation (FD) relation (equation (A.5)) in the appendix. In the case of a one-variable model (model 1), the matrix A consists of a single element λ, the stability parameter defined in equation (15). The steady state covariance matrix C reduces to a single element, the variance σ², which is determined as

$$\begin{equation} \sigma ^{2}=\frac{D}{2\mid \lambda \mid } \end{equation} \tag{ 17 }$$

with λ negative. Also, the lag-1 autocorrelation in the stationary state (equations (A.8) and (A.9) with τ = 1) is given by

$$\begin{equation} \rho (1)={\rm e}^{\lambda } . \end{equation} \tag{ 18 }$$

The variance σ² diverges and the lag-1 autocorrelation ρ(1) reaches its maximum value at the bifurcation point since the stability parameter λ is zero at this point. Rising variance and autocorrelation thus provide early signatures of impending regime shifts. The CSD close to the bifurcation point implies that the system's intrinsic rates of change are decreased so that the state of the system at time t resembles closely the state at time t − 1. This increased memory is measured by the lag-1 autocorrelation function. Also, because of a flatter basin of attraction close to the transition point, a given perturbation brings about a greater shift in the system's state, i.e., an increasing variance.

Two other early signatures of a regime shift which are not related to the CSD are skewness and flickering [1, 5]. The skewness γ of a probability distribution P(x) is a dimensionless measure of its degree of asymmetry and is defined as

$$\begin{equation} \gamma =\frac{\int (x-\mu )^{3}P(x)\,{\rm d}x}{\sigma ^{3}}, \end{equation} \tag{ 19 }$$

where μ and σ are respectively the mean and standard deviation of P(x). The variance σ², the coefficient of variation C_v and the third moment M₃ of P(x) are expressed as

$$\begin{eqnarray} \sigma ^{2}&=&\int (x-\mu )^{2}P(x)\,{\rm d}x \nonumber \\ M_{3}&=&\int (x-\mu )^{3}P(x)\,{\rm d}x \nonumber \\\ms C_{v}&=&\frac{\sigma }{\mu } . \end{eqnarray} \tag{ 20 }$$

The skewness of a probability distribution increases as the bifurcation point is approached. This is because of the asymmetry in fluctuations with a system exhibiting greater amplitude deviations in the direction of the state it is fated to switch to than in the other direction. The phenomenon of flickering is observed in the region of bistability with the system switching back and forth between the two attractor states before the bifurcation point is reached. We propose a quantitative measure of flickering which serves as an early signature of regime shift. In the region of bistability, the stability landscape has two minima corresponding to the two stable steady states. Noise-induced transitions take the system from one valley to the other. The mean first passage time (MFPT) refers to the average first exit time from a valley [26, 28]. Let T₁ and T₂ be the MFPTs for the exits respectively from valleys 1 and 2. The times are indicative of the amount of flickering present in the system. The MFPT T₂ becomes zero at the bifurcation point where the steady state 2 loses stability. The ratio $r=\frac{T_{1}}{T_{2}}$ measures the asymmetry in the exit times and diverges at the bifurcation point, at which T₂ → 0 is approached. The quantity r thus provides an early signature of the regime shift.

4.1. Model 1 (one variable)

We first present the results on early signatures in the case of model 1. In the region of bistability, the SSPD P_st(x) is bimodal with two distinct peaks. P_st(x) may be obtained via a solution of the FP equation (equation (12)) or by a numerical solution of the LE (equation (7)). Figure 4 shows an example of a bimodal distribution with two distinct peaks at x = x₁ and x₃. The inset of the figure shows the corresponding stochastic potential profile (equation (14)) with x₂ denoting the location of the hill, and x₁and x₃ denoting the minima of the left and right valleys, respectively. We find the normalized probability distributions, $P_{l}^{n}(x)$ and $P_{u}^{n}(x)$, corresponding to the low and high expression states, respectively, using the following procedure.

• (i)  
  From the stochastic potential profile ϕ_F(x), determine the position of the hill at x₂ and compute P_st(x₂). The point x₂ corresponds to the minimum of the probability distribution between its two peaks at x₁and x₃. The point x₂ is chosen to be the left cut-off point, $x_{L}^{u}$, for $P_{u}^{n}(x)$. The right cut-off point, $x_{R}^{u}$, of $P_{u}^{n}(x)$ is obtained as a solution of the equation

  $$\begin{equation} P_{\rm st}(x)=P_{\rm st}\big(x_{2}^{u}\big),\; x>x_{3} . \end{equation} \tag{ 21 }$$
• (ii)  
  If P_st(0) < P_st(x₂), then the lower cut-off point, $x_{L}^{l}$, of $P_{l}^{n}(x)$ is determined from the solution of the equation

  $$\begin{equation} P_{\rm st}(x)=P_{\rm st}\big(x_{R}^{l}\big),\; x<x_{1}, \end{equation} \tag{ 22 }$$

  where $x_{R}^{l}$, the right cut-off point of $P_{l}^{n}(x)$, is given by $x_{R}^{l}=x_{L}^{u}$. If P_st(0) > P_st(x₂), $x_{L}^{l}=0$.

Other cut-off procedures, e.g., that in [5], have been proposed to isolate the low and high expression probability distributions but the general results on the early signatures are qualitatively similar.

Zoom In Zoom Out Reset image size

In the case of model 1, we first consider only an additive noise in the LE (equation (7), g₁(x) = 0). The additive noise Γ(t) represents noise arising from an external perturbative influence or originating from some missing information because of rate equation approximations [31]. The SSPD, P_st(x), in the region of bistability is obtained from equation (12) putting g₁(x), g₁(x') = 0. Following the prescription already given, one determines the cut-off points of the distribution $P_{l}^{n}(x)$. The normalized distributions are obtained from

$$\begin{equation} P_{l}^{n}(x)=\frac{P_{{\rm st}}^{\prime }(x)}{\int _{x_{L}^{l}}^{x_{R}^{l}}P_{{\rm st}}(x)\,{\rm d}x}\;\; \mbox{for}\;\; x_{L}^{l}<x<x_{R}^{l} \end{equation} \tag{ 23 }$$

$$\begin{equation} P_{u}^{n}(x)=\frac{P_{{\rm st}}^{\prime }(x)}{\int _{x_{L}^{u}}^{x_{R}^{u}}P_{{\rm st}}(x)\,{\rm d}x}\;\; \mbox{for}\;\; x_{L}^{u}<x<x_{R}^{u}. \end{equation} \tag{ 24 }$$

In equations (23) and (24), $P_{{\rm st}}^{^{\prime }}(x)$ is not normalized. With a knowledge of the normalized probability distributions $P_{l}^{n}(x)$ and $P_{u}^{n}(x)$, one can compute the skewness γ, variance σ², third moment M₃ and the coefficient of variation C_v using equations (19) and (20). Figure 5 shows the plots of these quantities (solid lines) versus the bifurcation parameter k_a for the probability distribution $P_{u}^{n}(x)$, i.e., considering the system to be in the high expression state. The sudden regime shift in the deterministic case occurs at the lower bifurcation point k_{a, 1} = 1.334. The parameter values used for the computation are k_d = 1, k_s = 15, k_p = 0.03, J₀ = 0.01 and J₁ = 1. The strength d₂ of the additive noise is d₂ = 0.25. One finds that all the four quantities |γ|, σ², M₃ and C_v increase as the lower bifurcation point is approached thus providing early signatures of a regime shift. The quantities, however, reach their maxima before the deterministic bifurcation point is reached and then start decreasing. The quantities, though providing early signatures, cannot provide knowledge of the bifurcation point. We next include an additional multiplicative noise term in the LE (equation (7)) with

$$\begin{equation} g_{1}(x)=\frac{k_{a}^{\prime }}{k_{a}^{\prime }+k_{d}}. \end{equation} \tag{ 25 }$$

The multiplicative noise is associated with the protein synthesis rate constant J₁ in the active state, i.e., J₁ → J₁ + ε(t). The origin of multiplicative noise lies in the fact that the rate constants are expected to fluctuate in time due to the inherently stochastic nature of gene expression as well as due to stochastic influences such as fluctuations in the number of regulatory molecules and RNA polymerases. With both the additive and multiplicative noise terms present, the SSPDs $P_{l}^{n}(x)$ and $P_{u}^{n}(x)$ are computed following the procedure already described. The dotted curves in figure 5 show the variations of |γ|, σ², M₃ and C_v , associated with the probability distribution $P_{u}^{n}(x)$, as a function of the bifurcation parameter k_a. The parameter values are the same as before with d₁, the strength of the multiplicative noise term, having the value d₁ = 0.25. The cross-correlation coefficient λ_d in equation (8) is taken to be zero. Figure 6 shows how the skewness of $P_{u}^{n}(x)$ changes as a function of the additive noise strength d₂ for the parameter values k_a = 1.5, k_d = 1, k_s = 15, k_p = 0.03, J₀ = 0.01, J₁ = 1 and d₁ = 0 (no multiplicative noise). The value k_a = 1.5 is greater than the value k_{a, 1} = 1.344. The rising skewness is thus a signature of the noise-induced regime shift. In the case of the SSPD $P_{l}^{n}(x)$, one obtains early signatures of the upper bifurcation point k_{a, 2} = 5.176 similar to the ones shown in figures 5 and 6 in the case of $P_{u}^{n}(x)$, though the quantitative measures exhibit less prominent variation. We find that the early signatures of regime shifts are obtained in the cases, (i) where only additive noise is present and (ii) where additive as well as multiplicative types of noise are present.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Considering model 1, we next calculate the stability parameter λ (equation (15)), with f(x) given by the expression on the rhs of equation (4), in the region of bistability. Knowing λ, the return time $T_{R}\big(=\frac{1}{|\lambda |}\big)$ as well as the lag-1 autocorrelation ρ(1) and the variance σ² (equations (17) and (18)) are also determined. Figure 7 shows the variation of λ, T_R, ρ(1) and σ² as a function of the bifurcation parameter k_a and for the parameter values k_d = 1, k_s = 15, k_p = 0.03, J₀ = 0.01, J₁ = 1 and σ_d = 0.25. The stability parameter λ becomes zero at the lower (k_{a, 1} = 1.344) and the upper (k_{a, 2} = 5.176) bifurcation points as expected. The associated return time T_R diverges at the bifurcation points, a characteristic feature of the CSD. The lag-1 autocorrelation ρ(1) reaches the maximum value and the variance σ² diverges at the bifurcation points. These quantities are good indicators of regime shifts and carry distinct signatures of the bifurcation points. The return time T_R is known to satisfy a general scaling law $T_{R}\sim |B-B_{i}|^{-\frac{1}{2}}$ where B and B_i stand for the bifurcation parameter and point respectively [32]. In figure 8, the scaling is demonstrated for the bifurcation parameter k_a and the bifurcation points k_{a, 1} = 1.344 and k_{a, 2} = 5.176. The exponent is very close to $-\frac{1}{2}$ in each case.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 9 shows the variation of r = T₁/T₂ versus the bifurcation parameter k_a. The ratio is seen to diverge at the lower bifurcation point k_{a, 1} = 1.344. The parameter values are k_d = 1, k_s = 15, k_p = 0.03, J₀ = 0.01, J₁ = 1 and d₂ = 0.25 (only additive noise is considered). The time-series data shown in the inset of figure 9 are obtained via numerical solution of the LE (equation (7)) using the algorithm described in [36]. In the limit of large times, the steady state is assumed to be reached. The solid lines in the inset mark the stable expression states $x_{l}^{s}$ and $x_{h}^{s}$ and the dotted line corresponds to the unstable steady state x^u. The values are obtained from a solution of the deterministic rate equation. The MFPTs T₁ and T₂ are computed using the method outlined in [37]. Let us consider a bistable potential with the stable steady states at x₁ and x₃ (x₁ < x₃) which are separated by an unstable steady state at x₂, termed the barrier state (figure 4). The MFPT, T(x; a, b) is the average time of the first exit from the interval (a, b) and satisfies the equation [28, 35]

$$\begin{equation} -1=A(x)\frac{{\rm d}T(x)}{{\rm d}x}+\frac{1}{2}B(x)\frac{{\rm d}^{2}T(x)}{{\rm d}x^{2}}, \end{equation} \tag{ 26 }$$

where A(x) and B(x) appear in the associated FP equation (equation (9)). The MFPT T(x₁)( = T₁) for exit from the basin of attraction of the stable steady state at x₁ is obtained as a solution of equation (26) with the interval (a, b) = (0, x₂) and boundary conditions given by

$$\begin{equation} T^{\prime }(a;a,b)=0\;\; \mbox{and}\;\; T(b;a,b)=0. \end{equation} \tag{ 27 }$$

The prime denotes differentiation with respect to x, with reflecting and absorbing boundary conditions prevailing at a and b, respectively [28, 35]. Following the same procedure, the MFPT T(x₃)( = T₂) for exit from the basin with the stable steady state at x₃ can be calculated from equation (26). The interval now is (a, b) = (x₂, ∞) with x₂ and ∞ serving as absorbing and reflecting boundary points respectively, i.e.,

$$\begin{equation} T(a;a,b)=0\;\;\mbox{and}\;\;T^{\prime }(b;a,b)=0 . \end{equation} \tag{ 28 }$$

Zoom In Zoom Out Reset image size

4.2. Model 2 (two variables)

We next consider model 2, the dynamics of which are governed by the set of two equations (5) and (6). In [25], the phase diagram of the model in the b_k-b_s plane has been determined which has four different regions: (i) monostability with only one fixed point, (ii) bistability with three fixed points, two stable and one unstable, (iii) excitability involving three fixed points only one of which, corresponding to low ComK value, is stable. The competent fixed point (high ComK level) has the characteristic of an unstable spiral and the mid-ComK fixed point is a saddle point, and (iv) only one fixed point exists which is unstable. The system exhibits limit cycle oscillations between the mid-ComK and high-ComK levels. For the purpose of our study, a part of the b_k-b_s phase diagram has been recomputed and shown in figure 10. The diagram shows three different regions: monostable, bistable and excitable. The monostable region is again of two types: monostable low (low ComK level) and monostable high (high ComK level). The boundaries between the regions depicted by solid lines correspond to the saddle-node (SN) bifurcation [13]. Two vertical lines 1 and 2 are drawn in the phase diagram at the points b_{k, 1} = 0.0875 and b_{k, 2} = 0.15, respectively. Line 1 intersects the phase boundaries at the three points b_{s, 1} = 0.7799, b_{s, 2} = 0.7868 and b_{s, 3} = 0.8094. Line 2 has two points of intersection: b_{s, 1} = 0.6175 and b_{s, 2} = 0.7504. Figure 11 shows the plots of λ_max, the real part of the dominant eigenvalue of the Jacobian matrix computed from equations (5) and (6), and the return time $T_{R}=\frac{1}{|\lambda _{{\rm max}}|}$ versus b_s along line 1 ((a) and (b)) and along line 2 ((c) and (d)). For a specific steady state, the Jacobian matrix J has the following structure [13]:

$$\begin{equation} \mathbf {J}=\left( \begin{array}{@{}cc@{}}f_{11} & f_{12}\\ f_{21} & f_{22} \end{array} \right)_{{\rm SS}} . \end{equation} \tag{ 29 }$$

The suffix SS stands for steady state, i.e., the matrix elements are to be computed at the fixed point (x*, y*). A matrix element $f_{ij}=\frac{\partial f_{i}}{\partial x_{j}}$ (i = 1, 2, j = 1, 2), where

$$\begin{equation} f_{1}(x_{1},x_{2})=a_{k}+\frac{b_{k}\,x_{1}^{n}}{k_{0}^{n}+x_{1}^{n}}-\frac{x_{1}}{1+x_{1}+x_{2}} \end{equation} \tag{ 30 }$$

$$\begin{equation} f_{2}(x_{1},x_{2})=\frac{b_{s}}{1+(x_{1}/k_{1})^{p}}-\frac{x_{2}}{1+x_{1}+x_{2}} \end{equation} \tag{ 31 }$$

with x₁ = K, x₂ = S, n = 2 and p = 5. The parameter values are the same as in [25]: a_k = 0.005, k₀ = 0.2 and k₁ = 0.222. The eigenvalue of J are λ₁ and λ₂, and λ_max is the real part of the dominant eigenvalue. In figure 11, one notes that λ_max becomes zero at the SN bifurcation points as expected and the corresponding return time T_R diverges at a bifurcation point. The shaded regions in the figure denote the regions of bistability in which the two stable steady states correspond to the low and high ComK levels, respectively. The level of ComS is anticorrelated with that of ComK. The solid (dotted) lines in figure 11 are associated with the stable steady states representing low (high) ComK levels. Along line 2 and at the bifurcation point b_{s, 1} = 0.6176 (figures 11(c) and (d)), the stable steady state corresponding to the high ComK level loses stability (λ_max = 0, T_R diverges). At the boundary point b_{s, 2} = 0.7504, the state representing the low ComK level loses stability. A steady state is stable if the real parts of both λ₁ and λ₂ are negative. At the point b_{s, 1} = 0.7799 along line 1 (figures 11(a) and (b)), there is no loss of stability of a steady state and hence no CSD with diverging T_R is observed at this point. The low ComK level of the monostable region continues to remain stable as one enters the region of excitability at the point b_{s, 1} along line 1. At the point b_{s, 2} = 0.7868, the system enters the region of bistability from a region of excitability. The high ComK stable steady state loses its stability at this point with λ_max = 0 and a divergent T_R (dotted branch in figures 11(a) and (b)). The low ComK state loses stability at the point b_{s, 3} = 0.8094 when the system passes from a region of bistability to a region of monostable high ComK level. The eigenvalue λ_max = 0 and T_R diverges at this point. The main point to note from the results is that there is no CSD and diverging return time T_R at the transition point b_{s, 1} along line 1 between the regions of monostability and excitability. On the other hand, at the point b_{s, 1} along line 2, separating the regions of monostability and bistability, one of the expression states, namely, the high ComK state, loses stability with a divergent T_R.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We next consider the time evolution of the ComK–ComS system to be stochastic in nature. As in the case of model 1, the FD relation (equation (A.5)) can be used to calculate the stationary state variances of the fluctuations around the steady state and also the lag-1 autocorrelation function ρ₁₁(τ = 1) (equations (A.8) and (A.9)). The Jacobian matrix A has the form shown in equation (29) with the elements calculated from equations (30) and (31) using the relationship $f_{ij}=\frac{\partial f_{i}}{\partial x_{j}}$ (i = 1, 2, j = 1, 2). The stoichiometric matrix is given by

$$\begin{equation} \mathbf {S}=\left( \begin{array}{@{}cccc@{}}1 & -1 & 0 & 0\\ 0 & 0 & 1 & -1 \end{array}\right) . \end{equation} \tag{ 32 }$$

The first and second rows of S correspond to ComK and ComS molecular numbers respectively. The 'elementary complex reactions' [26, 27] considered for the computations are the composite reactions (equations (5) and (6) with x₁ = K and x₂ = S ),

$$\begin{eqnarray} &&X_{1} \mathop{\longrightarrow}^{a_{k}+\frac{b_{k}\,x_{1}^{n}}{k_{0}^{n}+x_{1}^{n}}} X_{1}+1 \nonumber\\ &&X_{1} \mathop{\longrightarrow}^{\frac{x_{1}}{1+x_{1}+x_{2}}}X_{1}-1 \nonumber\\ &&X_{2} \mathop{\longrightarrow}^{\frac{b_{s}}{1+\left(\frac{x_{1}}{k_{1}}\right)^{p}}} X_{2}+1 \nonumber\\ &&X_{2} \mathop{\longrightarrow}^{\frac{x_{2}}{1+x_{1}+x_{2}}}X_{2}-1. \end{eqnarray} \tag{ 33 }$$

The reaction propensity vector (equation (A.2)) is given by

$$\begin{eqnarray} \fl \mathbf {f}(\mathbf {x})= \left( \begin{array}{@{}cccc@{}}a_{k}+\frac{b_{k}x_{1}^{n}}{k_{0}^{n}+x_{1}^{n}} & 0 & 0 & 0 \\ 0 & \frac{x_{1}}{1+x_{1}+x_{2}} & 0 & 0 \\ 0 & 0 & \frac{b_{s}}{1+\left(\frac{x_{1}}{k_{1}}\right)^{p}} & 0 \\\ms 0 & 0 & 0 & \frac{x_{2}}{1+x_{1}+x_{2}} \end{array}\right) . \end{eqnarray} \tag{ 34 }$$

With knowledge of the stoichiometric matrix S and the reaction propensity vector f(x), the diffusion matrix D in the FD relation can be calculated using equation (A.7). Substituting the computed A and D matrices in the FD relation, the variances and the covariances are determined. Similarly, with the help of equations (A.8) and (A.9), the lag-1 autocorrelation function ρ₁₁(τ) = 〈δx₁(t + τ)δx₁(t)〉 with τ = 1 is calculated. Figures 12(a) and (b) exhibit the plots of the variance $\sigma ^{2}=\langle\delta x_{1}^{2}\rangle$, i.e., the variance of the ComK fluctuations, as a function of the bifurcation parameters b_s. Figures 12(c) and (d) exhibit the lag-1 autocorrelation function ρ₁₁(1), associated with the ComK fluctuations, as a function of the bifurcation parameter b_s. The parameter values used in the computations are the same as in the case of figure 11. In the cases of figures 12(a) and (b), the computations are carried out along the line b_{k, 1} = 0.0875 which traverses successively through the regimes of monostability, excitability, and bistability and monostability. In the cases of figures 12(c) and (d), the calculations are carried out along the line b_{k, 2} = 0.15. The shaded regions in figure 12 represent the regions of bistability with the solid (dotted) lines associated with low (high) ComK levels. From the plots one finds that along the line 2 (b_{k, 2} = 0.15), the high ComK level loses stability at the bifurcation point b_{s, 1} = 0.6176 (figures 12(c) and (d)). The sudden regime shift from the high to the low ComK state is signaled by a diverging variance as the bifurcation point is approached and the lag-1 autocorrelation function attaining its maximum value at the point. Similar signatures are obtained at the other bifurcation point b_{s, 2} = 0.7504, where the low ComK state loses stability. In figures 12(a) and (b), at the point b_{s, 1} = 0.7799 along line 1 (b_{k, 1} = 0.0875), as one enters a region of excitability from a region of monostability, the low ComK level of the monostable region continues to be stable. Since no state loses stability at the point, the variance σ² and the lag-1 autocorrelation function ρ₁₁(1) do not provide any signatures of regime change. At the point b_{s, 2} = 0.7868, the system enters the region of bistability from a region of excitability. The high ComK state loses its stability at this point signaled by a diverging variance and with ρ₁₁(1) attaining a maximum at this point. The low ComK level loses stability at the point b_{s, 3} = 0.8094 when the system transits from a region of bistability to a region of monostability with a high ComK level. Experimentally, the entry into a region of excitability/bistability from a region of monostability (low ComK level) is identified by the appearance of a bimodal distribution in the ComK levels as observed in single cell flow-cytometry measurements. In the case of excitability, however, there are no accompanying signatures in the measurable quantities such as variance and lag-1 autocorrelation function. In the case of bistability, sharp rises in σ² and ρ₁₁(1) indicates an impending regime shift.

Zoom In Zoom Out Reset image size

In this paper, we investigate the early signatures of sudden regime shifts occurring at the bifurcation points associated with a fold-bifurcation model. The basic concepts and methodology are applicable to a general class of bifurcation phenomena occurring in diverse systems. Our focus, however, is on obtaining the quantitative estimates of the early signatures of regime shifts in the gene expression dynamics of competence development in B. subtilis. While quantities such as λ_max and T_R are computed with a knowledge of the deterministic rate equations, stochastic formalisms based on the LE, the FP equation and the LNA have been employed to calculate the other quantities. The early indicators of regime shifts, which include the CSD, variance, lag-1 autocorrelation, skewness and the ratio r of MFPTs, are shown to exhibit distinctive variations as a bifurcation point is approached. Although the literature on the early signatures of regime shifts in ecosystems, financial markets and complex diseases is extensive [1–8], the issue has not been systematically addressed in the case of gene expression dynamics, a fundamental activity in the living cell. Our study is the first in this direction and illustrates how the early signatures provide knowledge in advance of an impending regime shift. Some of these signatures may also provide clues on the physical principles underlying gene expression dynamics. As already pointed out, λ_max and T_R do not exhibit any distinctive features (figures 11(a) and (b)) when the system passes from a region of excitability (bimodal distributions in protein levels) to a region of monostability (unimodal distribution of low ComK levels). This is in contrast to the case when the system passes from a region of bistability to one of monostability. The quantity λ_max, associated with the steady state which loses stability, becomes zero at the bifurcation point and the return time T_R diverges (figures 11(c) and (d)). Experiments detecting the CSD, as the bifurcation point is approached, are difficult to carry out. Observation of the CSD in recent experiments [33, 34] provides pointers for carrying out further such experiments. The CSD experiment, if designed in the appropriate manner, would be able to distinguish between the bistability versus excitability paradigm in the case of competence development in B. subtilis. Experiments based on flow-cytometry and time-lapse fluorescence microscopy are easier to carry out and provide estimates of the variance σ² and the lag-1 autocorrelation function ρ₁₁(1) [38, 39]. Similar to figure 11, no distinctive signatures are obtained in terms of σ² and ρ₁₁(1) as one crosses from a region of excitability to that of monostability (figures 12(a) and (b)), whereas prominent signatures indicate the passage from bistability to monostability (figures 12(c) and (d)). Thus, flow-cytometry and time-lapse fluorescence microscopy measurements for a range of values of the bifurcation parameter would be able to distinguish between bistability and excitability as the physical mechanism underlying competence developement.

One outcome of obtaining knowledge of the early signatures is that specific risk aversion strategies may be developed if the sudden transition is to a regime which is undesirable due to various considerations. Examples of such regime shifts include asthma attacks [7], epileptic seizures [8] and the sudden deterioration of complex diseases [18]. One may add the development of persistence of pathogens such as M. tuberculosis in the human lung granulomas to the list. Recent experiments [40, 41] on a sister species M. smegmatis provide evidence that a fraction of the mycobacterial population (the population of persisters) is able to survive under nutrient depletion. To achieve this, the mycobacteria adopt the strategy of generating phenotypic heterogeneity in the form of two distinct subpopulations. The concentration of a key regulatory protein Rel is high in one subpopulation (which develops into the persister subpopulation) and low in the other. In the subpopulation with a high Rel level, the stringent response pathway is initiated which help the mycobacteria to avoid death and adapt to nutrient depletion. The experiments [40, 41] show that the principle underlying the development of phenotypic heterogeneity is based on bistability and noise-induced transitions between the expression states corresponding to the low and high Rel levels. The problem of persisters is that they are not killed by antibiotic drugs and wait for an opportune moment to restart an infection [42]. Early signatures of a regime shift from the normal to the persistent state would help in developing measures preventing the switch to persistence. Similar studies could be carried out on other systems in which bifurcation phenomena are responsible for the development of heterogeneity, choice of cell fate [19] or regime shifts leading to new types of dynamical behavior.

MP acknowledges the support by UGC, India, see Lett. no. F.2-8/2002(SA-I) dated 23.11.2011. SG acknowledges the support by CSIR, India, under grant no. 09/015(0361)/2009-EMR-I.

In this appendix, we describe briefly the linear noise approximation (LNA) to the chemical master equation (CME) [26, 27] and introduce the notation for the relevant quantities. We consider an intracellular biochemical system with volume Ω and N different chemical components. The concentrations of the components are represented in the form of the vector x = (x₁, x₂, ..., x_N)^T, where T denotes the transpose. The chemical constituents take part in R elementary reactions. The state of the system is given by x which changes due to the occurrence of any one of the R reactions. We define the integers S_ij, i = 1, 2, ..., N , j = 1, 2, ..., R to be the elements of a stoichiometric matrix S. The number of molecules of the chemical components i changes from X_i to X_i + S_ij when the jth reaction takes place.

The deterministic dynamics of the system are described by the rate equations

$$\begin{equation} \frac{{\rm d}x_{i}}{{\rm d}t}= {\sum_{j=1}^{R} }\; S_{ij}\,f_{j}(x)\:\:(i=1,2,\ldots ,N). \end{equation} \tag{ A.1 }$$

In a compact notation

$$\begin{equation} \dot{\mathbf {x}}={ \mathbf {Sf(x)}},\:\:\mathbf {f}(\mathbf {x}) =(f_{1}(\mathbf {x}),\ldots ,f_{R}(\mathbf {x}))^{T}, \end{equation} \tag{ A.2 }$$

where f(x) defines the reaction propensity vector. In the steady state, $\dot {\mathbf {x}}=0$ with the state vector x_s determined from the condition f(x_s) = 0. Let δx denote a weak perturbation applied to the steady state, i.e., the new state vector x = x_s + δx. On Taylor expansion of the rate vector $\frac{{\rm d}\mathbf {x}}{{\rm d}t}$ about the steady state and retaining only terms linear in δx, one obtains

$$\begin{equation} \frac{{\rm d}}{{\rm d}t}\delta \mathbf {x}=\mathbf {A}\delta \mathbf {x}, \end{equation} \tag{ A.3 }$$

where A is the Jacobian matrix the elements of which are given by

$$\begin{equation} A_{ij}= {\sum _{{k=1}}^{R}} {S_{ik}} \:\:\frac{\partial f_{k}}{\partial x_{j}} . \end{equation} \tag{ A.4 }$$

The state x_s is stable if all the eigenvalues of A have negative real parts. A deterministic dynamical model, as described above, is appropriate for describing the dynamics of a system when the number of molecules, X_i  (i = 1, 2, ..., N) is large. In reality, the biomolecules participating in cellular reactions are mostly small in number so that a stochastic description of the dynamics is more valid.

The CME describes the rate of change of the probability distribution P(X₁, X₂, ..., X_N, t) of the numbers of the different chemical components. The CME is not exactly solvable in most cases and one has to take recourse to various approximate methods in order to solve the equation. The LNA provides an approximation to the CME via a large volume (Ω) expansion around the macroscopic steady state. Noise in the form of fluctuations around the steady state is expected to be small in the large Ω limit as the number of molecules scales with the volume. To the first order in the expansion, one obtains the set of deterministic rate equations, whereas in the second order, one obtains the linear FPE describing fluctuations about the steady state. The stationary solution of the linear FPE is given by a multivariate Gaussian distribution. One can further show that the covariance of the fluctuations about the deterministic steady state is given by the fluctuation–dissipation (FD) relation

$$\begin{equation} \mathbf {A}\mathbf {C}+(\mathbf {A}\mathbf {C})^{T}+\mathbf {D}=0, \end{equation} \tag{ A.5 }$$

where A is the Jacobian matrix (equation (A.4)), C = 〈δx δx^T〉 is the covariance matrix, the diagonal elements of which are the variances, and D is the diffusion matrix. The elements of the covariance matrix are

$$\begin{equation} C_{ij}=\langle \delta x_{i}\:\delta xj\rangle \;\;\;\;(i,j=1,2,\ldots ,N). \end{equation} \tag{ A.6 }$$

The matrix D has the form

$$\begin{equation} \mathbf {D}=\mathbf {S}\; \mbox{diag}(\mathbf {f}(\mathbf {x}))\;\mathbf {S}^{T}, \end{equation} \tag{ A.7 }$$

where diag(f(x)) is a diagonal matrix with the elements f_j(x), j = 1, 2, ..., R. The matrices A and D are computed at the stationary state x = x_s. Also, 〈δx_i〉 = 0, i = 1, 2, ..., N. Once A and D are determined, one can determine the elements of the covariance matrix (especially, the variances) from the FD relation (A.5). The time correlation matrix for δx is given by

$$\begin{equation*} \langle \delta \mathbf {x}(t+\tau )(\delta \mathbf {x}(t))^{T}\rangle =\exp (\mathbf {A}\tau )\; \mathbf {C}, \end{equation*}$$

with

$$\begin{equation} \langle \delta x_{i}(t+\tau )\delta x{}_{j}(t)\rangle ={\sum_{n} }{\sum _{m}} {\rm e}^{\lambda _{m}\tau }(T^{-1})_{mn}T_{im}C_{nj}. \end{equation} \tag{ A.8 }$$

In equation (A.8), λ_ms are the eigenvalues of the Jacobian matrix A and T is a matrix the columns of which represent the right eigenvectors of A.

The diagonal elements of the matrix are the autocovariances and τ defines the lag time. The lag-τ autocorrelation for the ith chemical component is

$$\begin{equation} \rho (\tau )=\frac{\langle \delta x_{i}(t+\tau )\delta x{}_{j}(t)\rangle }{\sqrt{var\; x_{i}(t+\tau )}\;\sqrt{var\; x_{i}(t)}} . \end{equation} \tag{ A.9 }$$