## Abstract

It has recently been shown that structural conditions on the reaction network, rather than a ‘fine-tuning’ of system parameters, often suffice to impart ‘absolute concentration robustness’ (ACR) on a wide class of biologically relevant, deterministically modelled mass-action systems. We show here that fundamentally different conclusions about the long-term behaviour of such systems are reached if the systems are instead modelled with stochastic dynamics and a discrete state space. Specifically, we characterize a large class of models that exhibit convergence to a positive robust equilibrium in the deterministic setting, whereas trajectories of the corresponding stochastic models are necessarily absorbed by a set of states that reside on the boundary of the state space, i.e. the system undergoes an extinction event. If the time to extinction is large relative to the relevant timescales of the system, the process will appear to settle down to a stationary distribution long before the inevitable extinction will occur. This quasi-stationary distribution is considered for two systems taken from the literature, and results consistent with ACR are recovered by showing that the quasi-stationary distribution of the robust species approaches a Poisson distribution.

## 1. Introduction

The interaction networks of chemical reaction systems of cellular processes are notoriously complex. Despite this, hidden within the complexity there are often underlying structures that, if properly quantified, give great insight into the dynamical or stationary behaviour of the system. In this vein, Shinar & Feinberg [1] have presented conditions on the structure of biochemical reaction networks that are sufficient to guarantee *absolute concentration robustness* (ACR) on a particular species of the network. When the dynamics of the system are modelled using ordinary differential equations with mass-action kinetics, a species is said to possess ACR if its concentration has the same value at every positive equilibrium concentration permitted by the system of equations, regardless of total molar concentrations. Such a property, which allows cells to respond in a uniform, predictable way given varying environments, is fundamental to many biological processes, including signal transduction cascades and gene regulatory networks [2].

In the two-component EnvZ/OmpR osmoregulatory signalling system, for example, it is important that the amount of phosphorylated OmpR, OmpR-P, which regulates the transcription of the porins OmpF and OmpC, is kept within tight bounds. It has been observed that while OmpR-P is sensitive to the availability of ADP and ATP, it is relatively insensitive to changes in the overall concentrations of the signalling proteins EnvZ and OmpR [3]. Concentration robustness has also been experimentally observed and studied in the two-component KdpD/KdpE [4], PhoQ/PhoP [5] and CpxA/CpxR [6] signalling systems, and in the IDHKP-IDH glyoxylate bypass regulation system [7,8]. Structural sources of robustness have also been identified in the bacterial chemotaxis pathway [9–11].

Consistent with these empirical results, it has been shown that deterministic mathematical models (i.e. ODE models) of the EnvZ/OmpR system yield equilibrium concentrations of OmpR-P that are stable and do not depend on the overall concentrations of either EnvZ or OmpR [12]. This phenomenon of equilibrium concentrations being independent of total molar concentrations was the basis of Shinar & Feinberg [1] where they presented their conditions on the structure of biochemical reaction networks which are sufficient to guarantee ACR on a particular species of the network. Shinar and Feinberg relate the capacity of a network to exhibit ACR to a structural parameter called the *deficiency*, which is well studied in *Chemical Reaction Network Theory* [13–15]. They do not consider stability of such equilibria directly, but they do apply their results to several mass-action models of biochemical networks for which stability is known, including the EnvZ/OmpR signal transduction network and the IDHKP-IDH glyoxylate bypass regulatory system. An interesting recent addition to this framework is the work by Karp *et al.* [16], where the authors carry out a linear analysis of formal expressions in a reaction network to find ACR and more general steady-state invariants.

In this work, we consider stochastically modelled systems satisfying essentially the same network conditions used by Shinar and Feinberg, and we show that strikingly different conclusions are reached pertaining to the long-term dynamics of the systems. For a wide class of biochemical reaction networks with stable ACR equilibria, we show that trajectories are necessarily absorbed by a set of states that reside on the boundary of the positive orthant. Hence, there is necessarily an irreversible ‘extinction’ event in the system. One immediate corollary to this is that the models admit no stationary distributions with mass near the equilibrium of the deterministic model. Our results therefore demonstrate fundamentally different long-term dynamics than those observed in the corresponding deterministic models. Stochastic modelling of chemical reactions is particularly relevant in models of intracellular dynamics because critical proteins may have a low copy number per cell.

