## Abstract

Absolute robustness allows biochemical networks to sustain a consistent steady-state output in the face of protein concentration variability from cell to cell. This property is structural and can be determined from the topology of the network alone regardless of rate parameters. An important question regarding these systems is the effect of discrete biochemical noise in the dynamical behaviour. In this paper, a variable freezing technique is developed to show that under mild hypotheses the corresponding stochastic system has a transiently robust behaviour. Specifically, after finite time the distribution of the output approximates a Poisson distribution, centred around the deterministic mean. The approximation becomes increasingly accurate, and it holds for increasingly long finite times, as the total protein concentrations grow to infinity. In particular, the stochastic system retains a transient, absolutely robust behaviour corresponding to the deterministic case. This result contrasts with the long-term dynamics of the stochastic system, which eventually must undergo an extinction event that eliminates robustness and is completely different from the deterministic dynamics. The transiently robust behaviour may be sufficient to carry out many forms of robust signal transduction and cellular decision-making in cellular organisms.

## 1. Introduction

Chemical reaction networks lie at the heart of many biological processes at the cell level such as signal transduction, development, gene expression and metabolism. These processes are often under strict regulation to prevent an unintended outcome given large amounts of noise. A typical protein may vary in concentration by tens of per cents from cell to cell even in seemingly identical populations [1–3]. An important scientific problem is to discover the mechanisms used by cells to produce a consistent phenotype in the face of this variability [4,5].

Work by Shinar & Feinberg [6] showed that under certain *structural* hypotheses on the underlying chemical network, the steady-state concentration of a particular species is independent of the total concentration of other proteins in the system. This property is called *absolute robustness* (or absolute concentration robustness), and it can be determined from the reaction network independently of rate constants. Experimental applications include some of the most common signalling pathways in bacteria, such as the Envz–Ompr system in *Escherichia coli* [7], and possibly other systems such as mammalian olfactory receptor development [8].

Protein copy numbers are often low inside each cell, i.e. it is not uncommon for key proteins in the system to have less than 10 copies per cell. The stochastic nature of the chemical events themselves becomes relevant in such systems, and it has been widely studied in a variety of contexts (e.g. [9–12]). A natural question is how the noise affects absolutely robust systems, and determining the stationary distribution of the absolutely robust variable. Anderson *et al*. [13] found that for many of these systems, when chemical noise is taken into account an extinction event takes place in the long term, destroying the absolute robustness property.

In this work, I consider the *transient behaviour* of absolutely robust systems under stochastic chemical noise. Preliminary simulations [13] have already hinted that the distribution of the absolutely robust variable after finite time resembles a Poisson distribution centred around the deterministic mean. Here this result will be mathematically demonstrated for a large class of systems, effectively establishing a stochastic counterpart to the work by Shinar & Feinberg. Most cellular signalling and decision-making processes take place under a limited timeframe, so that a transient distribution will often form the basis for these decisions when the final distribution takes a long time to be reached.

The concept of absolute robustness, as well as the relevant stochastic effects, can be best explained through a simple example. Consider the chemical reaction 1.1which is described using the system of differential equations 1.2

These reactions satisfy the mass conservation equation *A* + *B* = *N*, that is, the total protein amount *N* is conserved. In this way, *A* and *B* can be thought of as different forms of the same protein. At steady state, it holds that , and assuming , it follows that . Thus at steady state, *A* is independent of the total concentration *N* of the protein, a surprising result that does not hold for most reaction networks (figure 1*b*). This property is called the *absolute robustness* of *A* in this system. At steady state, the protein acts as a buffer and absorbs the excess protein. (We assume here that , otherwise *B* = 0 and the argument above breaks down.)

Gillespie simulations of this system [14] appear to show that its stochastic behaviour is analogous to that of the deterministic network. In figure 1*c*, this system was modelled for *t* = 10 and the distribution of *A* was plotted for two different total protein concentrations, choosing rate parameters so that . In both cases, the distribution approximates a Poisson distribution with *λ* = 10.

### 1.1. Two-component signalling

The behaviour observed above would be noteworthy but not particularly important if it was only true for that toy model. However, it also holds in larger generality for a family of networks, including a so-called *two-component system* of great importance in bacterial signalling [15–20]. Two-component systems typically consist of a sensory protein *X* that interacts with a regulatory protein *Y*. For instance, the EnvZ–OmpR system in *E. coli* regulates osmolarity, which can cause a damaging influx of water across the cell membrane if not properly controlled.

