Email updates

Keep up to date with the latest news and content from BMC Systems Biology and BioMed Central.

Open Access Research article

Bridging time scales in cellular decision making with a stochastic bistable switch

Steffen Waldherr*, Jingbo Wu and Frank Allgöwer

Author Affiliations

Institute for Systems Theory and Automatic Control, Universität Stuttgart, Pfaffenwaldring 9, Stuttgart, Germany

For all author emails, please log on.

BMC Systems Biology 2010, 4:108  doi:10.1186/1752-0509-4-108

The electronic version of this article is the complete one and can be found online at:

Received:1 March 2010
Accepted:9 August 2010
Published:9 August 2010

© 2010 Waldherr et al; licensee BioMed Central Ltd.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.



Cellular transformations which involve a significant phenotypical change of the cell's state use bistable biochemical switches as underlying decision systems. Some of these transformations act over a very long time scale on the cell population level, up to the entire lifespan of the organism.


In this work, we aim at linking cellular decisions taking place on a time scale of years to decades with the biochemical dynamics in signal transduction and gene regulation, occuring on a time scale of minutes to hours. We show that a stochastic bistable switch forms a viable biochemical mechanism to implement decision processes on long time scales. As a case study, the mechanism is applied to model the initiation of follicle growth in mammalian ovaries, where the physiological time scale of follicle pool depletion is on the order of the organism's lifespan. We construct a simple mathematical model for this process based on experimental evidence for the involved genetic mechanisms.


Despite the underlying stochasticity, the proposed mechanism turns out to yield reliable behavior in large populations of cells subject to the considered decision process. Our model explains how the physiological time constant may emerge from the intrinsic stochasticity of the underlying gene regulatory network. Apart from ovarian follicles, the proposed mechanism may also be of relevance for other physiological systems where cells take binary decisions over a long time scale.


The dynamics of biological systems span a wide range of temporal and spatial scales. The interactions among dynamical properties on different scales govern the overall behavior of the biological system, and thus form an important area of computational research in biology [1]. A particularly interesting question in this context is how the behavior on a slow time scale emerges mechanistically from the dynamics on fast time scales. For example, how do cell population dynamics in tissues, which may evolve on a time scale of months, years or even decades, originate from the dynamics of the underlying gene regulatory networks, with a time scale of just minutes to hours?

In this work, we aim at bridging the time scale from gene regulation to cellular transformation processes on the tissue or cell population level. We specifically consider cellular transformation processes based on a bistable biochemical switch. Such switches have two distinct stable stationary states, and the cell initiates a transformation when the switch changes from one stable state to the other one. Bistable switches have previously been used to model a large number of cellular transformation events, such as progression through cell cycle arrest in the maturation of Xenopus oocytes [2,3] or initiation of programmed cell death [4] and cellular differentiation [5] in higher organisms. Most models for these systems are constructed as deterministic models, and thus an external stimulus is required to induce changes in the switch's state. In addition, stochastic models for biochemical switches within a variety of biological processes have been formulated, for example the lac operon in E. coli [6,7], the genetic toggle switch [8], or a generic phosphorylation/dephosporylation cycle [9]. The typical questions that have been adressed by stochastic switch models are for example the steady state probability distribution of the different possible states of the switch [8], or the residence times in these states [9]. In the previously proposed stochastic models of bistable biochemical switches, cells are able to switch forth and back between the possible qualitative states of the switch. While this is appropriate if the switch serves to choose a cellular state based on environmental conditions, such as for example in the galactose utilization network in yeast [10], this feature should not be held up for transformation processes. In transformation processes, subsequent mechanisms, which are not included in the model description, are in place to ensure irreversibility once the switch changed its qualitative state from the initial condition. The most obvious example for such mechanisms is cell death, where the model of the biochemical switch does not hold anymore once the cell transitions to the "dead" state.

In this work, we consider irreversible transformation processes based on a stochastic switch model, which apparently do not require any external stimulus to be initiated, where the transition is based only on stochastic fluctuations. Despite the stochasticity, we see in this paper that the dynamics of the switch still follow reliable temporal characteristics. Reliable thereby means that in a large population of cells, the number of cells that have already initiated the transformation can be described deterministically with high accuracy. We propose a generic transformation process, where a phenotypical change in the state of a cell is initiated as soon as a bistable biochemical switch changes its internal state. In previous studies, random switching caused by internal fluctuations is usually attributed to pathological events [11]. In the mechanism proposed here, random switching has a regular physiological function.