Depending upon the total molecular abundances of the constituent species, it may be that such an extinction is a rare event on the relevant timescales of the system. In this case, and under the assumptions used throughout this work, the process will very likely seem to settle down to an equilibrium distribution long before the resulting instability will appear. This distribution is called a *quasi-stationary distribution*, and ACR-like results may still be obtained by consideration of this distribution. In fact, in two examples provided here, we observe that the quasi-stationary distribution of the absolutely robust species limits in a natural way to a Poisson distribution with mean value given by the concentration predicted by the deterministic model.

Section 2 provides a motivating example for our main results. In §3, we formally introduce the concept of a species exhibiting ACR in the deterministic modelling context, and present both the main theoretical result from [1] together with the main result being introduced here, theorem 2. A detailed proof of theorem 2, and more general results, can be found in the electronic supplementary material. In §4, we consider the quasi-stationary distributions for the example models considered here. We close with a brief discussion.

## 2. A motivating example

Consider the two-species activation/deactivation network
2.1where *A* is the active form of a protein, *B* is the inactive form, and *α* and *β* are positive rate constants [1]. Note that the inactive form *B* regulates both the activation and deactivation steps of the mechanism. The usual differential equations governing the time evolution of the molar concentrations of *A* and *B*, denoted *c _{A}* and

*c*here, are 2.2Setting the left-hand sides of the above equations to zero and solving yields the following values for the equilibrium concentrations: 2.3where is the total conserved protein concentration. Equation (2.3) shows that the deterministically modelled system (2.1) has ACR for the protein

_{B}*A*, as all positive equilibria, no matter the initial condition, must satisfy . It can also be easily checked that these equilibria are stable. The self-regulation of the mechanism (2.1), although quite simple, predicts the remarkable property that the concentration of the active protein is kept within tight bounds regardless of the value of the total protein concentration

*M*.

Now consider the usual stochastic model for the network (2.1), which treats the system as a continuous time Markov chain (CTMC) (see figure 1). Let **X*** _{A}*(

*t*) and

**X**

*(*

_{B}*t*) denote the individual counts of

*A*and

*B*, respectively, and model the reactions as discrete events which occur stochastically in time. Under usual assumptions on the rates, or propensities, of the reactions, the first reaction can only occur if

**X**

*(*

_{A}*t*) > 0 and

**X**

*(*

_{B}*t*) > 0, and the second only if

**X**

*(*

_{B}*t*) > 0, which implies no reaction may proceed if

*B*is depleted completely. For a more thorough introduction to the stochastic models for biochemical systems, see the electronic supplemental material or [17].

It is not hard to see that the long-term behaviour of this CTMC is different from that of the deterministic model (2.2). More specifically, it is possible for all of the inactive molecule to become active through the self-activation reaction *B → A*. This sends the chain to the state
2.4where . After this time, neither reaction may occur and so no active molecules *A* may be deactivated. Therefore, rather than having trajectories of the system spend most of their time near the value (2.3), over a sufficiently long time frame the inevitable outcome of the system is convergence to (2.4).

This disparity between the long-term predictions of the deterministic mass-action model and that of the stochastic model at first seems to be at odds with the established result that the stochastic chain is well approximated by the corresponding deterministic model when molecular counts are high. However, such results are valid only on *finite time intervals*, and therefore stand silent on the long-term behaviour of the models [18,19]. Isolated examples of models exhibiting such a fundamental difference between the long-term behaviour of the corresponding deterministic and stochastic models are well known in the biochemical literature [20–23].

## 3. Results pertaining to absolute concentration robustness

The main problem we consider here is: what structural conditions on the reaction network yield an absorption event similar to that of (2.1) for the corresponding stochastic system? We are particularly interested in networks for which the deterministic model predicts ACR, and in this context we will use the original theorem due to Shinar & Feinberg [1]. We briefly introduce some terminology from chemical reaction network theory, including the network parameters *n*, *ℓ* and *s*.

We let *n* denote the number of vertices of the reaction network. These vertices are the linear combinations of the species at either end of a reaction arrow and are called *complexes* in the chemical reaction network literature. We will use this word, although the reference [1] instead uses ‘nodes’ to avoid confusion with the biological meaning of the word ‘complex’. Note that we may naturally associate a complex with a non-negative vector *y* in which the *j*th component of *y* is the multiplicity of species *j* in that complex. For example, the complex *A* + *B* in system (2.1) has associated vector (1, 1), whereas the complex 2*B* has associated vector (0, 2).

