PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of interfaceThe Royal Society PublishingInterfaceAboutBrowse by SubjectAlertsFree Trial
 
J R Soc Interface. 2016 August; 13(121): 20160475.
PMCID: PMC5014071

Transient absolute robustness in stochastic biochemical networks

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.

Keywords: chemical reaction networks, absolute robustness, Poisson distribution, deficiency theory

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 [13]. 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. [912]). 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

equation image
1.1

which is described using the system of differential equations

equation image
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 An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i1.jpg, and assuming An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i2.jpg, it follows that An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i3.jpg. 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 An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i4.jpg acts as a buffer and absorbs the excess protein. (We assume here that An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i5.jpg, otherwise B = 0 and the argument above breaks down.)

Figure 1.
(a) The reactions of the so-called toy model of absolute robustness. (b) Steady states of the toy model for different total protein concentrations N. (c) Poisson distribution (red dotted), and histogram of Gillespie simulations for this model with k1 ...

Gillespie simulations of this system [14] 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 An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i6.jpg. 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 [1520]. 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.

Figure 2.
(a) The EnvZ–OmpR two-component pathway regulates osmolarity in E. coli. (b) Histograms of Yp in stochastic simulations of (c) at time t = 100 and two different total protein concentrations. Reaction An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i7.jpg has rate k3 = 5, all other reactions have ...

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 2c). EnvZ and OmpR are denoted by X and Y, respectively. EnvZ can be dephosphorylated, phosphorylated or simply bound to adenosine triphosphate (ATP) (An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i9.jpg, 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

equation image

Here An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i10.jpg is the deterministic steady state of Yp, that is, the distribution of Yp in the stochastic system hovers around the deterministic steady-state An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i11.jpg. In particular, An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i12.jpg 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.

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]:

equation image
2.1

Each node An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i13.jpg corresponds to the probability that A = i and An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i14.jpg at a given time. The node (N, 0) is an absorbing state, because the probability to leave it is An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i15.jpg. 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 An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i16.jpg, An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i17.jpg. In terms of these new parameters, the Markov process has the associated graph

equation image
2.2

Now, as An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i18.jpg this graph formally converges towards a reduced system

equation image
2.3

The reader might have encountered this graph before—it is the Markov process associated with a simple birth–death model

equation image
2.4

Moreover, a standard calculation shows that this reduced system at steady state has a Poisson distribution

equation image

with mean An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i19.jpg. 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 An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i20.jpg at time t for the process (2.2), and similarly An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i21.jpg 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 An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i22.jpg for all i. Since it also holds that An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i23.jpg, bringing both equations together

equation image

Because the solutions of (2.3) converge globally towards a Poisson distribution, after fixed finite time tt0 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 An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i24.jpg 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, An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i25.jpg, in such a way that as An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i26.jpg 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 2c). The graph of this network has An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i27.jpg nodes, known in this literature as complexes. The graph also has An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i28.jpg 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 An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i29.jpg. 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 [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 An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i30.jpg 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 An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i31.jpg, An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i32.jpg, R = 1, An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i33.jpg, and the non-terminal complexes A + B and B differ by A, therefore A must be absolutely robust. In the EnvZ–OmpR model, An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i34.jpg 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 An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i35.jpg and An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i36.jpg 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 An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i37.jpg 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 An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i38.jpg 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.

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,

equation image
4.1

which has r mass conservation relations

equation image
4.2

The variables An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i39.jpg are eliminated through a ‘variable freezing’ procedure as illustrated in the Introduction. Before doing this, note that only one variable An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i40.jpg appears in each mass conservation equation in (4.2) and the coefficient in front of An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i41.jpg 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 An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i42.jpg; however An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i43.jpg, An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i44.jpg. In the case r = 1, one can set d1 = 1 simply by rescaling N. We also ask that at least one of the variables An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i45.jpg 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 An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i46.jpg, the species An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i47.jpg becomes constant and disappears as a variable in the system, leading to the reduced network

equation image
4.3

where

equation image
4.4

To obtain the desired convergence result, fix qj and define instead An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i48.jpg. In practice, one can choose the qj and N so that the rate parameters kj have desired values, for instance, to fit experimental data.

Theorem 4.1. —

Suppose that a conservative chemical reaction network (4.1) satisfies An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i49.jpg and has two non-terminal nodes that differ by a single variable An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i50.jpg. Also suppose that after freezing An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i51.jpg, the irreducible network (4.3) has deficiency An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i52.jpg. Then the distribution of An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i53.jpg in the stochastic system, for finite t and sufficiently large N, is approximated by a Poisson distribution:

equation image
4.5

Moreover, the mean An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i54.jpg of the distribution is the deterministic steady state of An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i55.jpg in (4.1) in the limit as An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i56.jpg.

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 An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i57.jpg. 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.

Figure 3.
Cartoon depiction of the expected dynamics of the absolutely robust variable in a stochastic simulation, as a function of both time t and the total protein concentration parameter N.

Regarding the value An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i58.jpg, it is the steady state of An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i59.jpg in the reduced network (4.3). It is also the limit as An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i60.jpg of An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i61.jpg, the steady state of An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i62.jpg 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 An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i63.jpg. 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

equation image
5.1

Here yi(t) is the concentration of variable Ai at time t, An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i64.jpg, An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i65.jpg, and we use the shorthand notation An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i66.jpg, An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i67.jpg.

Also, the original system of reactions (4.1) can be reduced to n species in a standard way by simply writing An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i68.jpg in terms of the other variables,

equation image

The resulting network is described by the system of equations

equation image
5.2

Now, writing this network using the fixed parameters qj one obtains

equation image
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 An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i69.jpg. 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

equation image

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 An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i70.jpg, representing the number of molecules for each species. Denoting An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i71.jpg at time t for the reduced network (4.3), these probabilities evolve over time according to the ODE system

equation image
6.1

Here An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i72.jpg is the jth reaction vector, An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i73.jpg for any non-negative integers i, m and An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i74.jpg.

Moving on to the original network (4.1), one can once again eliminate all the variables An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i75.jpg by writing them in terms of the other variables, An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i76.jpg. In this case,

equation image

Using the notation for qj instead of kj, this can be written as

equation image
6.2

One can verify that for any non-negative integers N, i and s, it holds An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i77.jpg as An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i78.jpg. Applying this principle above, as An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i79.jpg 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

equation image
6.3

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 [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 An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i80.jpg, where An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i81.jpg 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. An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i82.jpg and An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i83.jpg for all An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i84.jpg. One can use the shorthand notation An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i85.jpg. Then a kth-order Taylor approximation can be used to accurately approximate x(t) over a finite range of values of t:

equation image

Recall that the matrix entry An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i86.jpg can be written as

equation image

where the sum is taken over all directed paths from I to J of length An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i87.jpg. One can (at least formally) write the master equation of the reduced system (6.1) as An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i88.jpg. For a fixed state I, approximate

equation image
6.4

As An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i89.jpg, the paths of length up to k in both systems coincide, and An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i90.jpg converges to An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i91.jpg. 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 An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i92.jpg, admits a positive equilibrium, and two non-terminal linkage classes differ by the species An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i93.jpg. Moreover, after freezing the variables An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i94.jpg it reduces to equation (4.3) which has deficiency An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i95.jpg 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 An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i96.jpg and weak reversibility, a result by Anderson et al. [27] ensures that the system converges globally towards a product-form Poisson stationary distribution

equation image
7.1

where An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i97.jpg 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

equation image

Then the absolutely robust variable An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i98.jpg has Poisson distribution with mean An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i99.jpg in this limit, leading to the main equation (4.5)

equation image

for the original system (4.1).

At this point, we know that the transient distribution of An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i100.jpg approaches a Poisson distribution with mean An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i101.jpg, but we still need to relate An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i102.jpg with a deterministic steady-state value for (4.1) (recall that An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i103.jpg is just the steady state for An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i104.jpg in the reduced system (4.3)).

For this, we will use some properties of deterministic systems with An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i105.jpg and weak reversibility. By the deficiency zero theorem [30,31], An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i106.jpg 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 An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i107.jpg, it follows that there are positive steady states An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i108.jpg of the original biochemical network (4.1) satisfying mass conservation equation (4.2), such that

equation image

That means that An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i109.jpg as An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i110.jpg. That is, for large N the steady states of both reaction networks are similar, and the Poisson-like distribution of An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i111.jpg is approximately centred around the steady state An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i112.jpg of An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i113.jpg. 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 An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i114.jpg, An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i115.jpg such that

equation image

where An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i116.jpg are the reaction vectors of the system. The reduced network (4.3) has very similar reaction vectors An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i117.jpg, and obviously it also holds An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i118.jpg, that is, the reduced network is also consistent. But any consistent network must have a positive equilibrium for some kinetic rates An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i119.jpg (for instance, letting An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i120.jpg and An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i121.jpg 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 4a. This network has the sole mass conservation law

equation image

which makes it a conservative network because it involves all species. It also has An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i122.jpg linkage classes, n + 3 complexes, rank R = n, and An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i123.jpg. Finally, the non-terminal complexes An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i124.jpg and An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i125.jpg 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 An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i126.jpg are equal to zero and A1 = N.

Figure 4.
(a,b) A generalization of the toy model using n+1 species, and its corresponding reduced network after freezing the variable An+1. (c) Network resulting from the EnvZ–OmpR model after freezing the variables Y and XT. (d) The same network after ...

By formally freezing the variable An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i127.jpg, we obtain the network from figure 4b, which has deficiency An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i128.jpg. The parameter scaling is An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i129.jpg, An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i130.jpg and An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i131.jpg for i = 2 … n.

The main theorem applies to conclude that for the original network in figure 4a

equation image

where An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i132.jpg is the steady state for A1 in the reduced network. One can also verify by inspection that An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i133.jpg 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

equation image

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, An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i134.jpg. 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 An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i135.jpg 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 An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i136.jpg, An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i137.jpg, one can apply the main theorem and conclude that

equation image

since An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i138.jpg 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 4e. This system satisfies the two mass conservation equations

equation image

and is therefore conservative, since every species is in one of these two equations. Also An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i139.jpg and the non-terminal complexes An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i140.jpg 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 An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i141.jpg (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

equation image

By freezing the variables An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i142.jpg and An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i143.jpg, the system is reduced to one with deficiency zero. By the theorem, the transient behaviour of I is approximately Poisson as An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i144.jpg. Since An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i145.jpg and An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i146.jpg, the growth of N is equivalent to that of Etot, Itot:

equation image

The mean λI approximates the deterministic steady state of 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 An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i147.jpg, and An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i148.jpg the total protein concentrations.

The idea is to find values of the qj, in such a way that when N = N0, it holds that An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i149.jpg and all parameters reach their experimentally determined values. The parameters kj have already been defined in (4.4) as An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i150.jpg, where

equation image

So let us define

equation image

This way when N = N0, it holds An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i151.jpg. 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 An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i152.jpg of An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i153.jpg 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 An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i154.jpg is the unique steady state of (4.3), and that there exists a sequence of positive steady states An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i155.jpg of (4.1) satisfying mass conservation (4.2), such that An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i156.jpg as An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i157.jpg. Also recall that An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i158.jpg is an absolutely robust species in (4.1).

Denote An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i159.jpg and An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i160.jpg in the special case N0 = 1. Note that

equation image

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 An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i161.jpg, for large N one could replace the mass conservation equation (4.2) with

equation image

without changing the steady-state concentration of An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i162.jpg. That is, An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i163.jpg. As An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i164.jpg,

equation image
9.1

which does not depend on N0.

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 An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i165.jpg. 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

equation image

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 [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 (An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i166.jpg) but also spontaneously promotes it (An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i167.jpg). 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, EIp lead to a the condition An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i168.jpg. 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.

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 xi(t) of the master equation (6.2) towards the solution yi(t) of (6.1) as An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i169.jpg.

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 An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i170.jpg corresponding to a Markov process with countably many variables is projected to a system An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i171.jpg 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) An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i172.jpg for every t > 0 and An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i173.jpg and (ii) An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i174.jpg, for all t > 0. See lemma 5.0.1 and theorem 5.0.2 of [29] 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 An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i175.jpg, excluding the absorbing state (N, 0). Let An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i176.jpg be the FSP of the reduced system's master equation (6.1) to SN. Assume that the system An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i177.jpg and the reduced system (6.1) both have the same initial condition b with finite support.

Proposition A.1. —

For every fixed t > 0 and state An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i178.jpg,

equation image

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 An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i179.jpg the set of states An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i180.jpg that have an outedge leading out of SN. The standard results from FSP theory guarantee that (i) An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i181.jpg for every An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i182.jpg and (ii) An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i183.jpg. Then for every An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i184.jpg,

equation image

We now bound the different terms in this expression: An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i185.jpg, An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i186.jpg, for An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i187.jpg and An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i188.jpg An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i189.jpg. Finally, define An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i190.jpg, the product Poisson distribution that is a steady state of the reduced system (6.3) by the deficiency assumptions [27]. Let An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i191.jpg be such that An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i192.jpg for all states J. Then a monotonicity property for (6.3) implies that An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i193.jpg for all s > 0 [38]. Putting all together we have

equation image
A 1

for all An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i194.jpg.

On the other hand, |I| is large for any An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i195.jpg, in the sense that there exists a constant An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i196.jpg such that An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i197.jpg for any such I. This is because there is at least one frozen variable An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i198.jpg that is off for every An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i199.jpg, and

equation image

Therefore, An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i200.jpg for An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i201.jpg.

If I is an absorbing state and An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i202.jpg is any state that has an edge pointing towards I, then An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i203.jpg for some j, and An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i204.jpg. Therefore as An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i205.jpg, |J| grows linearly for all An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i206.jpg. The right-hand side of (A 1) converges to zero since An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i207.jpg decreases exponentially as An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i208.jpg while An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i209.jpg only grows polynomially.[filled square]

One can now return to equation (6.3), the convergence of the solutions An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i210.jpg of the original master equation (6.2) towards the solution An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i211.jpg of (6.1). Suppose that an identical initial condition b = (bI) with finite support is used for both systems.

Proposition A.2. —

For every fixed t > 0 and state I,

equation image

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 An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i212.jpg as An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i213.jpg. Fix t > 0 and An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i214.jpg, and choose N1 sufficiently large that

equation image

Now we construct another FSP, that of the original master equation (6.2) with N = N2 > N1, restricted to the set An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i215.jpg, and call the solution of that projection wI(t). The intuition is that as An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i216.jpg, it holds An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i217.jpg, and at the same time one can compare wI(t) with An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i218.jpg using FSP techniques.

In fact, both the variables wI(t) and zI(t) are defined on the finite set An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i219.jpg, and as An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i220.jpg their infinitesimal generators converge. Suppose that N2 is sufficiently large that

equation image

Then

equation image

On the other hand,

equation image

Therefore for every An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i221.jpg and every sufficiently large N2,

equation image

where we are using the inequality An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i222.jpg for all An external file that holds a picture, illustration, etc.
Object name is rsif20160475-i223.jpg given by the FSP framework.[filled square]

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.

References

1. Alon U. 2007. An introduction to systems biology: design principles of biological circuits. Boca Raton, FL: Chapman & Hall/CRC.
2. Elowitz MB, Levine AJ, Siggia ED, Swain PS 2002. Stochastic gene expression in a single cell. Science 297, 1183–1186. (doi:10.1126/science.1070919) [PubMed]
3. Ozbudak EM, Thattai M, Kurtser I, Grossman AD, van Oudenaarden A 2002. Regulation of noise in the expression of a single gene. Nat. Genet. 31, 69–73. (doi:10.1038/ng869) [PubMed]
4. Altschuler SJ, Wu LF 2010. Cellular heterogeneity: do differences make a difference? Cell 141, 559–563. (doi:10.1016/j.cell.2010.04.033) [PMC free article] [PubMed]
5. Bluthgen N, Legewie S 2013. Robustness of signal transduction pathways. Cell. Mol. Life Sci. 70, 2259–2269. (doi:10.1007/s00018-012-1162-7) [PubMed]
6. Shinar G, Feinberg M 2010. Structural sources of robustness in biochemical reaction networks. Science 327, 1389–1391. (doi:10.1126/science.1183372) [PubMed]
7. Batchelor E, Goulian M 2003. Robustness and the cycle of phosphorylation and dephosphorylation in a two-component regulatory system. Proc. Natl Acad. Sci. USA 100, 691–696. (doi:10.1073/pnas.0234782100) [PubMed]
8. Tian X-J, Zhang H, Xing J 2015. Three-layer regulation leads to monoallelic olfactory receptor expression. (http://arxiv.org/abs/1505.05179)
9. Huang HJ, Velazquez JJL 2013. Bistable stochastic biochemical networks: large chemical networks and systems with many molecules. J. Math. Chem. 51, 2074–2103. (doi:10.1007/s10910-013-0200-5)
10. van Kampen NG. 2007. Stochastic processes in physics and chemistry, 3rd edn Amsterdam, The Netherlands: North-Holland Personal Library.
11. Lee CH. 2006. Stochastic analysis of biochemical reaction networks. Dissertation, University of Minnesota Mathematics Department.
12. Wilkinson DJ. 2012. Stochastic modelling for systems biology. Boca Raton, FL: CRC Press.
13. Anderson D, Enciso G, Johnston M 2014. Stochastic analysis of biochemical reaction networks with absolute concentration robustness. J. R. Soc. Interface 11, 20130943 (doi:10.1098/rsif.2013.0943) [PMC free article] [PubMed]
14. Gillespie DT. 1977. Exact stochastic simulation of coupled chemical reactions. J. Phys. Chem. 81, 2340–2361. (doi:10.1021/j100540a008)
15. Goulian M. 2010. Two-component signaling circuit structure and properties. Curr. Opin. Microbiol. 13, 184–189. (doi:10.1016/j.mib.2010.01.009) [PMC free article] [PubMed]
16. Laub MT, Goulian M 2007. Specificity in two-component signal transduction pathways. Annu. Rev. Genet. 41, 121–145. (doi:10.1146/annurev.genet.41.042007.170548) [PubMed]
17. Stock AM, Robinson VL, Goudreau PN 2000. Two-component signal transduction. Annu. Rev. Biochem. 69, 183–215. (doi:10.1146/annurev.biochem.69.1.183) [PubMed]
18. Maity AK, Bandyopadhyay A, Chaudhury P, Bank SK 2014. Role of functionality in two-component signal transduction: a stochastic study. Phys. Rev. E 89, 032713 (doi:10.1103/PhysRevE.89.032713) [PubMed]
19. Shinar G, Milo R, Rodriguez Martinez M, Alon U 2007. Input-output robustness in simple bacterial signaling systems. Proc. Natl Acad. Sci. USA 104, 19 931–19 935. (doi:10.1073/pnas.0706792104) [PubMed]
20. Straube R. 2014. Reciprocal regulation as a source of ultrasensitivity in two-component systems with bifunctional sensor kinase. PLoS Comp. Biol. 10, e1003614 (doi:10.1371/journal.pcbi.1003614) [PMC free article] [PubMed]
21. Zhu Y, Qin L, Yoshida T, Inouye M 2000. Phosphatase activity of histidine kinase EnvZ without kinase catalytic domain. Proc. Natl Acad. Sci. USA 97, 7808–7813. (doi:10.1073/pnas.97.14.7808) [PubMed]
22. Liggett TM. 2010. Continuous time Markov processes. Providence, RI: American Mathematical Society.
23. Heathcote R. 2000. The mathematics of infectious diseases. SIAM Rev. 42, 599–653. (doi:10.1137/S0036144500371907)
24. Kuske R, Greenwood P, Gordilla L 2007. Sustained oscillations via coherence resonance in SIR. J. Theor. Biol. 245, 459–469. (doi:10.1016/j.jtbi.2006.10.029) [PubMed]
25. Ingalls B. 2013. Mathematical modeling in systems biology: an introduction. Cambridge, MA: MIT Press.
26. Horn F, Jackson R 1972. General mass action kinetics. Arch. Rat. Mech. Anal. 47, 187–194. (doi:10.1007/BF00251225)
27. Anderson D, Craciun G, Kurtz TG 2011. Product-form stationary distributions for deficiency zero chemical reaction networks. Bull. Math. Biol. 72, 1947–1970. (doi:10.1007/s11538-010-9517-4) [PubMed]
28. Joshi B, Shiu A 2013. Atoms of multistationarity in chemical reaction networks. J. Math. Chem. 51, 153–178. (doi:10.1007/s10910-012-0072-0)
29. Munsky B. 2008. The finite state projection approach for the solution of the master equation and its applications to stochastic gene regulatory networks. Dissertation, UC Santa Barbara.
30. Feinberg M. 1979. Lectures on chemical reaction networks. Madison, WI: The Mathematical Research Center, University of Wisconsin; See http://www.chbmeng.ohio-state.edu/?feinberg/research/.
31. Gunawardena J. 2003. Chemical reaction network theory for in silico biologists. Lecture Notes, Bauer Center for Genomics Research Cambridge, MA: Harvard University.
32. Feinberg M. 1995. Multiple steady states for chemical reaction networks of deficiency one. Arch. Rat. Mech. Anal. 132, 371–406. (doi:10.1007/BF00375615)
33. Banaji M, , Pantea C.2016. Some results on injectivity and multistationarity in chemical reaction networks. SIAM J. App. Dyn. Syst.15, 807–869. ( doi:10.1137/15M1034441)
34. Joshi B, , Shiu A. Which small networks are multistationary? (http://arxiv.org/abs/1603.06441. )
35. Karp RL, Perez-Millan M, Dasgupta T, Dickenstein A, Gunawardena J 2012. Complex-linear invariants of biochemical networks. J. Theor. Biol. 311, 130–138. (doi:10.1016/j.jtbi.2012.07.004) [PubMed]
36. Dexter JP, Gunawardena J 2013. Dimerization and bifunctionality confer robustness to the isocitrate dehydrogenase regulatory system in Escherichia coli. J. Biol. Chem. 288, 5770–5778. (doi:10.1074/jbc.M112.339226) [PMC free article] [PubMed]
37. Munsky B, Khammash M 2008. The finite state projection approach for the analysis of stochastic noise in gene networks. IEEE Trans. Autom. Contr. 53, 201–214. (doi:10.1109/TAC.2007.911361)
38. Smith H. 2014. Monotone dynamical systems, Mathematical Surveys and Monographs 41 Providence, RI: American Mathematical Society.

Articles from Journal of the Royal Society Interface are provided here courtesy of The Royal Society