The EnvZ–OmpR osmoregulatory system has been well characterized and can be described as follows (figure 2*a*). EnvZ has a cytoplasmic osmosensory domain, and it phosphorylates OmpR through a histidine kinase domain. Phosphorylation activates OmpR, which forms a dimer and regulates the transcription of outer membrane porin proteins OmpC and OmpF. These two proteins form pores of two different sizes that regulate the passage of nutrients and other molecules through the cell membrane. Under low osmolarity conditions, the concentration of active OmpR is low, which leads to the expression of OmpF through OmpR binding to high affinity OmpF promoter sites. When osmolarity is high, active OmpR is in high concentration, leading to OmpC expression through OmpR binding to low affinity OmpC promoter sites, and also resulting in OmpF suppression through binding of OmpR to low affinity, OmpF inhibitory sites. In this way, low and high osmolarities lead to OmpF and OmpC expression, respectively.

Importantly, EnvZ not only phosphorylates OmpR, but it can also dephosphorylate it through a separate phosphatase domain [21]. Intuitively, one could think that if total EnvZ concentration were increased, then phosphorylated OmpR concentration might not dramatically change. This is because more EnvZ would more rapidly phosphorylate, but also more rapidly dephosphorylate, OmpR. This experiment was carried out by Batchelor & Goulian [7], who indeed observed a roughly ±30% change in OmpR-p, over two orders of magnitude change in total EnvZ concentration and under different osmolarity conditions. Surprisingly, OmpR-p concentration was also fairly constant over two orders of magnitude change of *total OmpR itself*. Thus, as total OmpR increases, the active OmpR concentration remains roughly constant while the rest of the protein is in inactive form.

Shinar & Feinberg [6] considered the following model of EnvZ–OmpR dynamics (figure 2*c*). EnvZ and OmpR are denoted by *X* and *Y*, respectively. EnvZ can be dephosphorylated, phosphorylated or simply bound to adenosine triphosphate (ATP) (, respectively). When phosphorylated, it has the ability to transfer this phosphorylation to OmpR through a classic Michaelis–Menten reaction. When bound to ATP, it acts as a phosphatase to dephosphorylate OmpR through a similar enzymatic mechanism. Standard mass action rates are used to define the ODE system for this model.

Shinar & Feinberg then demonstrated that OmpR-p is absolutely robust in this system, that is, *Y _{p}* has the same concentration at every positive steady state. In that sense, the output of this system is the same regardless of the total amount of the proteins EnvZ and OmpR.

In the current work, it is shown that the OmprR-p concentration has a transiently robust behaviour in the stochastic case. Specifically, the distribution after finite time *t* for large total protein concentrations *X*_{tot}, *Y*_{tot} approaches a Poisson distribution
Here is the deterministic steady state of *Y _{p}*, that is, the distribution of

*Y*in the stochastic system hovers around the deterministic steady-state . In particular, for reasonably large

_{p}*t*and large total protein concentrations.

This is consistent with Gillespie simulations as shown in figure 2*b*. The distribution of OmpR-p was calculated after a fixed time *t* = 100 and compared with a Poisson distribution centred around the deterministic mean for different total EnvZ and OmpR concentrations. A stationary distribution is found independent of the total protein concentrations.

The long-term picture for any single Gillespie simulation of this system is more complex, since an absorbing boundary event eventually must take place, which destroys the absolute robustness of the system. In this case, an absorbing state is ultimately reached when all EnvZ and OmpR protein is in the phosphorylated state *X _{p}* and

*Y*, respectively. In that case, none of the reactions in the model can take place, i.e. this is an absorbing steady state that is not robust to total protein concentrations.

_{p}However, simulations also indicate that it would likely take a long time to reach this absorbing boundary. Figure 2*d* shows this through a number of Gillespie simulations with different total protein concentrations. In each case, the system was simulated until it reached extinction, and this time was marked with a blue dot. The fact that the times lie on a rough straight line on a semilog plot is evidence that the average time it takes to reach this absorbing state is a roughly exponential function of total protein concentrations.

In the next section, the main ideas for the proof are outlined in the simple case of the toy model (1.1); the associated Markov process is compared with that of a reduced chemical reaction network that has a simpler dynamics. Section 3 introduces background notation and some existing results for absolute robustness. The main result is stated in detail in §4, and proved in §§5–7 as well as appendix A. Section 8 discusses several examples in detail, including an analysis of the EnvZ–OmpR model, followed by a discussion.

## 2. Toy model analysis

The ideas used in the proof of the main result can once again be best introduced using the toy model (1.1). The stochastic system consists of the following continuous-time Markov process [22]: 2.1