A striking example for the kind of transformation processes we aim to describe is involved in mammalian oocyte maturation. In mammalian females, all or almost all of the oocytes that will ovulate through the organism's life-span are already present at birth or shortly thereafter as a population of so-called primordial follicles. Throughout the organism's reproductive life, follicles undergo the primordial to primary transition, which marks the start of a development process that will eventually lead to either ovulation or removal of the oocyte through atresia [12,13]. In this way, there is a steady supply of mature follicles for ovulation, while the pool of primordial follicles is gradually depleted. The mechanisms through which the follicle transition is initiated are largely unknown, although a number of ovarian factors that may be relevant have been identified experimentally [14-16]. Importantly, the transition seems to be regulated locally in the ovary, and not through the endocrine system [17]. An astonishing observation in this process is that in one follicle, the transition may occur already a few months after generation of the primordial follicle pool, while another follicle may stay several decades (for organisms with a sufficiently long lifespan) in the resting stage before growth is initiated. From the medical side, a misregulation of this process is implicated in premature ovarian failure due to follicle depletion, which is a major reason for infertility in human females. By way of a case study, we apply the proposed transformation mechanism to the problem of growth initiation in ovarian follicles. Including also cell-cell interactions supported by experimental evidence, we obtain a physiologically plausible model for this process, showing very good agreement with human clinical data on a time scale of several decades.


Deterministic model of a bistable switch

The model of a bistable switch that we use is based on a positive feedback loop between two components. Consider a biochemical reaction network involving the two molecular species X and Y. Mathematically, the temporal evolution of the amounts of the two species is described with the ordinary differential equation

<a onClick="popup('','MathML',630,470);return false;" target="_blank" href="">View MathML</a>


where x and y denote the amounts of X and Y, respectively. The network is illustrated in Figure 1. The vector (x,y)T will be referred to as the microstate of the biochemical reaction system. Ultrasensitivity, which is required to achieve bistability [2], is generated by the Hill-type production rates v2 and v4. In the sequel, we will assume that the molecular species X and Y represent gene transcripts, and the amounts x and y indicate the respective transcript copy number. The nominal parameter values that we use are given in Table 1. For simplicity, we assume that the parameters are symmetric, i.e. V1 = V2, M1 = M2 and u1 = u2. The parameter values are within the physiological range for typical gene transcription processes. In particular, the degradation rate of 0.01<a onClick="popup('','MathML',630,470);return false;" target="_blank" href="">View MathML</a> corresponds to a gene transcript half-life time of about 70 minutes. Typical transcript half-life times in mammalian cells are in a range from tens of minutes to several hours [18], but can of course vary significantly depending on the gene and regulatory influences, with an estimated variation of 200 fold among different genes [19]. The minimal transcription rate of X is given by k1 and corresponds to 3.3 transcripts that are produced per hour. The transcription rate upon maximal activation is given by V1,2 and corresponds to 33 transcripts produced per hour. Upon maximal activation, this would yield a steady state mRNA copy number of 55 molecules per cell. The typical range of mRNA copy numbers in mammalian cells seems to be on the order of a few to hundreds [20,21].

Table 1. Nominal parameter values for the bistable switch model (1).

thumbnailFigure 1. Network schematic for the bistable switch model (1).

For two-dimensional systems, it is convenient to check bistability by considering nullclines in the state space [22]. With this graphical representation, it is also easy to evaluate how good the two stable states are actually separated [23]. The nullclines for the model given in (1), with nominal parameter values, are depicted in Figure 2A. From the figure, it is clear that there are three equilibrium points, labelled I, II and III. A stability analysis of the equilibrium points shows that the deterministic system described by (1) is bistable, and the corresponding reaction network implements a bistable switch. We construct a macrostate for this system by defining the two sets Ωoff, Ωon ⊂ ℝ2 corresponding to the switch being off or on, respectively. Ωoff contains the equilibrium point I, and Ωon contains III. For our model, we define

thumbnailFigure 2. Characterisation of the phase space in the bistable switch model (1). A: Phase space diagram for the deterministic model of the bistable switch. Black lines are nullclines for the variables x and y in the deterministic switch model (1), with their intersections corresponding to equilibria of the switch. I and III are stable equilibrium points, II is an unstable one. Trajectories converge to either I or III, depending on the initial condition, as shown for the sample trajectories plotted as light blue lines. B: Schematic illustration of the configuration space for the Markov process (5) describing the cell transformation process. Circular nodes below the dashed line correspond to possible configurations (X; Y)T of the switch, and the arrows between the nodes correspond to transitions in the configuration due to reactions. The configurations above the dashed line are collapsed into the on state, which is assumed to be irreversible due to subsequent transformation processes.

<a onClick="popup('','MathML',630,470);return false;" target="_blank" href="">View MathML</a>

with suitable parameter L. With model parameters as given in Table 1, a suitable choice which we will use in this work is L = 55.

Stochastic model of a bistable switch

The deterministic model of the bistable switch discussed in the previous section is suitable to describe the existence of two distinct macrostates, corresponding to stable equilibrium points in the model. However, to capture transitions between these macrostates which are caused by intrinsic fluctuations, a stochastic model has to be considered. In a stochastic setting, the amounts of molecular species may only take discrete values from the set <a onClick="popup('','MathML',630,470);return false;" target="_blank" href="">View MathML</a>. The stochastic state of the switch at time t is given by the discrete probability distribution p(X, Y, t), which for each microstate <a onClick="popup('','MathML',630,470);return false;" target="_blank" href="">View MathML</a> gives the probability that the switch is in the microstate (X, Y )T at time t:

<a onClick="popup('','MathML',630,470);return false;" target="_blank" href="">View MathML</a>


