Stochastic hybrid systems for studying biochemical processes

Abhyudai Singh, João P. Hespanha

Abstract

Many protein and mRNA species occur at low molecular counts within cells, and hence are subject to large stochastic fluctuations in copy numbers over time. Development of computationally tractable frameworks for modelling stochastic fluctuations in population counts is essential to understand how noise at the cellular level affects biological function and phenotype. We show that stochastic hybrid systems (SHSs) provide a convenient framework for modelling the time evolution of population counts of different chemical species involved in a set of biochemical reactions. We illustrate recently developed techniques that allow fast computations of the statistical moments of the population count, without having to run computationally expensive Monte Carlo simulations of the biochemical reactions. Finally, we review different examples from the literature that illustrate the benefits of using SHSs for modelling biochemical processes.

1. Introduction

Deterministic hybrid systems that integrate continuous dynamics with discrete events have been used to model a wide array of biological processes that exhibit switching behaviour (Belta et al. 2004; Ghosh & Tomlin 2004; Lincoln & Tiwari 2004; Tanaka et al. 2008). However, deterministic frameworks often fail to capture biochemical processes within living cells where low population counts of mRNAs and proteins can set the stage for significant stochastic effects (Raj & van Oudenaarden 2008). The inherent stochastic nature of cellular processes has motivated the use of stochastic hybrid systems (SHSs) for modelling biological phenomena at the single-cell level (Hu et al. 2004; Julius et al. 2008; Lygeros et al. 2008). SHSs (formally defined in §2) combine the generality of hybrid systems with stochastic processes.

Biochemical reactions inside cells are often modelled using a stochastic formulation, which takes into account the inherent randomness of thermal molecular motion (Gillespie 1976). In the stochastic formulation of chemical kinetics, reactions are treated as probabilistic events that change the population counts of individual chemical species based on the stoichiometry of the reactions. We show in §3 that the time evolution of the number of molecules of different chemical species involved in a set of chemical reactions can be modelled using an SHS. We use this SHS formalism to compute the time derivatives of the lower order statistical moments (for example, means, standard deviations, correlation etc.) of the population count. It turns out that the differential equations that describe the time evolution of the lower order statistical moments are generally ‘not closed’, in the sense that the right-hand sides of these equations depend on higher order moments. In §4, we illustrate recently developed techniques that ‘close’ these moment equations by expressing high order moments as nonlinear functions of lower order moments. The resulting closed differential equations provide quick and efficient computations of the statistical moments for a given chemical reaction network, without having to use Monte Carlo simulations that come at a significant computational cost.

Finally, in §5, we review different examples from the literature that illustrate how SHSs have been used to model uncertainty in different biochemical processes. These examples include: (i) modelling random gene activation/inactivation in gene regulatory networks, (ii) modelling the effects of perturbations on stochastic chemical kinetics, and (iii) modelling noise-induced transitions between different steady states of a multi-stable biological network.

2. Stochastic hybrid systems

The state space of an SHS is composed of a continuous component x that takes values in Euclidean space Embedded Image and a discrete component q that takes values in a finite set {q1,…,qN}. The continuous state evolves according to the ordinary differential equation (ODE) Embedded Image2.1 where the vector field f depends on the discrete state of the SHS. During the evolution of the above ODE, stochastic transitions or jumps may occur which change both the discrete state and the continuous state of the SHS. More specifically, these random transitions are characterized by a family of k reset maps Embedded Image2.2 and a corresponding family of k transition intensities Embedded Image2.3 In essence, between stochastic transitions, the discrete state remains constant, whereas the continuous state flows according to equation (2.1). At transition times, the continuous and discrete states are reset according to equation (2.2). The frequency with which different transitions occur is determined by the transition intensities (2.3). In particular, the probability that a particular transition will occur in an ‘elementary interval’ (t,t + dt], and therefore the corresponding reset map will be ‘activated’ is given by λi(q(t),x(t),t)dt. Equations (2.1)–(2.3) define an SHS, which is often conveniently represented by a directed graph as in figure 1. We refer interested readers to Hespanha (2005) for a mathematically precise characterization of an SHS and an algorithm to run Monte Carlo simulations of its trajectories.