Each node corresponds to the probability that *A* = *i* and at a given time. The node (*N*, 0) is an absorbing state, because the probability to leave it is . This means that at steady distribution, with probability 1 it holds that *A* = *N*, *B* = 0. Figure 1*d* estimates the mean time before extinction to be roughly exponential as a function of *N*. However before this state is reached, Gillespie simulations reach a transient distribution similar to a Poisson distribution with mean *k*_{2}/*k*_{1}.

To illustrate why this may be the case, consider the rescaled rate parameters , . In terms of these new parameters, the Markov process has the associated graph
2.2Now, as this graph formally converges towards a *reduced* system
2.3The reader might have encountered this graph before—it is the Markov process associated with a simple birth–death model
2.4Moreover, a standard calculation shows that this reduced system at steady state has a Poisson distribution
with mean . Now, the ODE systems associated with (2.2) and (2.3) cannot be expected to have the same steady-state behaviour, and in fact they behave very differently at steady state. However, as their associated vector fields are similar for large *N*, one can reasonably expect that after *finite* time *t* they have similar solutions. Define at time *t* for the process (2.2), and similarly at time *t* for (2.3). The convergence of (2.2) towards the reduced system (2.3) after fixed finite time *t* will be rigorously proved in appendix A, in the sense that for all *i*. Since it also holds that , bringing both equations together
Because the solutions of (2.3) converge globally towards a Poisson distribution, after fixed finite time *t* ≥ *t*_{0} the distribution of the model (2.2) for large *N* approximates a Poisson distribution centred around *λ*.

The question remains whether one can systematically associate with other absolutely robust systems a reduced network such as the birth–death process above. In order to do this, recall that the original toy model satisfies the mass conservation equation *A* + *B* = *N*. As *A* approaches *k*_{2}/*k*_{1} at steady state, then for large values of *N*. This motivates the step of formally setting *B* = *N*, i.e. *freezing* the value of *B*, which leads to its removal from the system since it is now constant. The new chemical reaction system after removal of *B* and its incorporation into the reaction rates is none other than (2.4).

It will be shown below that this variable freezing procedure can be done in large generality for other networks. In order to freeze a variable, it only needs be part of a mass conservation relation. Multiple variables can be frozen in this way, if there are multiple mass conservation relations available. Note that after freezing a variable, its associated mass conservation relation does not apply anymore.

Regarding the rate parameter scaling and the rate parameters *q _{i}*, in this case, one can carry out a rescaling of time so that the values

*k*are recovered and the same result holds. But in the more general analysis below, the rate parameters

_{i}*k*will be defined as functions of

_{i}*N*, , in such a way that as the network converges towards a fixed reduced system. A similar treatment of parameter values is sometimes carried out, for example, in stochastic ecological models and even the classic SIR model [23,24].

## 3. Previous absolute robustness results

To describe the more general results on absolute robustness, including the main result in this paper, we need a few definitions; consider again the Envz–OmpR osmolarity system (figure 2*c*). The graph of this network has nodes, known in this literature as *complexes*. The graph also has connected components, and the stoichiometry matrix has rank *R* = 5. The stoichiometry matrix contains the net number of molecules of each species that are produced or destroyed in each of the reactions [25]. The so-called *deficiency* of the system is defined as . A complex *Q* in the network is *terminal* if every random walk along the directed graph edges beginning at *Q* eventually returns to *Q*. Thus *X _{p}* + Y is

*non-terminal*, since such a random walk would eventually reach

*X*+

*Y*and not come back.

_{p}A steady state of the chemical reaction is *positive* if all its components have positive values. A variable *A* in the system is *absolutely robust* if all positive steady states have the same value. Assume mass action kinetics for the biochemical networks throughout for simplicity [26]. Also, assume for convenience that the reaction network admits a positive steady state—otherwise there would be nothing to prove about the system.

The main theorem by Shinar & Feinberg [6] for deterministic systems can be stated as follows: given a chemical reaction network such that and two different non-terminal complexes differ by exactly one species *A*, it must follow that *A* is absolutely robust. In the toy model, it holds , , *R* = 1, , and the non-terminal complexes *A* + *B* and *B* differ by *A*, therefore *A* must be absolutely robust. In the EnvZ–OmpR model, and the non-terminal complexes *XT* + *Y _{P}* and

*XT*differ by

*Y*, implying that

_{P}*Y*is absolutely robust.

_{P}This result is particularly striking in that dynamical information is concluded from structural network information alone, i.e. without knowledge of the rate constants. Of course, the rate constants must be known in order to calculate the absolutely robust value of *A* concentration, but the fact that it is robust does not require this knowledge.

A chemical reaction network is *conservative* if every species is included in some mass conservation relation. For example, the EnvZ–OmpR system is conservative because the mass conservation equations and include every single species in the system. Given a state of the system with concentrations for each species, a complex in the chemical reaction network is said to be *off* if at least one of the species in the complex has concentration equal to zero.

The main theorem proved by Anderson *et al*. [13] addresses the long-term behaviour of stochastic, absolutely robust networks: given a conservative chemical reaction network such that and two different non-terminal complexes differ by exactly one species *A*, an extinction event must eventually take place, such that every non-terminal complex turns off and stays off. This can be used to deduce much information about the stationary distribution. For example in the EnvZ–OmpR network, all the complexes labelled in blue must be off at stationary distribution according to the theorem. In particular, the species *X*, *XT*, *X _{p}*Y,

*XTY*must have concentration zero at steady state, which implies from the mass concentration equation for

_{p}*X*

_{tot}. Since the complex

*X*+ Y is also off, it follows that

_{p}*Y*= 0 because

*X*> 0. This in turn means that

_{p}*Y*

_{tot}=

*Y*at steady state, from the mass conservation equation for

_{p}*Y*

_{tot}. One can conclude that there is a single absorbing state for this system, namely that in which

*X*

_{tot}=

*X*,

_{p}*Y*

_{tot}=

*Y*and all other species have concentration zero. The question however remains as to what is the behaviour of such a system before this eventual absorption takes place.

_{p}## 4. Main result statement

Having described the relevant background, the following is the main theorem of this work. Consider a chemical reaction network with *n* + *r* species and *m* reactions,
4.1which has *r* mass conservation relations
4.2

The variables are eliminated through a ‘variable freezing’ procedure as illustrated in the Introduction. Before doing this, note that only one variable appears in each mass conservation equation in (4.2) and the coefficient in front of is equal to one. This could be done by carrying out linear operations on existing conservation relations as will be shown in the applications. No restrictions are made on the sign of ; however , . In the case *r* = 1, one can set *d*_{1} = 1 simply by rescaling *N*. We also ask that *at least one of the variables* *is turned off at steady state*. This can be easily verified by looking at the reaction network and using the main result in [13].

After formally setting , the species becomes constant and disappears as a variable in the system, leading to the reduced network 4.3where 4.4

To obtain the desired convergence result, fix *q _{j}* and define instead . In practice, one can choose the

*q*and

_{j}*N*so that the rate parameters

*k*have desired values, for instance, to fit experimental data.

_{j}### Theorem 4.1.

*Suppose that a conservative chemical reaction network* (*4.1*) *satisfies* *and has two non-terminal nodes that differ by a single variable* *. Also suppose that after freezing* *, the irreducible network (4.3) has deficiency* *. Then the distribution of* *in the stochastic system, for finite t and sufficiently large N, is approximated by a Poisson distribution:*
4.5

*Moreover, the mean* *of the distribution is the deterministic steady state of* *in* (*4.1*) *in the limit as* .

This result is proved by showing that the stochastic system associated with (4.1) has very similar dynamics after fixed finite time to that of (4.3). Since (4.3) has a Poisson stationary distribution for each of its variables, due to a result from chemical reaction network theory [27], then this is also the case for (4.1). The ‘irreducible’ assumption for (4.3) means that the Markov process associated with it has a strongly connected graph. This assumption is necessary to conclude that the solutions of (4.3) converge globally towards the unique stationary distribution.

The overall picture of the dynamical behaviour is described in figure 3. If *N* is fixed and the stochastic system associated with (4.1) evolves over time, the solution will eventually undergo an extinction event. However for fixed values of *t*, larger values of *N* eventually lead to a Poisson-like distribution for . The values of *t* for which this holds need only be larger than or equal to *t*_{0}, the time that it takes for the stochastic reduced network (4.3) to reasonably converge towards a Poisson distribution; that convergence is independent of *N*.

Regarding the value , it is the steady state of in the reduced network (4.3). It is also the limit as of , the steady state of for the original system (4.1). As such, it retains the absolute robustness properties of the original network.

Note the similarity with the assumptions of the existing results for deterministic and stochastic absolutely robust networks. The main assumptions made beyond those in the original result by Shinar and Feinberg are that the network is conservative and that after the variable freezing procedure the deficiency is zero, all of which can be verified easily without knowledge of parameter values. The irreducibility of the reduced network can also be readily verified.

The idea of freezing the value of a variable has been suggested previously in a deterministic context, notably in work by Joshi & Shiu [28].

## 5. Deterministic network convergence

This section is concerned with proving the convergence of the deterministic dynamics in the network (4.1) towards that of (4.3) after finite time for . It also constitutes a warm up for a similar analysis in the stochastic case. The reduced network (4.3) has the associated system of equations
5.1Here *y _{i}*(

*t*) is the concentration of variable

*A*at time

_{i}*t*, , , and we use the shorthand notation , .

Also, the original system of reactions (4.1) can be reduced to *n* species in a standard way by simply writing in terms of the other variables,
The resulting network is described by the system of equations
5.2Now, writing this network using the fixed parameters *q _{j}* one obtains
5.3

This allows one to compare the behaviour of the reduced ODE system (5.1) with that of the original ODE system (5.3) for . As *N* increases, it is clear that the equation on the right-hand side of (5.3) converges towards that of (5.1). This convergence holds pointwise and is also uniform on compact sets. One can conclude that
for every *t* > 0 and *i* = 1, … ,*n*, provided both systems start with the same initial condition and both systems are defined for that time interval.

## 6. Stochastic network convergence

A similar argument can be carried out to compare the master equation of both systems in the stochastic case. The state of the network at any time is given by a vector of non-negative integers , representing the number of molecules for each species. Denoting at time *t* for the reduced network (4.3), these probabilities evolve over time according to the ODE system
6.1

Here is the *j*th reaction vector, for any non-negative integers *i*, *m* and .

Moving on to the original network (4.1), one can once again eliminate all the variables by writing them in terms of the other variables, . In this case,
Using the notation for *q _{j}* instead of

*k*, this can be written as 6.2

_{j}One can verify that for any non-negative integers *N*, *i* and *s*, it holds as . Applying this principle above, as in the above equation the right-hand side converges pointwise towards that of (6.1).

An argument as described in the Introduction would imply that if the same initial condition is used for both systems (6.2) and (6.1), and if each is integrated over a finite interval [0, *t*], then for each *I* it holds that
6.3The problem is that one of the two systems of equations, that of the reduced network, has infinitely many state variables *I*, so that this convergence is a non-trivial problem. Equation (6.3) is rigorously proved in appendix A using a technique known as finite-state projection (FSP), originally developed by Munsky [29], Khammash and co-workers [24]. FSP is normally used to carry out numerical approximations of master equation solutions and to establish numerical error bounds. But it can also be helpful as a form of theoretical analysis.

An intuitive explanation of this convergence for finite time *t* is the following. Suppose that the master equation (6.2) is written in compact form as , where is the intensity of the transition from *I* to *J* for every two states *I*, *J*. Suppose that at time *t* = 0 the system is in state *J*, i.e. and for all . One can use the shorthand notation . Then a *k*th-order Taylor approximation can be used to accurately approximate *x*(*t*) over a finite range of values of *t*:
Recall that the matrix entry can be written as
where the sum is taken over all directed paths from *I* to *J* of length . One can (at least formally) write the master equation of the reduced system (6.1) as . For a fixed state *I*, approximate
6.4As , the paths of length up to *k* in both systems coincide, and converges to . Although the matrix *L* is infinite, only finitely many paths are considered in the equation above. As *k* is increased, the range of values of *t* for which the approximation applies increases but remains bounded, stressing the fact that this approximation is only valid for finite *t*. In other words, if *N* is fixed and *t* grows indefinitely, equation (6.4) eventually breaks down.

## 7. Proof of the theorem

One can now finish the proof of the main theorem, including equation (4.5) on the Poisson-like transient distribution of the absolutely robust variable. Recall that the original biochemical network (4.1) has deficiency , admits a positive equilibrium, and two non-terminal linkage classes differ by the species . Moreover, after freezing the variables it reduces to equation (4.3) which has deficiency and is irreducible as a Markov process.

There is a mature theory for stochastic chemical reaction networks with deficiency zero. But for this theory to apply it is important for the network to have a technical property known as weak reversibility. It will be left to the end of this section to prove that under the given hypotheses this is the case. Now, given and weak reversibility, a result by Anderson *et al*. [27] ensures that the system converges globally towards a product-form Poisson stationary distribution
7.1where is the unique steady state of the corresponding deterministic system (4.3). At stationary distribution each of the species *A _{i}* in the stochastic reaction network has marginal Poisson distribution and is uncorrelated from the other variables, a remarkable result for a nonlinear system. The product-form Poisson distribution is the global attractor of the master equation because this system is irreducible.

From equations (6.3) and (7.1), it follows that Then the absolutely robust variable has Poisson distribution with mean in this limit, leading to the main equation (4.5) for the original system (4.1).

At this point, we know that the transient distribution of approaches a Poisson distribution with mean , but we still need to relate with a deterministic steady-state value for (4.1) (recall that is just the steady state for in the *reduced* system (4.3)).

For this, we will use some properties of deterministic systems with and weak reversibility. By the deficiency zero theorem [30,31], is the unique steady state of (4.3), has only positive entries and is locally stable. This means that small changes to the reduced vector field (5.1) lead to small changes in the steady-state concentrations. Given the convergence of (5.3) towards the reduced system (5.1) as , it follows that there are positive steady states of *the original biochemical network* (4.1) satisfying mass conservation equation (4.2), such that
That means that as . That is, for large *N* the steady states of both reaction networks are similar, and the Poisson-like distribution of is approximately centred around the steady state of . This completes the proof of the theorem.

We conclude the section with a proof that the reduced chemical master equation (6.1) is weakly reversible. Since the original deterministic network (4.1) has a positive steady state, it follows [32] that it is *consistent*, that is, there exist , such that
where are the reaction vectors of the system. The reduced network (4.3) has very similar reaction vectors , and obviously it also holds , that is, the reduced network is also consistent. But any consistent network must have a positive equilibrium for some kinetic rates (for instance, letting and for all *i*). Since the reduced network has deficiency zero, it follows from the deficiency zero theorem that this network is weakly reversible [33,34].

## 8. Examples

In this section, we illustrate the main result in the context of several examples. Consider first a generalization of the toy model (1.1) that has *n* + 1 variables rather than two, given by the network in figure 4*a*. This network has the sole mass conservation law
which makes it a conservative network because it involves all species. It also has linkage classes, *n* + 3 complexes, rank *R* = *n*, and . Finally, the non-terminal complexes and differ by *A*_{1}. This means that by the Shinar and Feinberg theorem, *A*_{1} is absolutely robust, and by the Anderson, Enciso and Johnston theorem, eventually an extinction event takes place in which are equal to zero and *A*_{1} = *N*.

By formally freezing the variable , we obtain the network from figure 4*b*, which has deficiency . The parameter scaling is , and for *i* = 2 … *n*.

The main theorem applies to conclude that for the original network in figure 4*a*
where is the steady state for *A*_{1} in the reduced network. One can also verify by inspection that is also the steady state for the original network, for arbitrary *N*.

The second example is the EnvZ–OmpR model that has been used throughout this work. EnvZ, denoted by *X*, is a bifunctional sensory protein, and OmpR, denoted by *Y*, is its associated actuator. Phosphorylated EnvZ activates OmpR by phosphorylation, and EnvZ bound to ATP dephosphorylates OmpR (figure 2*c*). The mass conservation equations in this case are
This network is conservative since all species are involved in at least one of the two equations. By calculating the rank of the stoichiometry matrix, . Moreover, the complexes *XT* and *XT* + *Y _{p}* differ by

*Y*, hence

_{p}*Y*is absolutely robust in the deterministic system. As was discussed previously, the long-term behaviour of any stochastic simulation is necessarily

_{p}*X*=

_{p}*Y*

_{tot},

*Y*

_{p}*=*

*Y*

_{tot}.

By trial and error, one can find that freezing the variables *Y* and *XT* leads to the reduced network in figure 4*c*. At first, there seems to be a problem, since this network still has deficiency and is not weakly reversible. However notice that in this system the variable *X* is decoupled from the rest of the network, i.e. it does not affect any other variable. Without altering the dynamics of the deterministic or stochastic network, one can eliminate this variable and obtain the network in figure 4*d*, which indeed has deficiency zero and is weakly reversible. Finally, note that the mass conservation equations include exactly one frozen variable each, as required in the format (4.2), and that both of the frozen variables eventually turn off during the extinction event. Setting , , one can apply the main theorem and conclude that
since is equivalent to the growth of the total protein concentrations.

The next example is the so-called core absolutely robust module discussed by Shinar & Feinberg [6]. This is a simple model of a two-component signal transduction network, with a bifunctional sensory protein *E* and a regulatory protein *I*. Each of the two proteins may be phosphorylated or dephosphorylated, and they can also form dimer and trimer complexes as shown in figure 4*e*. This system satisfies the two mass conservation equations
and is therefore conservative, since every species is in one of these two equations. Also and the non-terminal complexes and *EI _{p}* differ by

*I*, so that

*I*is absolutely robust.

By trying out different combinations of species, one can see that freezing the species *E* and *EI _{p}* leads to a network with deficiency (figure 4

*f*). Notice that

*EI*forms a non-terminal complex, so it must turn off at stationary distribution for any stochastic simulation. Since the mass conservation equations are not in the correct format (exactly one frozen variable must appear in each equation), one can rearrange as By freezing the variables and , the system is reduced to one with deficiency zero. By the theorem, the transient behaviour of

_{p}*I*is approximately Poisson as . Since and , the growth of

*N*is equivalent to that of

*E*

_{tot},

*I*

_{tot}: The mean

*λ*approximates the deterministic steady state of

_{I}*I*in the core absolutely robust network as the total protein concentrations grow.

## 9. Fitting to experimental data

Suppose that a modeller is interested in describing the behaviour of a specific biochemical system that appears to be absolutely robust. After writing down the chemical reactions (4.1), the modeller estimates the value of the rate parameters and total protein concentrations on the right-hand side of (4.2) based on various forms of experimental data. Let us call the baseline rate parameters , and the total protein concentrations.

The idea is to find values of the *q _{j}*, in such a way that when

*N*=

*N*

_{0}, it holds that and all parameters reach their experimentally determined values. The parameters

*k*have already been defined in (4.4) as , where So let us define This way when

_{j}*N*=

*N*

_{0}, it holds . Since the main result shows that for fixed

*t*> 0 the distribution of

*A*is approximately Poisson for large

_{j}*N*, it is hoped that

*N*

_{0}is sufficiently large that this approximation holds.

As it turns out, the absolute robustness of the original network (4.1) provides a form of absolute robustness to the reduced network (4.3). Specifically, the steady-state concentration of in the reduced network is independent of *N*_{0}. This is surprising because the *q _{j}* depend on

*N*

_{0}, so in principle one would expect that the steady states of that network also depend on this parameter.

To see why this is the case, recall that is the unique steady state of (4.3), and that there exists a sequence of positive steady states of (4.1) satisfying mass conservation (4.2), such that as . Also recall that is an absolutely robust species in (4.1).

Denote and in the special case *N*_{0} = 1. Note that
That is, *k _{j}* would have the same value if

*N*

_{0}was replaced with 1 and

*N*was replaced with

*N*/

