|Home | About | Journals | Submit | Contact Us | Français|
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.
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  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 , and possibly other systems such as mammalian olfactory receptor development .
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.  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  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
which is described using the system of differential equations
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 1b). 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  appear to show that its stochastic behaviour is analogous to that of the deterministic network. In figure 1c, 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.
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 2a). 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 . 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 , 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  considered the following model of EnvZ–OmpR dynamics (figure 2c). 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, Yp 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 Xtot, Ytot approaches a Poisson distribution
Here is the deterministic steady state of Yp, that is, the distribution of Yp in the stochastic system hovers around the deterministic steady-state . In particular, for reasonably large t and large total protein concentrations.
This is consistent with Gillespie simulations as shown in figure 2b. 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 Xp and Yp, 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.
However, simulations also indicate that it would likely take a long time to reach this absorbing boundary. Figure 2d 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.
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 :
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 1d 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 k2/k1.
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
Now, as this graph formally converges towards a reduced system
The reader might have encountered this graph before—it is the Markov process associated with a simple birth–death model
Moreover, 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 ≥ t0 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 k2/k1 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 qi, in this case, one can carry out a rescaling of time so that the values ki are recovered and the same result holds. But in the more general analysis below, the rate parameters ki will be defined as functions of 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].
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 2c). 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 . 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 Xp + Y is non-terminal, since such a random walk would eventually reach X + Yp and not come back.
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 . 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  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 + YP and XT differ by YP, implying that YP is absolutely robust.
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.  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, XpY, XTYp must have concentration zero at steady state, which implies from the mass concentration equation for Xtot. Since the complex Xp + Y is also off, it follows that Y = 0 because Xp > 0. This in turn means that Ytot = Yp at steady state, from the mass conservation equation for Ytot. One can conclude that there is a single absorbing state for this system, namely that in which Xtot = Xp, Ytot = Yp 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.
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,
which has r mass conservation relations
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 d1 = 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 .
After formally setting , the species becomes constant and disappears as a variable in the system, leading to the reduced network
To obtain the desired convergence result, fix qj and define instead . In practice, one can choose the qj and N so that the rate parameters kj have desired values, for instance, to fit experimental data.
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:
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 , 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 t0, 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 .
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
Here yi(t) is the concentration of variable Ai at time 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
Now, writing this network using the fixed parameters qj one obtains
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.
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
Here is the jth 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 qj instead of kj, this can be written as
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
The 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 , Khammash and co-workers . 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 kth-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
As , 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.
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.  ensures that the system converges globally towards a product-form Poisson stationary distribution
where is the unique steady state of the corresponding deterministic system (4.3). At stationary distribution each of the species Ai 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  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].
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 4a. 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 A1. This means that by the Shinar and Feinberg theorem, A1 is absolutely robust, and by the Anderson, Enciso and Johnston theorem, eventually an extinction event takes place in which are equal to zero and A1 = N.
By formally freezing the variable , we obtain the network from figure 4b, 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 4a
where is the steady state for A1 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 2c). 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 + Yp differ by Yp, hence Yp is absolutely robust in the deterministic system. As was discussed previously, the long-term behaviour of any stochastic simulation is necessarily Xp = Ytot, Yp = Ytot.
By trial and error, one can find that freezing the variables Y and XT leads to the reduced network in figure 4c. 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 4d, 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 . 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 4e. 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 EIp 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 EIp leads to a network with deficiency (figure 4f). Notice that EIp 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 I is approximately Poisson as . Since and , the growth of N is equivalent to that of Etot, Itot:
The mean λI approximates the deterministic steady state of I in the core absolutely robust network as the total protein concentrations grow.
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 qj, in such a way that when N = N0, it holds that and all parameters reach their experimentally determined values. The parameters kj have already been defined in (4.4) as , where
So let us define
This way when N = N0, it holds . Since the main result shows that for fixed t > 0 the distribution of Aj is approximately Poisson for large N, it is hoped that N0 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 N0. This is surprising because the qj depend on N0, 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 N0 = 1. Note that
That is, kj would have the same value if N0 was replaced with 1 and N was replaced with N/N0. 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 ,
which does not depend on N0.
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 . 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 2c), the three proteins X, Xp, 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, Xp, 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 .
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 . 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, EIp 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.
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.
We address in detail equation (6.3), the convergence of the solution xi(t) of the master equation (6.2) towards the solution yi(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  for full proofs.
Define SN 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 SN corresponds to the states , excluding the absorbing state (N, 0). Let be the FSP of the reduced system's master equation (6.1) to SN. Assume that the system and the reduced system (6.1) both have the same initial condition b with finite support.
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 SN. 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 . Let be such that for all states J. Then a monotonicity property for (6.3) implies that for all s > 0 . Putting all together we have
for 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 = (bI) with finite support is used for both systems.
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 N1 sufficiently large that
Now we construct another FSP, that of the original master equation (6.2) with N = N2 > N1, restricted to the set , and call the solution of that projection wI(t). The intuition is that as , it holds , and at the same time one can compare wI(t) with using FSP techniques.
In fact, both the variables wI(t) and zI(t) are defined on the finite set , and as their infinitesimal generators converge. Suppose that N2 is sufficiently large that
On the other hand,
Therefore for every and every sufficiently large N2,
where we are using the inequality for all given by the FSP framework.
The author has no competing interests.
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.