We let *ℓ* denote the number of connected components of the reaction network and associate with every reaction a *reaction vector*, which determines the counts of the molecules gained and lost in one instance of that reaction. For example, for the reaction *y* → *y′*, the reaction vector is *y′* − *y*, where we have slightly abused notation by writing the vector associated with a complex in the place of the complex itself. We denote by *s* the dimension of the span of all of the reaction vectors. For the activation/deactivation network (2.1), there are four complexes, {*A* + *B*, 2*B*, *B*, *A*}, and two connected components, {*A* + *B*, 2*B*} and {*B*, *A*}. It follows that *n* = 4 and *ℓ* = 2. We also note that the reaction *A* + *B* → 2*B* has associated reaction vector (−1, 1), since the system loses one *A* molecule and gains a net of one *B* molecule owing to one instance of the reaction. The reaction vector for the reaction *B* → *A* is (1, −1) so that *s* = 1.

The *deficiency* of a network is defined to be . The deficiency is known to only take non-negative values and has been used to show a variety of steady-state results for mass-action systems, both deterministic and stochastic [13,24–30]. For the network (2.1), we have
so the deficiency is one.

We say two complexes are *strongly linked* if there is a directed path of reactions from the first to the second, and also a directed path from the second back to the first. A strong linkage class of a reaction network is a maximal subset of complexes that are strongly linked to each other. A strong linkage class is furthermore called *terminal* if no complex in the class reacts to a complex in another strong linkage class. Complexes that do not belong to a terminal strong linkage class are called *non-terminal*. For example, in the network (2.1), the complexes *A* + *B* and *B* are non-terminal, whereas the complexes 2*B* and *A* are terminal. See also the network in figure 2, where the terminal and non-terminal complexes are labelled in different colours. Finally, we say that two complexes ‘differ only in species *S*’ if the difference between them is a non-zero multiple of a single species *S*. For example, in (2.1), the complexes *A* + *B* and *B* are both non-terminal and also differ only in species *A*.

We can now state the main theorem of [1].

### Theorem 3.1.

*Consider a deterministic mass-action system that admits a positive steady state and suppose that the deficiency of the underlying reaction network is one. If, in the network, there are two non-terminal complexes that differ only in species S, then the system has ACR in S*.