*N*

_{0}. Moreover, by absolute robustness of , for large

*N*one could replace the mass conservation equation (4.2) with without changing the steady-state concentration of . That is, . As , 9.1which does not depend on

*N*

_{0}.

## 10. Discussion

This paper has described the stochastic behaviour of systems that are absolutely robust in the deterministic framework, and it has shown that a form of transient robustness is preserved before an eventual population collapse. This dynamics is analogous to the behaviour of stochastic ecological models involving, for example, a predator and prey in an island. While the dynamics can resemble that of the corresponding deterministic system after finite time, in the long run one of the species will die off after an unlikely sequence of events, causing a drastic qualitative change. Because this extinction appears to take place at a very long time scale (e.g. in exponential time as a function of the population), in many applications the relevant outcomes are determined away from equilibrium.

The main result here is framed in a very similar way as the original result from Shinar and Feinberg, involving networks with deficiency one and mass action kinetics. This was done in part for convenience, because it is now known that absolute robustness holds for many systems of higher deficiency [35]. Simulations appear to show that the conservative assumption is not necessary for the result to hold, which could be proved in future work. Similarly, no attempt was made to generalize the above results to non-mass action systems, e.g. involving different exponents in the kinetic reaction rates.

Importantly, the proof of the main result implies more than was stated, namely, the distribution of other variables in the original system is also approximately Poisson in the limit as . This has a number of consequences for the behaviour of these systems as predicted by methods such as fluctuation dissipation or the linear noise approximation, which will be studied in future work. Overall, it appears that variable freezing can be applied in a wide range of low-deficiency systems to prove that their behaviour transiently approximates the behaviour of zero deficiency systems, both in the deterministic and stochastic cases.