Figure 1.

Graphical representation of an SHS. Each vertex corresponds to a discrete mode and each edge to a transition between discrete modes. The vertices are labelled with the corresponding discrete mode and the vector fields that determine the evolution of the continuous state in that particular mode. The edges are labelled with the transition intensities and their corresponding reset maps.

Since the time evolution of the continuous state between stochastic transitions is deterministic, these SHSs have often been referred to in the literature as piecewise deterministic Markov processes (Davis 1993). However, equation (2.1) can be modified to allow the continuous state x to evolve according to a stochastic differential equation (SDE) rather than an ODE (Hespanha & Singh 2005). Another straightforward extension of the above SHS is to include deterministic transitions between discrete states, where a transition is triggered when a certain ‘guard condition’ is satisfied (Hespanha 2005). We show next that SHSs are convenient for modelling the temporal dynamics of population counts of different species involved in a set of chemical reactions.

3. Stochastic modelling of chemical reactions

We begin this section by reviewing the stochastic formulation of chemical kinetics.

(a) Stochastic formulation of chemical kinetics

Consider a spatially uniform mixture of n chemical species X1, X2,…,Xn in a fixed volume V involved in a system of k reactions R1, R2,…,Rk of the form Embedded Image where uij is the stoichiometry associated with the jth reactant of the ith reaction and vij is the stoichiometry associated with the jth product of the ith reaction. In the following, we denote by xi(t) the number of molecules of the species Xi at time t.

At high population counts, the time evolution of x = [x1,…,xn]T can be treated as a continuous and deterministic process governed by an ODE, often referred to in the literature as chemical rate equations (Wilkinson 1980). However, this deterministic framework fails within single cells where many species occur at very low molecular counts and change by discrete integer amounts whenever a reaction occurs inside the cell. The time evolution of such low-copy biochemical species is more accurately represented by a stochastic formulation of chemical kinetics which treats x(t) as a stochastic process (McQuarrie 1967).

In the stochastic formulation of chemical kinetics, each reaction is a probabilistic event and is assigned a probability that it will occur in the next ‘infinitesimal’ time interval (t,t + dt]. This probability is given by the propensity function of the reaction, which is a product of the following two terms.

  • — The number hi(x) of distinct molecular reactant combinations for the reaction Ri present in the volume V at time t.

  • — The probability ci dt that a particular reactant combination of Ri will actually react on (t,t + dt]. The constant ci for each chemical reaction depends on the physical properties of the reacting molecules and the temperature of the system and is typically experimentally determined.

Table 1 shows the form of the propensity function for different reaction types (Gillespie 1976). In summary, the stochastic formulation treats reactions as a set of stochastic channels, and whenever a particular channel ‘fires’ the molecular counts change based on the stoichiometry of that reaction. Moreover, the frequencies at which these channels ‘fire’ are determined by the propensity functions of the individual reactions.

View this table:
Table 1.

Propensity functions for different reaction types.

(b) Representing chemical reactions as a stochastic hybrid system

The evolution of the number of molecules x1,x2,…,xn can be generated by an SHS. Since the number of molecules takes values in the discrete set of integers, they can be regarded as part of either the discrete or the continuous state of the SHS. However, it will turn out to be more convenient to view them as part of a continuous state. In this case, the SHS has a single discrete mode, which we omit for simplicity. The continuous state of the SHS consists of the vector x = [x1,…,xn]T and has trivial dynamics, that is, Embedded Image3.1 Each of the reactions is represented using a reset map defined by the stoichiometry Embedded Image3.2 and a corresponding transition intensity Embedded Image3.3 given by the propensity function of the reaction. Thus, between reactions, the population count remains constant and whenever the ith reaction ‘fires’, the state x is reset according to equation (3.2); furthermore, the probability of the activation taking place in an ‘infinitesimal’ time interval (t,t + dt] is given by λi(x)dt.

