Modeling stochastic phenotype switching and bet-hedging in bacteria: stochastic nonlinear dynamics and critical state identification

Abstract

Fluctuating environments pose tremendous challenges to bacterial populations. It is observed in numerous bacterial species that individual cells can stochastically switch among multiple phenotypes for the population to survive in rapidly changing environments. This kind of phenotypic heterogeneity with stochastic phenotype switching is generally understood to be an adaptive bet-hedging strategy. Mathematical models are essential to gain a deeper insight into the principle behind bet-hedging and the pattern behind experimental data. Traditional deterministic models cannot provide a correct description of stochastic phenotype switching and bet-hedging, and traditional Markov chain models at the cellular level fail to explain their underlying molecular mechanisms. In this paper, we propose a nonlinear stochastic model of multistable bacterial systems at the molecular level. It turns out that our model not only provides a clear description of stochastic phenotype switching and bet-hedging within isogenic bacterial populations, but also provides a deeper insight into the analysis of multidimensional experimental data. Moreover, we use some deep mathematical theories to show that our stochastic model and traditional Markov chain models are essentially consistent and reflect the dynamic behavior of the bacterial system at two different time scales. In addition, we provide a quantitative characterization of the critical state of multistable bacterial systems and develop an effective data-driven method to identify the critical state without resorting to specific mathematical models.

Graphical abstract

Keywords

phenotypic heterogeneity / phenotypic variation / multistability / gene network / stochastic gene expression

INTRODUCTION

Bacteria in the wild exist in ever-changing environments and have to surmount the challenges posed by environmental fluctuations. Numerous experiments have confirmed that multiple distinct phenotypes can coexist within an isogenic bacterial population [ 1- 10]. This phenotypic heterogeneity in genetically identical cells has received increasing attention in recent years since it could help the bacteria to survive in rapidly changing environments. In recent years, there have been intense studies on phenotypic heterogeneity in response to changing environments [ 11, 12]. In the framework of traditional population genetics, a bacterial population enhances its fitness via genetic changes caused by mutation or recombination. However, extracellular conditions can change so rapidly that adaptation only by gene mutation or recombination would be too slow. One solution to this problem is to allow individual cells to stochastically switch among multiple phenotypes without genetic changes, a phenomenon widely known as stochastic phenotype switching [ 13- 18]. Generally, the multiple phenotypes within an isogenic bacterial population result from the multiple steady-state expression levels of a group of stress-related genes. Such kind of gene expression pattern with multiple steady-state expression levels are widely known as multistability [ 3, 7].

Phenotypic heterogeneity is a widespread phenomenon in the bacterial realm. Examples of phenotypic heterogeneity include the lactose utilization in Escherichia coli [ 19], the lysis-lysogeny switch of bacteriophage, competence development in Bacillus subtilis [ 20- 22], sporulation in Bacillus subtilis [ 23- 26], and persistence in Mycobacterium tuberculosis [ 27- 29]. The potential function of phenotypic heterogeneity with stochastic phenotype switching is generally understood to be a bet-hedging strategy [ 7, 12, 30, 31], a term originating from finance. In response to fluctuating environments, a heterogeneous bacterial population could optimize its fitness by altering the proportion of cells in each subpopulation via stochastic phenotype switching to achieve an optimal ‘investment portfolio’.

To study the evolution of heterogeneous bacterial populations, a number of Markov chain models have been proposed at the cellular level [ 1, 6, 14, 16, 17, 32- 34]. These models assumed a priori that the bacterial population has multiple distinct phenotypes. In these Markov chain models, each phenotype is modeled as a state and stochastic phenotypic switching is modeled as the state transition of the Markov chain. However, these models take phenotypic heterogeneity and stochastic phenotypic switching for granted and fail to account for their underlying molecular mechanisms.

Recent research has demonstrated that phenotypic heterogeneity within isogenic bacterial populations often results from the feedback circuitry of the gene regulatory network [ 35, 36]. To account for the molecular mechanism of phenotypic heterogeneity, a number of deterministic models have been proposed at the molecular level [ 19, 22, 27, 29, 37]. In these models, different steady-states of gene expression are described as different stable fixed points (attractors) of a deterministic system composed of several ordinary differential equations which are written down based on the regulatory relationships behind the gene network. However, deterministic models cannot provide a correct description of many important experimental phenomena, such as stochastic phenotype switching and bet-hedging. In every deterministic model, if the expression level of an individual cell lies in an attraction basin at a particular time, it will never leave this attraction basin and thus phenotype switching will never occur.

Although deterministic models can give rise to multiple attractors and attraction basins, they do not allow transitions among different attraction basins. One solution to this problem is to consider stochastic effects, which allow the system to switch among different attraction basins and thus drive stochastic phenotype switching. This fact is analogous to the simulated annealing techniques in optimization problems, in which noise is indispensable to make the search escape from the trap of local minimum points and reach the global minimum point. To better understand the role that stochastic effects play in bistable systems, Qian and coworkers studied the relations between deterministic and stochastic nonlinear dynamics in great detail [ 38- 43]. However, their models are somehow so abstract and oversimplified that they cannot be directly applied to practical problems with experimental data and observations.

Stochastic effects are extremely important not simply because they are indispensable for the model to generate phenotype switching, but because gene expression is an inherently stochastic process. Recent developments of single-cell and single-molecule experiments have shown that many important cellular processes, such as transcription, translation, replication, and gene regulation, are inherently stochastic [ 44- 53]. Due to the stochastic effects, the expression levels of the stress-related genes in a multistable system will have a multimodal distribution.

In this paper, we propose a unified nonlinear stochastic model of multistable bacterial systems at the molecular level based on a core double-positive-feedback gene network. By studying its stochastic nonlinear dynamics, we show that our model not only provides a clear description of phenotypic heterogeneity, stochastic phenotype switching, and bet-hedging within isogenic bacterial populations, but also provides a deeper insight into the analysis of multidimensional experimental data, such as gene expression data and the data of more comprehensive indicators like the forward scatter (FSC) and the side scatter (SSC) measured by flow cytometry.

Next, we use the mathematical tool of large deviation theory established by Freidlin and Wentzell [ 54] to show that every multistable dynamical system under a small random perturbation can be approximated by a Markov chain with multiple states, each corresponding to an attraction basin of the multistable system. In this way, our stochastic model at the molecular level can be reduced to a Markov chain model at the cellular level. This justifies the wide applications of previous Markov chain models of population evolution, in particular, the Markov chain model proposed by Lander and coworkers [ 33] about the dynamics of the phenotypic proportions in human breast cancer cell lines.

In addition, we point out a widespread misunderstanding on the analysis of gene expression data, inspired by our recent work about antibiotic resistance in Escherichia coli. Previous work tended to think that phenotypic heterogeneity can be identified by the multistable expression of a single pivotal gene (reviewed in [ 2, 7]). However, phenotypic heterogeneity in bacterial populations often results from the interaction of a group of stress-related genes. We use simulation results to show that in many cases, the expression data of a group of genes give rise to an apparent multimodal distribution, however, we cannot observe the multistable expression if we only focus on the expression data of a single gene. This suggests the traditional method to identify phenotypic heterogeneity by measuring the expression of a single pivotal gene is sometimes ineffective.