As this paper neared completion the author was made aware of an unpublished manuscript by D. Anderson, D. Cappelletti and T. Kurtz on the topic of stochastic absolute robustness. While there is overlap with that work, there are significant differences in the assumptions and particularly in the methodology of both forms of analysis, so that the two studies are broadly complementary.

Recall that the main result was applied to a specific model of two-component signalling involving the osmolarity regulatory proteins EnvZ and OmpR. To what extent could this result apply to other networks? The key for the application of the result is the existence of a bifunctional enzyme that can carry out a modification (such as phosphorylation) and separately undo it through another reaction domain. This feature is very common in two-component systems, in that the sensory protein is usually a bifunctional histidine kinase/phosphatase. The fact that the same protein is carrying out two separate, similar reactions introduces a redundancy into the system, which decreases the rank of the stoichiometry matrix of the network, and thereby increases the deficiency. For example, in the EnvZ–OmpR model (figure 2*c*), the three proteins *X*, *X _{p}*,

*XT*can turn into each other and therefore are equivalent in the network for the purposes of rank calculation. One also has the overall reactions where

*X*

_{*}stands for one of the equivalent proteins

*X*,

*X*,

_{p}*XT*. There is a linear dependence between these two reactions which reduces the rank of the reaction matrix, and therefore increases the deficiency from zero to one. This relationship between enzyme bifunctionality and deficiency likely also applies to other models of two-component signalling, as well as to other models involving a bifunctional enzyme. Notice that if one replaced the separate phosphorylation and dephosphorylation reactions by a single, reversible enzymatic reaction, this argument would break down and the deficiency of the system would not increase. This could be an important reason for the rather mysterious existence of a separate phosphatase domain in EnvZ and other histidine kinases [21].

