## 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 and a *discrete component* **q** that takes values in a finite set {*q*_{1},…,*q*_{N}}. The continuous state evolves according to the ordinary differential equation (ODE)
2.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*
2.2
and a corresponding family of *k* *transition intensities*
2.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* + d*t*], and therefore the corresponding reset map will be ‘activated’ is given by *λ*_{i}(**q**(*t*),**x**(*t*),*t*)d*t*. 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.

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 *X*_{1}, *X*_{2},…,*X*_{n} in a fixed volume *V* involved in a system of *k* reactions *R*_{1}, *R*_{2},…,*R*_{k} of the form
where *u*_{ij} is the stoichiometry associated with the *j*th reactant of the *i*th reaction and *v*_{ij} is the stoichiometry associated with the *j*th product of the *i*th reaction. In the following, we denote by *x*_{i}(*t*) the number of molecules of the species *X*_{i} at time *t*.

At high population counts, the time evolution of **x** = [**x**_{1},…,**x**_{n}]^{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* + d*t*]. This probability is given by the propensity function of the reaction, which is a product of the following two terms.

— The number

*h*_{i}(**x**) of distinct molecular reactant combinations for the reaction*R*_{i}present in the volume*V*at time*t*.— The probability

*c*_{i}d*t*that a particular reactant combination of*R*_{i}will actually react on (*t*,*t*+ d*t*]. The constant*c*_{i}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.

### (b) Representing chemical reactions as a stochastic hybrid system

The evolution of the number of molecules **x**_{1},**x**_{2},…,**x**_{n} 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** = [**x**_{1},…,**x**_{n}]^{T} and has trivial dynamics, that is,
3.1
Each of the reactions is represented using a reset map defined by the stoichiometry
3.2
and a corresponding transition intensity
3.3
given by the propensity function of the reaction. Thus, between reactions, the population count remains constant and whenever the *i*th 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* + d*t*] is given by *λ*_{i}(**x**)d*t*.

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).

## 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 {*m*_{1},*m*_{2},…,*m*_{n}} of *n* non-negative integers, we define the *uncentred* moment of **x** = [**x**_{1},…,**x**_{n}]^{T} to be
4.1
where *E* stands for the expected value. We refer to the sum as the *order of the moment*. For example, consider a system of reactions with two species with population counts **x**_{1} and **x**_{2}. Then, the first-order moments are given by
4.2
the second-order moments are given by
4.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):
4.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 whenever a reaction occurs and the frequency with which reaction occurs, summed up over all chemical reactions. Since the propensity functions *c*_{i}*h*_{i}(**x**) are polynomials in the population count **x** (table 1), it follows from equation (4.4) that the time derivative of an uncentred moment 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. (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 *m*th-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 *c*_{i}*h*_{i}(**x**). This implies from table 1 that we only have reactions of the form or *X*_{i} → *X*_{j}. In this case, it follows from equation (4.4) that the time derivative of an *m*th-order moment 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
4.5
for some appropriate constant vector 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:
4.6

The first two reactions represent mRNA transcription from DNA at a rate *c*_{1} and mRNA degradation at a constant rate *c*_{2}. The last two reactions correspond to protein translation from the mRNA and protein degradation at rates *c*_{3} and *c*_{4}, respectively. Let **x**_{1} and **x**_{2} denote the population count of the mRNA and the protein, respectively. Then, the time evolution of **x** = [**x**_{1},**x**_{2}]^{T} can be represented by an SHS with trivial continuous dynamics, four reset maps
4.7a
4.7b
and corresponding transition intensities *λ*_{1}(**x**) = *c*_{1}, *λ*_{2}(**x**) = *c*_{2}**x**_{1}, *λ*_{3}(**x**) = *c*_{3}**x**_{1} and *λ*_{4}(**x**) = *c*_{4}**x**_{2}. 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:
4.8
A steady-state analysis of the above moment dynamics shows that the steady-state variance *σ*^{2} of protein levels is given by and scales linearly with the steady-state value of the mean protein level . 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
4.9
for an appropriate constant vector , constant matrices *A* and *B* and a (time-varying) vector containing moments of order **M** + 1 and higher. For example, consider the following reactions:
4.10
where species *X* is produced as a monomer at a constant rate *c*_{1} and dimerizes to form *Y* . The dimer *Y* then decays at a constant rate *c*_{3}. Let **x**_{1} and **x**_{2} denote the population count of species *X* and *Y* , respectively. Then, the time evolution of **x** = [**x**_{1},**x**_{2}]^{T} can be represented by an SHS with trivial continuous dynamics, three reset maps
4.11
and corresponding transition intensities *λ*_{1}(**x**) = *c*_{1}, *λ*_{2}(**x**) = (*c*_{2}/2)**x**_{1}(**x**_{1} − 1) and *λ*_{3}(**x**) = *c*_{3}**x**_{2}. 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
4.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 as nonlinear functions of lower order moments in *μ*, as in . 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
4.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 *t*_{0} and set of initial conditions (Hespanha & Singh 2005). In particular, this derivative-matching approach attempts to determine nonlinear functions *φ* for which
4.14
holds for deterministic initial conditions **x**(*t*_{0}) = **x**_{0} 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 *t*_{0}, 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* = *t*_{0} 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 and . 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.

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
4.15
scales as ||**x**_{0}||^{−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 **x**_{0} 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 4*a*). 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 4*b*). 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.

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 2009*a*,*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 2^{N} 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 4*b*, 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.* 2008*a*) and nutrient stress response in *E. coli* (Cinquemani *et al.* 2008*b*). 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.

### (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.* (2008*b*). 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.

## Footnotes

One contribution of 12 to a Theme Issue ‘Theory of hybrid dynamical systems and its applications to biological and medical systems’.

- © 2010 The Royal Society