To describe the temporal evolution of the probability distribution, we use the chemical master equation (CME) [24]. The reaction network for the bistable switch is not described with elementary reactions only, and thus it is not possible to construct the CME according to its rigorous derivation [25]. However, a theoretical investigation by Rao and Arkin [26] has shown that as an approximation, the propensity functions for state transitions can be taken from the according deterministic reaction rate laws. Thus, for the bistable switch described above, we can formulate the CME

<a onClick="popup('','MathML',630,470);return false;" target="_blank" href="">View MathML</a>


for <a onClick="popup('','MathML',630,470);return false;" target="_blank" href="">View MathML</a>, where the reaction propensities vi, i = 1,..., 5, are the same expressions as in the deterministic model (1).

In the stochastic model (3), we aim to identify the qualitative states on and off as in the deterministic model. For many biochemical systems, the stable equilibrium states in the deterministic description correspond to peaks in the probability distribution p(X, Y, t) [10], although there are also cases where this is not true, for example systems where extinction of molecular species is possible [27]. For the stochastic switch model (3), simulations suggest that we indeed obtain two peaks in the probability distribution close to the stable equilibrium points of the deterministic model (1) (see Figure 3).

thumbnailFigure 3. Steady state probability distribution for the stochastic bistable switch. 500 realizations of the stochastic reaction network model (3) were generated using the Gillespie algorithm in the stochastic simulation software Dizzy [41,42]. Each realization was for a simulated time of 300 years, and the steady state probability distribution was generated from the samples after discarding a transient phase of 50 years simulated time, using a total of about 5 · 107 data points.

In the stochastic description, we can compute the probabilities that the switch is in any of its two macrostates directly from a solution of the CME. Define poff (t) and pon(t) as the probabilities that the switch is off and on, respectively. Given a solution of the CME, these can be computed by summing up the probabilities that the system is in the corresponding microstates, i.e. <a onClick="popup('','MathML',630,470);return false;" target="_blank" href="">View MathML</a>, and equivalently for poff(t).

A transformation process modelled with a stochastic switch

Cellular transformation processes are often based on a bistable biochemical or genetic switch. In the initial state of the cell, the switch would be in the off state. Switching to the on state implies a significant change in the amount of an involved signaling molecule, e.g. a transcription factor. If the on state is maintained for some time, this change would result in a larger phenotypical change of the cell, e.g. through significant changes in gene expression. The mechanisms that induce this change are not included in the stochastic switch model, but from a signaling perspective downstream of it.

Most transformation processes rely on specific external stimuli, and the cell will initiate the transformation upon encountering the required stimulus. There are however examples where such a stimulus is not strictly required, and this is the case that we are dealing with in this paper. Moreover, we will focus on the behavior of cell populations, studying the problem how the temporal dynamics of the transformation process evolve in a pool of many cells.

The basic mechanism that actually triggers the bistable switch in our model without an external stimulus are the intrinsic fluctuations of concentrations in any biochemical reaction network, that are due to the stochastic nature of chemical reactions. As a rare event, these fluctuations may become so large that the microstate of the system crosses the separatrix between the domains of attraction in the deterministic system. As a consequence, the microstate around the other stable equilibrium point will become strongly attractive, and the switch will change its macrostate to on with a high probability. In this paper, we assume that the transformation is irreversible, which fits well to the process of follicle growth initiation. Also other processes such as programmed cell death are irreversible.

The described transformation process is easily modelled as a continuous-time Markov process. If the switch is in the macrostate off, then we directly use the microstates and transition probabilities of the underlying biochemical reaction network to model the transformation process. To account for the irreversibility of the transformation, all microstates <a onClick="popup('','MathML',630,470);return false;" target="_blank" href="">View MathML</a> are collapsed to one state of the Markov process, labeled with "on" in Figure 2B, which is an absorbing state. The transitions of other microstates to the absorbing state are governed by the propensity functions for the corresponding transitions in the underlying biochemical network. The resulting state space for the Markov process model of the transformation process is shown in Figure 2B.

In our model of the stochastic switch, the macrostate off is defined by a compact region in state space. As a consequence, the Markov model of the considered transformation process has a finite state space, and can therefore be treated computationally with standard approaches. Let P(t) ∈ ℝn denote the complete probability state vector of the system,

<a onClick="popup('','MathML',630,470);return false;" target="_blank" href="">View MathML</a>


The master equation can be written as the linear ordinary differential equation

<a onClick="popup('','MathML',630,470);return false;" target="_blank" href="">View MathML</a>


where A ∈ ℝn × n is the state transition matrix. The matrix A can be computed directly from the values of the reaction propensity functions in each microstate [28]. The differential equation (5) can be solved using standard tools for numerical integration. For the results described in this paper, we used the ode15s function in MATLAB (The MathWorks, Natick, MA) to obtain a numerical solution of (5).

Results and Discussion

A hypothetical mechanism for oocyte maturation

In this section, we suggest a biochemical mechanism that offers a molecular explanation for the large depletion times of several decades in the human oocyte pool. The model is based on experimental evidence obtained in a very informative series of studies by Skinner and colleagues (see [13] for a review), where the influence of several ovarian factors on the primordial to primary transition as well as some interactions between them have been studied. Because a positive feedback loop is necessary for a bistable switch [29], we have specifically searched for such an interconnection.