Finally, we use our stochastic model to provide an answer to the important question of identifying the critical state of multistable bacterial systems. In our stochastic model, there is a saddle lying on the boundary of two adjacent attraction basins which characterizes a critical state between two steady-states of gene expression. The critical state is not targeted in previous work since it is rarely observed in experiments and cannot be estimated by simple statistical analysis of gene expression data. However, the identification of the critical state has drawn increasing attention in recent years since it is closely related to the early diagnosis of complex diseases [ 55, 56]. In this paper, we develop an effective method to identify the critical state of multistable bacterial systems by using the time-course data of gene expression without resorting to specific mathematical models.

METHODS

The parameters used to draw the figures in the main text are chosen as α = 1, β = 2, γ = 0.8, K = 0.05, δ = 0, and n = 6. The inducer concentration a and the two noise levels, ϵ and η, are chosen appropriately when drawing different figures.

RESULTS

Model

In natural bacterial systems, phenotypic heterogeneity and stochastic phenotype switching always originate from the feedback circuitry of the regulatory network which governs a group of stress-related genes. To better understand the general principles behind phenotypic heterogeneity, we illustrate the gene regulatory networks that govern some best-understood multistable systems in bacteria (Figure. 1A-D). A crucial similarity shared by these examples is that the wiring of the gene regulatory network forms a double-positive-feedback loop. To establish a unified model of these bacterial systems, we focus on the core double-positive-feedback gene network depicted in Figure. 1E, where, protein X is the product of a pivotal stress-related gene, i.e., gene X, protein Y is a transcription factor which activates the expression of gene X, and A is an inducer whose concentration reflects the fluctuations in extracellular environmental conditions, such as fluctuations in temperature, pH, and concentrations of nutrients and toxins [ 23, 27].

We use lowercase letters x, y, and a to denote the concentrations of X, Y, and A, respectively. Since gene expression is an inherently stochastic process, the dynamics of x and y can be described by the following two-dimensional system of stochastic differential equations (SDE):

(1) \{\begin{array}{l}\dot{x}=-\alpha (x-F(a,y))+\sqrt{\mathrm{2}ϵ}{\xi }_x,\hfill \\ \dot{y}=-\beta (y-G(x))+\sqrt{\mathrm{2}\eta }{\xi }_y,\hfill \end{array}

where F(a, y) describes the activation of protein X by inducer A and protein Y, G(x) describes the activation of protein Y by protein X, and α and β are two parameters characterizing the response speeds of proteins X and Y, respectively. In addition, {\xi }_x and {\xi }_y are two independent standard white noises satisfying \langle {\xi }_x(t)\rangle =\langle {\xi }_y(t)\rangle =\mathrm{0} and \langle {\xi }_x(t){\xi }_x(t^{\prime })\rangle =\langle {\xi }_y(t){\xi }_y(t^{\prime })\rangle =\delta (t^{\prime }-t) . Since the fluctuations in the levels of proteins X and Y can be different, we use two noise levels ϵ and η to describe their stochastic fluctuations. We emphasize here that the noise in gene regulatory networks generally comes from a great number of sources and may not be subsumed into white noises. Fortunately, the specific noise distributions will hardly affect the main results of this paper. To make our discussion friendly to both theoretical and experimental biologists, we would like to use white noises to describe noise in gene expression in this paper. For more reasonability about this assumption, please see Discussion.

If we ignore the stochastic effects, then the stochastic system (1) reduces to the following deterministic system as the two noise levels, ϵ and η, tend to zero:

(2) \dot{x}=-\alpha (x-F(a,y)),\dot{y}=-\beta (y-G(x)).

We note that the fixed points of this deterministic system are the solutions to the following equations: x-F(a,G(x))=\mathrm{0} and y=G(x) . In natural bacterial systems, the most common expression of F(a,G(x)) has the following form (see Supplementary Information): F(a,G(x))={\displaystyle \frac{\gamma x^n}{K+x^n}}+(\mu a+\delta ), where the Hill function \gamma x^n/(K+x^n) with n > 1 describes the activation of gene X by protein Y, the term μa describes the activation of protein X by inducer A, and the term δ describes a basal expression of gene X independent of the activation of protein Y.

It turns out that the deterministic system (2) has one or three fixed points under different inducer concentrations (Figure. 2A). To be specific, the inducer concentration a has two threshold levels, a₀ and a₁. If a < a₀ or a > a₁, the system has only one fixed point (Figure. 2A). If a < a₀, the only attractor (x_L, y_L) describes the phenotype of low-expressing cells in which gene X is inactivated. If a > a₁, the only attractor (x_H, y_H) describes the phenotype of high- expressing cells in which gene X is activated. If a₀ < a < a₁, however, the system has three fixed points, including two attractors, (x_L, y_L) and (x_H, y_H), and a saddle (x_M; y_M) (Figure. 2A, B), where the two attractors describe the phenotypes of low- and high-expressing cells, respectively, whereas the saddle is a critical state between the two steady-states of gene expression. Mathematically, each attractor of a deterministic system has an attraction basin, and two adjacent attraction basins are separated by a boundary (Figure. 2B).

Phenotypic heterogeneity and bet-hedging

In recent decade, single-cell and single-molecule experiments have made significant progresses and shown that gene expression is an inherently stochastic process. Although deterministic models, such as the deterministic system (2), can give rise to multiple attractors and attraction basins, they cannot provide a correct description of stochastic phenotype switching and bet-hedging. These two facts show that the reduction from the stochastic system (1) to the deterministic system (2) is inappropriate.

In the following discussion, we focus on the stochastic system (1) in which the two noise levels, ϵ and η, are strictly positive. Although the system does not satisfy detailed balance, we can still obtain an approximate steady-state probability distribution p_s(a, x, y) of the system (see Supplementary Information), which has the following form:

(4) p_s(a,x,y)={\displaystyle \frac{\mathrm{1}}{Z}}\exp \left\{-{\displaystyle \frac{\mathrm{1}}{ϵ}}U(a,x,y)\right\},

where Z is a normalization constant and U(a, x, y) is an approximate global potential, also called landscape, of the system defined as

(5) U(a,x,y)={\displaystyle \frac{ϵ\beta }{\mathrm{2}\eta }}{(y-G(x))}^{\mathrm{2}}+\alpha {\displaystyle {\int }_{\mathrm{0}}^x(u-F(a,G(u)))du.}

We make a crucial observation that the fixed points of the deterministic system (2) are exactly the solutions to the equation \partial xU(a,x,y)={\partial }_{y}U(a,x,y)=\mathrm{0} . This shows that the attractors of the deterministic system are the local minimum points of the approximate global potential U(a, x, y) and thus are the local maximum points of the steady-state distribution p_s(a, x, y). From (Figure. 2C-E), we see that the steady-state distribution of the levels of proteins X and Y is controlled by the inducer concentration a. If a < a₀, the steady-state distribution has only the left peak, suggesting that the bacterial population contains almost exclusively low-expressing cells under favorable conditions (Figure. 2C). If a₀ < a < a₁, the steady-state distribution has both the left and right peaks, each corresponding to a phenotype. With the increase of the inducer concentration, a larger fraction of cells will switch from the low- to the high-expressing subpopulation to maximize survival (Figure. 2D)). If a > a₁, the left peak of the steady-state distribution disappears, suggesting that the bacterial population contains almost almost exclusively high-expressing cells under unfavorable conditions (Figure. 2E). The above discussion clearly shows that how the bet-hedging strategy could help the bacterial population better adapt to rapidly changing environmental conditions.