It is surprising that although the abstract result does not require protein bifunctionality, most of the biologically motivated examples of this theorem in the literature involve bifunctional enzymes (see also [19,36]). The simplest example of absolute robustness, the toy model described in the Introduction, can also be thought of as involving a bifunctional enzyme *B*, which catalyses the degradation of *A* () but also spontaneously promotes it (). Another interesting example of enzyme ‘bifunctionality’ is the histone demethylase LSD1 involved in olfactory receptor development [8]. This protein eliminates methylation sites in histones which are associated with olfactory receptor expression. However, some of these methylation sites promote gene expression, while others inhibit it. In this way, LSD1 can carry out both gene activation and gene inhibition. Absolute robustness may play a strong role in this heavily regulated system.

The specific form for the mass conservation equations (4.2) introduces certain restrictions on the total protein concentrations. For example, in the core absolutely robust model example the chosen frozen variables *E*, *EI _{p}* lead to a the condition . However, simulations give no indication that such restriction is necessary. One intuitive way to generalize the algorithm is to iteratively freeze one variable after another, and to apply the main result in sequence, first with the very last network, which has only one mass conservation equation left, then with the second to last, and so forth. Each step could then be used to guarantee an approximately Poisson distribution for the network in the previous iteration. Such a procedure could be proved in the future.

