Product-Form Stationary Distributions for Deficiency Zero Networks with Non-mass Action Kinetics.

In many applications, for example when computing statistics of fast subsystems in a multiscale setting, we wish to find the stationary distributions of systems of continuous-time Markov chains. Here we present a class of models that appears naturally in certain averaging approaches whose stationary distributions can be computed explicitly. In particular, we study continuous-time Markov chain models for biochemical interaction systems with non-mass action kinetics whose network satisfies a certain constraint. Analogous with previous related results, the distributions can be written in product form.

Introduction

Biological interaction systems are typically modeled in one of three ways. If the counts of the constituent species are high, then their concentrations are often modeled via a system of ordinary differential equations with state space , where is the number of species. If the counts are moderate (perhaps order or ), then they may be approximated by some form of continuous diffusion process (Gillespie 2000; Van Kampen 2007). If the counts are low, then the system is typically modeled stochastically as a continuous-time Markov chain in (Anderson and Kurtz 2011, 2015). We often want to understand the stationary behavior of the model under consideration. For deterministic models, understanding the stationary behavior usually entails characterizing the stable fixed points of the system, whereas for stochastic models, we require the calculation of the stationary distribution.

Stationary distributions are also useful in a multiscale setting, where the stationary statistics of a fast subsystem can be utilized in the approximation of the dynamics of the slow variables, which are typically of most interest. When an analytical form for the stationary distribution of the fast subsystem is not known, numerical approximations can be used. However, these computations are often expensive and part of an “inner loop,” typically making this calculation the rate-limiting step of the analysis. In the context of biochemical reaction networks, quasi-equilibrium (QE)-based approximations lead to fast subsystems which preserve mass action kinetics (Goutsias 2005; Janssen 1989; Thomas et al. 2012; Cao et al. 2005; Weinan et al. 2005). However, more recent improvements in stochastic averaging can lead to fast subsystems with non-mass action kinetics, and this observation was the motivation for the present work (Cotter 2016; Cotter and Erban 2016; Cotter et al. 2011).

One class of interaction networks that has been quite successfully analyzed, and that appears ubiquitously as fast subsystems, are those that are weakly reversible and have a deficiency of zero (see Appendix). For this class of models, and under the assumption of mass action kinetics, the fixed points of the deterministic models (Anderson 2011, 2008; Craciun 2015; Feinberg 1979, 1987; Gunawardena 2003) and the stationary distributions for the stochastic models, have been fully characterized (Anderson et al. 2010, 2015; Cappelletti and Wiuf 2016; Van Kampen 1976). In fact, it is the study of this class of networks that is largely responsible for the development of the field of chemical reaction network theory (Feinberg 1972, 1979; Gunawardena 2003; Horn 1972), a branch of applied mathematics in which the dynamical properties of the mathematical model are related to the structural properties of the interaction network.

In this article, we return to stochastically modeled interaction networks that are weakly reversible and have a deficiency of zero, though we consider propensity functions (also called intensity functions or rate functions) that are more general than mass action kinetics. However, we add a certain condition to the rates within the reaction network (see Assumption 1 below). Following Anderson et al. (2010), which was motivated by the work of Kelly (1979) who discovered the product-form stationary distribution of certain stochastically modeled queuing networks, we provide the form of the stationary distribution for this class of models. In particular, and in similarity with the main results of Anderson et al. (2010), in Theorem 2 we show that the distribution is of product form and that the key parameter of the distribution is a complex-balanced equilibrium value of an associated deterministically modeled system with mass action kinetics. The result of Anderson et al. (2010) can now be viewed as a special case of this new result.

The paper proceeds as follows: In Sect. 2, we introduce the formal mathematical model of interest. In Sect. 3, we provide the main theorem of this article, which characterizes the stationary distribution for the class of models of interest. In Sect. 4, we provide a series of examples which demonstrate the usefulness of the main result. Some brief concluding remarks are given in Sect. 5.

We assume throughout that the reader is familiar with terminology from chemical reaction network theory. However, we provide in Appendix all necessary terminology and results from this field that are used in the present work.

Mathematical Model

We consider a system with d chemical species, , undergoing reactions which alter the state of the system. For concreteness, we suppose there are distinct reaction channels. For the kth reaction channel, we denote by the vectors representing the number of molecules of each species consumed and created in one instance of that reaction, respectively. Note that is the net change in the system due to one instance of the kth reaction. We associate each such and with a linear combination of the species in which the coefficient of is , the ith element of . For example, if for a system consisting of two species, then we associate with the linear combination . Under this association, each and is termed a complex of the system. We denote any reaction by the notation , where is the source, or reactant, complex and is the product complex. We note that each complex may appear as both a source complex and a product complex in the system. The set of all complexes will be denoted by .

