Reaction networks are systems in which the populations of a finite number of species evolve through predefined interactions. Such networks are found as modeling tools in many biological disciplines such as biochemistry, ecology, epidemiology, immunology, systems biology and synthetic biology. It is now well-established that, for small population sizes, stochastic models for biochemical reaction networks are necessary to capture randomness in the interactions. The tools for analyzing such models, however, still lag far behind their deterministic counterparts. In this paper, we bridge this gap by developing a constructive framework for examining the long-term behavior and stability properties of the reaction dynamics in a stochastic setting. In particular, we address the problems of determining ergodicity of the reaction dynamics, which is analogous to having a globally attracting fixed point for deterministic dynamics. We also examine when the statistical moments of the underlying process remain bounded with time and when they converge to their steady state values. The framework we develop relies on a blend of ideas from probability theory, linear algebra and optimization theory. We demonstrate that the stability properties of a wide class of biological networks can be assessed from our sufficient theoretical conditions that can be recast as efficient and scalable linear programs, well-known for their tractability. It is notably shown that the computational complexity is often linear in the number of species. We illustrate the validity, the efficiency and the wide applicability of our results on several reaction networks arising in biochemistry, systems biology, epidemiology and ecology. The biological implications of the results as well as an example of a non-ergodic biological network are also discussed.
Let us now formally describe reaction networks. Motivated by the literature on chemical kinetics, we refer to the network participants as molecules which may belong to one of species. There are reactions in the network and for any, the stoichiometric vector denotes the change in the number of molecules in each of the species due to the -th reaction. Consider the deterministic model for the reaction network described above. In this setting, the state of the system is described by a vector of concentrations of the species which we denote by. The concentration of a species is simply its molecular count divided by the system volume. Let be the flux associated with the -th reaction (see [8]). To ensure positivity of the system, we require that whenever and. If the initial state is, then the evolution of concentrations is given by which satisfies the Reaction Rate Equations (RRE) of the form (1) We are interested in the long-term behavior and stability of our reaction dynamics. More precisely, we would like to check if the following conditions are satisfied. DC1 For any, there is a compact set such that for all. DC2 There exists a compact set such that for any, we have for large values of. DC3 There is a such that for any we have as. The first condition, DC1, says that for any, the entire trajectory stays within some compact set. We would expect this to be true for most realistic systems. Hence a violation of this property may suggest a flaw in the deterministic model. The second condition, DC2, says that there is an attractor set for the dynamics, where all the trajectories eventually lie, irrespective of their starting point. The last condition, DC3, says that there is a globally attracting fixed point for the deterministic model. Using techniques from the theory of dynamical systems [30], [46], one can verify these conditions, without the need of simulating the deterministic model. There is also a general theory to check condition DC3 for reaction networks satisfying mass-action kinetics (see [47]–[50]). Broadly speaking, these three conditions present different ways of saying that the reaction dynamics is “well-behaved”. Our goal in this paper is to develop a theoretical and computational framework for verifying conditions similar to DC1, DC2 and DC3 for stochastic models of reaction networks. Consider the stochastic model corresponding to the reaction network described above. In this setting, the firing of reactions are discrete events and the state of the system refers to the vector of molecular counts of the species. When the state is, the -th reaction fires after a random time which is exponentially distributed with rate. The functions are known as the propensity functions in the literature. To ensure positivity of the system, we require that if, then, where is the set of non-negative integers. The dynamics can be represented by the Markov process where is the initial state. Note that if, then is the number of molecules of at time. It is important to select a suitable state space for the Markov process representing the reaction dynamics. We choose to be a non-empty subset of satisfying the following properties: Observe that part (A) ensures that if then for all and hence can be taken to be the state space of all the Markov processes describing the stochastic reaction network with an initial state in. Part (B) implies that the reaction dynamics cannot be contained in a proper subset of. The role of this assumption will become clear in the next section, when we discuss the issue of state space irreducibility. Note that in certain cases, such as the pure-birth network, a suitable state space satisfying the above criteria cannot be found. There also exist cases where the above criteria restricts the choice of state space. For example, for the pure-death network, the only possible choice for state space is. Finally we remark that if the reactions in a network satisfy a conservation relation then the state space must be chosen with an initial condition in mind. For example, for the network, the sum of molecular counts of and is preserved by the reactions. Hence if we wish to study the stochastic dynamics with the initial sum as, then the correct choice for state space is. Let denote the space of probability distributions over, endowed with the weak topology which is metrized by the Prohorov metric (see [51]). For any let denote the following probability (2) Defining, for any, we can view as an element in. In fact, is the distribution at time of the Markov process. The dynamics of is given by the Chemical Master Equation (CME) which has the following form: (3) where if and for all. Theoretically, one can find for any and, by solving this system. However this system consists of as many equations as the number of elements in. Hence an explicit solution is only possible when is finite, which only happens in very restrictive cases where all the reactions preserve some conservation relation. Typically, is infinite and solving this system analytically or even numerically is nearly impossible, except in some restrictive cases (see [16]). From now on, we assume that is infinite. The above discussion shows that at the level of distributions, we can view the stochastic dynamics as the deterministic dynamics, which satisfies the CME. However, the major difficulty in analyzing this deterministic dynamics is that it occurs over an infinite dimensional space. Nevertheless we can recast the conditions DC1, DC2 and DC3 in the stochastic setting as below. SC1 For any, there is a compact set such that for all. SC2 There exists a compact set such that for any we have for large values of. SC3 There is a such that for any we have as. Each of the above conditions give an important insight about the long-term behavior and stability of the stochastic dynamics. The first condition, SC1, says that for every we can find a finite set such that each puts at least of its mass in. In other words, the probability that the state of the underlying Markov process at any time is inside is greater than. We would expect this to be true for most realistic models. If condition SC2 holds then the evolution of distributions have a compact attractor set in, where all the trajectories eventually lie irrespective of their starting point. This suggests that in the long run, the family of processes, spend most of their time on the same set of states. The last condition SC3 says that the evolution of distributions have a globally attracting fixed point. If this holds, then the Markov process representing the reaction dynamics is ergodic with as the unique stationary distribution. For understanding the long-term behavior of a stochastic process, ergodicity is a desirable property to have. In the long-run, the proportion of time spent by any trajectory of an ergodic process, in any subset of the state space is equal to the stationary probability of that subset (see (12) ). In other words, information about the stationary distribution can be obtained by observing just one trajectory for a sufficiently long time. Such a result can have important applications. For example, consider a culture with a large number of identical cells with each cell having the same reaction network. If we can show that this intracellular network is ergodic, then by observing the long-term reaction dynamics in a single cell, using for example. time-lapse microscopy, we can obtain statistical information about all the cells at stationarity. Conversely, ergodicity allows us to obtain the stationary distribution of a single-cell by observing the distribution over the population, using for example flow cytometry. In this paper we develop a general framework for checking conditions SC1, SC2 and SC3. However, the scope of our paper is broader than that. As mentioned in the introduction, we obtain easily computable bounds for the statistical moments of the underlying Markov process and investigate when these moments converge with time. We also present conditions for the distribution of the process to be light-tailed.
In this section we discuss the main results of our paper. In particular, we explain how conditions SC1, SC2 and SC3 can be verified without having to simulate the trajectories of the Markov process representing the reaction dynamics. Intuitively, these conditions can only hold if the Markov process has a low probability of hitting states that have a very large size. In our case, the states are vectors in and so we can measure their size by using any norm on. The central theme of this paper is to demonstrate that for many networks, long-term behavior can be easily analyzed by choosing the right norm for measuring the state sizes. This right norm has the form (4) where is a positive vector in satisfying the following condition. In this section, we formally present the main results of our paper. Their proofs are given in the Supplementary Material S1.
Let us then consider a unimolecular reaction network which involves species that interact through reaction channels of the form: (15) where, , and. The reaction rates, and are positive real numbers. In accordance with (3), the reactions are indexed from to, and corresponding propensities and stoichiometries are denoted by and, respectively. Similar results are now presented for stochastic bimolecular reaction networks which, in addition to the unimolecular reactions (15), also involve bimolecular reactions of the form: (18) where, , , and. The reaction rates and are positive real numbers. In this section we illustrate that stochastic and deterministic models of the same reaction network may exhibit very different qualitative behaviors. Therefore assessing ergodicity or the convergence of moments of a stochastic model from the stability properties of the corresponding deterministic model is, in general, incorrect. To support this claim, we consider two reaction networks. The goal of this section is to compute a compact set that is attractive for the first-order moment of using the optimization problems (17) or (20). Due to the moment closure problem [54], analytical expressions for the steady-state values of the moments of bimolecular reaction networks are not available, and hence this is an important class of networks to analyze. Consider the following bimolecular reaction network (29) representing a dimerization process, i. e. dimerizes to. It is easily seen that this network is irreducible since any point in the state-space can be reached from any other point in a finite number of reactions having nonzero propensities. Choosing in, e. g. , yields that and, hence the network is exponentially ergodic, and all the moments are bounded and converging. On solving the optimization problem (20) with numerical values, , and, we get that which coincides with the theoretical value. One can regard to be an attractive compact set in which the first-order moments of eventually lie. To validate this calculation, Monte-Carlo simulations were performed which yield (30) showing the correctness of the attractive compact set. To further illustrate this result, several trajectories of and for different initial conditions are plotted in Figure 3. We now discuss how the computation of an attractive compact set for the first-order moments can be used to assess whether a closure method leads to a result that is consistent with the stochastic dynamics. The idea is to check whether the closed system converges towards a value which lies within the compact set. Let us consider the reaction network (29) and close the first-order moments equations by neglecting the second order cumulant, i. e. neglecting the variance. By doing so, we get the model (31) where and are the approximate first-order moments of the system. The unique positive equilibrium point for this model is given by (32) where. With the same parameter values as before, we find that and and therefore for, showing that the state of the closed system converges to the boundary of the compact set. Note that SSA also predicts that the trajectories of the first-order moments converge to the boundary of this set. However the actual equilibrium values for the first-order moments of the stochastic dynamics are and, which differ from the ones obtained with the closure method. This discrepancy is expected since the variance has been neglected. This example shows how attractive compact sets for the moments can be used as a test for the moment-closure methods by checking whether the closed system predicts trajectories that converge inside those compact sets. However, note that in the current state, these compact sets can only be used to obtain a lower bound on the closure-error whenever the trajectories of the closed dynamics converge to a point outside the compact set. In such a case, the lower bound on the closure-error is simply given by the distance between the equilibrium point of the closed-system (33) where is the attractive (convex) compact set and is the equilibrium point of the closed dynamics. Let us consider the feedback loop network of Figure 4 represented by the reaction network (34) where is mRNA and is the corresponding protein. The dimer acts back on the gene expression through an arbitrary bounded nonnegative function. We have the following result: Result 16 For any positive values of the rate parameters and any bounded nonnegative function, the feedback loop with dimerization (34) is ergodic and all the moments are bounded and globally converging. Let us consider the stochastic switch of [63] described by the unimolecular stochastic reaction network (35) Above and represent mRNAs and proteins of gene, respectively. The functions and are arbitrary bounded nonnegative functions. We have the following result: Result 17 For any positive values of the rate parameters and any bounded nonnegative functions and, the stochastic switch (35) is ergodic and all the moments are bounded and globally converging. We consider here the stochastic repressilator of Figure 5 (see also [42]) involving genes. The reaction network corresponding to this -gene repressilator is given by (36) where, . Above, and are the mRNA and protein corresponding to gene. We have the following result: Result 18 For any positive values of the rate parameters and, the stochastic -gene repressilator (36) is ergodic and all the moments are bounded and globally converging. We consider here the following SIR-model, similar to the one in [64], defined as (37) where birth and death reactions represent people entering and leaving the process, respectively. The only bimolecular reaction is the contamination reaction which turns one susceptible person into an infectious one. The two last reactions represent how infectious people are recovering and how recovered people become susceptible again. We then have the following result: Result 19 For any positive values of the rate parameters, the SIR-model (37) is ergodic and all the moments are bounded and globally converging. Let us consider the circadian oscillator of [65], depicted in Figure 6, which is a network involving 9 species and 18 reactions. Applying the developed theory on this model, we obtain the following result: Result 20 For any positive values of the rate parameters, the circadian clock model of [65] is ergodic and all the moments are bounded and globally converging. Using, for instance, the values of [65] and solving for the optimization problem (20) using linprog and Yalmip [66], we find that and. Typical trajectories for the proteins A, R and C are depicted in Figure 7 where we can observe the expected oscillatory behavior. When averaging the populations of the proteins A, R and C over a population of 2000 cells, we obtain the sample-average trajectories depicted in Figure 8. Convergence to stationary values is easily seen. Moreover, from the ergodicity property, we can even state that these fixed points for the sample-averages are globally attracting and that they coincide with the asymptotic time-average (dashed lines). The steady-state average values for the proteins A, R and C are given by 222. 1797,534. 8853 and 549. 7195, respectively. Let us consider one of the oscillatory p53 models of [67], which is described by the reactions (38) where is the number of p53 molecules, the number of precursor of Mdm2 molecules and the number of molecules of Mdm2. The function implements a nonlinear feedback on the degradation rate of p53. We have the following result: Result 21 For any positive values of the rate parameters, the oscillatory p53 model (38) is ergodic and all the moments are bounded and globally converging. We consider here the stochastic reaction network (39) which is an open analogue of the deterministic Lotka-Volterra system of [68]. The first set of reactions represent immigration, the second one reproduction, the third one competition due to overpopulation and the last one deaths/migrations. We obtain then the following result, which is a stochastic analogue of the results in [69] obtained in the deterministic setting: Theorem 22 Let us define and assume that one of the following conditions hold: Then, the stochastic reaction network (39) is ergodic and all the moments up to order are bounded and globally converging. In order to illustrate that the method can be applied to systems with more general mass-action kinetics, we consider the stochastic version of the well-known Schlögl model [70]: (40) where is the main molecule in the network. The above model is derived in the supplementary material S1 where we have assumed that the other molecular populations do not vary over time. Note that in the present form the model has an infinite state-space and involves a single trimolecular reaction. We then have the following result. Theorem 23 For any positive values of the rate parameters and any positive values for and, the Markov process describing the Schlögl model (40) is exponentially ergodic. Note, however, that we cannot say anything on the stability of the moments (besides the fact that the first order-moment converges) since the condition DD2 does not hold here due to the presence of a cubic term. Note that extending the condition DD2 to handle more general cases, such as this one, might be possible.
The central theme of this paper is to verify the ergodicity and moment boundedness of reaction networks in the stochastic setting. Note that even though we mainly consider mass-action kinetics in this paper, the framework also applies to more general kinetics described, for instance, by Hill functions (see the examples on the repressilator and the stochastic switch) and more general mass-action kinetics. These results have several interesting and important biological implications. For example, the ergodicity of a network shows that population-level information could be obtained by observing a single trajectory for a long time. Such an insight can be used to leverage different experimental techniques for a given application. For example, consider a clonal cell population with each cell having a gene-expression network that is ergodic. Then the stationary distribution (at the population level) of the species involved in this network can be ascertained by observing a single cell over time. In other words, to obtain stationary distributions one can either collect samples over time from a single cell (e. g. using time-lapse microscopy) or one can take a snapshot of the entire cell population at some fixed time (e. g. using flow-cytometry). Due to ergodicity, both these approaches will yield the same information. Hence, far from being a technical condition, ergodicity can have far reaching experimental implications. As a property of a network, ergodicity also sheds important light on the long range behaviors that can be exhibited by that network. One may expect that most endogenous biochemical networks to be ergodic in order to achieve robustness with respect to variability in initial conditions and kinetic parameters, thus ensuring proper biological functions in spite of environmental disturbances. As also mentioned in the introduction, ergodicity is a non-trivial property which needs to be carefully established and cannot be generically assumed. To illustrate this, let us consider a simplified version of the model of carcinogenesis considered in [71] which is given by (41) where, . When, the trajectories of the species grow unbounded, as shown in Figure 9, emphasizing then non-ergodicity of the model for this choice of parameters. The ideas we use for analysis can also be applied for rationally designing circuits in synthetic biology, where it is important that the network be (structurally) ergodic in order to ensure that the dynamics has the desired behavior irrespective of the initial conditions. Such a design is crucial because the initial conditions are usually unknown or difficult to control at certain times, e. g. after cell division or after the transfection of plasmids in the cell. Our results on boundedness and convergence of statistical moments enable verification of the suitability of a stochastic model and to characterize the properties of its steady-state distributions, even if such a distribution is not explicitly computable. One application of this is to provide justifications and insights for using moment closure techniques which have been extensively used to study stochastic chemical reaction networks. Some of these techniques [72], [73] are based on manipulations of the moment generating function of the underlying stochastic process. The existence of this moment generating function is implicitly assumed in such techniques but it may not always hold, thereby jeopardizing the validity of the technique. In this article, we show that under certain conditions, the distribution of the stochastic process is uniformly light-tailed, which proves that the moment generating function exists for all time. Certain moment closure techniques (see [74], [75]) prescribe ways to approximate higher order moments as a function of lower order moments. Such an approximation is, however, only reasonable if the higher order moments are bounded over time. This can be easily assessed with our approach and one can even quantify the error by explicitly computing the moment bounds as described in this article. Finally, the techniques developed here will prove invaluable for designing synthetic biological control systems and circuits whose objective is to steer the moments of the network of interest to a specific steady-state value. Until now, no theory has provided guidance for such a design. The specifics are outside the scope of this article and will be pursued elsewhere.
In many biological disciplines, computational modeling of interaction networks is the key for understanding biological phenomena. Such networks are traditionally studied using deterministic models. However, it has been recently recognized that when the populations are small in size, the inherent random effects become significant and to incorporate them, a stochastic modeling paradigm is necessary. Hence, stochastic models of reaction networks have been broadly adopted and extensively used. Such models, for instance, form a cornerstone for studying heterogeneity in clonal cell populations. In biological applications, one is often interested in knowing the long-term behavior and stability properties of reaction networks even with incomplete knowledge of the model parameters. However for stochastic models, no analytical tools are known for this purpose, forcing many researchers to use a simulation-based approach, which is highly unsatisfactory. To address this issue, we develop a theoretical and computational framework for determining the long-term behavior and stability properties for stochastic reaction networks. Our approach is based on a mixture of ideas from probability theory, linear algebra and optimization theory. We illustrate the broad applicability of our results by considering examples from various biological areas. The biological implications of our results are discussed as well.
ecology and environmental sciences
genetic networks
synthetic biology
metabolic networks
signaling networks
systems ecology
mathematics
algebra
network analysis
genome analysis
stochastic processes
computer and information sciences
regulatory networks
probability theory
systems biology
biochemistry
biochemical simulations
linear algebra
ecology
gene regulatory networks
genetics
biology and life sciences
physical sciences
genomics
computational biology
A Scalable Computational Framework for Establishing Long-Term Behavior of Stochastic Reaction Networks