We have already seen the system (2.1) has a positive steady state if *M* > *β*/*α*, that the underlying network has a deficiency of one, and that the non-terminal complexes *A* + *B* and *B* differ only in the species *A*. Thus, theorem 3.1 could be used to guarantee that the system exhibits ACR in *A* even if the equilibrium could not have been calculated explicitly.

Note that theorem 3.1 stands silent on whether or not the equilibria of the deterministically modelled systems are stable, either locally or globally. In fact, there do exist ACR models for which the set of equilibria is unstable. For example,
satisfies the requirements of theorem 3.1, but the equilibrium is *unstable.* The equilibria for the models considered in this article, however, are known to be globally stable.

### 3.1. Stochastic differences in robustness

The network (2.1) demonstrates that biochemical reaction networks satisfying the hypotheses of theorem 3.1, and in which the ACR equilibria are known to be stable, may fail to exhibit similar stability when modelled stochastically. This disparity in the long-term dynamics is not restricted to only a few systems. In this section, we provide a theorem proved in the electronic supplemental material which demonstrates that a large class of biochemical reaction networks satisfying the assumptions of theorem 3.1, and thereby exhibit ACR when modelled deterministically, have absorbing boundary states. In that case, when the deterministic ACR equilibria are stable, the stochastic model exhibits fundamentally different long-term dynamics.

In order to state the result corresponding to theorem 3.1 for stochastically modelled systems, we need a few more basic definitions. A chemical reaction network is said to be *conservative* if there is a vector *w* with strictly positive components for which *w* · (*y*′ − *y*) = 0 for all reactions *y → y*′. Note that in this case, there is necessarily a conserved quantity *M* > 0 so that for the given species set {*S*_{1}, *S*_{2}, *…* , *S _{m}*} we have that
is invariant to the dynamics of the system. For example, we have already seen that is a conserved quantity for the system (2.1).

Enumerating the reactions arbitrarily, we denote by *λ _{k}*(

*x*) the rate, or propensity, function of the reaction

*y*′

_{k}→ y*, and note that it is reasonable to assume that*

_{k}*λ*(

_{k}*x*) > 0 if and only if

*x*≥

_{i}*yk*, which says a reaction may only occur if there are a sufficient number of molecules in the system. We say any family of rate functions satisfying this simple condition is

_{i}*stoichiometrically admissible*. For example, the rate function

*λ*(

_{k}*x*) =

*x*would be stoichiometrically admissible for the reaction

_{B}*B → A*, but not for the reaction

*A*+

*B →*2

*B*. A common choice for the propensities λ

*(*

_{k}*x*) is

*stochastic mass-action*where ,

*κ*

_{k}is the rate constant and

*V*is the volume of the reaction vessel. Note that under the assumption of mass-action kinetics, the propensity function is proportional to the number of ways in which one can choose the molecules necessary for the reaction to occur. Also, note that stochastic mass-action kinetics is stoichiometrically admissible.

We say a complex *y _{k}* is

*turned off*at a particular state value

*x*if

*x*<

_{i}*y*for at least one

_{ki}*i*; otherwise, we say the complex is

*turned on*. Note that

*λ*(

_{k}*x*) = 0 for all

*x*at which

*y*is turned off. Next, we say that a complex

_{k}*y*is

*dominated*by the complex

*y*′ (denoted ) if for all

*i*= 1,

*…*,

*m*. In particular, whenever two complexes differ in only one species, one necessarily dominates the other. The intuition here is that if

*y*′ is ever turned off, then so is

*y*.

Finally, we remind the reader that a state **X** of a CTMC is *recurrent* if the chain satisfying **X**(0) = **X** returns to **X** with probability one, and that a recurrent state is *positive recurrent* if the expected value of the return time is finite. It is also a basic fact that stationary distributions only give mass to *positive recurrent* states, showing that the long-term dynamics are restricted to those states.

The following is the main theoretical result of the paper and should be compared with theorem 3.1. It, along with more general results allowing for *higher deficiency*, is proved in the electronic supplemental material.

### Theorem 3.2.

*Consider a reaction network which is conservative, has a deficiency of one, and for which the deterministically modelled mass-action system admits a positive equilibrium for some choice of rate constants. Suppose that, in the network, there are two non-terminal complexes, y _{1} and y_{2} say, for which*

*. Then, for any choice of stoichiometrically admissible kinetics, all non-terminal complexes of the network are turned off at each positive recurrent state of the stochastically modelled system.*

Hence, a trajectory of the stochastically modelled system will, with a probability of one, be absorbed by a set of states for which all of the non-terminal complexes are turned off. In particular, the propensity of reactions out of the complexes *y*_{1} and *y*_{2} will be zero.

For instance, in the simple system (2.1), we have , and theorem 3.2 states that the complexes *A* + *B* and *B* are turned off at any positive recurrent state. This is easily verified since the only positive recurrent state is given by (2.4). In the electronic supplementary material, we also provide a theorem characterizing the limiting behaviour of the system *after* the absorption event. In particular, we show that the reduced network is weakly reversible and has a deficiency of zero, and so admits a stationary distribution which is a product of Poisson random variables [25]. We also provide in the electronic supplementary material an example in which the conclusions of the theorem do not hold if only the conservation requirement is dropped.

### 3.2. EnvZ/OmpR signalling system

In order to demonstrate theorem 3.2 on a more complicated example, we now consider a model of the two-component EnvZ/OmpR signalling system in *Escherichia coli* [3], which was also studied in [1]. The histidine kinase EnvZ is sensitive to extracellular osmolarity, the input of the system, and in its active phosphorylated form EnvZ-P is able to phosphorylate the response regulator OmpR into the active form OmpR-P. OmpR-P in turn signals the transcription of the porins OmpF and OmpC. The level of OmpR-P can therefore be thought of as the output of the system. EnvZ is also known to play a role in the regulation of the level of OmpR-P through dephosphorylation [31].

Multiple mechanisms have been proposed for the EnvZ/OmpR system [3,12,32]. We will consider the mechanism given in figure 2 which was proposed in [1,12]. Here, it is imagined that ADP (*D*) and ATP (*T*) interact with EnvZ (*X*) to produce bound complexes, but that only ATP can successfully transfer the phosphate group (*P*) to EnvZ to form EnvZ-P (*X _{p}*). EnvZ-P may then transfer the phosphate group to OmpR (

*Y*) to form OmpR-P (

*Y*) while the modified EnvZ–ADP complex regulates the dephosphorylation of OmpR-P. Assuming [

_{p}*D*], [

*T*] and [

*P*] are of sufficient quantity to be relatively unchanged by the course of the reaction, we may incorporate them into the rate constants, yielding the network contained in figure 2.

Note that the non-terminal complexes *XD* and *XD* + *Y _{p}* differ only in the species

*Y*. As the network also has deficiency one (see the electronic supplemental material), by theorem 3.1 we may conclude that the deterministically modelled mass-action system exhibits ACR in

_{p}*Y*. It is shown in the electronic supplemental materials of [1] that the ACR value is 3.1

_{p}The corresponding equilibrium is stable and, consequently, the deterministic model predicts that the active form of the response regulator, OmpR-P, is robust to the overall level of the signalling proteins EnvZ and OmpR. That is to say, the prediction is that the system will exhibit a similar response regardless of differences in these internal characteristics.

We now consider the stochastic model for the network in figure 2. It can be seen that, in addition to satisfying the assumptions of theorem 3.1, the network has the conservation relations 3.2

As the sum of these two relations has support on all species, it follows that the network is conservative. By theorem 3.2, the stochastic model converges in finite time to a state (or set of states) for which all non-terminal complexes are turned off.

Note that, as the species *X*, *XD*, *XT*, *X _{p}Y* and

*XDY*are also non-terminal complexes, theorem 3.2 guarantees that each will be zero after the inevitable absorption event. As

_{p}*X*

_{tot}is conserved, we may also conclude that

*X*=

_{p}*X*

_{tot}at this time. Finally, as

*X*+

_{p}*Y*is also non-terminal, we may also conclude that

*Y*= 0. In this way, the unique sink of the stochastic chain is seen to be 3.3

## 4. Time until absorption and quasi-stationary distributions

A disparity between the long-term behaviour of deterministic and stochastic models of some reaction networks is well known in the literature [20–23]. The activation/deactivation network (2.1) is a canonical example that mimics the well-studied stochastic susceptible–infected–susceptible (SIS) epidemic model [33–36], where *A* and *B* correspond to the number of healthy and infected individuals, respectively. In that setting, the state (2.4) corresponds to a true ‘extinction state’. A similar extinction effect is also achieved in many population biology models where stochastic effects may irreversibly drive a population to zero [37].

Although the chain associated with (2.1) will inevitably converge to the extinction state (2.4), such an event may be exceedingly rare on biologically reasonable timescales. In such situations, the absorbing state is not of practical concern, and, in fact, the process may seem to settle to a stationary distribution. This distribution is called the *quasi-stationary probability distribution*, and it is useful in analysing the transient behaviour of a model before the absorption event. We suspect that stable ACR-like behaviour may still be attained in the present context by consideration of this distribution.

Suppose we denote the absorbing states of a process **X**(*t*) by *∂A*. If we define , where **x** is an arbitrary state in the state space, and , then the transition probabilities conditioned upon non-extinction, *q*** _{x}**(

*t*), are given for

*t*≥ 0 by 4.1

The limiting vector , if it exists, is the *quasi-limiting* distribution of the process, which equals the *quasi-stationary* distribution in the present context. See [38] for proofs of the facts that such a distribution exists, and is unique, in the setting of our theorem 3.2. See [37] for a recent survey article on quasi-stationary distributions in the context of population processes.

Reconsider the activation/deactivation system (2.1). We begin by reimagining the chain as a birth–death process following **X*** _{B}*(

*t*) with the extinction state

**X**

*= 0 corresponding to (2.4). The first reaction corresponds to a ‘birth’ since*

_{B}*B*is increased by one while the second reaction corresponds to a ‘death’ since

*B*is decreased by one. The chain for

**X**

*can then be determined by the conservation relation . Under mass-action kinetics, the corresponding birth and death propensities,*

_{A}*λ*(

*i*) and

*μ*(

*i*), respectively, are given by 4.2where

*i*= 0,

*…*,

*M*corresponds to the state with

*i*molecules of

*B*. Note that

*λ*(0) =

*λ*(

*M*) =

*μ*(0) = 0 so that

**X**

*= 0 is an absorbing state and*

_{B}**X**

*=*

_{B}*M*is a reflecting state. Up to rescaling of the rate constants, this chain is identical to the stochastic SIS epidemic model considered in [33–35].

We are interested in whether the quasi-stationary distribution displays a form of robustness in **X*** _{A}* with respect to changes in overall molecularity

*M*. Furthermore, we are interested in the case of large

*M*, as that is when the time to extinction is large. In the electronic supplemental material, we prove that as

*M → ∞*, the quasi-stationary distribution for

**X**

*approaches the Poisson distribution 4.3which has previously been shown to be the limit of one commonly used approximation to the quasi-stationary distribution [39]. In figure 3, we graphically demonstrate the convergence by providing realizations of the quasi-stationary distribution for different values of*

_{A}*β*and

*M*.

Now consider the EnvZ/OmpR signalling system from figure 2. As in the simple activation/deactivation network, the expected time before entering the state (3.3) may be very large. We therefore consider the quasi-stationary distribution of the process. This distribution cannot be computed using a simple analytic formula, so we approximate it numerically and point the reader to the electronic supplementary material for the details of our computational methods. We use the parameter values given in the table in figure 2*b*, which have been chosen to be in close agreement with the values in [32]. It was found in [41] that a typical *E. coli* cell has roughly 100 total molecules of EnvZ and 3500 molecules of OmpR. Therefore, for our simulations we chose a ratio of *Y*_{tot} : *X*_{tot} of 35 : 1. Based on the deterministic model, the anticipated mean of the ACR species *Y _{p}* in the quasi-stationary distribution is 25.

The results of the simulations are shown in figure 4. In figure 4*a* we see that, for low molecularity of EnvZ and OmpR (1 and 35, respectively), the chains converge to the boundary state (3.3), even while the quasi-stationary distribution becomes apparent. In figure 4*b*, we see that the quasi-stationary distribution of *Y _{p}* appears to converge to a Poisson distribution centred around the deterministic steady state as the total molecularity grows.

## 5. Discussion

Robustness and stability in the face of varying environments are of fundamental importance to the proper functioning of many biological processes, and understanding this behaviour is one arena where mathematics can play a role in elucidating biological phenomena. In this paper, together with its electronic supplemental material, we have outlined a class of structural conditions that are sufficient to guarantee convergence of trajectories of stochastically modelled systems to an absorbing boundary set. Notably, these conditions overlap significantly with a set of structural conditions that are known to confer ACR on the corresponding deterministic models. For such ACR models with stable equilibria, our results present a significant disparity in the predictions for the limiting behaviour of the systems.

This work highlights several points. First, it is surprising that the long-term behaviour of such a large class of systems considered in [1] is fundamentally different for the stochastic and deterministic models. Second, deterministic models are typically unable to capture trapping phenomena such as those described here unless they are artificially modified, for instance by adding degradation terms. We have shown here that for a wide range of models no such ad hoc modifications are necessary for stochastically modelled systems. Third, in regions where the time to absorption for the stochastic model is large relative to the timescale of the system, the proper object of study when considering the stochastic analogue of ACR behaviour is the quasi-stationary distribution, as opposed to the stationary distribution.

This work suggests a number of promising avenues for future research. First, finding more general conditions for which the conclusions of theorem 3.2 hold will be a focus. In particular, weakening the requirement that the system possesses a conservation relation will allow the results to be applicable to more models arising in ecology and population processes, which typically do not satisfy such an assumption, though do often possess low numbers of the constituent species [42]. Second, the question of when equilibria exhibiting ACR in the deterministic modelling context are stable or unstable is, to the best of the authors' knowledge, currently open and has to date received surprisingly little consideration in the literature. Third, we have observed in the two examples considered here the recurrence of the Poisson distribution when analysing a certain limiting behaviour of their quasi-stationary distributions. It is a suspicion of the authors that this phenomenon applies more generally to other systems satisfying ACR, and this will be investigated.

Results of the type presented here are not only of theoretical interest. In particular, they are exactly the types of results required in order to automate the multi-scale reduction methods for stochastic models of biochemical processes currently being explored in the probability literature [43]. For example, in order to determine the behaviour of species operating on a timescale that is slower than other species, it is necessary to first understand the long-term dynamics of the species on the *fast timescale* in order to perform either the necessary stochastic averaging, or to recognize that some species will have gone extinct.

This work is part of a growing research field in which mathematical methods are developed in order to rise above the bewildering complexity of biochemical processes. Algebraic methods, for example, have been successfully employed in a number of areas related to the equilibria of mass-action systems, where the steady states of such systems form a real algebraic variety [26,44,45]. We believe there is significant progress to be made towards the understanding of stochastically modelled systems through analysis of the underlying network, and we hope to uncover the important substructures hidden in the complexity of biochemical networks that inform system behaviour, both on short and long time frames.

## Funding statement

D.F.A. was supported by NSF grant nos. DMS-1009275 and DMS-1318832. M.D.J. was supported by NSF grant no. DMS-1009275 and NIH grant no. R01-GM086881. G.A.E. was supported by NSF grant nos. DMS-1122478 and DMS-1129008.

## Acknowledgements

We gratefully acknowledge the American Institute of Mathematics (AIM) for hosting a workshop at which this research was initiated.

- Received October 14, 2013.
- Accepted January 17, 2014.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.