Primary ovarian follicles are composed by three main cell types: a single oocyte as the main component, and granulosa and theca cells surrounding the oocyte [13]. Experimental evidence suggests a positive feedback circuit involving two ovarian factors that are relevant in the primordial to primary transition: the factor KIT ligand (KITL) is produced by granulosa cells and stimulates both the oocyte and surrounding theca cells to promote follicle development. Moreover, KITL stimulates the production of both keratinocyte growth factor (KGF) and hepatocyte growth factor (HGF) in the surrounding theca cells. KGF and HGF themself stimulate the production of KITL in the granulosa cells, thus providing a positive feedback loop [30]. Moreover, the oocyte of primordial and developing follicles produces basic fibroblast growth factor (bFGF), which acts on surrounding granulosa cells and has been shown to increase the expression of KITL [16].

These pieces of experimental evidence thus support the hypothetical mechanism that is shown in Figure 4. Our simplistic mathematical model presented in (1) and Figure 1 can be used to describe this mechanism, where the variable x represents granulosa-derived KITL activity and y represents theca-derived KGF and HGF activity. The reaction v1 describes the influence on KITL expression of oocyte-derived bFGF, which is here assumed to be constant. The reactions v2 and v4 arise from the positive feedback interconnection, whereas v3 and v5 describe a constitutive degradation of KITL, KGF and HGF.

thumbnailFigure 4. Hypothetical biochemical network for the primordial to primary transition in ovarian follicles.

The stochastic switch generates reliable long-term behavior

The differential equation (5) that governs the initiation probabilities of the irreversible transformation process is a linear ordinary differential equation, so in principle it can be solved analytically. Due to the size of the system (n = 1653 in this example), this is however not feasible. Yet, we can characterize the probability that a given cell has initiated the transformation process by the explicit formula

<a onClick="popup('','MathML',630,470);return false;" target="_blank" href="">View MathML</a>


where c1 > 0, 0 <λ1 < Re (λi) for i = 2, ..., n, and the ci(t) are polynomials in t. The mathematical derivation of (6) is provided in the appendix.

From considering the general form of the analytical solution given in (6), we obtain two important conclusions about the stochastic transformation process. First, we observe that the probability for a given cell to initiate the transformation tends to 1 as time increases. Second, because λ1 is the dominant decay rate, for larger times t ≫ 0 the probability of not having initiated the transformation can be approximated by <a onClick="popup('','MathML',630,470);return false;" target="_blank" href="">View MathML</a>, a simple exponential decay. For the biochemical parameter values given in Table 1, the numerical solution for pon(t) is shown in Figure 5A. For these parameter values, which are in the physiological range for the considered biological processes, we indeed get to a time scale of years to decades in the probability of the transformation event, with a half-life time of about 5.9 years. Let us now move to the population level, and consider a pool of cells, each of them being subject to the considered transformation process with a bistable switch. In the first step, we make the simplistic assumption that no interactions among the cells are taking place, so individual transformations are probabilistically independent events. The number of remaining cells Nr(t) can easily be characterized by a binomial distribution as

<a onClick="popup('','MathML',630,470);return false;" target="_blank" href="">View MathML</a>


where N0 is the initial number of cells in the pool. The properties of the binomial distribution give the expected number of cells remaining in the pool as

<a onClick="popup('','MathML',630,470);return false;" target="_blank" href="">View MathML</a>


The probability distribution P (Nr(t) = N) for the population size in the transformation process considered in this paper is shown in Figure 5B as a function of both cell number N and time t. The number of initial cells N0 = 106 was chosen from the reported range of ovarian follicles, 7 · 105 to 2 · 106 in human females at birth [31]. For each point in time, the distribution has a very sharp peak, which indicates that the average value E[Nr(t)] is a reliable prediction for the number of cells that have already undergone the transformation at a given time.

A relevant characteristic of the considered process is the time at which the initial cell population is depleted, i.e. when nearly all cells have undergone the transformation. To make this notion precise, we introduce the depletion number Nd. The depletion time Td is defined as the smallest time t such that Nr(t) ≤ Nd, i.e. only Nd cells are remaining in the initial population. For the process of follicle growth initiation, we use Nd = 103, which has been considered to mark the onset of menopause [32].

thumbnailFigure 5. Dynamical characteristics of the stochastic bistable switch on the single cell level and the population level. A: Probability of transformation event pon(t) B: Population size probability distribution over time. C: Probability density function of the depletion time Td.

The cumulative probability distribution function for the depletion time Td is computed from the distribution obtained in (7) as

<a onClick="popup('','MathML',630,470);return false;" target="_blank" href="">View MathML</a>


The probability density function for the depletion time is computed by taking the derivative of the cumulative probability distribution function (9). The resulting probability density function for the depletion time in follicle growth initiation is shown in Figure 5C. From the density function, the expected value and the standard deviation are obtained as E[Td] = 59.1 years and <a onClick="popup('','MathML',630,470);return false;" target="_blank" href="">View MathML</a> years, respectively.

