Abstract
Background
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.
Results
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.
Conclusions
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.
Background
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 lifespan are already present at birth or shortly thereafter as a population of socalled 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 [1416]. 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 cellcell 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.
Methods
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
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 Hilltype production rates v_{2 }and v_{4}. 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. V_{1 }= V_{2}, M_{1 }= M_{2 }and u_{1 }= u_{2}. The parameter values are within the physiological range for typical gene transcription processes. In particular, the degradation rate of 0.01 corresponds to a gene transcript halflife time of about 70 minutes. Typical transcript halflife 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 k_{1 }and corresponds to 3.3 transcripts that are produced per hour. The transcription rate upon maximal activation is given by V_{1,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).
Figure 1. Network schematic for the bistable switch model (1).
For twodimensional 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
Figure 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.
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 . The stochastic state of the switch at time t is given by the discrete probability distribution p(X, Y, t), which for each microstate gives the probability that the switch is in the microstate (X, Y )^{T }at time t:
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
for , where the reaction propensities v_{i}, 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).
Figure 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 · 10^{7 }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 p_{off }(t) and p_{on}(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. , and equivalently for p_{off}(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 continuoustime 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 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,
The master equation can be written as the linear ordinary differential equation
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 granulosaderived KITL activity and y represents thecaderived KGF and HGF activity. The reaction v_{1 }describes the influence on KITL expression of oocytederived bFGF, which is here assumed to be constant. The reactions v_{2 }and v_{4 }arise from the positive feedback interconnection, whereas v_{3 }and v_{5 }describe a constitutive degradation of KITL, KGF and HGF.
Figure 4. Hypothetical biochemical network for the primordial to primary transition in ovarian follicles.
The stochastic switch generates reliable longterm 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
where c_{1 }> 0, 0 <λ_{1 }< Re (λ_{i}) for i = 2, ..., n, and the c_{i}(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 simple exponential decay. For the biochemical parameter values given in Table 1, the numerical solution for p_{on}(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 halflife 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 N_{r}(t) can easily be characterized by a binomial distribution as
where N_{0 }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
The probability distribution P (N_{r}(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 N_{0 }= 10^{6 }was chosen from the reported range of ovarian follicles, 7 · 10^{5 }to 2 · 10^{6 }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[N_{r}(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 N_{d}. The depletion time T_{d }is defined as the smallest time t such that N_{r}(t) ≤ N_{d}, i.e. only N_{d }cells are remaining in the initial population. For the process of follicle growth initiation, we use N_{d }= 10^{3}, which has been considered to mark the onset of menopause [32].
Figure 5. Dynamical characteristics of the stochastic bistable switch on the single cell level and the population level. A: Probability of transformation event p_{on}(t) B: Population size probability distribution over time. C: Probability density function of the depletion time T_{d}.
The cumulative probability distribution function for the depletion time T_{d }is computed from the distribution obtained in (7) as
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[T_{d}] = 59.1 years and years, respectively.
The expected value for T_{d }can also be computed by solving . Using (6), it can thus be approximated by , 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 T_{d }= 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 . 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 AntiMü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 k_{1 }in the original model given in (1) from a constant parameter to the expression
where k_{1,max }is the maximal production rate of KITL, n_{2 }is the number of growing, AMH producing follicles, and K_{n }is an additional parameter. While follicle development is a complex process with many intermediate stages [31], in this analysis we use a simple twostate population model, where n_{1 }denotes the number of primordial follicles, and n_{2 }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 k_{1 }depending on n_{2 }as defined in (10), we obtain a dependency of λ_{1 }on n_{2}. 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 delaydifferential equations
where λ_{1}(n_{2}) is the decay rate computed from the transition matrix A(n_{2}) in (5), with k_{1}(n_{2}) as in (10). Using the parameters in Table 3, the population model given by (11) now predicts a depletion time of T_{d }= 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 K_{n }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.
Conclusions
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 continuoustime Markov process. We show that the proposed stochastic mechanism generates reliable longterm 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 longterm 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 laterstage 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 p_{on}(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
where A_{rev }∈ ℝ^{(n  1) × (n  1) }describes the interactions among the nonabsorbing states, and a_{abs }∈ ℝ^{1 × (n  1) }describes the transition propensities to the absorbing state.
Let us first derive some essential properties of the matrix A_{rev}. Since A is a stochastic matrix, we have
for i = 1, ..., n  1 i.e. A_{rev }is diagonally dominant. Thus, Gersgorin's theorem [39] asserts that all eigenvalues of A_{rev }have a nonpositive real part. Even more, since a_{abs }is nonzero, (13) holds with a strict inequality for at least one i. Thus, by Theorem 10.7.2 in [39], all eigenvalues of A_{rev }have negative real part. By the properties of the considered biochemical network, A_{rev }is irreducible, and its offdiagonal elements are nonnegative. From Corollary 4.3.2 in [40], it follows that A_{rev }has an eigenvalue λ_{1 }∈ ℝ with algebraic multiplicity 1 and a strictly positive corresponding eigenvector v_{1 }such that Re λ <λ_{1 }for all λ ≠ λ_{1 }in the spectrum of A_{rev}.
Denoting P_{rev }= (P_{1}, ..., P_{n1})^{T }we have . From the previously derived properties of the matrix A_{rev}, the general solution of this differential equation is given by
where are polynomials in t and is a constant coefficient, depending on the initial condition P_{rev}(0). The condition P_{rev}(t) ≥ 0 for all t implies that . For a nonnegative initial condition P_{rev}(0) with at least one positive element, we have . The transition probability p_{on}(t) is computed as
Acknowledgements
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.
References

Martins ML, Ferreira SC Jr, Vilela MJ: Multiscale models for biological systems.
Curr Opin Colloid Interface Sci 2010, 15:1823. Publisher Full Text

Ferrell JE, Xiong W: Bistability in cell signaling: How to make continuous processes discontinuous, and reversible processes irreversible.
Chaos 2001, 11:227236. PubMed Abstract  Publisher Full Text

Ferrell JE, Machleder EM: The biochemical basis of an allornone cell fate switch in Xenopus oocytes.
Science 1998, 280(5365):895898. PubMed Abstract  Publisher Full Text

Eissing T, Conzelmann H, Gilles ED, Allgöwer F, Bullinger E, Scheurich P: Bistability Analyses of a Caspase Activation Model for Receptorinduced Apoptosis.
J Biol Chem 2004, 279(35):3689297. PubMed Abstract  Publisher Full Text

Chickarmane V, Enver T, Peterson C: Computational modeling of the hematopoietic erythroidmyeloid switch reveals insights into cooperativity, priming, and irreversibility.
PLoS Comput Biol 2009, 5:e1000268. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

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:73049. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Kaufmann BB, Yang Q, Mettetal JT, van Oudenaarden A: Heritable Stochastic Switching Revealed by SingleCell Genealogy.
PLoS Biology 2007, 5(9):e239. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Tian T, Burrage K: Stochastic models for regulatory networks of the genetic toggle switch.
Proc Natl Acad Sci 2006, 103(22):83728377. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

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

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

Isaacs FJ, Hasty J, Cantor CR, Collins JJ: Prediction and measurement of an autoregulatory genetic module.
Proc Natl Acad Sci 2003, 100(13):77147719. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Fortune JE, Cushman RA, Wahl CM, Kito S: The primordial to primary follicle transition.
Mol Cell Endocrinol 2000, 163(12):5360. PubMed Abstract  Publisher Full Text

Skinner MK: Regulation of primordial follicle assembly and development.
Hum Reprod Update 2005, 11(5):461471. PubMed Abstract  Publisher Full Text

Parrott JA, Skinner MK: Kitligand/stem cell factor induces primordial follicle development and initiates folliculogenesis.
Endocrinol 1999, 140(9):42624271. Publisher Full Text

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):215218. PubMed Abstract  Publisher Full Text

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(12):1925. PubMed Abstract  Publisher Full Text

BrawTal R: The initiation of follicle growth: the oocyte or the somatic cells?
Mol Cell Endocrinol 2002, 187(12):1118. PubMed Abstract  Publisher Full Text

Ross J: mRNA stability in mammalian cells.
Microbiol Rev 1995, 59(3):423450. PubMed Abstract  PubMed Central Full Text

Hargrove JL, Hulsey MG, Beale EG: The kinetics of mammalian gene expression.
Bioessays 1991, 13(12):667674. PubMed Abstract  Publisher Full Text

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

Warren L, Bryder D, Weissman IL, Quake SR: Transcription factor profiling in individual hematopoietic progenitors by digital RTPCR.
Proc Natl Acad Sci 2006, 103(47):1780717812. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

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:591601. PubMed Abstract  Publisher Full Text

Cherry JL, Adler FR: How to make a Biological Switch.
J Theor Biol 2000, 203(2):117133. PubMed Abstract  Publisher Full Text

van Kampen NG: Stochastic processes in physics and chemistry. NorthHolland Amsterdam; 1981.

Gillespie DT: A rigorous derivation of the chemical master equation.
Physica A: Statist Theor Phys 1992, 188(13):404425. Publisher Full Text

Rao CV, Arkin AP: Stochastic chemical kinetics and the quasisteadystate assumption: Application to the Gillespie algorithm.
J Chem Phys 2003, 118:49995010. Publisher Full Text

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

Kaufman M, Soule C, Thomas R: A new necessary condition on interaction graphs for multistationarity.
J Theor Biol 2007, 248(4):675685. PubMed Abstract  Publisher Full Text

Parrott JA, Skinner MK: Thecal cellgranulosa 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):22402245. Publisher Full Text

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

Faddy MJ, Gosden RG: A model conforming the decline in follicle numbers to the age of menopause in women.
Hum Reprod 1996, 11(7):14841486. PubMed Abstract  Publisher Full Text

Gougeon A, Ecochard R, Thalabard JC: Agerelated changes of the population of human ovarian follicles: increase in the disappearance rate of nongrowing and earlygrowing follicles in aging women.
Biol Reprod 1994, 50(3):653663. PubMed Abstract  Publisher Full Text

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):20142018. PubMed Abstract  Publisher Full Text

Nilsson E, Rogers N, Skinner MK: Actions of antiMullerian hormone on the ovarian transcriptome to inhibit primordial to primary follicle transition.
Reproduction 2007, 134(2):209221. PubMed Abstract  Publisher Full Text

Visser JA, Themmen APN: AntiMüllerian hormone and folliculogenesis.
Mol Cell Endocrinol 2005, 234(12):8186. PubMed Abstract  Publisher Full Text

Faddy MJ, Gosden RG: A mathematical model of follicle dynamics in the human ovary.
Hum Reprod 1995, 10(4):770775. PubMed Abstract  Publisher Full Text

SilvaButtkus 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:456461. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

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

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.

Gillespie DT: Exact stochastic simulation of coupled chemical reactions.
J Phys Chem 1977, 81(25):23402361. Publisher Full Text

Ramsey S, Orrell D, Bolouri H: Dizzy: Stochastic Simulation of LargeScale Genetic Regulatory Networks.
J Bioinform Comput Biol 2005, 3(2):415436. PubMed Abstract  Publisher Full Text