Definition 1

Let , and denote the sets of species, complexes, and reactions, respectively. The triple is called a chemical reaction network.

Throughout, we assume that for each .

Deterministic Model

The usual deterministic model for a chemical reaction network assumes that the vector of concentrations for the species satisfies a differential equation of the formwhere is the (state-dependent) rate of the kth reaction channel. If we assume that each function satisfies deterministic mass action kinetics, then

Definition 2

An equilibrium value of (1) is said to be complex balanced if for each ,where the sum on the left, respectively right, is over those reactions with as source complex, respectively product complex. In the special case of mass action kinetics, c is complex balanced if and only if

Stochastic Model, Previous Results, and Assumptions

The usual stochastic model for a reaction network treats the system as a continuous-time Markov chain for which the rate of transition from state to state is , where is a suitably chosen intensity function (the intensity functions are also termed propensity functions in the literature). This stochastic process can be characterized in a variety of useful ways (Anderson and Kurtz 2011, 2015). For example, it is the stochastic process with state space and infinitesimal generatorwhere . Kolmogorov’s forward equation for this model, termed the chemical master equation in the biology literature, is where is the probability the process is in state at time , given an initial distribution of , and is the adjoint of A. Note that (3) implies that a stationary distribution for the model, , must satisfyNote also that (4) implies that is in the null space of . Thus, and as is well known, finding reduces to solving , with , where A is reinterpreted as a generator matrix.

Of particular interest to us are the rate functions, . Under the assumption of stochastic mass action kinetics, we havewhere are the reaction rate constants. Here and throughout, we interpret any product of the form to be equal to one. In Anderson et al. (2010), the authors considered a class of stochastically modeled reaction networks with mass action kinetics (those with a deficiency of zero and are weakly reversible, see Appendix) and characterized their stationary distributions as products of Poisson distributions. However, they did not just consider mass action kinetics in Anderson et al. (2010), but any kinetics satisfying the functional formso long as for each i. In particular, they proved the following.

Theorem 1

(Anderson et al. 2010) Let be a reaction network, and let be a choice of positive rate constants. Suppose that modeled deterministically with mass action kinetics and rate constants , the system is complex balanced with complex-balanced equilibrium . Then, the stochastically modeled system with intensity functions (6) admits the invariant measureThe measure can be normalized to a stationary distribution so long as it is summable.

The connection between weakly reversible, deficiency zero networks, and complex-balanced equilibria is given in Appendix.

In this paper, we generalize Theorem 1 by showing how to find the stationary distributions for models with rates that do not seem to satisfy the form (6). However, we add an assumption pertaining to the form of the rates within the reaction network, which we describe now. We begin by partitioning the set of species into two sets, . We will say that if . Otherwise, we say that , noting that . We now assume that the intensity functions for the reaction network are of the formwhere and if and only if .

Assumption 1

A stochastically modeled reaction network satisfies this assumption if it satisfies the partition described above and has intensity functions of the form (8).

We provide an example to clarify the notation.

Example 1

Consider the stochastically modeled system with reaction networkwhere the intensity functions are placed next to the reaction arrows. Here , with , and . The assumption (8) then supposes that for appropriate functions , we haveFor example, valid choices include , , and , in which case the form for the stationary distribution does not follow immediately from Theorem 1.

Main Result

Here we state and prove our main result.

Theorem 2

Let be a reaction network satisfying Assumption 1, and let be a choice of positive rate constants. Suppose that modeled deterministically with mass action kinetics and rate constants the system is complex balanced with complex-balanced equilibrium . Then the stochastically modeled system with intensity functions (8) admits the invariant measureThe measure can be normalized to a stationary distribution so long as it is summable.

Note that the theorem applies to models with reaction networks satisfying Assumption 1 and that are weakly reversible and have a deficiency of zero. See Appendix.

Proof

First note that if , then Theorem 2 is the same as Theorem 1 and there is nothing to show. Thus, we suppose .

The proof proceeds in the following manner. First, for each we will demonstrate the existence of a function for whichNext, we will apply Theorem 1 and prove that the resulting distribution is indeed given by (10).

Let . We begin by settingFor an integer, we may define recursively via the formulaNote that is a well-defined function since for each by assumption. It is clear that (11) is satisfied with this choice of .

For , we may now writeHence, we may apply Theorem 1 and conclude thatis an invariant measure for the system, where c is a complex-balanced fixed point for the deterministic system. It remains to show thatFirst note that if , then both sides of (13) are equal to one. For the time being, assume that . Under this assumption, the right-hand side of (13) iswhere in the final equality we used that for (and that ), and the left-hand side iswhere we again used that when . Hence, (13) is verified when .