The expected value for Td can also be computed by solving <a onClick="popup('','MathML',630,470);return false;" target="_blank" href="">View MathML</a>. Using (6), it can thus be approximated by <a onClick="popup('','MathML',630,470);return false;" target="_blank" href="">View MathML</a>, where λ1 is the dominant decay rate of the process.

Next, we compare the computed statistical characteristics of the follicle depletion process to medical data. Explicit follicle counts are only sparsely available. However, the available pieces of data indicate that fluctuations in actual follicle numbers are larger than predicted by our model [33]. Concerning the depletion time, a recent medical study suggests an average age of 51.1 years for the onset of menopause, with a standard deviation of 3.8 years [34]. Our model predicts a depletion time of Td = 59.1 years, which is reasonably close to the experimentally observed depletion time. However, the standard deviation of 0.27 years in our model is significantly less than observed from medical data. In summary, even though our model is based on a highly stochastic process, the analysis reveals that it leads to much more reliable temporal characteristics than observed in the real system. This indicates that stochastic effects alone may not be sufficient to explain the heterogeneity observed in the follicle depletion process.

An alternative explanation would be by heterogeneous parameter values among individual organisms. This explanation is also supported by statistical analyses of medical data [34], where it is suggested that the onset of menopause is largely based on genetic factors, which would be related to parameter values in our model. To investigate this possibility, we have computed the expected depletion times for different parameter values. The computation was based on the eigenvalues of the transition matrix A and the approximation <a onClick="popup('','MathML',630,470);return false;" target="_blank" href="">View MathML</a>. The results are given in Table 2. From these results, we note that even small parameter variations in the model of the bistable switch lead to very large variations in the expected depletion time. This is not realistic for a biological system, and in the following section we explore mechanisms to increase the robustness of the depletion time with respect to parameter variations.

Table 2. Expected depletion times (years) in the single cell model (5).

Increased robustness by interactions on the population level

In the last section, we have characterized the properties of the transformation process based on a bistable switch, with the depletion time of a pool of cells subject to the transformation as characteristic output of the model. We have shown that the proposed model produces reliable depletion times, in the sense of a small standard deviation, for fixed values of the biochemical parameters. However, we have also observed that the average depletion time in the basic model is quite sensitive to variations in the biochemical parameters. Clearly, this large sensitivity is not acceptable in a model that should be a meaningful representation of the primordial to primary follicle transition. In this section, we propose an additional mechanism that reduces the sensitivity of the average depletion time significantly.

The additional mechanism is based on the experimental observation that follicles in later stages of development actively suppress the primordial to primary transition in resting follicles [13]. The inhibition of follicle growth initiation is mediated by the Anti-Mülerian hormone (AMH), which is produced by growing follicles and interferes with stimulatory signals by KITL, bFGF, and KGF [35]. Although AMH is known to signal via SMAD proteins [36], the molecular mechanisms of follicle growth inhibition by AMH seem to be unknown. To include the inhibitory effect into the simplistic switch model (1), we assume that the rate of KITL production in primordial follicles is reduced with an increasing number of growing follicles. This is achieved by changing k1 in the original model given in (1) from a constant parameter to the expression

<a onClick="popup('','MathML',630,470);return false;" target="_blank" href="">View MathML</a>


where k1,max is the maximal production rate of KITL, n2 is the number of growing, AMH producing follicles, and Kn is an additional parameter. While follicle development is a complex process with many intermediate stages [31], in this analysis we use a simple two-state population model, where n1 denotes the number of primordial follicles, and n2 the number of growing follicles. The assumptions of the model are that primordial follicles initiate growth with a rate as determined by λ1 in (6). Due to k1 depending on n2 as defined in (10), we obtain a dependency of λ1 on n2. Growing follicles are assumed to stay in this stage for a constant amount of time t, after which they leave the pool either through ovulation or atresia. From these specifications, one can derive a model given by the system of delay-differential equations

<a onClick="popup('','MathML',630,470);return false;" target="_blank" href="">View MathML</a>


where λ1(n2) is the decay rate computed from the transition matrix A(n2) in (5), with k1(n2) as in (10). Using the parameters in Table 3, the population model given by (11) now predicts a depletion time of Td = 50.0 years, which is almost equal to the depletion time suggested by the medical study [34]. The development of the ovarian follicle pool over time, as predicted by the model in (11), is shown in Figure 6. The prediction is compared to clinical data of follicle numbers at different ages taken from [37]. Although the parameters have only been adjusted to the depletion time, the predicted time course is reasonable close to the clinical data. In particular, the proposed population model (11) intrinsically captures the previously observed increase in the follicle depletion rate at an age of approximatively 38 years [37]. In order to investigate the sensitivity of the extended model to variations in the biochemical parameters, we have computed again the expected depletion times for different parameter values. The results are given in Table 4. The variation in the depletion time is significantly reduced, compared to the model (5), where follicle interactions are neglected. It should also be pointed out that the depletion time is quite insensitive towards variations in the two parameters Kn and τ which were newly introduced in the population model. This result illustrates that the robustness of the depletion time with respect to parameter variations may be substantially increased by adding interactions among individual follicles to the proposed model of the transformation process.