## Competing interests

The author has no competing interests.

## Funding

The author would like to thank the Isaac Newton Institute for Mathematical Sciences, Cambridge, for support and hospitality during the programme Stochastic Dynamical Systems in Biology, where work on this paper was undertaken. This work was supported by EPSRC grant no. EP/K032208/1. This material is based upon work supported by the National Science Foundation under grant nos. DMS-1122478 and 1616233.

## Acknowledgements

I would like to thank Simon Cotter, Michael Cranston, Oleg Igoshin, Badal Joshi, Ilya Kachkovskiy, Jay Newby and Abhyudai Singh for many helpful suggestions. Credit is due to Badal Joshi in particular for the suggested use of consistent networks in this context.

## Appendix A

We address in detail equation (6.3), the convergence of the solution *x _{i}*(

*t*) of the master equation (6.2) towards the solution

*y*(

_{i}*t*) of (6.1) as .

The technique known as finite state projection, or FSP for short, was developed originally by Mustafa Khammash and Brian Munsky to truncate the master equation of a chemical reaction and establish error bounds compared with the original solution (e.g. [29,37]). In an FSP analysis, a linear system corresponding to a Markov process with countably many variables is projected to a system which is identical to the original system except that all states outside of a finite region of interest *S* are clumped together into a single absorbing state *G*. All the reaction edges and intensities that are internal to *S* are unchanged, but all fluxes originally leading out of *S* are now aimed towards *G*, and all fluxes leading into *S* from outside *S* are eliminated from the system altogether. The basic result that we use regarding such projections is that if both systems are initialized at time *t* = 0 with the same initial condition *b* with support in *S*, then (i) for every *t* > 0 and and (ii) , for all *t* > 0. See lemma 5.0.1 and theorem 5.0.2 of [29] for full proofs.