In many cases, biochemical reactions can be divided into subsystems of fast and slow reactions. For example, the binding and unbinding of a transcription factor to a promoter typically occurs at much faster time scales than the process of transcription, which creates mRNAs from DNA. When such differences in time scale exist between reactions, slow reactions or reactions containing low-copy molecular species should be modelled using stochastic transitions and resets, as in the stochastic formulation of chemical kinetics. However, fast reactions or reactions containing high-copy molecular species can be modelled using ODEs (or in some cases SDEs) resulting in a reduced approximate SHS where the dynamics of the continuous state is no longer trivial (figure 2). These reduced SHSs provide much faster simulation times than the original SHS, with only a marginal decrease in accuracy (Neogi 2004; Salis & Kaznessis 2005; Chen et al. 2009) and have been used to model a wide array of biological processes including lactose regulation in Escherichia coli (Julius et al. 2008), the human immunodeficiency virus transactivation network (Griffith et al. 2006) and synthetic gene networks (Bortolussi & Policriti 2008).

Figure 2.

SHS representation of chemical reactions, where reactions are modelled as stochastic transitions that reset the population count x based on the stoichiometry of the reaction (a). Note that this SHS representation of the chemical reaction network is a special case of the SHS given by equations (2.1)–(2.3) as it has only a single discrete mode, which we omit for simplicity. To reduce the computational costs of simulating these SHSs, fast reactions or reactions containing high-copy molecular species (solid arrows) are often modelled using differential equations (b).

4. Computing moment dynamics of a chemically reacting system

The stochastic formulation of chemical reactions permits the computation of the probability density function of the population count x(t), which is often done through various Monte Carlo techniques at a significant computational cost (Gillespie 1976, 2001). Since one is often interested in computing only the first- and second-order moments for the number of molecules of the different species involved, much time and effort can be saved by applying approximate methods to directly compute these low-order moments (for example, Van Kampen’s (2001) linear noise approximation), without actually having to solve for the probability density function. In this section, we describe a procedure to compute the time evolution of the statistical moments of x(t) for an arbitrary set of chemical reactions.

Given a vector {m1,m2,…,mn} of n non-negative integers, we define the uncentred moment of x = [x1,…,xn]T to be Embedded Image4.1 where E stands for the expected value. We refer to the sum Embedded Image as the order of the moment. For example, consider a system of reactions with two species with population counts x1 and x2. Then, the first-order moments are given by Embedded Image4.2 the second-order moments are given by Embedded Image4.3 and so on. The SHS formalism for chemical reactions introduced in the previous section allows a straightforward derivation of the moment dynamics using Dynkin’s equation (Davis 1993). Applying the Dynkin equation to the SHS given by equations (3.1)–(3.3) gives the following time derivative of an uncentred moment of the population count (Singh & Hespanha 2006): Embedded Image4.4 The right-hand side of the above equation can be interpreted as the expected value of the product of the change in the monomial Embedded Image whenever a reaction occurs and the frequency with which reaction occurs, summed up over all chemical reactions. Since the propensity functions cihi(x) are polynomials in the population count x (table 1), it follows from equation (4.4) that the time derivative of an uncentred moment Embedded Image will be a linear combination of uncentred moments of x(t). Note that this result will also hold for SHSs where the continuous dynamics is non-trivial, i.e. Embedded Image (as in figure 2), as long as the vector field f is a polynomial in x (see Hespanha & Singh 2005). Thus, for chemically reacting systems, uncentred moments of the population count evolve according to a linear system of equations. However, these moment equations may not always be ‘closed’ in the sense that the time derivative of an mth-order moment may depend on moments of order higher than m. We discuss techniques used for solving such a system of equations below, but we first consider a class of chemical reactions where the moment dynamics is always ‘closed’.

(a) Linear system of chemical reactions