Table 3. Nominal parameter values for the population model (11).

Table 4. Expected depletion times (years) in the population model (11).

thumbnailFigure 6. Evolution of follicle number. Model predictions from (11) (line) vs. clinical data from [37] (crosses).


In this paper, we deal with a fundamental question in the development of multicellular organisms: How can biochemical reactions and genetic mechanism acting on the scale of minutes in individual cells generate dynamics with characteristic times of years to decades on the tissue level? As a possible mechanism to achieve this, we propose a generic transformation process based on a bistable stochastic switch. From the underlying genetic interactions and biochemical reactions, the process can be modelled as a continuous-time Markov process. We show that the proposed stochastic mechanism generates reliable long-term behavior on the population level, with cells undergoing the transformation with an exponentially decaying rate. Thereby, the decay rate is equal to the dominant eigenvalue of the transition matrix describing the underlying biochemical network. Due to bistability of the considered switch, this dominant eigenvalue corresponds to very slow dynamics, thus leading to the very long timescale as observed in the simulations. We pose the hypothesis that a biological instance of this mechanism is present in the development of ovarian follicles. To describe this process, we constructed a simple model of a bistable switch in the primordial to primary transition for ovarian follicles. The model is based on experimentally determined factors and their interactions in the different cell types constituting the ovarian follicles. Although it is not assured that a bistable switch in ovarian follicles will indeed be based on the factors that we have used here, the basic mechanism would work equivalently well with other factors.

Despite its simplicity, our model explains well how the long-term characteristics of follicle development may reliably be generated by biochemical reactions occurring on much shorter time scales. Keeping the model simple serves two purposes: first, it shows that the dynamics of follicle growth initiation can be generated by a quite simple mechanism. Clearly, additional pathways and regulatory feedback interactions that we have not included in this model can be expected to be present in the system. These may serve to increase robustness of the network, or to provide additional inputs to control the transition rate, e.g. for the endocrine system. Second, the simplicity of the model allows us to solve the chemical master equation for the network numerically, and thus to obtain a good quantitative description of the model behavior.

As a possible shortcoming of the basic model on the single cell level, we observe an unrealistic large sensitivity of the follicle depletion time with respect to parameter variations. By adding the experimentally supported inhibition of follicle growth initiation by later-stage growing follicles, the sensitivity of the depletion time could be reduced significantly. Apart from the inhibition included in the model, other interactions among individual follicles seem to play a role in the primordial to primary transition [38]. We envision that the inclusion of more regulatory interactions may further decrease the sensitivity of the depletion time with respect to parameter variations to a physiologically plausible level.

Authors' contributions

SW conceived of the study, participated in its design, helped with the computational implementation and data analysis, and drafted the manuscript. JW carried out the computational implementation and helped with the data analysis. FA participated in the design of the study, and helped to draft the manuscript. All authors read and approved the final manuscript.

Appendix: Computation of the transition probability

In this section, we prove that the probability that a given cell has undergone the considered transformation process is given by pon(t) as in (6). The proof is based on considering the solution of the underlying CME (5).

Since the last microstate is an absorbing state of the Markov process, (5) can be written as

<a onClick="popup('','MathML',630,470);return false;" target="_blank" href="">View MathML</a>


where Arev ∈ ℝ(n - 1) × (n - 1) describes the interactions among the non-absorbing states, and aabs ∈ ℝ1 × (n - 1) describes the transition propensities to the absorbing state.

Let us first derive some essential properties of the matrix Arev. Since A is a stochastic matrix, we have

<a onClick="popup('','MathML',630,470);return false;" target="_blank" href="">View MathML</a>


for i = 1, ..., n - 1 i.e. Arev is diagonally dominant. Thus, Gersgorin's theorem [39] asserts that all eigenvalues of Arev have a non-positive real part. Even more, since aabs is non-zero, (13) holds with a strict inequality for at least one i. Thus, by Theorem 10.7.2 in [39], all eigenvalues of Arev have negative real part. By the properties of the considered biochemical network, Arev is irreducible, and its off-diagonal elements are non-negative. From Corollary 4.3.2 in [40], it follows that Arev has an eigenvalue λ1 ∈ ℝ with algebraic multiplicity 1 and a strictly positive corresponding eigenvector v1 such that Re λ <λ1 for all λ λ1 in the spectrum of Arev.

Denoting Prev = (P1, ..., Pn-1)T we have <a onClick="popup('','MathML',630,470);return false;" target="_blank" href="">View MathML</a>. From the previously derived properties of the matrix Arev, the general solution of this differential equation is given by

<a onClick="popup('','MathML',630,470);return false;" target="_blank" href="">View MathML</a>


where <a onClick="popup('','MathML',630,470);return false;" target="_blank" href="">View MathML</a> are polynomials in t and <a onClick="popup('','MathML',630,470);return false;" target="_blank" href="">View MathML</a> is a constant coefficient, depending on the initial condition Prev(0). The condition Prev(t) ≥ 0 for all t implies that <a onClick="popup('','MathML',630,470);return false;" target="_blank" href="">View MathML</a>. For a non-negative initial condition Prev(0) with at least one positive element, we have <a onClick="popup('','MathML',630,470);return false;" target="_blank" href="">View MathML</a>. The transition probability pon(t) is computed as