Define *S _{N}* to be the finite set of non-absorbing states in the master equation (6.2) of the original system. So for instance, in the toy model of the Introduction

*S*corresponds to the states , excluding the absorbing state (

_{N}*N*, 0). Let be the FSP of the reduced system's master equation (6.1) to

*S*. Assume that the system and the reduced system (6.1) both have the same initial condition

_{N}*b*with finite support.

### Proposition A.1.

*For every fixed t > 0 and state* *,*

*Proof.* It will be used here that the reduced network associated with the master equation (6.1) has deficiency zero and is weakly reversible. Denote by the set of states that have an outedge leading out of *S _{N}*. The standard results from FSP theory guarantee that (i) for every and (ii) . Then for every ,

We now bound the different terms in this expression: , , for and . Finally, define , the product Poisson distribution that is a steady state of the reduced system (6.3) by the deficiency assumptions [27]. Let be such that for all states *J*. Then a monotonicity property for (6.3) implies that for all *s* > 0 [38]. Putting all together we have
A 1for all .

On the other hand, |*I*| is large for any , in the sense that there exists a constant such that for any such *I*. This is because there is at least one frozen variable that is off for every , and

Therefore, for .

If *I* is an absorbing state and is any state that has an edge pointing towards *I*, then for some *j*, and . Therefore as , |*J*| grows linearly for all . The right-hand side of (A 1) converges to zero since decreases exponentially as while only grows polynomially.▪

One can now return to equation (6.3), the convergence of the solutions of the original master equation (6.2) towards the solution of (6.1). Suppose that an identical initial condition *b* = (*b _{I}*) with finite support is used for both systems.

### Proposition A.2.

*For every fixed t > 0 and state I,*

*Proof.* Once again it will be important that the reduced network (6.1) has deficiency zero and is weakly reversible. Notice that in the proof of proposition A.1 a stronger result was proved, namely that as . Fix *t* > 0 and , and choose *N*_{1} sufficiently large that

Now we construct another FSP, that of the original master equation (6.2) with *N* = *N*_{2} > *N*_{1}, restricted to the set , and call the solution of that projection *w _{I}*(

*t*). The intuition is that as , it holds , and at the same time one can compare

*w*(

_{I}*t*) with using FSP techniques.

In fact, both the variables *w _{I}*(

*t*) and

*z*(

_{I}*t*) are defined on the finite set , and as their infinitesimal generators converge. Suppose that

*N*

_{2}is sufficiently large that Then

On the other hand,

Therefore for every and every sufficiently large *N*_{2},
where we are using the inequality for all given by the FSP framework.▪

- Received June 15, 2016.
- Accepted August 5, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.