Consider a system of chemical reactions where all reactions have linear propensity functions cihi(x). This implies from table 1 that we only have reactions of the form Embedded Image or XiXj. In this case, it follows from equation (4.4) that the time derivative of an mth-order moment Embedded Image is a linear combination of moments of x of order up to m. Hence, if we construct a vector μ consisting of the first M order moments of x, then its time evolution is given by Embedded Image4.5 for some appropriate constant vector Embedded Image and constant matrix A. Thus, for a chemically reacting system with linear propensity functions, the moments of the population count can always be computed by solving equation (4.5). We illustrate this point with a stochastic model of gene expression, which is given by the following set of chemical reactions: Embedded Image4.6

The first two reactions represent mRNA transcription from DNA at a rate c1 and mRNA degradation at a constant rate c2. The last two reactions correspond to protein translation from the mRNA and protein degradation at rates c3 and c4, respectively. Let x1 and x2 denote the population count of the mRNA and the protein, respectively. Then, the time evolution of x = [x1,x2]T can be represented by an SHS with trivial continuous dynamics, four reset maps Embedded Image4.7a Embedded Image4.7b and corresponding transition intensities λ1(x) = c1, λ2(x) = c2x1, λ3(x) = c3x1 and λ4(x) = c4x2. From equation (4.4), the time evolutions of the first- and second-order moments of the mRNA and the protein count are given by the following system of linear equations: Embedded Image4.8 A steady-state analysis of the above moment dynamics shows that the steady-state variance σ2 of protein levels is given by Embedded Image and scales linearly with the steady-state value of the mean protein level Embedded Image. This linear scaling of variance with the mean protein levels obtained from the stochastic model is consistent with experimental measurements of gene expression noise in both eukaryotes (Bar-Even et al. 2006; Newman et al. 2006) and prokaryotes (Ozbudak et al. 2002).

(b) Nonlinear system of chemical reactions

We next consider the scenario where the system of chemical reactions contains at least one reaction with a nonlinear propensity function. Then, the time derivative of the vector μ consisting of the first M-order moments of x is given by Embedded Image4.9 for an appropriate constant vector Embedded Image, constant matrices A and B and a (time-varying) vector Embedded Image containing moments of order M + 1 and higher. For example, consider the following reactions: Embedded Image4.10 where species X is produced as a monomer at a constant rate c1 and dimerizes to form Y . The dimer Y then decays at a constant rate c3. Let x1 and x2 denote the population count of species X and Y , respectively. Then, the time evolution of x = [x1,x2]T can be represented by an SHS with trivial continuous dynamics, three reset maps Embedded Image4.11 and corresponding transition intensities λ1(x) = c1, λ2(x) = (c2/2)x1(x1 − 1) and λ3(x) = c3x2. Using Dynkin’s equation for the above SHS, the time evolution of the first- and second-order moments of the population count is given by equation (4.9) where Embedded Image4.12 and is dependent on the third-order moments of the population count. This example illustrates the general principle that nonlinear propensity functions result in a ‘non-closed’ system of moment equations, where the dynamics of the lower order moments depends on higher order moments. For analysis purposes, the time evolution of the vector μ is often made to be closed by approximating the higher order moments Embedded Image as nonlinear functions of lower order moments in μ, as in Embedded Image. This procedure is referred to in the literature as moment closure (Nasell 2003; Gillespie 2009) and results in a nonlinear approximated moment dynamics given by Embedded Image4.13 where the state of this closed system ν(t) can be viewed as an approximation for μ(t).

Recent work has proposed a novel moment closure technique based on derivative matching, where the moment closure is done by matching time derivatives of the exact (not closed) moment equations with that of the approximate (closed) moment equations at some initial time t0 and set of initial conditions (Hespanha & Singh 2005). In particular, this derivative-matching approach attempts to determine nonlinear functions φ for which Embedded Image4.14 holds for deterministic initial conditions x(t0) = x0 with probability one. The main rationale for doing so is that, if a sufficiently large number of derivatives of μ(t) and ν(t) match point-wise at an initial time t0, then from a Taylor series argument the trajectories of μ(t) and ν(t) will remain close at least locally in time. Singh & Hespanha (2006, 2007) provide explicit formulas to construct a class of functions φ for which equation (4.14) holds approximately for all i ≥ 1, i.e. all time derivatives of ν(t) and μ(t) match at t = t0 with small errors. Table 2 provides the nonlinear approximations for all possible third-order moments as a function of the first- and second-order moments based on the above derivative-matching moment closure technique. We close the dynamics of the first- and second-order moments for the reaction set of equation (4.10) by using the corresponding nonlinear approximations for the third-order moments Embedded Image and Embedded Image. Numerical solutions of the resulting approximated moment dynamics are shown in figure 3. The procedure described here to generate approximated moment dynamics can be fully automated. The software StochDynTools (Hespanha 2006) is available to compute the approximated moment dynamics starting from a simple ASCII description of the chemical reactions involved.