Generally, the steady-state gene expression levels in a multistable bacterial system have a multimodal distribution, which can be viewed as the superposition of multiple monomodal distributions, each concentrated within an attraction basin. The attractors are the locally mostprobable states and thus are most likely to be observed in experiments. Based on the stochastic system (1), we simulate the time course of the levels of proteins X and Y in a single cell (Figure. 2B). The simulation result shows that the gene expression data are generally distributed around the attractors and are rarely distributed around the boundaries of the attraction basins. These facts clearly show that each phenotype of a bacterial population cannot be simply described as an attractor of the deterministic model, but should be understood as a monomodal distribution concentrated within an attraction basin.

A widespread misunderstanding on the analysis of gene expression data

Phenotypic heterogeneity in isogenic bacterial populations often results from the interaction of a group of stress-related genes and biochemical species. Previous work tended to think that phenotypic heterogeneity can be identified by the multistable expression of a single pivotal gene (reviewed in [ 2, 7]). In experiments, however, it often happens that the histogram of the steady-state expression data of a single gene does not display a multimodal distribution, and only when a subpopulation of cells are sorted out to start from some extreme initial conditions, the multimodal distribution can be observed at certain times before reaching the steady-state. Thus it is rather difficult to decide whether the bacterial population has multiple phenotypes or not. To explain these experimental phenomena, we point out that phenotypic heterogeneity often results from the interaction of a group of stress-related genes and may not be observed if we only focus on the expression data of a single gene. In our recent study about antibiotic resistance in Escherichia coli (unpublished work), we found that the expression data of the hydrolase gene only lead to a monomodal distribution, whereas the expression data of a group of stress-related genes lead to an apparent multimodal distribution.

We now use our stochastic model to account for this interesting phenomenon. Based on the stochastic system (1), we simulate the steady-state levels of proteins X and Y in 500,000 virtual cells under a set of model parameters (Figure. 2F). From the simulation result, we see that the two-dimensional gene expression data are distributed around two attractors in the phase plane and lead to a bimodal distribution, which can be viewed as the superposition of two monomodal distributions. Although these two monomodal distributions are concentrated within two different attraction basins in the phase plane, there is an obvious overlap between their marginal distributions, whose superposition, which represents the steady-state distribution of the level of protein Y, has only one peak (Figure. 2G). This suggests that the traditional idea to identify phenotypic heterogeneity by the multistable expression of a single pivotal gene is sometimes ineffective. The above discussion also shows that our stochastic model can help us gain a deeper insight into the pattern behind multidimensional experimental data.

From the molecular level to the cellular level

We have seen that the stochastic model (1) provides a clear description of phenotypic heterogeneity and bet-hedging within isogenic bacterial populations at the molecular level. However, more widely used models in early work are Markov chain models at the cellular level. These models assume a priori that the bacterial population has multiple distinct phenotypes, each of which corresponds to a state of the Markov chain and can switch to other phenotypes with certain transition rates. This raises the question of whether the two kinds of models at two different levels, our stochastic model at the molecular level and the Markov chain models at the cellular level, are consistent in some way or not. We now use the mathematical tool of the large deviation theory established by Freidlin and Wentzell [ 54] to answer this question.

The Freidlin-Wentzell theory is mainly concerned about the dynamic behavior of a multistable dynamical system under random perturbations, such as the stochastic system (1), when the noise level is not too large. The conclusions of the Freidlin-Wentzell theory are not very intuitive at first sight and the proofs of them are rather tedious. Readers who are interested in the mathematical aspects of the Freidlin-Wentzell theory may refer to [ 54, 57]. To make readers understand this useful mathematical tool, we would like to list the major results of the Freidlin-Wentzell theory as follows.

Basic Result 1. No matter how small the noise level is, the accumulation of the stochastic forces will make the system escape from the trap of an attraction basin and enter another attraction basin. Before the system escapes from an attraction basin, it will spend most of the time staying around the attractor and spend little time staying around the boundary of the attraction basin. These facts can be seen from our numerical simulation in Figure. 2B and Figure. 3A.

Basic Result 2. Each point x in an attraction basin has a local potential V(x) called the quasi-potential. If the system has a global potential U(x), as in Equation. (5), then the quasi-potential V(x) can be calculated explicitly as

(6) V(x)=\mathrm{2}(U(x)-U(x_{\mathrm{0}})),

where x₀ is the attractor in the attraction basin. When the system escapes from one attraction basin to another, it must cross the boundary around a specific point y₀ where the quasi-potential attains its minimum. In most cases, the minimum point y₀ of the quasi-potential on the boundary is exactly the saddle of the system. These facts can be seen from our numerical simulation in Figure. 2B.

Basic Result 3. The time needed for the system to escape from an attraction basin is referred to as the escape time. The escape time T from an attraction basin approximately follows an exponential distribution. The mean escape time \langle T\rangle , which is approximately the time constant of the exponential distribution, has the form of

(7) \langle T\rangle \doteq k\exp \left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}ϵ}}{V}_{\mathrm{0}}\right),

where k is a positive constant, ϵ is the noise level, and V₀ = V (y₀) is the minimum value of the quasi-potential on the boundary. By Basic Result 2, if the system has a global potential U(x), then V₀ / 2 = U (y₀) - U (x₀), which represents the potential difference between the minimum point y₀ of the potential on the boundary and the attractor x₀. These facts can be seen from our mathematical derivations in Supplementary Information.

It is a well-known result that the time needed for a Markov chain to make a state transition follows the exponential distribution. The Freidlin-Wentzell theory tells us that the escape time from each attraction basin approximately follows an exponential distribution. This shows that if we combine each attraction basin as a state, then the stochastic system with multiple attractors can be approximated by a Markov chain with multiple states at the time scale of exp(1 / ϵ). When the noise level ϵ is small, exp(1 / ϵ) becomes very large. This shows that the approximate Markov chain reflects the long-term dynamic behavior of the stochastic system.

According to the above discussion, each stochastic model of a multistable system at the molecular level can be reduced to a Markov chain model at the cellular level. These two kinds of models at two different levels are essentially consistent and reflect the dynamic behavior of the system at two different time scales. For instance, the stochastic model (1) proposed in this paper has two attraction basins when a₀ < a < a₁, and thus can be approximated by a Markov chain model with two states, each corresponding to a phenotype. If the feedback architecture of the gene network becomes more complicated, then the stochastic system may display three or more attraction basins and thus can be approximated by a Markov chain model with three or

more states. For examples of Markov chain models with two, three, or four states, readers may refer to [ 16, 17, 33].

Stochastic phenotype switching

To survive in rapidly changing environments, a heterogenous bacterial population may allow individual cells to stochastically switch among multiple phenotypes, ensuring that some cells are always prepared for an unforeseen environmental fluctuation. This kind of phenotype switching is stochastic and temporary: An individual cell may switch to an alternative state at a random time and switch back again after some random time. Even without a significant change in environmental conditions, stochastic phenotype switching still exists. Stochastic phenotype switching has been observed in a wide range of bacterial species. As an example, upon encountering nutrient limitation, a minority of Bacillus subtilis cells transiently enter the competent state with the capability for DNA uptake from the environment before returning to vegetative growth [ 20].