<a onClick="popup('','MathML',630,470);return false;" target="_blank" href="">View MathML</a>


where <a onClick="popup('','MathML',630,470);return false;" target="_blank" href="">View MathML</a> and <a onClick="popup('','MathML',630,470);return false;" target="_blank" href="">View MathML</a>.


SW and FA acknowledge support by the German Research Foundation (DFG) through the Cluster of Excellence in Simulation Technology (EXC 310) at the University of Stuttgart.


  1. Martins ML, Ferreira SC Jr, Vilela MJ: Multiscale models for biological systems.

    Curr Opin Colloid Interface Sci 2010, 15:18-23. Publisher Full Text OpenURL

  2. Ferrell JE, Xiong W: Bistability in cell signaling: How to make continuous processes discontinuous, and reversible processes irreversible.

    Chaos 2001, 11:227-236. PubMed Abstract | Publisher Full Text OpenURL

  3. Ferrell JE, Machleder EM: The biochemical basis of an all-or-none cell fate switch in Xenopus oocytes.

    Science 1998, 280(5365):895-898. PubMed Abstract | Publisher Full Text OpenURL

  4. Eissing T, Conzelmann H, Gilles ED, Allgöwer F, Bullinger E, Scheurich P: Bistability Analyses of a Caspase Activation Model for Receptor-induced Apoptosis.

    J Biol Chem 2004, 279(35):36892-97. PubMed Abstract | Publisher Full Text OpenURL

  5. Chickarmane V, Enver T, Peterson C: Computational modeling of the hematopoietic erythroid-myeloid switch reveals insights into cooperativity, priming, and irreversibility.

    PLoS Comput Biol 2009, 5:e1000268. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  6. Mettetal JT, Muzzey D, Pedraza JM, Ozbudak EM, van Oudenaarden A: Predicting stochastic gene expression dynamics in single cells.

    Proc Natl Acad Sci 2006, 103:7304-9. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  7. Kaufmann BB, Yang Q, Mettetal JT, van Oudenaarden A: Heritable Stochastic Switching Revealed by Single-Cell Genealogy.

    PLoS Biology 2007, 5(9):e239. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  8. Tian T, Burrage K: Stochastic models for regulatory networks of the genetic toggle switch.

    Proc Natl Acad Sci 2006, 103(22):8372-8377. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  9. Krishnamurthy S, Smith E, Krakauer D, Fontana W: The stochastic behavior of a molecular switching circuit with feedback.

    Biol Direct 2007, 2:13. PubMed Abstract | BioMed Central Full Text | PubMed Central Full Text OpenURL

  10. Song C, Phenix H, Abedi V, Scott M, Ingalls BP, Kaern M, Perkins TJ: Estimating the stochastic bifurcation structure of cellular networks.

    PLoS Comput Biol 2010, 6(3):e1000699. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  11. Isaacs FJ, Hasty J, Cantor CR, Collins JJ: Prediction and measurement of an autoregulatory genetic module.

    Proc Natl Acad Sci 2003, 100(13):7714-7719. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  12. Fortune JE, Cushman RA, Wahl CM, Kito S: The primordial to primary follicle transition.

    Mol Cell Endocrinol 2000, 163(1-2):53-60. PubMed Abstract | Publisher Full Text OpenURL

  13. Skinner MK: Regulation of primordial follicle assembly and development.

    Hum Reprod Update 2005, 11(5):461-471. PubMed Abstract | Publisher Full Text OpenURL

  14. Parrott JA, Skinner MK: Kit-ligand/stem cell factor induces primordial follicle development and initiates folliculogenesis.

    Endocrinol 1999, 140(9):4262-4271. Publisher Full Text OpenURL

  15. Castrillon DH, Miao L, Kollipara R, Horner JW, DePinho RA: Suppression of ovarian follicle activation in mice by the transcription factor Foxo3a.

    Science 2003, 301(5630):215-218. PubMed Abstract | Publisher Full Text OpenURL

  16. Nilsson EE, Skinner MK: Kit ligand and basic fibroblast growth factor interactions in the induction of ovarian primordial to primary follicle transition.

    Mol Cell Endocrinol 2004, 214(1-2):19-25. PubMed Abstract | Publisher Full Text OpenURL

  17. Braw-Tal R: The initiation of follicle growth: the oocyte or the somatic cells?

    Mol Cell Endocrinol 2002, 187(1-2):11-18. PubMed Abstract | Publisher Full Text OpenURL

  18. Ross J: mRNA stability in mammalian cells.

    Microbiol Rev 1995, 59(3):423-450. PubMed Abstract | PubMed Central Full Text OpenURL

  19. Hargrove JL, Hulsey MG, Beale EG: The kinetics of mammalian gene expression.

    Bioessays 1991, 13(12):667-674. PubMed Abstract | Publisher Full Text OpenURL

  20. Raj A, Peskin CS, Tranchina D, Vargas DY, Tyagi S: Stochastic mRNA Synthesis in Mammalian Cells.

    PLoS Biol 2006, 4(10):e309. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  21. Warren L, Bryder D, Weissman IL, Quake SR: Transcription factor profiling in individual hematopoietic progenitors by digital RT-PCR.

    Proc Natl Acad Sci 2006, 103(47):17807-17812. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  22. Eissing T, Waldherr S, Allgöwer F, Scheurich P, Bullinger E: Steady state and (bi-)stability evaluation of simple protease signalling networks.

    BioSystems 2007, 90:591-601. PubMed Abstract | Publisher Full Text OpenURL

  23. Cherry JL, Adler FR: How to make a Biological Switch.

    J Theor Biol 2000, 203(2):117-133. PubMed Abstract | Publisher Full Text OpenURL

  24. van Kampen NG: Stochastic processes in physics and chemistry. North-Holland Amsterdam; 1981. OpenURL

  25. Gillespie DT: A rigorous derivation of the chemical master equation.

    Physica A: Statist Theor Phys 1992, 188(1-3):404-425. Publisher Full Text OpenURL

  26. Rao CV, Arkin AP: Stochastic chemical kinetics and the quasi-steady-state assumption: Application to the Gillespie algorithm.

    J Chem Phys 2003, 118:4999-5010. Publisher Full Text OpenURL

  27. Nuño J, Tarazona P: Lifetimes of small catalytic networks.

    Bull Math Biol 1994, 56(5):875-898. OpenURL

  28. Munsky B, Khammash M: The finite state projection algorithm for the solution of the chemical master equation.

    J Chem Phys 2006, 124(4):044104. PubMed Abstract | Publisher Full Text OpenURL

  29. Kaufman M, Soule C, Thomas R: A new necessary condition on interaction graphs for multistationarity.

    J Theor Biol 2007, 248(4):675-685. PubMed Abstract | Publisher Full Text OpenURL

  30. Parrott JA, Skinner MK: Thecal cell-granulosa cell interactions involve a positive feedback loop among keratinocyte growth factor, hepatocyte growth factor, and Kit ligand during ovarian follicular development.

    Endocrinol 1998, 139(5):2240-2245. Publisher Full Text OpenURL

  31. Yeh J, Adashi EY: The ovarian life cycle. In Reproductive Endocrinology. Edited by Yen SSC, Jaffe RB. Saunders Philadelphia; 1999:153-190. OpenURL

  32. Faddy MJ, Gosden RG: A model conforming the decline in follicle numbers to the age of menopause in women.

    Hum Reprod 1996, 11(7):1484-1486. PubMed Abstract | Publisher Full Text OpenURL

  33. Gougeon A, Ecochard R, Thalabard JC: Age-related changes of the population of human ovarian follicles: increase in the disappearance rate of non-growing and early-growing follicles in aging women.

    Biol Reprod 1994, 50(3):653-663. PubMed Abstract | Publisher Full Text OpenURL

  34. de Bruin JP, Bovenhuis H, van Noord PA, Pearson PL, van Arendonk JA, te Velde ER, Kuurman WW, Dorland M: The role of genetic factors in age at natural menopause.

    Hum Reprod 2001, 16(9):2014-2018. PubMed Abstract | Publisher Full Text OpenURL

  35. Nilsson E, Rogers N, Skinner MK: Actions of anti-Mullerian hormone on the ovarian transcriptome to inhibit primordial to primary follicle transition.

    Reproduction 2007, 134(2):209-221. PubMed Abstract | Publisher Full Text OpenURL

  36. Visser JA, Themmen APN: Anti-Müllerian hormone and folliculogenesis.

    Mol Cell Endocrinol 2005, 234(1-2):81-86. PubMed Abstract | Publisher Full Text OpenURL

  37. Faddy MJ, Gosden RG: A mathematical model of follicle dynamics in the human ovary.

    Hum Reprod 1995, 10(4):770-775. PubMed Abstract | Publisher Full Text OpenURL

  38. Silva-Buttkus PD, Marcelli G, Franks S, Stark J, Hardy K: Inferring biological mechanisms from spatial analysis: Prediction of a local inhibitor in the ovary.

    Proc Natl Acad Sci 2009, 106:456-461. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  39. Lancaster P, Tismenetsky M: The Theory of Matrices. San Diego: Academic Press; 1985. OpenURL

  40. Smith HL: Monotone Dynamical Systems. An Introduction to the Theory of Competitive and Cooperative Systems, Volume 41 of. Mathematical Surveys and Monographs. Providence, Rhode Island: American Mathematical Society; 1995. OpenURL

  41. Gillespie DT: Exact stochastic simulation of coupled chemical reactions.

    J Phys Chem 1977, 81(25):2340-2361. Publisher Full Text OpenURL

  42. Ramsey S, Orrell D, Bolouri H: Dizzy: Stochastic Simulation of Large-Scale Genetic Regulatory Networks.

    J Bioinform Comput Biol 2005, 3(2):415-436. PubMed Abstract | Publisher Full Text OpenURL