View this table:
Table 2.

Moment closure approximation for third-order moments.

Figure 3.

(a,b) The time evolution of the means (solid lines)±one s.d. (dashed lines) obtained from the approximated moment dynamics (4.13) corresponding to M = 2 (dark blue) and M = 3 (red) for the reaction set (4.10). A sample Monte Carlo run is shown in light blue. The means ± one s.d. at the final time are X = 2.094±1.021 (M = 2, dark blue); X = 2.036±1.107 (M = 3, red); Y = 4.592±1.932 (M = 2, dark blue); Y = 4.597±1.936 (M = 3, red). (c,d) The distribution of the molecular counts of species X and Y at the final time obtained from 20 000 Monte Carlo simulations done using the stochastic simulation algorithm (Gillespie 1976). The means ± one s.d. obtained from Monte Carlo simulations are X = 1.956±1.188 and Y = 4.583±1.937. Note that as one increases M, there is an improvement in the moment estimates for species X. However, there is no significant change in the moment estimates of Y , and the moment trajectories corresponding to M = 2 and M = 3 lie on top of each other. Reaction parameters were taken as c1 = 100 molecules h−1, c2 = 30 molecules−1 h−1 and c3 = 20 h−1.

A striking feature of the above moment closure technique is that its accuracy can be arbitrarily increased by reducing the error in matching the time derivative between the exact and the approximate moment dynamics. More specifically, if we close the dynamics of the first M-order moments, then the derivative matching error given by Embedded Image4.15 scales as ||x0||−M (Singh & Hespanha 2006). Thus, by increasing M, which corresponds to including higher order moments in the vector μ, the approximated moment dynamics (4.13) provides more accurate approximations to the exact moment dynamics (4.9), as long as the elements of x0 are larger than one.

Another striking feature of the above moment closure technique is that the nonlinear functions φ, which express high order moments as functions of lower order moments, are consistent with lognormal distributions. Singh & Hespanha (in press) perform a systematic comparison of the derivative matching moment closure technique to an alternative procedure, where moment closure is performed by setting the third and higher order cumulants of x equal to zero (Goutsias 2007; Gomez-Uribe & Verghese 2007; Lee et al. 2009). As for a Gaussian distribution all cumulants of order three and higher are equal to zero, the zero-cumulant moment closure is consistent with Gaussian distributions. Singh & Hespanha (in press) show that at low population counts the derivative-matching moment-closure technique provides more accurate estimates of the lower order moments than the zero-cumulant moment closure. Intuitively, this occurs because at low population regimes molecular counts have highly skewed distributions that are much closer to a lognormal distribution than to a Gaussian distribution. In fact, using the zero-cumulant moment-closure technique for low population species often results in unstable approximated moment dynamics with unbounded solutions (Nasell 2003).

In summary, moment closure allows one to approximate the moment dynamics of a chemically reacting system by an ODE, which provides fast and efficient computation of stochasticity in a given reaction network.

5. Modelling gene regulatory networks and other biological processes

In this section, we briefly review how SHSs have been used to model uncertainties arising from different sources in biochemical processes.

(a) Modelling gene regulatory networks