To better understand the principle behind stochastic phenotype switching, we simulate the time course of the level of protein X in an individual cell based on the stochastic system (1) (Figure. 3A). The simulation result shows that the cell switches between the low- and high-expression states at certain random times. The Freidlin-Wentzell theory shows that the escape times T_L and T_H from the low- and high-expression states approximately follow exponential distributions. The mean escape times, \langle T_L\rangle and \langle T_H\rangle , have the following form (see Supplementary Information and Basic Result 3):

(8) \begin{array}{l}\langle T_L\rangle \doteq {\displaystyle \frac{\mathrm{2}\pi }{\sqrt{{\kappa }_L{\kappa }_M}}}\exp \left({\displaystyle \frac{\mathrm{1}}{ϵ}}\mathrm{\Delta }{U}_L\right),\hfill \\ \langle T_H\rangle \doteq {\displaystyle \frac{\mathrm{2}\pi }{\sqrt{{\kappa }_H{\kappa }_M}}}\exp \left({\displaystyle \frac{\mathrm{1}}{ϵ}}\mathrm{\Delta }{U}_H\right),\hfill \end{array}

where \mathrm{\Delta }{U}_L and \mathrm{\Delta }{U}_H are the potential differences between the saddle and the two attractors (Figure. 3C) and {\kappa }_L , {\kappa }_M , and {\kappa }_H are the curvatures of the one-dimensional effective potential U(a, x,G(x)) at x_L, x_M, and x_H, respectively. At the time scale of exp (1 / ϵ), the dynamic behavior of the stochastic system (1) can be approximated by a Markov chain model with two states (Figure. 3B).

According to Equation. (8), the mean escape time is an exponential function of the potential barrier. The higher the potential barrier, the longer time is needed for a cell to make a state transition. Equation (8) also shows that the mean escape times have the time scale of exp (1 / ϵ). This suggests that stochastic phenotype switching is a long-term dynamic behavior of the system. When the noise level ϵ is small, the escape time may be longer than the time of cell division. For a fraction of cells, stochastic phenotype switching may not occur in a single cell cycle, and thus they can pass their phenotypic state to the next generation.

If there were no stochastic effects, the phenotype of each individual cell would never change. Specifically, if the low-expressing cells are sorted out at a particular time, then all cells will stay in the low-expression state forever and phenotype switching is impossible (Figure. 3C). If the stochastic effects cannot be ignored, however, the accumulation of stochastic forces will drive individual cells to surmount the potential barrier and make a state transition (Figure. 3D). This suggests that stochasticity in gene expression is the driving force to generate stochastic phenotype switching.

Importance of the critical state

We have seen that at certain ranges of the inducer concentration, the two attractors of the stochastic system (1) are separated by a boundary, forming two attraction basins. The saddle of the system lies exactly on the boundary of the two attraction basins (Figure. 2B) and thus characterizes a critical state between the two steady-states of gene expression. This saddle is not targeted in early work since it is rarely observed in experiments and cannot be estimated by simple statistical analysis of gene expression data. However, the identification of the critical state has drawn increasing attention in recent years due to the following three reasons.

First, the saddle represents the critical level of gene expression. Recent studies on complex diseases show that any disease progression can be divided into a normal state, a pre-disease state, and a disease state [ 55, 56], similar to the low-expression state, the critical state, and the high-expression state described in this paper. Once the expression levels of the disease-related genes in a person is close to the saddle, we have good reasons to believe that this person is in a pre-disease state and is at high risk of disease progression. This suggests that the identification of the critical state is closely related to the early diagnosis of complex diseases.

Second, the saddle is the most important point on the boundary of two adjacent attraction basins. According to the Freidlin-Wentzell theory, when the system escapes from one attraction basin to another, it must cross the boundary around a specific point where the quasi-potential attains its minimum. To see this, we simulate the time course of the levels of proteins X and Y in an individual cell based on the stochastic model (1) (Figure. 2B). The simulation result shows that the protein levels of the cell stay around the attractors at most times and cross the boundary of the two attraction basins around the saddle, which is exactly the minimum point of the potential U(a, x, y) on the boundary.

Third, the saddle characterizes a critical state of the transition between multiple attraction basins. To accomplish stochastic phenotype switching, the system needs first to climb up the potential from one attractor to the saddle, and then to fall down the potential from the saddle to another attractor. Before reaching the saddle, the accumulation of stochastic forces will drive the system to climb up the potential against the potential gradient. This process in general will take rather a long time. Once the system crosses the saddle, it will reach another attractor along the potential gradient in a very short time. This shows that the dynamic features of a multistable system before and after reaching the saddle are totally different. To be more precise, let T_u denote the time needed for the system to climb up the potential from one attractor to the saddle and let T_d denote the time needed for the system to fall down the potential from the saddle to another attractor. The Freidlin-Wentzell theory shows that the ratio of T_u to T_d has the time scale of exp(1 / ϵ). This suggests that the process of climbing up the potential is much longer than that of falling down the potential, which is consistent with the old saying: Diseases come on horseback, but go away on foot.

Identification of the critical state

The critical state of a multistable system is important in many ways. Thus it is natural to ask whether we can identify the critical state in an effective way according to the noisy data of gene expression. Recently, Chen et al. [ 58] have developed a method of identifying the leading network in complex diseases by evaluating a kind of network entropy and Gore et al. [ 59] have used generic statistical indicators to provide early warning signals for catastrophic collapse in the budding yeast population. Inspired by their ideas, in this section, we shall develop an effective method to identify the critical state of a multistable system by using the time-course data of gene expression. Moreover, we shall validate the effectiveness of our method through both theoretical derivations and numerical simulations.

We have seen that at certain ranges of the inducer concentration, the pivotal gene, gene X, has two steady-state expression levels, x_L and x_H, and a critical expression level x_M. We assume that we have measured the level of protein X in each individual cell within a bacterial population at two times, t and t + h, where the interval h of two successive measurements is chosen to have the time scale of 1 / ϵ, which is much shorter than the time scale exp (1 / ϵ) of stochastic phenotype switching. Intuitively, if the protein level in an individual cell is around x_L or x_H at time t, then the protein level at time t + h should be also around x_L or x_H since the interval h is much shorter than the time needed for stochastic phenotype switching (Figure. 4A). However, if the protein level in an individual cell is around x_M at time t, then the protein level at time t + h will become rather scattered since the critical state is rather unstable (Figure. 4A).

Let x(t) denote the level of protein X at time t, whose value can differ significantly between two individual cells. The above discussion illuminates us to define the variance D(x) at x as the variance of x(t + h) conditioned on the information of x(t) = x. More precisely, we define the variance D(x) at x as

(9) D(x)=\mathrm{V}\mathrm{a}\mathrm{r}\mathrm{i}\mathrm{a}\mathrm{n}\mathrm{c}\mathrm{e}(x(t+h)|x(t)=x),

where Variance (Z)=\langle {(Z-\langle Z\rangle )}^{\mathrm{2}}\rangle is the variance of the random variable Z. According to the above discussion, the variance D(x) around x_L or x_H should be small since the distribution of x(t + h) will be rather concentrated if x(t) is around x_L or x_H, whereas the variance D(x) around x_M should be large since the distribution of x(t + h) will be rather scattered if x(t) is around x_M. This suggests that we may detect the critical state by seeking the maximum point of the variance function D(x).

We next validate the above intuitive discussion from the theoretical point of view. In fact, the theoretical expression of the variance function D(x) has the following form (see Supplementary Information):