We will now prove that (13) holds in general by induction. We suppose that (13) holds for all , where , and will show it to hold at . Using (12), the left-hand side of (13) evaluated at iswhere the final equality is an application of (12) with in place of the variable z. Continuing, we havewhere the final equalities are straightforward. Note that (14) is the right-hand side of (13) evaluated at , and so the proof is complete.

Examples

Example 1: Motivating Example

First, we consider a motivating example arising from model reduction, through constrained averaging (Cotter 2016; Cotter and Erban 2016; Cotter et al. 2011), of the following system:where the intensity functions are placed next to the reaction arrows. Note that the intensities of all the reactions follow mass action kinetics. We consider this system in a parameter regime where the reversible dimerization reactions are occurring more frequently than the production of and the degradation of . Both and are changed by the fast reactions, but the quantity is invariant with respect to the fast reactions, and as such is the slow variable in this system. We wish to reduce the dynamics of this system to a model only concerned with the possible changes in S:where and are the effective rates of the system.

Using the QE approximation (QEA), and , where is the stationary distribution for the systemunder the assumption that . Since the system (17) satisfies the necessary conditions of the results of Anderson et al. (2010) (weak reversibility and deficiency of zero), the invariant distribution is known exactly.

In comparison, the constrained approach requires us to find the invariant distribution of the following system:subject to . Readers interested in seeing how this is derived should refer to Cotter (2016). This network is weakly reversible and has a deficiency of zero. However, the form of the rates in this system does not satisfy the conditions specified in Anderson et al. (2010). In the context of constrained averaging, this lack of a closed form for the stationary distribution would result in the need for some form of approximation of the stationary distribution. There are two common methods utilized for performing this approximation. One possibility would be to perform exhaustive stochastic simulation of the system (18). Another option involves finding the distribution by finding the null space of the adjoint of the generator (see the discussion in and around (5)). However, as the state space of (18) will typically be huge, the latter method often involves truncating the state space and approximating the actual distribution with that of the stationary distribution of the truncated system (Cotter 2016). Both approaches will lead to approximation errors and varying amounts of computational cost. However, note that the system (18) does satisfy Assumption 1, with and . We denote the rate of dimerization by and its reverse by . Therefore,with defined in the final equality. The form of the rate of the reverse reaction is much simpler and is given bywhich defines .

By Theorem 2, we can write down the stationary distribution of this system. The complex-balanced equilibrium of the associated deterministically modeled system with is given by . Then by Theorem 2, and by recalling that all states in the domain satisfy , the stationary distribution for is given bywhere is a normalization constant and s is the conserved quantity. Note that the indicator function in (in the denominator) has disappeared since it is always equal to one over the domain of the product.

We can compare (19) with the distribution of (17), which arises from the QEA, and also with the distribution of the full system (15) conditioned on (which can be approximated by finding the null space of the adjoint of the generator of the full system on the truncated domain). First, we consider the QEA approximation. The invariant distribution of the fast subsystem (17) can be found using Theorem 1 and is given bywhere is the complex-balanced equilibrium for this system satisfying and is a normalizing constant.

Since the full system (15) does not have a deficiency of zero, we are not able to find its invariant distribution directly. However, by truncating the state space appropriately, we are able to approximate the full distribution by constructing the generator on this truncated state space and finding the null space of the adjoint.

Once we have approximated the null space of the truncated generator, we can find the approximation of by taking the probabilities of all states with and renormalizing. In what follows, we truncated the domain of the generator to .

We consider the system (15) with parameters given by:Note that it is not obvious from these rates that the reactions with rates and are in fact the slow reactions in this system. The invariant density is largely concentrated in a small region centered close to the point . By using the approximation of the invariant density that we have computed on the truncated domain, we can compute the expected ratio between occurrences of the fast reactions with rates and with the slow reactions with rates and . For this choice of parameters, the expected proportion of the total reactions which are fast reactions (dimerization/disassociation) is . This indicates a difference in timescales between these reactions, but the difference is not particularly stark, and as such we would expect there to be significant error in any approximation relying on the QEA.

Figure 1 shows the three approximations of the distribution for the system (15) with parameters given by (21). The constrained and QEA approximations are computed using (19) and (20), respectively, with the normalizing constants computed numerically. As would be expected in this parameter regime, the constrained approximation is far more accurate than the QEA.

Fig. 1

(Color figure online) Approximations of the distribution for the system (15) with parameters given by (21) using constrained averaging, QEA averaging, and through approximation of the invariant distribution of the full system on .

We can quantify the accuracy of each of the approximations by computing the relative differences with the distribution computed using the full generator on the truncated domain. This relative difference was for the QEA, in comparison with for the constrained approximation. This demonstrates the improvement in approximation that can be achieved by using constrained averaging, and which motivates the need for results like Theorem 2 which take non-mass action kinetics into account.

Example 2: Dimerization

We begin by considering a model consisting of only proteins, denoted P, and dimers, denoted D. We suppose there are two mechanisms by which the proteins are dimerized: random interactions between the protein molecules, and via a catalyst. The rate at which random interactions lead to the formation of dimers can be taken to be of mass action form. Assuming the concentrations are such that the catalyst is acting at capacity, the rate of formation of the dimers due to the catalyst can be faithfully modeled as a constant (so long as there are at least two proteins to make the reaction happen). Thus, letting and denote the numbers of proteins and dimers in the model, respectively, the reaction system can be represented aswhere , and are given parameters. Note that after an obvious change of variables the stationary distribution of this particular model is provided in (19).

The reactions in (22) are typically a subset of the reactions in a larger system. For example, the actual model of interest may bewhere M represents an mRNA molecule. This is a standard model for dimer production. Depending upon the relevant timescales in the system, we may want to take . If the reactions are appreciably faster than those of (23), then an obvious path simulation strategy presents itself in the mold of the slow-scale SSA (Cao et al. 2005) (i.e., stochastic averaging):Being able to analytically calculate the stationary distribution in step 1 allows us to bypass the need to numerically approximate the stationary distribution, as is commonly done (Weinan et al. 2005; Weinan and Vanden-Eijnden 2007).

for the current values of P and D, find the stationary distribution of the fast subsystem analytically,

determine the effective rates of the reduced model where the expectations are with respect to the distribution found in step 1.

simulate forward in time using the stochastic simulation algorithm (Gillespie 1976) or the next reaction method (Anderson 2007; Gibson and Bruck 2000), and return to step 1.

Example 3

In Sects. 4.1 and 4.2, we applied Theorem 2 on the common motif . In this example, we present another common motif, , for which Theorem 2 is also useful. As opposed to the specific models we considered in the previous examples, here we present a more general framework in which the specific form of the propensity function for the reaction is arbitrary. This situation is common when undertaking certain types of averaging arguments (Cotter 2016).

Let us suppose that the effective dynamics of a slow variable in a larger system can be modeled aswhere and are positive constants and where is a conditional expectation of the fast variables , conditioned on . The conditional expectation could in general be highly nonlinear, and not of the form (6) required by Theorem 1. Supposing that if , Theorem 2 says that, up to a normalization constant, the invariant distribution of S is given byIn practice, the normalization constant could be approximated by summing over an appropriate domain.

Example 4

This example will demonstrate the difficulties that can arise when a reaction is added that involves a species in the set with multiplicity not equal to . In particular, consider the systemwherewhere . Note that we have a non-mass action kinetics rate for the reaction . There is another reaction involving , but the amount of molecules involved in this reaction is not a multiple of 2. This means that we cannot apply the result of Theorem 2 to this system unless satisfies the recurrence relation detailed in the proof of Theorem 2. In particular, it must satisfy , orThis recurrence relation defines a unique function for each . In general, the function can oscillate wildly. Let us consider, for example, when is given as above. In this case,

Figure 2 demonstrates how and the amplitude of the oscillations grow with , for the case . It is clear that this function does not represent any physical reaction rate arising from chemistry, but we can still write down the invariant distribution of the system (25).

Fig. 2

(Color figure online) as given in (25) with .

The complex balance equilibrium of the associated mass action kinetics system to (25) is given by . Therefore, the stationary distribution is given by:where is a normalizing constant and is given by (26) with . Note that the value of a here dictates the oddness or evenness of the quantity , which is preserved by each of the reactions.

Figure 3 shows the invariant distribution (27) of the system (25) with an even initial value of , and for the following parameters:The normalization constant was approximated by summing the values of all states .

Fig. 3

(Color figure online) Stationary distribution (27) of the system (25) with parameters given by (28), with the assumption that the initial value of is even. The normalization constant was approximated by summing the value of all states .

Conclusions

In this paper, we provided the stationary distributions for the class of stochastically modeled reaction networks with non-mass action kinetics that satisfy our Assumption 1 and that admit a complex-balanced equilibrium when modeled deterministically with mass action kinetics. Similarly to the results of Anderson et al. (2010), we showed that the stationary distributions are of product form.

We motivated the need for such results through consideration of modern averaging techniques. In particular, Theorem 2 significantly reduces the computational cost of finding the invariant distribution of most fast subsystems in a multiscale setting, therefore making accurate approximations very cheap to compute. This in turn opens up more possibilities, for example approximation of likelihoods via multiscale reductions in the context of parameter inference for biochemical networks.