Gene regulatory networks consist of a collection of genes that regulate the transcriptional activity of each other through their expressed proteins. SHSs have frequently been used to model the uncertainties associated with activation/inactivation of a gene in response to binding/unbinding of protein complexes to its promoter. We illustrate this with the simplest possible network: an autoregulatory gene network where a protein inhibits or activates its own gene expression (figure 4a). Zeiser et al. (2009) models autoregulatory gene networks as an SHS with two discrete states, which represent a gene in an ‘ON’ or ‘OFF’ state (figure 4b). The protein count represents the continuous state of the SHS and evolves according to a linear differential equation with production and degradation in the ‘ON’ state, and only degradation in the ‘OFF’ state. For simplicity, the authors combine transcription and translation into a single protein production term and do not explicitly consider the mRNA dynamics. Random gene activation/inactivation is modelled through stochastic transitions between the discrete states with frequencies that are dependent on the protein count. In this SHS model, noise in protein levels only comes from stochastic promoter transitions between different transcriptional states, which is likely to be true when these transitions occur at much slower time scales than those of the protein production and decay.

Figure 4.

(a) An autoregulatory gene network where a protein regulates its own expression. (b) An SHS model of the autoregulatory gene network with two discrete states representing gene ‘ON’ and ‘OFF’ states. The protein count x exponentially increases to c1/c2 in the gene ‘ON’ state and decreases exponentially to zero in the gene ‘OFF’ state. Transitions between transcriptionally ‘ON’ and ‘OFF’ states occur with transition intensities that are dependent on the protein count. As there are no changes in the protein population when the gene turns ‘ON’ or ‘OFF’, the reset maps for this SHS are the identity map. (c) A singe realization of the protein count which shows increase and decrease in protein levels as the SHS transitions between ‘ON’ and ‘OFF’ states (Zeiser et al. 2009).

Mathematical models of autoregulatory gene networks have been extensively used to study negative feedback regulation, where a protein inhibits its own expression. Such negative feedback loops occur commonly in many cellular genes (Alon 2007) and have been hypothesized as a mechanism to reduce stochastic fluctuations in protein levels (Becskei & Serrano 2000). A negative feedback can be easily implemented in the above SHS model by assuming that the promoter is more likely to transition to the OFF state if the protein count increases within the cell. Analysis of these models not only predicts conditions under which feedback will provide the best suppression of gene expression noise but also determines the fundamental limits of noise suppression possible through negative autoregulation (Singh & Hespanha 2009a,b). Counterintuitively, these models also show that, in some cases, introducing a negative feedback may actually increase gene expression noise rather than decreasing it (Stekel & Jenkins 2008; Zeiser et al. 2009).

The autoregulatory gene network model presented above can be easily extended to a network of N genes. Assuming each gene can be either ‘ON’ or ‘OFF’, then the SHS will have 2N discrete states with each discrete state corresponding to some set of genes being in the ‘ON’ state and others being in the ‘OFF’ state. As in figure 4b, the protein population either exponentially grows or decays depending on the transcriptional status of its gene. SHS models of gene networks, in which genes stochastically transition between transcriptional states and protein counts evolve according to linear differential equations, have proven to be very useful for both parameter identification and modelling of subtilin production in Bacillus subtilis (Cinquemani et al. 2008a) and nutrient stress response in E. coli (Cinquemani et al. 2008b). An important caveat of using these models is that many mRNA species occur at very low molecular counts within cells (Bar-Even et al. 2006). Thus, by modelling transcription as a completely deterministic process, one may fail to capture the stochasticity in protein levels due to thermal fluctuations in the corresponding mRNA counts. However, this limitation can be obviated by incorporating mRNA production and degradation as stochastic events, while still modelling protein translation and degradation as ODEs.

(b) Modelling perturbations of stochastic chemical kinetics

The SHS framework for chemical reacting systems introduced above provides a convenient model to investigate perturbations to the system caused by discrete events (figure 5). This point is best illustrated by the recent work of Riley et al. (2009) that investigates the effect of drug treatment on the sugar cataract development process. This process can be expressed as a set of biochemical reactions, and hence modelled as an SHS with a single discrete state (figure 2). To study the effects of drug treatment, the authors introduce a new discrete state, which models the time evolution of the biochemical reactions in the presence of the drug. The transitions between the two discrete states of the SHS are determined by the criterion used to add or remove the drug from the system. Using stochastic reachability analysis methods for SHSs, the authors were able to compute the probability that a patient will develop a cataract with no drug treatment and with an intermittent drug treatment.

Figure 5.

Modelling perturbations of chemical kinetics using SHSs. The discrete mode q1 represents a set of chemical reactions where some reactions are modelled deterministically as differential equations while others are modelled stochastically through transitional intensities and reset maps. The discrete mode q2 represents the same set of reactions under a perturbation (such as additions of a drug) which can alter the vector field f and the transition intensities/reset maps. Drug addition and removal are modelled via stochastic transitions between the discrete modes which are governed by the transition intensities λ3(q1,x) and λ4(q2,x) and their corresponding reset maps. SHSs such as the one illustrated here have recently been used to study disease progression under intermittent drug treatment (see Riley et al. 2009).

(c) Modelling complex dynamics with multi-stability

SHSs have also been used to model random transitions between different stable states of a multi-stable biological network. These transitions occur due to noise in protein levels, which causes the system to move from the region of attraction of one stable state to another. Noise-induced transitions between alternative stable states play a key role in mediating cell-fate decisions in stem cells (Losick & Desplan 2008) and viruses (Singh & Weinberger 2009). A well-known example of this is the lactose regulation system in E. coli, where, in response to external lactose, a colony of identical cells bifurcates into two distinct populations: either fully induced with high lactose metabolizing enzyme levels or uninduced with no enzymes (Novick & Weiner 1957). This bifurcation at the population level occurs owing to an underlying bi-stability in the lactose metabolic network, with stochastic fluctuations in regulatory molecules causing single cells to converge to either one of the two stable steady states. Starting from a full SHS model of the lactose metabolic network, Julius et al. (2008) built a reduced model that can quickly predict the fraction of induced and uninduced cells in response to a given concentration of lactose. This reduced model allowed the design of control feedback laws that can robustly steer a colony of cells to a desired fraction of induced cell, using external lactose as a control input. Perhaps, one day, similar reduced SHS models of multi-stable networks underlying cell-fate decisions in stem cells can be controlled for possible therapeutic benefit.

6. Conclusions

Intracellular processes are driven by reactant molecules randomly diffusing and colliding within the cell and are thus inherently stochastic. We showed that the time evolution of the number of molecules of different species involved in an arbitrary system of biochemical reactions can be modelled as an SHS. The SHS framework is useful as it provides the flexibility of modelling fast reactions deterministically, using differential equations, while slow reactions can evolve stochastically using stochastic transitions that reset the population count based on the stoichiometry of the reaction (figure 2).

The SHS framework also provides a method to compute the moment dynamics for a given set of chemical reactions. We illustrated a novel moment-closure scheme based on derivative matching, which allows one to approximate the time evolution of the statistical moments by a nonlinear ODE. These equations not only provide quick computations of the means, standard deviation, correlation etc., but can also be used for investigating how stochasticity in molecular counts is affected by the parameters of the reaction network. Comparisons with alternative techniques showed that this moment-closure scheme performs best at low population counts, where stochastic effects are most likely to be manifested.

Given the recent evidence that proteins involved in various cellular pathways exhibit stochastic fluctuations in their copy numbers (Bar-Even et al. 2006), SHSs will likely find increased use to answer fundamental questions in noise biology. For example, what mechanisms ensure that, in spite of noise in protein levels, signal transmission and information processing occur with sufficiently high fidelity inside cells? How is noise at the cellular level exploited to create variability at the population level for dealing with environmental uncertainties? Finally, SHS models will also likely play an important role in parameter identification as has been illustrated by Cinquemani et al. (2008b). Models that take into account the inherent stochastic nature of biochemical reactions will likely provide better parameter estimates from biological data than purely deterministic models that consider all variation as measurement noise.

Acknowledgements

J.P.H. acknowledges funding from the Institute for Collaborative Biotechnologies through grant DAAD19-03-D-0004 from the US Army Research Office and by the National Science Foundation under grant numbers ECCS-0725485 and ECCS-0835847.

    References

    View Abstract