(10) D(x)\doteq {(x_H-x_L)}^{\mathrm{2}}{p}_L(x)p_H(x)+ϵ\left({\displaystyle \frac{p_L(x)}{{\kappa }_L}}+{\displaystyle \frac{p_H(x)}{{\kappa }_H}}\right),

where p_L(x) and p_H(x) are two functions (see Supplementary Information for specific expressions) satisfying p_L(x) + p_H(x) = 1. We denote the maximum point of the variance function D(x) by x_max. By Equation (10), we can easily see that x_max satisfie

(11) p_H(x_{max})={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}+{\displaystyle \frac{({\kappa }_L-{\kappa }_H)ϵ}{\mathrm{2}{(x_H-x_L)}^{\mathrm{2}}{\kappa }_L{\kappa }_H}}.

In order to find the location of the maximum point x_max, we depict the graph of the function p_H(x) in Figure. 4B, from which we see that the function p_H(x) is sigmoidal with a critical transition around x_M. With the decrease of the noise level ϵ, the slope of the function p_H(x) at x_M tends to infinity. This fact, together with Equation (11), indicates that when the noise level ϵ is small, the maximum point x_max of the variance function D(x) is very close to x_M. This suggests that the variance function D(x) may provide a clear signal for the position of critical state.

A natural and important question is that whether we can estimate the variance function D(x) according to the noisy data of gene expression. The answer is of course affirmative. Recent experimental techniques such as fluorescent labeling and microfluidic devices allow us to measure the expression level of the pivotal gene in each individual cell within a bacterial population at a series of times t₁, …, t_m with interval h. We assume that the experimental data of the n-th cell have the form of x(n, t₁), …, x(n, t_m), where x(n, t_i) represents the expression level of the n-th cell at time t_i. For each time t_i, if x(n, t_i) is in a given small neighborhood of x, we pick out its next datum x(n, t_i + 1), whereas if x(n, t_i) is not in the given small neighborhood of x, we throw away its next datum. By evaluating the sample variance of those data which are picked out, we can obtain a good estimation of the variance D(x) at x.

To validate the effectiveness of our method, we simulate the time course of the level of protein X in 4,800 virtual cells based on the stochastic model (1) and estimate the variance D(x) at several discrete levels of x (red circles in Figure. 4C), from which we see that the simulation result coincides perfectly with the theoretical result (blue line in Figure. 4C). Moreover, we see that the variance function D(x) changes drastically around x_M, whereas there is no significant change in D(x) around x_L and x_H (Figure. 4C). This suggests that our method is effective in detecting the critical state of a multistable system by using the noisy data of gene expression, even if no detailed mathematical model is available.

Before leaving this section, we would like to point out that our method of identifying the critical state is not only effective in simple multistable systems with stable fixed points, but also effective in complicated multistable systems with periodic orbits or stable limit cycles, in which case the time-course data of the system will oscillate. To see this, we consider the following two-dimensional system of stochastic differential equations:

(12) \{\begin{array}{l}\dot{x}=y+\sqrt{\mathrm{2}ϵ}\xi x,\hfill \\ \dot{y}=-x^{\mathrm{3}}+x+\sqrt{\mathrm{2}\eta }\xi y,\hfill \end{array}

where ϵ and η are two noise levels and ξ_x and ξ_y are two independent standard white noises. If ϵ = η = 0, then the stochastic system is reduced to a deterministic system whose phase portrait is illustrated in Figure. 4D. The origin (0,0) (orange dot) is the only saddle and thus is the only critical point of the system. Interestingly, the stable and unstable manifolds at the critical point match up exactly, forming two homoclinic orbits. Moreover, the system has many periodic orbits inside and outside the two homoclinic orbits.

To see whether our method can be applied to identify the critical point of the stochastic system (12), we simulate the dynamics of the system at time h = 3:5 for an ensemble of 60 samples under two different initial values (red and orange dots) (Figure. 4D). The simulation result shows that the samples at time h (green dots) are rather concentrated if the system starts from the point ( - 1,0.45) (red dot) on a periodic orbit, whereas the samples at time h (blue dots) become rather scattered if the system starts from the critical point (orange dot). This suggests that the variance at the critical point is much larger than the variances at the points lying on the periodic orbits. Accordingly, we see that our method of identifying the critical state by seeking the maximum point of the variance function may be applied to various kinds of complicated multistable systems, not restricted to simple multistable systems with stable fixed points, such as the stochastic system (1).

DISCUSSION

Comparison with previous work

Multistability in biologic systems have been widely studied in recent two decades and it has become an important recurring theme in cell signaling. In this paper, we study a class of multistable biologic systems, i.e., heterogeneous bacterial populations with stochastic phenotype switching and also provide a data-driven method to identify the critical state of multistable bacterial systems. In recent literature, an influential paper about multistability is the paper of Angeli, Ferrell, and Sontag [ 60], who studied a large class of biologic systems with positive feedbacks and also provided a method of detecting bistability, bifurcations, and hysteresis in these systems. The themes of their paper is closely related to ours. Thus we find it necessary to discuss the similarities and differences between their paper and ours.

First, both the paper of Angeli et al. and ours studied multistability in biologic systems with a positive-feedback network. Recent studies show that multistability always arises in biologic systems that contain a positive-feedback loop and it has been proved that the existence of at least one positive-feedback loop is a necessary condition for the existence of multiple steady-states [ 61- 63]. In our paper, we show that phenotypic heterogeneity in bacterial populations often results from the underlying double-positive-feedback gene networks, which is consistent with the above theoretical results.

Second, the paper of Angeli et al. studied the deterministic nonlinear dynamics of multistable biologic systems and our paper studies the stochastic nonlinear dynamics of heterogenous bacterial populations. In fact, deterministic models may lead to phenotypic heterogeneity, but they cannot explain the widely observed phenomena of stochastic phenotype switching and bet-hedging within isogenic bacterial populations. In our paper, we use our stochastic model to provide a clear description of stochastic phenotype switching and bet-hedging within heterogenous bacterial populations and study the role that stochastic effects play in generating these important experimental phenomena.

Third, in the paper of Angeli et al., the relationship between the deterministic model of multistable systems and the widely used Markov chain model of population evolution is not clear. In our paper, we use the Freidlin-Wentzell theory to show that our stochastic model at the molecular level can be approximated by a Markov chain model at the cellular level, which reflects the long-term dynamics of multistable bacterial systems, when the noise level is small. It turns out that only by considering the stochastic dynamics of multistable bacterial systems can we unify the models at the two different levels.

Fourth, the paper of Angeli et al. provides a possible method to detect multistability and bifurcations in a large class of multistable systems satisfying the so-called “monotonicity” and “open-loop steady-state response” assumptions. Their method strongly depends on the properties of the system when the feedback is blocked. To detect multistability and bifurcations, the response data of the open-loop, feedback-blocked system to input stimuli must be obtained. This requirement is obvious too strong for realistic biologic systems. In our paper, however, we provide an effective data-driven method to identify the critical state of multistable biologic systems. In our method, only the time-course data of gene expression in individual cells are needed, even if no detailed mathematical model is available.

Reasonability of our stochastic model

In our stochastic model, we have used white noises with two different noise levels to describe the stochastic fluctuations in proteins X and Y, respectively. In this section, we shall explain the reasonability of this assumption. The content of this section is somewhat technical. Readers who are unfamiliar with knowledge of stochastic processes can skip this part.

In fact, the most precise model of the gene network in living cells is the chemical master equation (CME) [ 64]. Mathematically, the CME is the equation satisfied by the probability distribution of a Markov jump process which describes the copy number fluctuations of all participating macromolecules in the gene network. However, the dimension of the CME model is often so high that the theoretical expressions of many important quantities related to the system cannot be explicitly calculated. This makes the CME model difficult to be directly applied to practical problems with experimental data and observations.

To simplify the CME model, the mathematician Kurtz [ 65- 67] proved in his pioneering work that every CME model can be approximated reasonably well by the so-called chemical Langevin equation when the volume of the system is large. Mathematically, the chemical Langevin Equation [ 68] is the equation satisfied by the probability distribution of an SDE model which describes the concentration fluctuations of all participating macromolecules in the gene network. To be more precise, we assume that there are n macromolecules involved in the gene network whose concentrations are denoted by x₁, x₂, …, x_n. The general form of an SDE model which governs x₁, x₂; …, x_n is given by

(13) \dot{x}=b_i(x)+\sqrt{\mathrm{2}{ϵ}_i}{\displaystyle \sum _{j=\mathrm{1}}^n{\sigma }_{ij}(x){\xi }_j},i=\mathrm{1},\mathrm{2},\cdots ,n,

where x = (x₁, x₂, … , x_n), b(x) = (b₁(x), b₂(x), …, b_n(x)) is the drift coefficient of the SDE, σ(x) = (σ_ij(x))_n×n is the diffusion coefficient of the SDE, and ξ₁,ξ₂, …, ξ_n are n independent standard white noises. Mathematically, it can be proved that when the matrix σ(x) is bounded in both below and above, almost all the major properties of the SDE model related to the Freidlin-Wentzell theory will change little if σ(x) is regarded as the identity matrix [ 54, 57], in which case the SDE model (13) can be simplified as

(14) {\dot{x}}_i=b_i(x)+\sqrt{\mathrm{2}{ϵ}_i}{\xi }_i,\hspace{1em}i=\mathrm{1},\mathrm{2},\cdots ,n,

where the noise of the stochastic system is described by white noises. In fact, the SDE has been applied to model the concentration fluctuations of macromolecules in gene networks in previous work [ 20]. Some reviews on this topic can be found in [ 69- 71].

In the core double-positive-feedback gene network depicted in Figure. 1E, there are only two participating macromolecules, proteins X and Y, if the inducer concentration is regarded as a parameter. In this case, the simplified SDE model (14) is exactly our two-dimensional stochastic model (1). To make biologists, especially experimental biologists, understand the main results of this paper, we would like to use white noises to describe stochasticity in gene expression and use the simplified SDE (14) to model multistable bacterial systems. Under this simplification, our stochastic model (1) can be simply understood as the random perturbation of the traditional deterministic model.

Strengths and deficiencies of the SDE model

Just as George Box’s famous saying said: all models are wrong, but some are useful. Our model is no exception. In this paper, we use the SDE to model the concentration fluctuations of the participating macromolecules in the gene network of heterogenous bacterial populations. However, the SDE model proposed in this paper has some deficiencies. In fact, in living bacterial systems, some macromolecules, such as mRNA and protein molecules, may exist at very low copy numbers [ 52], in which case the concept of concentration makes no sense. Thus the SDE model cannot provide a good approximation of the CME model when the copy numbers of the participating macromolecules are very small.

However, compared with the CME model, the SDE model proposed in this paper has many strengths. First, the dimension of the CME model is often too high to be directly applied to practical problems with experimental data and observations. If the gene network of the bacterial system contains n macromolecules and the maximal possible copy numbers of these n macromolecules are N₁, N₂, …, N_n. Then the dimension of the CME model will be N₁N₂ … N_n. However, the SDE model compresses the dimension of the system to a large extent. If the gene network of the bacterial system contains n macromolecules, then the dimension of the SDE model is only n. This makes the system easy to be analyzed and makes many important quantities related to the system easy to be explicitly calculated.

Second, the mathematical theory of the SDE model is well developed, whereas that of the CME model is poorly developed due to its high complexity. In this paper, we use the Freidlin-Wentzell theory to explain stochastic phenotype switching and bet-hedging within isogenic bacterial populations and the consistence between our stochastic model and the traditional Markov chain models of population evolution. We believe that an analog of the Freidlin-Wentzell theory must exist in the CME model. However, such a theory in the CME model is not well developed up till now.

Third, the CME model is so abstract that it is not easy for biologists to understand and follow. In fact, without some deep mathematical knowledge, it is very difficult to understand the relationship between the CME model and the traditional deterministic model. However, the SDE model can be simply viewed as the random perturbation of the deterministic model. This makes the SDE model easy to be understood for both theoretical and experimental biologists.

Due to these reasons, we choose to use the SDE to model heterogenous bacterial populations, instead of the CME. We believe that the SDE, or equivalently, the chemical Langevin equation, which does not lose much information of the complicated CME, is a useful tool in the modeling of gene regulatory networks with inherent noises and in the analysis of the stochastic nonlinear dynamics of biologic systems.

Potential of our stochastic model

Given the small size of a cell and the small copy numbers of participating macromolecules, cellular processes during gene expression are inherently stochastic. In this paper, we establish a unified nonlinear stochastic model of multistable bacterial systems at the molecular level based on a core double-positive feedback gene network (Figure. 1E) and provide a clear description of phenotypic heterogeneity, stochastic phenotype switching, and bet-hedging within isogenic bacterial populations. Although we have used the expression levels of two stress-related genes to establish our model, we point out here that the variables in our model are not necessarily the expression levels of the stress-related genes, but can include more comprehensive indicators measured by flow cytometry and other techniques, such as FSC (roughly proportional to cell size), SSC (roughly proportional to cell granularity and complexity), and the ATP concentration (positively correlated with total cell energy).

In our recent study about antibiotic resistance of Escherichia coli (unpublished work), we found that the one-dimensional expression data of the hydrolase gene only lead to a monomodal distribution, but the multidimensional expression data of a group of stress-related genes lead to an apparent multimodal distribution. This phenomenon is described in this paper. Our simulation result shows that although the expression data of a group of genes are distributed within multiple attraction basins in the phase plane, their marginal distribution may overlap to a large extent so that we may not be able to observe a multimodal distribution if we only focus on the expression data of a specific gene. In our future work, we shall further apply the general framework discussed in this paper to the specific problem of antibiotic resistance of Escherichia coli.

Biologic systems with multistability are ubiquitous in nature. Some fundamental cellular processes, such as decision-making processes in cell cycle progression [ 72], cell fate determination [ 73- 76], and apoptosis [ 77, 78], display multistable features. In addition, multistability is also involved in disease progression, which can be divided into a normal state, a pre-disease state, and a disease state [ 55, 56]. Although various multistable biologic systems have different feedback regulatory networks, the mathematical structures behind them are quite similar. We hope that the stochastic approach discussed in this paper can give enlightenment to the understanding of biologic systems with multistability and to the analysis of the related new phenomena and new questions.

SUPPLEMENTARY MATERIALS

The supplementary materials can be found online with this article at DOI 10.1007/s40484-014-0035-5.

References

[1] Kussell, E. and Leibler, S. (2005) Phenotypic diversity, population growth, and information in fluctuating environments. Science, 309, 2075-2078

[2] Smits, W. K., Kuipers, O. P. and Veening, J. W. (2006) Phenotypic variation in bacteria: the role of feedback regulation. Nat. Rev. Microbiol., 4, 259-271

[3] Dubnau, D. and Losick, R. (2006) Bistability in bacteria. Mol. Microbiol., 61, 564-572

[4] Avery, S. V. (2006) Microbial cell individuality and the underlying sources of heterogeneity. Nat. Rev. Microbiol., 4, 577-587

[5] Dhar, N. and McKinney, J. D. (2007) Microbial phenotypic heterogeneity and antibiotic tolerance. Curr. Opin. Microbiol., 10, 30-38

[6] Lu, T., Shen, T., Bennett, M. R., Wolynes, P. G. and Hasty, J. (2007) Phenotypic variability of growing cellular populations. Proc. Natl. Acad. Sci. USA, 104, 18982-18987

[7] Veening, J. W., Smits, W. K. and Kuipers, O. P. (2008) Bistability, epigenetics, and bet-hedging in bacteria. Annu. Rev. Microbiol., 62, 193-210

[8] Fraser, D. and Kaern, M. (2009) A chance at survival: gene expression noise and phenotypic diversification strategies. Mol. Microbiol., 71, 1333-1340

[9] Jablonka, E. and Raz, G. (2009) Transgenerational epigenetic inheritance: prevalence, mechanisms, and implications for the study of heredity and evolution. Q. Rev. Biol., 84, 131-176

[10] Snijder, B. and Pelkmans, L. (2011) Origins of regulated cell-to-cell variability. Nat. Rev. Mol. Cell Biol., 12, 119-125

[11] Mao, J., Blanchard, A. E. and Lu, T. (2014) Slow and steady wins the race: A bacterial exploitative competition strategy in fluctuating environments. ACS Synth. Biol.

[12] Rulands, S., Jahn, D. and Frey, E. (2014) Specialization and bet hedging in heterogeneous populations. Phys. Rev. Lett., 113, 108102

[13] Acar, M., Mettetal, J. T. and van Oudenaarden, A. (2008) Stochastic switching as a survival strategy in fluctuating environments. Nat. Genet., 40, 471-475

[14] Salathé, M., Van Cleve, J. and Feldman, M. W. (2009) Evolution of stochastic switching rates in asymmetric fitness landscapes. Genetics, 182, 1159-1164

[15] Leisner, M., Stingl, K., Frey, E. and Maier, B. (2008) Stochastic switching to competence. Curr. Opin. Microbiol., 11, 553-559

[16] Gaál, B., Pitchford, J. W. and Wood, A. J. (2010) Exact results for the evolution of stochastic switching in variable asymmetric environments. Genetics, 184, 1113-1119

[17] Libby E, Rainey PB (2011) Exclusion rules, bottlenecks and the evolution of stochastic phenotype switching. P. Roy. Soc. B-Biol. Sci., 278: 3574-3583

[18] Rainey, P. B., Beaumont, H. J., Ferguson, G. C., Gallie, J., Kost, C., Libby, E. and Zhang, X. X. (2011) The evolutionary emergence of stochastic phenotype switching in bacteria. Microb. Cell Fact., 10, S14

[19] Ozbudak, E. M., Thattai, M., Lim, H. N., Shraiman, B. I. and Van Oudenaarden, A. (2004) Multistability in the lactose utilization network of Escherichia coli. Nature, 427, 737-740

[20] Süel, G. M., Garcia-Ojalvo, J., Liberman, L. M. and Elowitz, M. B. (2006) An excitable gene regulatory circuit induces transient cellular differentiation. Nature, 440, 545-550

[21] Tsang, J. and Van Oudenaarden, A. (2006) Exciting fluctuations: monitoring competence induction dynamics at the single-cell level. Mol. Syst. Biol., 2, 0025

[22] Schultz, D., Ben Jacob, E., Onuchic, J. N. and Wolynes, P. G. (2007) Molecular level stochastic model for competence cycles in Bacillus subtilis. Proc. Natl. Acad. Sci. USA, 104, 17582-17587

[23] Sonenshein, A.L., Hoch, J.A. and Losick R. (2002) Bacillus subtilis and its closest relatives: from genes to cells Washington: American Society for Microbiology.

[24] Errington, J. (2003) Regulation of endospore formation in Bacillus subtilis. Nat. Rev. Microbiol., 1, 117-126

[25] Morohashi, M., Ohashi, Y., Tani, S., Ishii, K., Itaya, M., Nanamiya, H., Kawamura, F., Tomita, M. and Soga, T. (2007) Model-based definition of population heterogeneity and its effects on metabolism in sporulating Bacillus subtilis. J. Biochem., 142, 183-191

[26] de Jong, I. G., Veening, J. W. and Kuipers, O. P. (2010) Heterochronic phosphorelay gene expression as a source of heterogeneity in Bacillus subtilis spore formation. J. Bacteriol., 192, 2053-2067

[27] Sureka K., Ghosh, B., Dasugpta, A., Basu J., Kundu, M. and Bose Z., (2008) Positive feedback and noise activate the stringent response regulator rel in mycobacteria. PLoS One3: e1771.

[28] Gefen, O. and Balaban, N. Q. (2009) The importance of being persistent: heterogeneity of bacterial populations under antibiotic stress. FEMS Microbiol. Rev., 33, 704-717

[29] Ghosh, S., Sureka, K., Ghosh, B., Bose, I., Basu, J. and Kundu, M. (2011) Phenotypic heterogeneity in mycobacterial stringent response. BMC Syst. Biol., 5, 18

[30] Veening, J. W., Stewart, E. J., Berngruber, T. W., Taddei, F., Kuipers, O. P. and Hamoen, L. W. (2008) Bet-hedging and epigenetic inheritance in bacterial cell development. Proc. Natl. Acad. Sci. USA., 105, 4393-4398

[31] Beaumont, H. J., Gallie, J., Kost, C., Ferguson, G. C. and Rainey, P. B. (2009) Experimental evolution of bet hedging. Nature, 462, 90-93

[32] Wolf, D. M., Vazirani, V. V. and Arkin, A. P. (2005) Diversity in times of adversity: probabilistic strategies in microbial survival games. J. Theor. Biol., 234, 227-253

[33] Gupta, P. B., Fillmore, C. M., Jiang, G., Shapira, S. D., Tao, K., Kuperwasser, C. and Lander, E. S. (2011) Stochastic state transitions give rise to phenotypic equilibrium in populations of cancer cells. Cell, 146, 633-644

[34] Zhou D, Wu D, Li Z, Qian M, Zhang MQ (2013) Population dynamics of cancer cells with cell state conversions. Quant. Biol. 1, 1-8.

[35] Karmakar, R. and Bose, I. (2007) Positive feedback, stochasticity and genetic competence. Phys. Biol., 4, 29-37

[36] Mitrophanov, A. Y. and Groisman, E. A. (2008) Positive feedback in cellular control systems. BioEssays, 30, 542-555

[37] Mantzaris, N. V. (2007) From single-cell genetic architecture to cell population dynamics: quantitatively decomposing the effects of different population heterogeneity sources for a genetic network with positive feedback architecture. Biophys. J., 92, 4271-4288

[38] Vellela, M. and Qian, H. (2009) Stochastic dynamics and non-equilibrium thermodynamics of a bistable chemical system: the Schlögl model revisited. J. R. Soc. Interface, 6, 925-940

[39] Qian, H. and Bishop, L. M. (2010) The chemical master equation approach to nonequilibrium steady-state of open biochemical systems: linear single-molecule enzyme kinetics and nonlinear biochemical reaction networks. Int. J. Mol. Sci., 11, 3472-3500

[40] Qian, H. (2011) Nonlinear stochastic dynamics of mesoscopic homogeneous biochemical reaction systemsłan analytical theory. Nonlinearity, 24, R19-R49

[41] Ge, H. and Qian, H. (2011) Non-equilibrium phase transition in mesoscopic biochemical systems: from stochastic to nonlinear dynamics and beyond. J. R. Soc. Interface, 8, 107-116

[42] Qian, H. (2012) Cooperativity in cellular biochemical processes: noise-enhanced sensitivity, fluctuating enzyme, bistability with nonlinear feedback, and other mechanisms for sigmoidal responses. Annu. Rev. Biophys., 41, 179-204

[43] Qian, H. and Ge, H. (2012) Mesoscopic biochemical basis of isogenetic inheritance and canalization: stochasticity, nonlinearity, and emergent landscape. Mol Cell Biomech, 9, 1-30

[44] McAdams, H. H. and Arkin, A. (1997) Stochastic mechanisms in gene expression. Proc. Natl. Acad. Sci. USA., 94, 814-819

[45] Elowitz, M. B., Levine, A. J., Siggia, E. D. and Swain, P. S. (2002) Stochastic gene expression in a single cell. Science, 297, 1183-1186

[46] Ozbudak, E. M., Thattai, M., Kurtser, I., Grossman, A. D. and van Oudenaarden, A. (2002) Regulation of noise in the expression of a single gene. Nat. Genet., 31, 69-73

[47] Paulsson, J. (2004) Summing up the noise in gene networks. Nature, 427, 415-418

[48] Kaern, M., Elston, T. C., Blake, W. J. and Collins, J. J. (2005) Stochasticity in gene expression: from theories to phenotypes. Nat. Rev. Genet., 6, 451-464

[49] Raser, J. M. and O’Shea, E. K. (2005) Noise in gene expression: origins, consequences, and control. Science, 309, 2010-2013

[50] Cai, L., Friedman, N. and Xie, X. S. (2006) Stochastic protein expression in individual cells at the single molecule level. Nature, 440, 358-362

[51] Yu, J., Xiao, J., Ren, X., Lao, K. and Xie, X. S. (2006) Probing gene expression in live cells, one protein molecule at a time. Science, 311, 1600-1603

[52] Xie, X. S., Choi, P. J., Li, G. W., Lee, N. K. and Lia, G. (2008) Single-molecule approach to molecular biology in living bacterial cells. Annu. Rev. Biophys., 37, 417-444

[53] Sanchez, A., Choubey, S. and Kondev, J. (2013) Regulation of noise in gene expression. Annu. Rev. Biophys., 42, 469-491

[54] Freidlin, M. I., Szücs, J. and Wentzell, A. D. (2012) Random perturbations of dynamical systems, New York: Springer-Verlag

[55] Kim, D., Rath, O., Kolch, W. and Cho, K. H. (2007) A hidden oncogenic positive feedback loop caused by crosstalk between Wnt and ERK pathways. Oncogene, 26, 4571-4579

[56] Kellershohn, N. and Laurent, M. (2001) Prion diseases: dynamics of the infection and properties of the bistable transition. Biophys. J., 81, 2517-2529

[57] Olivieri, E. and Vares, M. E. (2005) Large deviations and metastability. UK: Cambridge University Press.

[58] Liu R, Li M.Y., Liu Z.P., Wu J.R., Chen L.N. and Aihara K. (2012) Identifying critical transitions and their leading biomolecular networks in complex diseases. Sci. Rep.,

[59] Dai, L., Vorselen, D., Korolev, K. S. and Gore, J. (2012) Generic indicators for loss of resilience before a tipping point leading to population collapse. Science, 336, 1175-1177

[60] Angeli, D., Ferrell, J. E. Jr and Sontag, E. D. (2004) Detection of multistability, bifurcations, and hysteresis in a large class of biological positive-feedback systems. Proc. Natl. Acad. Sci. USA., 101, 1822-1827

[61] Gouzé, J. L. (1998) Positive and negative circuits in dynamical systems. J. Biol. Syst., 6, 11-15

[62] Snoussi, E. H. (1998) Necessary conditions for multistationarity and stable periodicity. J. Biol. Syst., 6, 3-9

[63] Cinquin, O. and Demongeot, J. (2002) Positive and negative feedback: striking a balance between necessary antagonists. J. Theor. Biol., 216, 229-241

[64] Gillespie, D. T. (1992) A rigorous derivation of the chemical master equation. Physica A, 188, 404-425

[65] Kurtz, T. G. (1970) Solutions of ordinary differential equations as limits of pure jump markov processes. J. Appl. Probab., 7, 49-58

[66] Kurtz, T. G. (1971) Limit theorems for sequences of jump markov processes approximating ordinary differential processes. J. Appl. Probab., 8, 344-356

[67] Kurtz, T. G. (1972) The relationship between stochastic and deterministic models for chemical reactions. J. Chem. Phys., 57, 2976-2978

[68] Gillespie, D. T. (2000) The chemical langevin equation. J. Chem. Phys., 113, 297-306

[69] Rao, C. V., Wolf, D. M. and Arkin, A. P. (2002) Control, exploitation and tolerance of intracellular noise. Nature, 420, 231-237

[70] Wilkinson, D. J. (2009) Stochastic modelling for quantitative description of heterogeneous biological systems. Nat. Rev. Genet., 10, 122-133

[71] Meister A, Du C, Li YH, Wong WH (2014) Modeling stochastic noise in gene regulatory systems. Quant. Biol.: 1-29

[72] Yao, G., Lee, T. J., Mori, S., Nevins, J. R. and You, L. (2008) A bistable Rb-E2F switch underlies the restriction point. Nat. Cell Biol., 10, 476-482

[73] Ferrell, J. E. Jr and Machleder, E. M. (1998) The biochemical basis of an all-or-none cell fate switch in Xenopus oocytes. Science, 280, 895-898

[74] Xiong, W. and Ferrell, J. E. Jr. (2003) A positive-feedback-based bistable ‘memory module’ that governs a cell fate decision. Nature, 426, 460-465

[75] Corson, F. and Siggia, E. D. (2012) Geometry, epistasis, and developmental patterning. Proc. Natl. Acad. Sci. USA., 109, 5568-5575

[76] Yan, L., Yang, M., Guo, H., Yang, L., Wu, J., Li, R., Liu, P., Lian, Y., Zheng, X., Yan, J., (2013) Single-cell RNA-Seq profiling of human preimplantation embryos and embryonic stem cells. Nat. Struct. Mol. Biol., 20, 1131-1139

[77] Eissing, T., Conzelmann, H., Gilles, E. D., Allgöwer, F., Bullinger, E. and Scheurich, P. (2004) Bistability analyses of a caspase activation model for receptor-induced apoptosis. J. Biol. Chem., 279, 36892-36897

[78] Spencer, S. L., Gaudet, S., Albeck, J. G., Burke, J. M. and Sorger, P. K. (2009) Non-genetic origins of cell-to-cell variability in TRAIL-induced apoptosis. Nature, 459, 428-432

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg