Model reduction for stochastic chemical systems with abundant species

Smith, Stephen; Cianci, Claudia; Grima, Ramon

doi:10.1063/1.4936394

Biochemical processes typically involve many chemical species, some in abundance and some in low molecule numbers. We first identify the rate constant limits under which the concentrations of a given set of species will tend to infinity (the abundant species) while the concentrations of all other species remains constant (the non-abundant species). Subsequently, we prove that, in this limit, the fluctuations in the molecule numbers of non-abundant species are accurately described by a hybrid stochastic description consisting of a chemical master equation coupled to deterministic rate equations. This is a reduced description when compared to the conventional chemical master equation which describes the fluctuations in both abundant and non-abundant species. We show that the reduced master equation can be solved exactly for a number of biochemical networks involving gene expression and enzyme catalysis, whose conventional chemical master equation description is analytically impenetrable. We use the linear noise approximation to obtain approximate expressions for the difference between the variance of fluctuations in the non-abundant species as predicted by the hybrid approach and by the conventional chemical master equation. Furthermore, we show that surprisingly, irrespective of any separation in the mean molecule numbers of various species, the conventional and hybrid master equations exactly agree for a class of chemical systems.

I. INTRODUCTION

The chemical master equation (CME) is the accepted stochastic description of chemical reaction systems.¹ Since intrinsic noise roughly scales as the inverse square root of the mean number of molecules,¹ it follows that the CME provides a more accurate description than deterministic rate equations (REs), when species are in low concentrations. However, exact solution of the CME has proved to be impossible for all but the simplest systems (see for example Refs. 2–5), and Monte Carlo simulations using the stochastic simulation algorithm (SSA)⁶ are also time-consuming in many cases of interest. One way to bypass these issues is to use a hybrid model which treats some parts of the system using the SSA and the rest using a simulation method which is computationally more efficient. A common example of such hybrid modelling utilises time scale separation whereby some reactions occur on a fast time scale and are modelled using continuous approaches such as REs or chemical Langevin equations, while the rest of reactions occur on slow time scales and are modelled using the standard SSA.^7–9 Other methods which enable a considerable improvement over the SSA when time scale separation is present are: the nested SSA,¹⁰ a coarse-grained equation-free approach,¹¹ the constrained multiscale algorithm,^12,13 an approach based on finite-state projection,^14,15 the slow-scale SSA¹⁶ and the slow-scale linear-noise approximation.^17,18

In this paper, we consider a different type of hybrid model, one which uses a separation in the abundance of species (abundance separation) rather than time scale separation. In particular, we no longer split reactions into fast and slow but rather categorise species based on how abundant they are. These methods utilise a continuous approach to model the abundant species and a discrete approach to model the less abundant species. While less popular than time scale separation, some hybrid algorithms have been developed to take advantage of this idea (see, for example, Refs. 19–22). Stochastic simulations verify that these hybrid models can capture important features of the fully stochastic model. In particular, the model by Hellander and Lötstedt¹⁹ has been shown by Jahnke²¹ to be exact for monomolecular systems, i.e., the marginal distributions of non-abundant species in the hybrid model are exactly the same as the same obtained from the full stochastic model. More sophisticated (and computationally expensive) approaches have been postulated^21,22 to deal with systems which are not well described by the Hellander and Lötstedt hybrid model.

The advantages of methods using abundance separation over time scale separation are that: (i) the time scales of reactions are often unknown while the abundances are readily measurable; (ii) there is evidence suggesting that abundance separation is at least as significant, if not more, than time scale separation for biochemical networks inside cells. For example it has been shown that the mean number of proteins per E. Coli cell is roughly a thousand times that of the mean number of mRNA molecules per cell²³—in contrast, the ratio of protein to mRNA lifetime in E. Coli is expected to be above 1 but in the single digits.³ For mammalian cells, the same has been found; the median number of protein per cell is roughly 3000 times that of median number of mRNA per cell—in contrast, the ratio of median protein to median mRNA lifetime is about 5.²⁴ Clearly in these cases, abundance separation is significant while time scale separation is weak, and thus a method which takes advantage of the former appears to be ideal as a means to infer information about the stochastic dynamics of mRNA and of other proteins present in low copy numbers.

In this paper, we postulate a novel simple hybrid model based on abundance separation consisting of a CME for the non-abundant species coupled to REs for all species. Subsequently, we identify the rate constant limits under which the concentrations of a given set of species will tend to infinity (the abundant species) while the concentrations of all other species remain constant (the non-abundant species). This limit we shall refer to as the abundance or abundant limit. We show that, in this limit, the marginal distributions of the non-abundant species given by the hybrid model converge to the same distributions given by the CME of the full system. This fact is particularly useful when the hybrid model can be solved analytically, which is the case in several examples that we study. We illustrate the accuracy of our hybrid model by comparing the exact stochastic simulations of its reduced CME with exact stochastic simulations of the full CME. We also show that there are several chemical systems for which our hybrid model is exact even without abundance separation. In Section IV we offer an error analysis which provides an easy means to estimate the error incurred by the use of the hybrid model when the ratios of abundant to non-abundant species concentrations are finitely large. We conclude with a summary and discussion in Section V.

II. A REDUCED CHEMICAL MASTER EQUATION

In this section, we first propose a reduced CME which constitutes our hybrid model, and subsequently rigorously prove that it converges to the CME of the full system (that describing all species) in the abundance limit.

A. A heuristic reduction of the CME

The CME for a system of N chemical species which interact through R reactions has the form

\frac{d}{d t} P (n, t) = Ω \sum_{j = 1}^{R} (\prod_{i = 1}^{N} E_{i}^{- S_{i j}} - 1) {\hat{f}}_{j} (Ω, n) P (n, t),

(1)

where Ω is the volume in which the reactions occur, $E_{i}^{x}$ is the step operator which replaces n_i with n_i + x, the entries of the state vector n = (n₁, …, n_N) are the number of molecules of each species, P(n, t) is the probability of the system being in state n at time t, and the stoichiometric matrix elements S_ij are given by the net change in the number of molecules of species i when the jth reaction occurs. The probability that a reaction j occurs in a time interval [t, t + dt) is given by $Ω \hat{f} (n, Ω) d t$ ⁠, where ${\hat{f}}_{j}$ is a function of elements of the state vector and reaction rates. The REs are defined by

\frac{d ϕ_{i}}{d t} = \sum_{j = 1}^{R} S_{i j} f_{j} (\vec{ϕ}),

(2)

where $\vec{ϕ} = (ϕ_{1}, \dots, ϕ_{N})$ are the concentrations of the N chemical species, and $f_{j} (\vec{ϕ}) = lim_{Ω \to \infty} {\hat{f}}_{j} (Ω, Ω \vec{ϕ})$ ⁠.

We wish to reduce the number of species in this CME from N to M with M < N. Without loss of generality, we will keep species 1 to M and remove species M + 1 to N, which we consider to be “abundant.” We will do this by first summing the CME over n_M+1, …, n_N to leave us with an equation for the time evolution of the exact marginal distribution P^∗(n′, t), where n′ = (n₁, …, n_M), then subsequently we use a mean-field assumption to obtain a time evolution equation for the approximate marginal distribution $\tilde{P} (n^{'}, t)$ ⁠.

We will be considering the different possible forms that ${\hat{f}}_{j}$ can take assuming elementary reactions, specifically up to bimolecular reactions (reactions involving more than three molecules are rare in a biological setting). We will first investigate what happens to the CME if we sum over, say, n_h, using the notation that n_−h is the state vector n without the hth entry, in other words, $n_{- h} = (n_{1}, \dots, n_{h - 1}, n_{h + 1}, \dots, n_{N})$ ⁠. In what follows, we will use $X = (X_{1}, \dots, X_{N})$ to refer to the vector of chemical species, and $Y = (Y_{1}, \dots, Y_{N})$ to refer to the vector of random variables which give the number of molecules of each species. The state vector $n = (n_{1}, \dots, n_{N})$ therefore refers to a particular realisation of the random vector Y.

If reaction j does not feature X_h amongst its reactants, then ${\hat{f}}_{j}$ has no n_h dependence and the corresponding term in the CME remains unchanged.

If reaction j is a unimolecular reaction of the type X_h → … then ${\hat{f}}_{j} (Ω, n) = k_{j} n_{h} Ω^{- 1}$ ⁠, and we will have

\sum_{n_{h} = 0}^{\infty} k_{j} n_{h} Ω^{- 1} P (n, t) = k_{j} Ω^{- 1} 〈 Y_{h} | Y_{1} = n_{1}, Y_{2} = n_{2}, \dots 〉 P (n_{- h}, t) .

(3)

If reaction j is a bimolecular reaction of the type X_h + X_h → … then ${\hat{f}}_{j} (Ω, n) = k_{j} n_{h} (n_{h} - 1) Ω^{- 2}$ and we will have

\sum_{n_{h} = 0}^{\infty} k_{j} n_{h} (n_{h} - 1) Ω^{- 2} P (n, t) = k_{j} Ω^{- 2} 〈 Y_{h} (Y_{h} - 1) | Y_{1} = n_{1}, Y_{2} = n_{2}, \dots 〉 P (n_{- h}, t) .

(4)

Finally, if reaction j is a bimolecular reaction of the type X_h + X_i → … then ${\hat{f}}_{j} (Ω, n) = k_{j} n_{h} n_{i} Ω^{- 2}$ and we will have

\sum_{n_{h} = 0}^{\infty} k_{j} n_{h} n_{i} Ω^{- 2} P (n, t) = k_{j} Ω^{- 2} n_{i} 〈 Y_{h} | Y_{1} = n_{1}, Y_{2} = n_{2}, \dots 〉 P (n_{- h}, t) .

(5)

These results follow from the definition of conditional expectation,

〈 f (X) | Y = y 〉 = \sum_{x} f (x) \frac{P (X = x, Y = y)}{P (Y = y)} .

(6)

Given the above results, we can now see what will happen to each ${\hat{f}}_{j}$ if we sum over all h = M + 1, …, N. The result of this summation leads to new propensities involving conditional expectations which we call $f_{j}^{⋆}$ ⁠. It then follows that the exact marginal CME is given by

\frac{d}{d t} P^{⋆} (n^{'}, t) = Ω \sum_{j = 1}^{R} (\prod_{i = 1}^{M} E_{i}^{- S_{i j}} - 1) f_{j}^{⋆} (Ω, n) P^{⋆} (n^{'}, t),

(7)

where n′ = (n₁, …, n_M). In theory, we have simplified the CME while keeping it exact, but we should be careful because the dependence of the conditional expectations on n′ is currently unknown.

We proceed by making the heuristic mean-field assumption that these conditional expectations can be replaced by deterministic concentrations ϕ_i as defined earlier in Eq. (2), for example,

〈 Y_{h} | Y_{1} = n_{1}, \dots, Y_{M} = n_{M} 〉 \approx Ω ϕ_{h},

(8)

where we have approximated away all conditional dependence. We can correspondingly update the exact effective propensities $f_{j}^{⋆}$ with the approximate effective propensities ${\tilde{f}}_{j}$ ⁠. A general recipe for converting ${\hat{f}}_{j}$ to ${\tilde{f}}_{j}$ is given in Table I.

TABLE I.

The recipe for converting standard propensities ${\hat{f}}_{j}$ to reduced propensities ${\tilde{f}}_{j}$ ⁠. The subscripts i, r ≤ M refer to non-abundant species, while h, v > M refer to abundant species.

Propensity ${\hat{f}}_{j} (Ω, n)$	Reduced propensity $\tilde{f_{j}} (Ω, n^{'}, \vec{ϕ} (t))$
k_j	k_j
k_jn_iΩ⁻¹	k_jn_iΩ⁻¹
k_jn_in_rΩ⁻²	k_jn_in_rΩ⁻²
k_jn_i(n_i − 1)Ω⁻²	k_jn_i(n_i − 1)Ω⁻²
k_jn_in_hΩ⁻²	k_jn_iΩ⁻¹ϕ_h(t)
k_jn_hΩ⁻¹	k_jϕ_h(t)
k_jn_h(n_h − 1)Ω⁻²	$k_{j} ϕ_{h}^{2} (t)$
k_jn_hn_vΩ⁻²	k_jϕ_h(t) ϕ_v(t)

The approximate marginal CME is now

\frac{d}{d t} \tilde{P} (n^{'}, t) = Ω \sum_{j = 1}^{R} (\prod_{i = 1}^{M} E_{i}^{- S_{i j}} - 1) {\tilde{f}}_{j} (Ω, n^{'}, \vec{ϕ} (t)) \tilde{P} (n^{'}, t) .

(9)

In the rest of the paper, we will refer to Eq. (1) as the full CME and Eq. (9) as the reduced CME.

An alternative way of summarising our reduction method is that it consists of approximating a general chemical system,

s_{1 j} X_{1} + \dots + s_{N j} X_{N} \overset{k_{j}}{⟶} r_{1 j} X_{1} + \dots + r_{N j} X_{N}, j = 1, \dots, R,

(10)

by the reduced chemical system,

s_{1 j} X_{1} + \dots + s_{M j} X_{M} \overset{k_{j} ϕ_{M + 1}^{s_{M + 1, j}} \dots ϕ_{N}^{s_{N, j}}}{⟶} r_{1 j} X_{1} + \dots + r_{M j} X_{M}, j = 1, \dots, R,

(11)

when species X_M+1, …, X_N are abundant.

B. The abundant limit

We wish to show that the approximate marginal CME given in Eq. (9) is a good approximation to the exact marginal CME given in Eq. (7) when species X_M+1, …, X_N are abundant. To do this, we will need to define an abundance limit. Precisely, we want to know which parameters we should tweak in order that some species concentrations should go to infinity, while others stay constant. We will assume, without loss of generality, that we want to take the abundant limit of species X_N. For systems with multiple abundant species, we can just repeat the below process for each one in turn.

The convention we use for numbering reactions is given in Table II. We will introduce the functions a(i) and b(i) so that we can say that the bimolecular reaction with rate k_i has species X_a(i) and X_b(i) as its reactants, where a(i)≠b(i). We will have N input reactions with rate k_i, i = − 1, …, − N which lead to the production of each species, monomolecular reactions with rates k_i, i = 1, …, N, bimolecular reactions between different species with rates k_i, i = N + 1, …, L (for some L ∈ ℕ), and bimolecular reactions between the same species with rates k_i, i = L + 1, …, L + N.

TABLE II.

A convention for numbering reactions. Reactions −1, …, − N are input reactions, reactions 1, …, N are unimolecular, reactions N + 1, …, L are bimolecular and between different species and reactions L + 1, …, L + N are bimolecular and between the same species.

Reaction index	Reaction type	Reaction rate
i = − 1, …, − N	0̸ → …	k_i
i = 1, …, N	X_i → …	k_iϕ_i
i = N + 1, …, L	X_a(i) + X_b(i) → …	k_iϕ_a(i)ϕ_b(i)
i = L + 1, …, L + N	2X_i−L → …	$k_{i} ϕ_{i - L}^{2}$

Reaction index	Reaction type	Reaction rate
i = − 1, …, − N	0̸ → …	k_i
i = 1, …, N	X_i → …	k_iϕ_i
i = N + 1, …, L	X_a(i) + X_b(i) → …	k_iϕ_a(i)ϕ_b(i)
i = L + 1, …, L + N	2X_i−L → …	$k_{i} ϕ_{i - L}^{2}$

Now the rate equation for the concentration of X_r is

\frac{d}{d t} ϕ_{r} = \sum_{i = - N}^{- 1} S_{r i} k_{i} + \sum_{i = 1}^{N} S_{r i} k_{i} ϕ_{i} + \sum_{i = N + 1}^{L} k_{i} S_{r i} ϕ_{a (i)} ϕ_{b (i)} + \sum_{i = L + 1}^{L + N} S_{r i} k_{i} ϕ_{i - L}^{2},

(12)

where S_ri is the net change in the number of molecules of species X_r when reaction i occurs.

An intuitive means to obtain an abundant species X_N is to make the rate constants of the reactions which remove this species, to be very small, whilst the rest of the rate constants remain at their constant value. In particular, if exactly one molecule of X_N is a reactant, then we let $k_{j} \propto \frac{1}{x}$ ⁠; if two molecules of X_N are reactants, then we let $k_{j} \propto \frac{1}{x^{2}}$ ⁠, where x → ∞. This means that $k_{L + N} \propto \frac{1}{x^{2}}$ ⁠, $k_{N} \propto \frac{1}{x}$ ⁠, $k_{j} \propto \frac{1}{x}$ for j such that a(j) or b(j) equal N and j = N + 1, …, L. In what follows we shall refer to these rate constant limits as the abundance or abundant limit. Plugging in the aforementioned rates constant scalings and the trial solution,

(13)

ϕ_{i} = c_{i}, i \neq N,

ϕ_{N} = c_{N} x,

in Eq. , where c_i are constants independent of x, and considering steady-state by setting the time derivative to zero, one obtains a set of N simultaneous equations in the N constants c_i (i = 1, …, N). Importantly, the coefficients of these simultaneous quadratic equations are independent of x which implies that if these equations can be solved for c_i then the solution is independent of x, as they should indeed be, given the form of the trial solution above. Thus it follows that provided the simultaneous equations can be solved, the steady-state solution of the REs is given by Eq. . Note this implies that in the abundance limit, the ratio of the abundant to the non-abundant species concentrations, ϕ_N/ϕ_i (i≠N), scales as x where x → ∞, whilst the concentration of the non-abundant species remains constant.

Note also that we have here implicitly assumed that there is no chemical conservation law which involves an abundant species (chemical conservation laws which involve only non-abundant species are however allowed). This is since this law necessitates a finite upper bound on the concentrations which is contrary to the manner in which we have here defined the abundant limit.

For systems with multiple abundant species, the above recipe implies that $k_{j} \propto \frac{1}{x^{q}}$ where q is the total number of reactant molecules of abundant species involved in reaction j. For example, for the reaction X_h + X_v → … where both species X_h and X_v are abundant, the rate constant of the reaction scales as $\frac{1}{x^{2}}$ ⁠.

The limits here derived assume a steady-state rate equation description for all species. This derivation is here presented to simplify the presentation and since it is very intuitive. However , as we show in Sec. II C 3, the limits elucidated here, also constitute abundance limits for the time-dependent stochastic description.

C. Proof of N species abundant convergence

We will now use the limits derived in Section II B to prove that the approximate marginal distribution $\tilde{P} (n^{'}, t)$ governed by heuristic marginal CME (9) converges to the exact marginal distribution P^⋆(n′, t) governed by exact marginal CME (7) in the abundance limit.

1. Taylor expansion of exact marginal distribution

A full N species CME with R reactions has the form of Eq. (1). We will expand the solution P(n, t) as a Taylor series in time about t = 0. We assume deterministic initial conditions, $P (n, t = 0) = δ_{n_{1}}^{n_{1}^{0}} \dots δ_{n_{N}}^{n_{N}^{0}}$ ⁠, where $n_{i}^{0}$ denotes the initial value of n_i. We can write the Taylor expansion as

P (n, t) = \sum_{k = 0}^{\infty} P^{(k)} (n, 0) \frac{t^{k}}{k!},

(14)

where P^(k) is the kth time derivative of P. Since the full CME is a coupled set of first-order ordinary differential equations with constant coefficients, the Taylor series above is guaranteed to have an infinite radius of convergence by Fuchs’s theorem.²⁵

From this series, we can compute the marginal distribution,

P^{⋆} (n^{'}, t) = \sum_{n_{M + 1} = 0}^{\infty} \dots \sum_{n_{N} = 0}^{\infty} P (n, t) = \sum_{n_{M + 1} = 0}^{\infty} \dots \sum_{n_{N} = 0}^{\infty} \sum_{k = 0}^{\infty} P^{(k)} (n, 0) \frac{t^{k}}{k!} = \sum_{k = 0}^{\infty} P^{⋆ (k)} (n^{'}, 0) \frac{t^{k}}{k!} .

(15)

We already know $P (n, 0) = δ_{n_{1}}^{n_{1}^{0}} \dots δ_{n_{N}}^{n_{N}^{0}}$ ⁠, so our first problem is the second term of the expansion, which is the first time derivative. This is given by CME (1) itself,

\dot{P} (n, 0) = Ω \sum_{j = 1}^{R} (E_{1}^{- S_{1 j}} \dots E_{N}^{- S_{N j}} - 1) {\hat{f}}_{j} (Ω, n) P (n, 0) = Ω \sum_{j = 1}^{R} (E_{1}^{- S_{1 j}} \dots E_{N}^{- S_{N j}} - 1) {\hat{f}}_{j} (Ω, n) δ_{n_{1}}^{n_{1}^{0}} \dots δ_{n_{N}}^{n_{N}^{0}},

(16)

and thus we can compute the kth derivative,

P^{(k)} (n, 0) = Ω^{k} {(\sum_{j = 1}^{R} [(E_{1}^{- S_{1 j}} \dots E_{N}^{- S_{N j}} - 1) {\hat{f}}_{j} (Ω, n)])}^{k} δ_{n_{1}}^{n_{1}^{0}} \dots δ_{n_{N}}^{n_{N}^{0}},

(17)

where we are careful to note that $P^{(k)} (n, 0) ≢ {(\dot{P} (n, 0))}^{k}$ ⁠, since the E operators do not commute with the propensities ${\hat{f}}_{j}$ ⁠. Now, we can compute the marginal distribution P^⋆, which is made much simpler by the presence of the Kronecker-deltas,

P^{⋆ (k)} (n^{'}, 0) = \sum_{n_{M + 1} = 0}^{\infty} \dots \sum_{n_{N} = 0}^{\infty} P^{(k)} (n, 0) = Ω^{k} {(\sum_{j = 1}^{R} [(E_{1}^{- S_{1 j}} \dots E_{M}^{- S_{M j}} E_{M + 1^{'}}^{- S_{M + 1, j}} \dots E_{N^{'}}^{- S_{N j}} - 1) {\hat{f}}_{j} (Ω, n^{'}, n_{M + 1}^{0}, \dots, n_{N}^{0})])}^{k} δ_{n_{1}}^{n_{1}^{0}} \dots δ_{n_{M}}^{n_{M}^{0}},

(18)

where $E_{i^{'}}^{x}$ now acts on the initial conditions $n_{i}^{0}$ rather than the variable n_i.

2. Taylor expansion of the approximate reduced distribution

The approximate marginal distribution $\tilde{P} (n^{'}, t)$ is defined by the reduced CME,

\dot{\tilde{P}} (n^{'}, t) = Ω \sum_{j = 1}^{R} (E_{1}^{- S_{1 j}} \dots E_{M}^{- S_{M j}} - 1) {\tilde{f}}_{j} (Ω, n^{'}, \vec{ϕ} (t)) \tilde{P} (n^{'}, t) .

(19)

Its Taylor expansion about t = 0 is

\tilde{P} (n^{'}, t) = \sum_{k = 0}^{\infty} {\tilde{P}}^{(k)} (n^{'}, 0) \frac{t^{k}}{k!} .

Fuchs’s theorem guarantees that a first-order ordinary differential equation with time-dependent coefficients will have a radius of convergence at least as large as the minimum of the radius of convergence of the time-dependent parameters.²⁵ The reduced CME is a set of coupled first-order equations with time-dependent coefficients determined by the solution of the REs. Hence, if the REs admit a Taylor series solution with an infinite radius of convergence then the Taylor series of the reduced CME also does.

We already know $\tilde{P} (n^{'}, 0) = δ_{n_{1}}^{n_{1}^{0}} \dots δ_{n_{M}}^{n_{M}^{0}}$ ⁠, so our first problem is the second term of the expansion. Again, this is given by approximate CME (19),

\dot{\tilde{P}} (n^{'}, 0) = Ω \sum_{j = 1}^{R} (E_{1}^{- S_{1 j}} \dots E_{M}^{- S_{M j}} - 1) {\tilde{f}}_{j} (Ω, n^{'}, \vec{ϕ} (0)) \tilde{P} (n^{'}, 0) = Ω \sum_{j = 1}^{R} (E_{1}^{- S_{1 j}} \dots E_{M}^{- S_{M j}} - 1) {\tilde{f}}_{j} (Ω, n^{'}, \vec{ϕ} (0)) δ_{n_{1}}^{n_{1}^{0}} \dots δ_{n_{M}}^{n_{M}^{0}},

(20)

where we note that the propensities ${\tilde{f}}_{j}$ will in general depend on t as well as n′. This will cause complications, for instance, the second derivative has the form:

{\tilde{P}}^{(2)} (n^{'}, 0) = Ω \sum_{j_{1} = 1}^{R} (E_{1}^{- S_{1 j}} \dots E_{M}^{- S_{M j}} - 1) [{\dot{\tilde{f}}}_{j} (Ω, n^{'}, \vec{ϕ} (0)) \tilde{P} (n^{'}, 0) + {\tilde{f}}_{j} (Ω, n^{'}, \vec{ϕ} (0)) {\tilde{P}}^{(1)} (n^{'}, 0)] = Ω \sum_{j = 1}^{R} (E_{1}^{- S_{1 j}} \dots E_{M}^{- S_{M j}} - 1) {\dot{\tilde{f}}}_{j} (Ω, n^{'}, \vec{ϕ} (0)) δ_{n_{1}}^{n_{1}^{0}} \dots δ_{n_{M}}^{n_{M}^{0}} + Ω^{2} \sum_{j_{1} = 1}^{\infty} \sum_{j_{2} = 1}^{\infty} [(E_{1}^{- S_{1 j_{1}}} \dots E_{M}^{- S_{M j_{1}}} - 1) {\tilde{f}}_{j_{1}} (Ω, n^{'}, \vec{ϕ} (0)) (E_{1}^{- S_{1 j_{2}}} \dots E_{M}^{- S_{M j_{2}}} - 1) {\tilde{f}}_{j_{2}} (Ω, n^{'}, \vec{ϕ} (0)) δ_{n_{1}}^{n_{1}^{0}} \dots δ_{n_{M}}^{n_{M}^{0}}] .

(21)

We now have an extra term in our sum which depends on the time derivative of the ${\tilde{f}}_{j}$ ⁠, and if we take higher order Taylor coefficients, we get higher time derivatives of the ${\tilde{f}}_{j}$ ⁠. We get, using the notation used in Subsection II C 1,

{\tilde{P}}^{(k)} (n^{'}, 0) = Ω^{k} {(\sum_{j = 1}^{R} [(E_{1}^{- S_{1 j}} \dots E_{M}^{- S_{M j}} - 1) {\tilde{f}}_{j} (Ω, n^{'}, \vec{ϕ} (0))])}^{k} δ_{n_{1}}^{n_{1}^{0}} \dots δ_{n_{M}}^{n_{M}^{0}} + terms proportional to the time derivatives of the {\tilde{f}}_{j} .

(22)

3. Convergence of full and reduced Taylor series

The absolute difference between the kth terms of the two Taylor CMEs is given by the difference between Eqs. (18) and (22),

Ω^{k} [{(\sum_{j = 1}^{R} [(E_{1}^{- S_{1 j}} \dots E_{M}^{- S_{M j}} - 1) {\tilde{f}}_{j} (Ω, n^{'}, \vec{ϕ} (0))])}^{k} - {(\sum_{j = 1}^{R} [(E_{1}^{- S_{1 j}} \dots E_{M}^{- S_{M j}} E_{M + 1^{'}}^{- S_{M + 1, j}} \dots E_{N^{'}}^{- S_{N j}} - 1) {\hat{f}}_{j} (Ω, n^{'}, n_{M + 1}^{0}, \dots, n_{N}^{0})])}^{k}] δ_{n_{1}}^{n_{1}^{0}} \dots δ_{n_{M}}^{n_{M}^{0}} + terms proportional to the time derivatives of the {\tilde{f}}_{j} .

(23)

This will tend to zero in the abundant limit if we can prove two claims,

For any α_i ∈ ℤ, and j ∈ 1, …, R,

$[{\tilde{f}}_{j} (Ω, n_{1} + α_{1}, \dots, n_{M} + α_{M}, \vec{ϕ} (0)) - {\hat{f}}_{j} (Ω, n_{1} + α_{1}, \dots, n_{N}^{0} + α_{N})] \to 0,$
(24)
in the abundant limit.
The time derivatives of the ${\tilde{f}}_{j}$ tend to zero in the abundant limit.

To prove (I), we must recall the different possible forms that ${\hat{f}}_{j}$ and ${\tilde{f}}_{j}$ can take, which are given in Table III. As before, we say that the indices i, r ≤ M correspond to non-abundant species, while the indices h, v > M correspond to abundant species.

TABLE III.

The different forms that ${\tilde{f}}_{j}$ and ${\hat{f}}_{j}$ can take in our Taylor-expanded distributions. The α_i correspond to integer changes in species number introduced by the $E_{i}^{x}$ operators. Note that indices i, r ≤ M correspond to non-abundant species, while the indices h, v > M correspond to abundant species.

Reaction type	${\tilde{f}}_{j}$	${\hat{f}}_{j}$
0̸ → …	k_j	k_j
X_i → …	k_j(n_i + α_i)Ω⁻¹	k_j(n_i + α_i)Ω⁻¹
X_h → …	k_jϕ_h(0)	$k_{j} (n_{h}^{0} + α_{h}) Ω^{- 1}$
X_i + X_r → …	k_j(n_i + α_i)(n_r + α_r)Ω⁻²	k_j(n_i + α_i)(n_r + α_r)Ω⁻²
X_i + X_i → …	k_j(n_i + α_i)(n_i + α_i − 1)Ω⁻²	k_j(n_i + α_i)(n_i + α_i − 1)Ω⁻²
X_h + X_v → …	k_jϕ_h(0) ϕ_v(0)	$k_{j} (n_{h}^{0} + α_{h}) (n_{v}^{0} + α_{v}) Ω^{- 2}$
X_h + X_h → …	$k_{j} ϕ_{h}^{2} (0)$	$k_{j} (n_{h}^{0} + α_{h}) (n_{h}^{0} + α_{h} - 1) Ω^{- 2}$
X_i + X_h → …	k_j(n_i + α_i) ϕ_h(0)Ω⁻¹	$k_{j} (n_{i} + α_{i}) (n_{h}^{0} + α_{h}) Ω^{- 2}$

Convergence is trivial for reactions 0̸ → …, X_i → …, X_i + X_r → …, and X_i + X_i → …, since ${\tilde{f}}_{j}$ and ${\hat{f}}_{j}$ agree. For the other reactions, we have to decide how the initial conditions should scale in the abundant limit. We will suppose that the initial concentration for species X_M + 1, …, X_N should tend to infinity O(x) (since they are the abundant species), while the initial concentration for species 1, …, M should stay constant.

For the X_h → … reaction, we are interested in

k_{j} ϕ_{h} (0) - k_{j} (n_{h}^{0} + α_{h}) Ω^{- 1} .

(25)

However by definition, $ϕ_{h} (0) = \frac{n_{h}^{0}}{Ω}$ ⁠, so the absolute error becomes $\frac{k_{j} α_{h}}{Ω}$ which tends to zero since k_j → 0 in accordance with our abundant limit elucidated in Section II B.

For the X_h + X_v → … reaction, we are interested in

k_{j} ϕ_{h} (0) ϕ_{v} (0) - k_{j} (n_{h} + α_{h}) (n_{v} + α_{v}) Ω^{- 2},

(26)

which simplifies to

k_{j} (\frac{n_{h}^{0} α_{v}}{Ω^{2}} + \frac{n_{v}^{0} α_{h}}{Ω^{2}} + \frac{α_{h} α_{v}}{Ω^{2}}),

(27)

and which again tends to zero, since in the abundant limit $k_{j} \propto \frac{1}{x^{2}}$ while $n_{h}^{0}$ and $n_{v}^{0}$ are proportional to x and x → ∞.

For the X_h + X_h → … reaction, we are interested in

k_{j} ϕ_{h} {(0)}^{2} - k_{j} (n_{h}^{0} + α_{h}) (n_{h}^{0} + α_{h} - 1) Ω^{- 2},

(28)

which simplifies to

k_{j} (\frac{2 n_{h}^{0} α_{h}}{Ω^{2}} - \frac{n_{h}^{0}}{Ω^{2}} + \frac{α_{h}^{2} - α_{h}}{Ω^{2}}),

(29)

and which again tends to zero since $k_{j} \propto \frac{1}{x^{2}}$ while $n_{h}^{0} \propto x \to \infty$ ⁠.

Finally, for the X_i + X_h → … reaction we consider

k_{j} (n_{i} + α_{i}) ϕ_{h} (0) Ω^{- 1} - k_{j} (n_{i} + α_{i}) (n_{h}^{0} + α_{h}) Ω^{- 2},

(30)

which we write as

α_{h} k_{j} (n_{i} + α_{i}) Ω^{- 2},

(31)

and which also tends to zero since k_j → 0 in the abundance limit. We have therefore proved claim (I).

To prove claim (II), we will use the convention for reactions introduced in Table II. We will say that k₀ corresponds to the rate of a null reaction, k₁, …, k_N correspond to the rates of monomolecular reactions, k_N+1, …, k_L correspond to the rates of bimolecular reactions involving distinct species with reaction k_j involving species X_a(j) and X_b(j), and k_L+1, …, k_L+N correspond to homodimerisation reactions.

We will prove claim (II) for the reaction X_h + X_v → 0, with reduced propensity ${\tilde{f}}_{j} (\vec{ϕ} (0)) = k_{j} ϕ_{i} (0) ϕ_{r} (0)$ ⁠. Note that generally we have ${\tilde{f}}_{j} (Ω, n^{'}, \vec{ϕ} (0))$ but because the reaction involves two abundant species there is no explicit dependence on Ω and n′ and hence in what follows we use the less cumbersome notation ${\tilde{f}}_{j} (\vec{ϕ} (0))$ ⁠. Using Leibniz’s rule,

{\tilde{f}}_{j}^{(k)} (\vec{ϕ} (0)) = k_{j} \sum_{l = 0}^{k} (\binom{k}{l}) ϕ_{h}^{(k - l)} (0) ϕ_{v}^{(l)} (0) .

(32)

For this reaction, $k_{j} \propto \frac{1}{x^{2}}$ and also ϕ_h(0), ϕ_v(0) ∝ x → ∞. If we can prove that each derivative $ϕ_{i}^{(k)} (0)$ is bounded in the abundant limit, k > 0, then the above expression for ${\tilde{f}}_{j}^{(k)}$ will tend to zero in that limit, owing to the limiting prefactor k_j. First, the definition of the first derivative of ϕ_h is given by the rate equation

\frac{d ϕ_{h} (0)}{d t} = \sum_{j = - 1}^{- N} S_{h j} k_{j} + \sum_{j = 1}^{N} S_{h j} k_{j} ϕ_{j} (0) + \sum_{j = N + 1}^{L} k_{j} S_{h j} ϕ_{a (j)} (0) ϕ_{b (j)} (0) + \sum_{j = L + 1}^{L + N} k_{j} S_{h j} ϕ_{j - L} {(0)}^{2} .

(33)

Wherever we have a concentration which tends to infinity, (ϕ_j(0), j = M + 1, …, N) it is by definition cancelled out by the corresponding parameter k_j. So this expression remains constant, and therefore bounded, in the abundant limit. For the purposes of mathematical induction, assume now that all derivatives up to $ϕ_{h}^{(k)} (0)$ remain bounded in that limit. Then, using Leibniz’s rule,

ϕ_{h}^{(k + 1)} (0) = \sum_{j = 1}^{N} S_{h j} k_{j} ϕ_{j}^{(k)} (0) + \sum_{j = N + 1}^{L} S_{h j} k_{j} \sum_{q = 0}^{k} (\binom{k}{q}) ϕ_{a (j)}^{(k - q)} (0) ϕ_{b (j)}^{(q)} (0) + \sum_{j = L + 1}^{L + N} S_{h j} k_{j} \sum_{q = 0}^{k} (\binom{k}{q}) ϕ_{j - L}^{(k - q)} (0) ϕ_{j - L}^{(q)} (0) .

(34)

The only terms here that have a possibility of being unbounded are those involving zeroth derivatives of the ϕ_j, specifically,

\sum_{j = N + 1}^{L} S_{h j} k_{j} (ϕ_{a (j)}^{(k)} (0) ϕ_{b (j)} (0) + ϕ_{a (j)} (0) ϕ_{b (j)}^{(k)} (0)) + 2 \sum_{j = L + 1}^{L + N} S_{h j} k_{j} ϕ_{j - L}^{(k)} (0) ϕ_{j - L} (0) .

(35)

But these ϕ_j(0) will only go to infinity if j > M, and in that case the reaction rate k_j will go to zero at least fast enough to counteract the limiting concentrations. So therefore we have shown that each time derivative of the concentrations is bounded in the abundant limit, and therefore each time derivative of the ${\tilde{f}}_{j}$ goes to zero in that limit. The proofs for the other reactions are very similar. Therefore we have proved claim (II), and consequently we have proved the convergence of the full marginal CME and the reduced CME for all times, in the abundant limit.

III. THE ACCURACY OF THE REDUCTION FOR FINITELY LARGE ABUNDANT CONCENTRATIONS AND SPECIAL CASES

The abundant limit, as stated previously, is the limit that the rate constants of the reactions removing the abundant species go to zero (in a particular manner) which ensures that the ratios of the abundant to non-abundant concentrations go to infinity. It is in this limit that we have proved that the difference between the reduced and full marginal CME goes to zero. Generally, we are interested in the case where the ratio of the abundant to non-abundant concentrations is finitely large, not infinite. In this case, the reduced CME is approximate. Given two identical copies of a chemical system, one placed in a small volume and the other in a much larger volume, and given they have the same finite large ratio of abundant to non-abundant concentrations, we expect the reduced CME to be a better approximation for the system confined in the larger volume. The reason is that the number of molecules of abundant and non-abundant species is larger in the system confined in the larger volume and hence the REs in this case provide a good approximation to the dynamics of the abundant species,^26,27 the key principle behind our reduced CME. Another equivalent point of view is that in the larger volume the reactions occur on faster time scales than the reactions in low volumes, due to the larger number of interacting molecules and hence the dynamics of the abundant species in the larger volume are more amenable to being modelled by a continuous approach like the chemical Langevin equation or REs.²⁸

Of course, as Gillespie pointed out in Ref. 28, it is difficult to tell how large should the compartment volume be so that a macroscopic approach becomes a feasible approximation. Hence, in this section, we explore the accuracy of the reduced CME, for finite ratios of abundant to non-abundant concentrations, by means of stochastic simulations. In particular we will use the SSA⁶ to sample the probability distribution of full CME (1) and the Extrande algorithm²⁹ to sample the probability distribution of reduced CME (9) for the case where the rate constants of the reactions removing the abundant species scale as $k_{j} \propto \frac{1}{x^{q}}$ (where q is the total number of reactant molecules of abundant species involved in reaction j) and x is a finite real number. Note that the reduced CME cannot be sampled using the SSA since the propensities are generally time-dependent and hence the use of Extrande. An alternative to the use of the latter algorithm would be to use a method involving the numerical integration of reaction propensities.³⁰ We also show in this section that curiously, for some chemical systems, the exact and approximate marginals are identical even without taking abundance limit.

A. Stochastic simulations

1. Homodimerisation

We will investigate an open homodimerisation reaction studied in Ref. 26,

0̸ ⇌_{k_{3}}^{k_{1}} X_{1}, X_{1} + X_{1} \overset{k_{2}}{⟶} X_{2}, X_{2} \overset{k_{4}}{⟶} 0̸ .

(36)

Species X₁ is produced with rate k₁, and is either consumed with rate k₃ or else forms a dimer X₂ with rate k₂. The dimer X₂ is then consumed with rate k₄. We consider the case where X₁ is abundant and X₂ is not. The RE for the concentration of X₁ is

{\dot{ϕ}}_{1} = k_{1} - 2 k_{2} ϕ_{1}^{2} - k_{3} ϕ_{1},

(37)

which has the solution

ϕ_{1} (t) = - \frac{α tanh (\frac{(β - t) α}{2}) + k_{3}}{4 k_{2}},

(38)

where

α = \sqrt{8 k_{1} k_{2} + k_{3}^{2}}, β = \frac{2}{α} arctanh (\frac{4 k_{2} ϕ_{1} (0) + k_{3}}{- α}) .

(39)

Therefore, approximate reduced CME (9) for non-abundant species X₂ is given by

\dot{\tilde{P}} (n_{2}, t) = Ω k_{2} ϕ_{1}^{2} (t) (\tilde{P} (n_{2} - 1, t) - \tilde{P} (n_{2}, t)) + k_{4} ((n_{2} + 1) \tilde{P} (n_{2} + 1, t) - n_{2} \tilde{P} (n_{2}, t)) .

(40)

The approximate equilibrium distribution of this system is therefore Poisson $(\frac{Ω k_{2} ϕ_{1}^{2}}{k_{4}})$ ⁠. In Figs. 1 and 2, we compare the exact marginal distribution of X₂ from the full CME and the approximate marginal distributions of the reduced CME. These are obtained by means of an ensemble average of SSA and Extrande trajectories, respectively. In Fig. 1, we plot the distributions for four different relative abundances of X₁ and X₂ at the same time point t = 0.01. The abundance is adjusted by choosing the rate constants to scale as $k_{2} \propto \frac{1}{x^{2}}$ and $k_{3} \propto \frac{1}{x}$ (in accordance with the limits delineated in Section II B) and varying x over a certain finite range (see caption of Fig. 1). The distance between the two distributions clearly becomes smaller as the ratio of the abundant to non-abundant species concentrations increases, in line with the proof in Sec. II; the two distributions are practically indistinguishable when this ratio is of the order of 100. Nevertheless, we find that some salient features of the two distributions are fairly similar (namely, the position of mode and the width of distribution) over a large range of the ratio of the abundant to non-abundant species concentrations (0.4–263). In Fig. 2, we show that the high accuracy of the reduced CME also extends to predicting the whole time-evolution of the distribution. Both of these figures indicate that the results of simulations using the reduced CME bear a significant closeness to those obtained using the full CME under a wide range of abundances and hence point towards the utility of the reduced CME as a low dimensional approximation of the full CME.

FIG. 1.

View large Download slide

Exact marginal distribution of species X₂ in open homodimerisation reaction (36) (red line) compared with the approximate marginal distribution (blue circles) at time t = 0.01, and for different relative abundances $\frac{ϕ_{1} (t)}{ϕ_{2} (t)}$ ⁠. Parameters are k₁ = 10³, $k_{2} = \frac{100}{x^{2}}$ ⁠, $k_{3} = \frac{10}{x}$ ⁠, k₄ = 10, ϕ₁(0) = 4x, ϕ₂(0) = 0, Ω = 1. The parameter x equals 1, 10, 100, 1000 in (a)-(d), respectively.

FIG. 2.

View large Download slide

Time development of the distribution for species X₂ in open homodimerisation reaction (36) when abundance is high (ϕ₁(0) = 4000, ϕ₂(0) = 0). There is excellent agreement between exact and approximate marginals at all times. Parameter values are k₁ = 10³, k₂ = 10⁻⁴, k₃ = 10⁻², k₄ = 10, Ω = 1.

2. Genetic feedback loop

We next investigate a negative genetic feedback loop studied in Ref. 2,

D_{on} \overset{v_{0}}{⟶} D_{on} + P, P \overset{d_{0}}{⟶} 0̸, D_{off} ⇌_{d_{1}}^{v_{1}} D_{on} + P .

(41)

The “on” promoter D_on produces proteins P with rate v₀. These proteins bind to the promoter with rate d₁ to generate the inactive “off” promoter D_off, which can then unbind back into active promoter and protein with rate v₁. Furthermore, the protein P is consumed by a unimolecular reaction with rate d₀. We consider the case where P is abundant.

The REs for this system are given by

(42)

{\dot{ϕ}}_{1} = - v_{1} ϕ_{1} + d_{1} (\frac{1}{Ω} - ϕ_{1}) ϕ_{2},

{\dot{ϕ}}_{2} = v_{1} ϕ_{1} - d_{1} (\frac{1}{Ω} - ϕ_{1}) ϕ_{2} + v_{0} (\frac{1}{Ω} - ϕ_{1}) - d_{0} ϕ_{2},

where ϕ₁ is the concentration of D_off, ϕ₂ is the concentration of P, Ω is the volume, and $\frac{1}{Ω}$ is the total (fixed) gene concentration (equivalent to one gene). Unlike the previous example, these equations do not admit an analytical solution, so we must solve them numerically, and then use the numerical solution in the reduced CME,

\dot{\tilde{P}} (n_{1}, t) = v_{1} ((n_{1} + 1) \tilde{P} (n_{1} + 1, t) - n \tilde{P} (n_{1}, t)) + d_{1} ϕ_{2} (t) ((2 - n_{1}) \tilde{P} (n_{1} - 1, t) - (1 - n_{1}) \tilde{P} (n_{1}, t)),

(43)

where n₁ is the number of molecules of D_off. Since at any one time, the gene is either on or off, the distribution of D_off is Bernoulli. In particular, since the parameter of the Bernoulli distribution, θ(t), is equal to $\tilde{P} (1, t)$ ⁠, we can say

D_{off} \sim Bernoulli (θ (t)),

(44)

where $\dot{θ} (t)$ is obtained by setting n₁ = 1 in the above reduced CME,

\dot{θ} = - v_{1} θ (t) + d_{1} ϕ_{2} (t) (1 - θ (t)) .

(45)

Our reduction method therefore provides us with an “exact” solution at all times in this case, since we do not need to perform any stochastic simulations sampling the reduced CME, but rather just numerically solve the ordinary differential equation above.

Figure 3 shows how the approximate expression for θ given in Eq. (45) compares with the true θ (obtained by computing 〈n₁〉 from an ensemble average of SSA trajectories of the full CME) for different relative abundances at t = 50. The relative abundance is controlled by choosing the rate constants to scale as $d_{0}, d_{1} \propto \frac{1}{x}$ (in accordance with the limits delineated in Section II B) and varying x over a certain finite range (see caption of Fig. 3). For the parameter set chosen, we find that the approximation for θ using the reduced CME is in good qualitative agreement with that calculated from the full CME when the ratio of abundant to non-abundant concentrations varies over the range 10³–10⁷. In particular both the full and reduced CME predict that the probability of the gene being in the off state increases monotonically, in a step-like manner, as the protein concentration increases at constant gene concentration (consistently with a negative feedback loop). For low relative abundances (ratios less than a few hundreds), the approximate θ is almost double the true value implying that the reduced CME in this case over-estimates the strength of the negative feedback.

FIG. 3.

View large Download slide

The parameter θ = 〈n₁〉 for genetic feedback loop (41). The exact value of the parameter is obtained from an ensemble average of SSA trajectories of the full CME (blue line), while the approximate value is that obtained from our reduced CME via numerical solution of Eq. (45) (red circles). The time is fixed to t = 50. Parameter values are $Ω = 1, v_{0} = 1, v_{1} = 0.1, d_{0} = \frac{1}{x}, d_{1} = \frac{0.01}{x}, 1 < x < 1 0^{4}$ ⁠.

3. Metabolic network

We consider an arbitrarily large, sequential enzyme reaction network which has been previously associated with metabolism.^31,32 The network consists of N + 1 enzymes, each converting a substrate into a product which then serves as the substrate for the next enzyme in the cascade. The first substrate is produced by a zeroth-order reaction, and the final substrate is converted into a product which we ignore. We seek the approximate distribution of the enzymes, which we expect to be exact in the limit where substrates are abundant. The full chemical system is described by the scheme,

(46)

0̸ \overset{k_{in}}{⟶} S_{0},

E_{0} + S_{0} ⇌_{k_{- 1}^{0}}^{k_{1}^{0}} C_{0} \overset{k_{2}^{0}}{⟶} E_{0} + S_{1},

E_{1} + S_{1} ⇌_{k_{- 1}^{1}}^{k_{1}^{1}} C_{1} \overset{k_{2}^{1}}{⟶} E_{1} + S_{2},

⋮

E_{N} + S_{N} ⇌_{k_{- 1}^{N}}^{k_{1}^{N}} C_{N} \overset{k_{2}^{N}}{⟶} E_{N} .

The conservation law n_{E_i} + n_{C_i} = D_i is implied for each enzyme, where D_i is a given positive integer which represents the total number of both free and bound enzymes of type i, which is constant in time. According to our method, encapsulated by Eq. , the reduced chemical system takes the simpler form as follows:

E_{i} ⇌_{k_{- 1}^{i} + k_{2}^{i}}^{k_{1}^{i} ϕ_{S_{i} (t)}} C_{i}, i = 0, \dots, N .

(47)

It is hence clear that the molecule numbers of each enzyme species are binomially distributed in steady-state conditions,

\tilde{P} (n_{E_{i}}) = (\binom{D_{i}}{n_{E_{i}}}) {(\frac{k_{- 1}^{i} + k_{2}^{i}}{k_{1}^{i} ϕ_{S_{i}} + k_{- 1}^{i} + k_{2}^{i}})}^{n_{E_{i}}} {(\frac{k_{1}^{i} ϕ_{S_{i}}}{k_{1}^{i} ϕ_{S_{i}} + k_{- 1}^{i} + k_{2}^{i}})}^{D_{i} - n_{E_{i}}},

(48)

where ϕ_{S_i} is the steady-state solution of the REs of the full system, Eq. (46), given by

ϕ_{S_{i}} = \frac{k_{in} (k_{2}^{i} + k_{- 1}^{i})}{k_{1}^{i} (\frac{k_{2}^{i} D_{i}}{Ω} - k_{in})} .

(49)

For a time-dependent description, the reduced CME corresponding to reduced chemical system (47) cannot be exactly solved and stochastic simulations are required. In Fig. 4, we plot both the approximate and exact distributions (using Extrande for the reduced system and the SSA for the full system) of the enzyme E₁ at a fixed time for different abundances of substrate S₁. It is clear that the approximation improves as the substrate becomes more abundant than enzyme and is essentially exact in the bottom right panel where the relative abundance is $Ω \frac{ϕ_{S_{1}}}{D_{1}} = 10$ ⁠. It is also remarkable that the approximation is good even when there is essentially no clear separation in abundance, i.e., $\frac{Ω ϕ_{S_{1}}}{D_{1}} = 1$ ⁠. Indeed even though the approximation suffers quantitatively when the relative abundance is not high, yet it captures the main distinctive qualitative feature, namely, that the distribution changes from positive to negative skewness as a function of the relative abundance (the switch happens at a relative abundance between 1 and 10).

FIG. 4.

View large Download slide

Probability distribution of the number of molecules of species E₁ in the metabolic network model (46), for different substrate to enzyme abundance ratios: 0.01, 0.1, 1, and 10, respectively. This system consists of 11 distinct enzyme and 11 substrate species, though here we only look at the distribution of E₁ relative to the abundance of S₁, i.e., Ωϕ_S₁(t)/D₁. There is clear convergence between the approximate and exact distributions as the substrate becomes more abundant than the enzyme. Parameters are Ω = 1, t = 100, N = 10, D_i = 10, k_in = 18, $k_{- 1}^{i} = 1$ ⁠, $k_{2}^{i} = 2$ ⁠, i = 0, …, N.

For this system, we have the added benefit that the distribution of the number of molecules of each enzyme E_i is independent in the approximate description. This means that if we are interested in the distribution of the number of molecules of a given enzyme, say, E₁, then we only need to simulate the three reactions involving that particular enzyme, rather than the 3N + 4 reactions of the full system. There is therefore a marked reduction in computational time for our reduced SSA, particularly for large N, as shown in Fig. 5, where the approximate SSA is roughly 3 times faster than the exact SSA when N = 20. We note, however, that the computational time of the approximate method does increase slightly with N, owing to the need to solve a coupled system of 2N + 2 rate equations.

FIG. 5.

View large Download slide

Computational time taken to compute an individual trajectory of length 100 time units for species E₁ in the metabolic network model with a total of 2N + 2 species, using the Exrande algorithm for reduced (approximate) chemical system (47) and the SSA for full (exact) chemical system (46). Parameters are Ω = 1, D_i = 10, k_in = 18, $k_{1}^{i} = 0.02$ ⁠, $k_{- 1}^{i} = 1$ ⁠, $k_{2}^{i} = 2 \forall i$ ⁠. Simulations were performed using MATLAB on a computer with a 3 GHz Intel Core 2 Duo processor and 4 GB RAM. Note that as the total number of species increases, stochastic simulations of the reduced system (approximate) become significantly more computationally efficient than the SSA for the full system (exact).

4. Genetic oscillator with transcriptional feedback

We consider an arbitrarily large gene regulatory network which has been previously studied as a model of a circadian oscillator.^33,34 The mechanism is as follows. A protein P₁ is translated by mRNA, M, which is itself transcripted by a gene in the on-state, D_on. Subsequently, the protein P₁ generates P₂, and P₂ generates P₃, etc., until a final protein P_N is generated. The latter can bind to D_on to deactivate it as D_off, which can reversibly unbind into P_N and D_on. We seek the approximate distribution of the number of molecules of D_on and M, which we expect to be accurate when the proteins are abundant. The full chemical system is described by the scheme,

(50)

D_{on} \overset{v_{0}}{⟶} D_{on} + M, M \overset{d_{0}}{⟶} 0̸, M \overset{k_{1}}{⟶} M + P_{1},

P_{1} \overset{k_{2}}{⟶} P_{2} \overset{k_{3}}{⟶} \dots \overset{k_{N}}{⟶} P_{N} \overset{k_{N + 1}}{⟶} 0̸, D_{on} + P_{N} ⇌_{v_{1}}^{d_{1}} D_{off} .

According to our method, encapsulated by Eq. , the reduced chemical system takes the simpler form as follows:

D_{on} \overset{v_{0}}{⟶} D_{on} + M, M \overset{d_{0}}{⟶} 0̸, D_{on} ⇌_{v_{1}}^{d_{1} ϕ_{P_{N}} (t)} D_{off} .

(51)

We note that the distribution of D_on is independent of M, and is therefore simply $Bernoulli (\frac{v_{1}}{d_{1} ϕ_{P_{N}} (t) + v_{1}})$ ⁠. The steady-state distribution of M can be straightforwardly obtained using the method in Ref. 2 or else a fast implementation of the finite-state projection algorithm is equally effective.¹⁵ For a time-dependent description, however, we must use stochastic simulations to determine the accuracy of our method. In Fig. 6, we plot the time development of the distribution of the number of mRNA molecules, M, for parameters such that the steady state is characterised by a fixed large abundance of proteins, in particular, $\frac{ϕ_{P_{i}}}{ϕ_{M}} = 100, \forall i$ ⁠, which is a physically realistic ratio for some cells.²³ Remarkably the approximate distribution provides an excellent match to the exact distribution for all times, reproducing even the transition from unimodality to bimodality and back to unimodality as a function of time.

FIG. 6.

View large Download slide

Distribution of the number of mRNA, M, in genetic oscillator model (50), for parameter values which give the steady-state abundance $\frac{ϕ_{P_{N}}}{ϕ_{M}} = 100$ and at different times: 0.01, 0.06, 0.11, and 1, respectively. The exact and approximate distributions agree for all times. Remarkably, the approximate system also shows the bimodal behaviour characteristic of the full system, as shown by the bimodal distribution at t = 0.11. Parameters are N = 10, Ω = 1, v₀ = 100, d₀ = 1, v₁ = 1, d₁ = 0.01, k₁ = 1, k₂ = ⋯ = k_N+1 = 0.01.

In Fig. 7, we plot the computational time taken to simulate an individual trajectory of length 10 time units with the SSA for full system (50) and Extrande for approximate system (51). For the approximate system, only 4 reactions must be simulated, while the full system has N + 5 reactions. On the other hand, the approximate system requires the time-dependent solution of an N + 2 dimensional system of REs. This trade-off implies that for N = 1, 2, the full system is slightly faster, but for any N > 3, the approximate system is faster.

FIG. 7.

View large Download slide

Computational time taken to compute an individual trajectory of length 10 time units for the mRNA species in the genetic oscillator model using the Extrande algorithm for reduced (approximate) chemical system (51) and the SSA for full (exact) chemical system (50). The computational time is plotted as a function of the number of distinct protein species (N). Parameters are Ω = 1, v₀ = 100, d₀ = 1, v₁ = 1, d₁ = 0.01, k₁ = 1, k₂ = ⋯ = k_N+1 = 0.01. Simulations were performed using MATLAB on a computer with a 3 GHz Intel Core 2 Duo processor and 4 GB RAM. Note that as the total number of protein species increases, stochastic simulations of the reduced system (approximate) becomes significantly more computationally efficient than the SSA for the full system (exact).

B. Exact reductions

We have already shown that our approximation is exact in the limit of infinite concentration of the abundant species, but as we now show, surprisingly, there is also a wide class of systems where the method is exact, regardless of abundance separation.

1. Systems in detailed balance

We now show that for systems in detailed balance,¹ where the species we remove from the system (what we would previously call “abundant” species) are not involved in a chemical conservation law, the approximation is exact. Note that detailed balance is a property of some systems in steady-state conditions, and hence necessarily, the exactness mentioned does not apply to finite times, rather it applies only in the limit of infinitely long times.

Consider a detailed balance system of R reversible reactions, and N chemical species, X₁, …, X_N. Let us denote reaction j as

s_{1 j} X_{1} + \dots + s_{N j} X_{N} ⇌_{k_{j}^{'}}^{k_{j}} r_{1 j} X_{1} + \dots + r_{N j} X_{N},

(52)

where we allow a conservation law on the species X₁, …, X_M of the form

\sum_{i = 1}^{M} α_{i} n_{i} = k,

(53)

where α_i and k are time-independent constants. Application of the method described by van Kampen in Ref. 35 leads to an explicit expression for the steady-state solution of the CME of this chemical system which is given by a constrained multivariate Poisson distribution of the form

P (n) = C \frac{{(Ω ϕ_{1})}^{n_{1}} \dots {(Ω ϕ_{N})}^{n_{N}}}{n_{1}! \dots n_{N}!} δ (α_{1} n_{1} + \dots + α_{M} n_{M}, k),

(54)

where ϕ_i are the steady-state rate equation solutions, δ(…, …) is a Kronecker delta, and C is a normalisation constant. The marginal distribution is obtained by summing over n_M+1, …, n_N. The exact marginal is therefore,

P^{⋆} (n^{'}) = \sum_{n_{M + 1} = 0}^{\infty} \dots \sum_{n_{N} = 0}^{\infty} P (n) = C^{'} \frac{{(Ω ϕ_{1})}^{n_{1}} \dots {(Ω ϕ_{M})}^{n_{M}}}{n_{1}! \dots n_{M}!} δ (α_{1} n_{1} + \dots + α_{M} n_{M}, k),

(55)

where C′ is a normalisation constant.

Now the approximate reduction introduced in Section II is equivalent to approximating chemical system (52) by reduced chemical system,

s_{1 j} X_{1} + \dots + s_{M j} X_{M} ⇌_{k_{j}^{'} ϕ_{M + 1}^{r_{M + 1, j}} \dots ϕ_{N}^{r_{N, j}}}^{k_{j} ϕ_{M + 1}^{s_{M + 1, j}} \dots ϕ_{N}^{s_{N, j}}} r_{1 j} X_{1} + \dots + r_{M j} X_{M},

(56)

with the same conservation law as above. Application of the method in Ref. 35 to the reduced CME describing the above system immediately leads to a steady-state solution which is exactly the same as Eq. (55) since the RE solution of the reduced system is the same as that of the full system. Hence the approximation is exact for this class of chemical systems in detailed balance.

2. Open Michaelis-Menten reaction with one enzyme molecule

In the following example, we show that the approximation can be exact in steady-state conditions without taking any abundant limits, even if the system is not in detailed balance.

The open Michaelis-Menten reaction is given by

0̸ \overset{k_{in}}{⟶} S, S + E ⇌_{k_{1}^{'}}^{k_{1}} C \overset{k_{2}}{⟶} E + P,

(57)

where substrate molecules S are input into the system, they reversibly bind with enzyme E to form a complex C which in turn irreversibly decays into the original enzyme E and product molecules P.

We will consider the case with the conservation law n_E + n_C = 1, that is where there is just one enzyme molecule in the compartment. The CME describing the above reaction system is

\dot{P} (n_{S}, n_{E}) = k_{in} Ω (P (n_{S} - 1, n_{E}) - P (n_{S}, n_{E})) + \frac{k_{1}}{Ω} ((n_{S} + 1) (n_{E} + 1) P (n_{S} + 1, n_{E} + 1) - n_{S} n_{E} P (n_{S}, n_{E})) + k_{1}^{'} ((1 - n_{E} + 1) P (n_{S} - 1, n_{E} - 1) - (1 - n_{E}) P (n_{S}, n_{E})) + k_{2} ((1 - n_{E} + 1) P (n_{S}, n_{E} - 1) - (1 - n_{E}) P (n_{S}, n_{E})) .

(58)

This equation has been solved exactly in steady-state conditions in Appendix G of Schnoerr et al.³⁶ In particular therein it was shown that the average enzyme molecule number in steady-state condition is given by $〈 n_{E} 〉 = 1 - \frac{k_{in} Ω}{k_{2}}$ ⁠. This together with the fact that a single enzyme molecule, at any given time, can be in only one of two states, implies that the steady-state marginal distribution of enzyme number fluctuations is

P^{⋆} (n_{E}) \sim Bernoulli (1 - \frac{k_{in} Ω}{k_{2}}) .

(59)

Next, we show that our reduction gives exactly the same distribution, regardless of the abundance of substrate. The reduced CME describing enzyme fluctuations is given by

\dot{\tilde{P}} (n_{E}) = k_{1} ϕ_{S} ((n_{E} + 1) \tilde{P} (n_{E} + 1) - n \tilde{P} (n_{E})) + (k_{1}^{'} + k_{2}) ((2 - n_{E}) \tilde{P} (n_{E} - 1) - (1 - n) \tilde{P} (n_{E})) .

(60)

In steady state, setting n_E = 0 gives us

k_{1} ϕ_{S} \tilde{P} (1) - (k_{1}^{'} + k_{2}) \tilde{P} (0) = 0 .

(61)

Therefore, with the condition $\tilde{P} (0) + \tilde{P} (1) = 1$ ⁠, we find that

〈 n_{E} 〉 = \tilde{P} (1) = \frac{k_{1}^{'} + k_{2}}{k_{1} ϕ_{S} + k_{1}^{'} + k_{2}} .

(62)

The REs for this system are

(63)

{\dot{ϕ}}_{S} = k_{in} + k_{1}^{'} (\frac{1}{Ω} - ϕ_{E}) - k_{1} ϕ_{E} ϕ_{S},

{\dot{ϕ}}_{E} = - k_{1} ϕ_{E} ϕ_{S} + (k_{1}^{'} + k_{2}) (\frac{1}{Ω} - ϕ_{E}),

which possess a steady state solution

ϕ_{S} = \frac{k_{in}}{k_{1}} (\frac{k_{1}^{'} + k_{2}}{k_{2} - k_{in}}) .

(64)

Substituting Eq. (64) into Eq. (62), we find

〈 n_{E} 〉 = \tilde{P} (1) = 1 - \frac{k_{in} Ω}{k_{2}} = Ω ϕ_{E} .

(65)

As by arguments before, the steady-state distribution is Bernoulli and hence it follows that

\tilde{P} (n_{E}) \sim Bernoulli (1 - \frac{k_{in} Ω}{k_{2}}),

(66)

which is equal to exact solution Eq. (59).

By similar arguments, it can be easily deduced that the marginal distribution of any species which exists in two states and for which the average number of molecules predicted by the REs is the same as the CME, is exactly predicted by the reduced CME. The second criterion on average molecule numbers is bound to generally be the limiting one since it is typically not the case that the REs exactly agree with the mean concentrations calculated from the CME (see, for example, Ref. 26). For example for genetic feedback loop (41) the marginal distributions of the gene in the on or off state cannot be exactly predicted by the reduced CME because as shown in Ref. 2, the average number of genes in each state (equivalently the fraction of time spent in each state) predicted by the CME does not equal that of the REs.

IV. ESTIMATING THE APPROXIMATION ERROR OF THE HYBRID MODEL

As we have shown for most systems, the reduction is exact only in the limit of infinite concentrations of certain species, and the reduction is therefore an approximation if concentrations are finitely large.

We now investigate the use of the Linear Noise Approximation (LNA) to obtain an estimate of the error made by the use of the reduced CME. By comparing this estimate with that obtained from stochastic simulations of both the full and approximate systems, we demonstrate that the LNA’s estimate is accurate for a wide range of parameters and systems. Since the LNA is obtained by solving a system of coupled ordinary differential equations, our results suggest the use of the LNA as a computationally efficient means of estimating the error which bypasses lengthy stochastic simulations using stochastic simulations of the full and reduced CMEs.

The LNA is an approximation which assumes the fluctuations in each chemical species are normally distributed. More precisely, it is the leading order approximation of the system-size expansion of the CME¹ in the limit of large volumes. The general formulation of the LNA is as follows (see, for example, Ref. 37 for more details).

Consider a system of N chemical species, with R reactions, where the jth reaction is

s_{1 j} X_{1} + \dots + s_{N j} X_{N} \overset{k_{j}}{⟶} r_{1 j} X_{1} + \dots + r_{N j} X_{N} .

(67)

The REs for the system are then given by

\dot{\vec{ϕ}} = S \vec{f},

(68)

where we remind the reader that S is the stochiometric matrix with elements S_ij = r_ij − s_ij and $\vec{f}$ is the macroscopic propensity vector $\vec{f}$ defined as

f_{j} = k_{j} ϕ_{1}^{s_{1 j}} \dots ϕ_{N}^{s_{N j}} .

(69)

The Jacobian matrix J is the derivative of Eq. (68) with respect to $\vec{ϕ}$ ⁠,

J = S \frac{\partial \vec{f}}{\partial \vec{ϕ}} .

(70)

The diffusion matrix D is given by the matrix product,

D = S diag (\vec{f}) S^{T} .

(71)

The time-evolution of the second moments of the fluctuations is then approximately given by the Lyapunov differential equation,

\frac{d C}{d t} = J C + C J^{T} + D,

(72)

where ΩC_ij is the LNA estimate for the CME’s prediction of the covariance in the number fluctuations of species X_i and X_j.

Now, proposed reduction approximates reaction scheme (67) by

s_{1 j} X_{1} + \dots + s_{M j} X_{M} \overset{k_{j} ϕ_{M + 1}^{s_{M + 1, j}} \dots ϕ_{N}^{s_{N j}}}{⟶} r_{1 j} X_{1} + \dots + r_{M j} X_{M}, j = 1, \dots, R,

(73)

when species X_M+1, …, X_N are the abundant species. Note that, for the reduced system with M species, the abundant concentrations no longer function as concentrations, but instead as parameters like the reaction rates k_j. The REs remain unchanged, however. The Jacobian and diffusion matrices for this system, $\tilde{J}$ and $\tilde{D}$ ⁠, are hence the upper left M × M blocks of J and D previously defined in Eqs. (70) and (71), respectively. Thus, the LNA leads to an estimate of the reduced CME’s prediction of the covariance of fluctuations, $Ω \tilde{C}$ ⁠, which is the solution of the Lyapunov equation,

\frac{d \tilde{C}}{d t} = \tilde{J} \tilde{C} + \tilde{C} {\tilde{J}}^{T} + \tilde{D} .

(74)

Hence, it follows that the LNA’s estimate of the absolute relative difference in the variance predictions of the full and reduced CME’s for species i is given by

R_{i} = \frac{| C_{i i} - {\tilde{C}}_{i i} |}{C_{i i}} .

(75)

Of course, one can also calculate this quantity as a function of time, by solving the Lyapunov equations of the full and reduced CMEs numerically; however, in what follows we shall assume steady-state conditions to simplify the presentation.

Though its generally impossible to obtain a closed-form simple analytical solution to the LNA equations, one can show that the error R_i is approximately proportional to the inverse of the ratio of the abundant to non-abundant species concentrations in the abundant limit. The proof is as follows.

Referring to Table II, if species X_N is abundant, then the abundant limit consists of the rate constant scalings: $k_{L + N} \propto \frac{1}{x^{2}}$ ⁠, $k_{N} \propto \frac{1}{x}$ ⁠, $k_{j} \propto \frac{1}{x}$ for j such that a(j) or b(j) equal N (j denotes a bimolecular reaction which involves X_N and another species) and the steady-state concentration scalings: ϕ_i = c_i for i≠N, where c_i are constants independent of x and ϕ_N ∝ x. It is easy to verify using this limit and the REs given by Eq. (12) that the Jacobian matrix can be written as J = J₀ + yJ₁, where y = 1/x and J_i are matrices to be determined from the REs. In particular, J₀ is the Jacobian of the REs with the terms describing the removal of the abundant species set to zero. The same scaling form for the Jacobian is obtained for any number of abundant species.

On the other hand, the diffusion matrix D is unchanged under the abundance limit. This is since by Eq. (71), the elements of the D are linear functions of the macroscopic rate functions f_j (see Eq. (69)) which are unchanged by the abundance limit since each limit of k_j tending to zero will be counterbalanced by the opposite limit of a concentration of an abundant species tending to infinity.

Hence, in the abundance limit, Lyapunov Eq. (72) can be written as

\frac{d C}{d t} = (J_{0} + J_{1} y) C + C (J_{0}^{T} + J_{1}^{T} y) + D .

(76)

The form of this equation suggests a solution of the type C = C⁽⁰⁾ + C⁽¹⁾y + C⁽²⁾y² + O(y³). Indeed plugging this ansatz in the above equation, one transforms it into a coupled set of equations for the matrices C⁽ⁱ⁾ which can be solved iteratively, i.e., the equation for C⁽ⁱ⁾ depends on C^(j) where j < i except for C⁽⁰⁾ which is a function of J₀ and D only. Now the abundance limit is the limit y = 1/x → 0 and hence the relative error in the variance can be written as

R_{i} = \frac{| C_{i i} - C_{i i}^{(0)} |}{C_{i i}} = \frac{| C_{i i}^{(1)} |}{C_{i i}^{(0)}} y + O (y^{2}) .

(77)

Now in the abundance limit, the ratios of abundant to non-abundant concentrations are proportional to x = 1/y and hence it follows that in this limit, the relative error R_i is proportional to the inverse of these ratios.

Next we demonstrate the accuracy of the LNA’s estimate of the relative error in the predictions of the reduced CME, i.e., Eq. (75). This is done by comparison of the LNA estimate with the relative error directly computed from the SSA of the full CME and the steady-state analytical solution of the reduced CME for three examples of biochemical relevance.

A. Open Michaelis-Menten reaction with multiple enzyme molecules

We consider open Michaelis Menten system (57) with multiple number of enzyme molecules, i.e., (n_E + n_C) = E_T, where E_T is the total number of enzyme molecules. We consider the case where the substrate is much more abundant than the enzyme. Computing the LNA of the full CME and of the reduced CME for the non-abundant enzyme species, we find that the relative error in the variance of the enzyme number fluctuations, as given by Eq. (75), is

R = \frac{(1 - a) \frac{E_{T}}{Ω} + K_{M}}{{(1 - a)}^{2} \frac{E_{T}}{Ω} + K_{M}} - 1,

(78)

where $a = \frac{k_{in} Ω}{k_{2} E_{T}}$ and $K_{M} = \frac{k_{1}^{'} + k_{2}}{k_{1}}$ is the Michaelis-Menten constant.

The steady-state substrate concentration solution of the REs for this system is

ϕ_{S} = \frac{K_{M} a}{(1 - a)} .

(79)

We define $L = \frac{Ω ϕ_{S}}{E_{T}}$ as a measure of the relative abundance of substrate S. It then follows that R can be written as

R = \frac{a^{2}}{a (1 - a) + L} .

(80)

The condition 0 < a < 1 is a requirement for the existence of a steady state, and R is a monotonically increasing function of a, so the maximum possible value of R is at a = 1, in other words,

R_{\max} = \frac{1}{L} = \frac{E_{T}}{Ω ϕ_{S}} .

(81)

That is, if the substrate concentration is ten times the total enzyme concentration, then the percentage relative error in the reduced CME’s estimate of the variance of enzyme number fluctuations will be less than ten percent.

The reduced CME can in this case be exactly solved in steady-state conditions and one obtains a binomial distribution with parameters E_T and 1 − a describing the fluctuations in enzyme molecule numbers; indeed for the case E_T = 1, the binomial distribution reduces to the Bernoulli distribution found earlier for the open Michaelis Menten system with one enzyme molecule (see Eq. (59)). In Fig. 8, we use the variance calculated from this solution together with the variance calculated from time-averages of SSA (for the full CME) to compute the true error in the reduced CME’s variance of enzyme number fluctuations for the open Michaelis-Menten system. This is done for two different volumes, Ω = 1 and Ω = 10³. The true error is also compared in the same figure with the LNA estimate given by Eq. (80). The relative concentrations of substrate and enzyme are controlled by setting the rate constant k₁ proportional to 1/x and varying x (in accordance with the abundance limits discussed earlier; see caption of Fig. 4 for details). The LNA estimates are reasonably good for both volumes but practically indistinguishable from the true error for the larger volume of Ω = 10³. This is to be expected since the LNA becomes exact in the limit of large volumes.

FIG. 8.

View large Download slide

Relative error in the reduced CME’s prediction of the variance of fluctuations of non-abundant species for open Michaelis-Menten reaction (57) in steady-state conditions. The true relative error in the variance of non-abundant enzyme number fluctuations at Ω = 1 (red circles) and Ω = 10³ (green crosses) is computed using a time-average of SSA (for the full CME) and the analytical steady-state distribution solution of the reduced CME and compared with the LNA estimate given by Eq. (80) (blue line). For this example, the LNA gives a reasonable estimate of the true error at Ω = 1 and an excellent estimate at Ω = 10³. L is the relative abundance defined as the substrate concentration divided by the total enzyme concentration (according to the REs). Parameter values are $k_{in} = 10, k_{1} = \frac{1}{x}, k_{1}^{'} = 1 0^{2}, k_{2} = 1, E_{T} = 20, 10 < x < 1 0^{4}$ ⁠.

B. Open homodimerisation reaction

Next, we use the LNA to estimate the errors in the reduced CME description for homodimerisation example (36). We consider the case in which species X₁ is abundant compared to species X₂. Choosing the scalings k₂ = c₁/x² and k₃ = c₂/x (where c_i are proportionality constants), it follows by the considerations of Section II B that the steady-state concentration of X₁ is proportional to x while that of X₂ is a constant; hence, by varying x we have a convenient way to control the ratio of the two concentrations. In particular, in steady-state, the solution of the REs is given by

(82)

ϕ_{1} = \frac{\sqrt{c_{2}^{2} + 8 c_{1} k_{1}} - c_{2}}{4 c_{1}} x,

ϕ_{2} = \frac{c_{2}^{2} + 4 c_{1} k_{1} - c_{2} \sqrt{c_{2}^{2} + 8 c_{1} k_{1}}}{8 c_{1} k_{4}} .

The LNA relative error in the variance of the fluctuations of the non-abundant species X₂ as given by Eq. is

R = \frac{λ_{0}}{λ_{1} + x λ_{2}},

(83)

where $λ = (c_{2} + \sqrt{c_{2}^{2} + 8 c_{1} k_{1}}) k_{4} / 2 c_{1} k_{1}$ ⁠, $λ_{0} = 4 c_{1}^{2} ϕ_{2} λ^{2} (k_{1} - ϕ_{2} λ (c_{2} + 4 c_{1} ϕ_{2} λ))$ ⁠, $λ_{1} = c_{2}^{2} k_{4} + 8 c_{1} c_{2} ϕ_{2} k_{4} λ + c_{1} c_{2}^{2} ϕ_{2} λ^{2} + 4 c_{1}^{2} ϕ_{2} (k_{1} + 4 ϕ_{2} k_{4}) λ^{2} + 4 c_{1}^{2} c_{2} ϕ_{2}^{2} λ^{3}$ are constants. Note that the ratio of the abundant to the non-abundant species concentrations is proportional to x. Hence, by Eq. (83), it follows that the relative error has a maximum equal to λ₀/λ₁ and decreases monotonically as X₁ becomes more abundant relative to X₂. Next we test the accuracy of the LNA estimate.

The reduced CME for this system can be exactly solved in steady-state conditions and one obtains a Poisson distribution for the fluctuations in the number of molecules of X₂ with parameter $Ω k_{2} ϕ_{1}^{2} / k_{4}$ ⁠. In Fig. 9, we use the variance calculated from this solution together with the variance calculated from time-averages of SSA (for the full CME) to compute the true error in the reduced CME’s variance of number of X₂ fluctuations. This is done for two different volumes, Ω = 1 and Ω = 10³. The true error is also compared in the same figure with the LNA estimate given by Eq. (83). As for the previous example, the LNA accuracy is good across a wide range of volumes and becomes particularly accurate in the limit of large volumes. It is also noteworthy that the LNA estimate of the relative error is good over a wide range of relative abundances; in particular it even provides an accurate value (about 0.3) for the maximum relative error which occurs in the limit of small relative abundance of X₁ compared to X₂.

FIG. 9.

View large Download slide

Relative error in the reduced CME’s estimate of the variance of fluctuations of non-abundant species for open homodimerisation reaction (36) in steady-state conditions. The true relative error in the variance of the fluctuations of the non-abundant species X₂ at Ω = 1 (red circles) and Ω = 10³ (green crosses) is computed using the time-averages of SSA for the full CME and the analytical steady-state distribution of the reduced CME and compared with the LNA estimate given by Eq. (83) (blue line). For this example, the LNA gives a good estimate of the true error for Ω = 1 and an excellent estimate for Ω = 10³. Parameter values are k₁ = 10³, k₄ = 10, $k_{2} = \frac{100}{x^{2}}$ ⁠, $k_{3} = \frac{10}{x}$ ⁠, where 1 < x < 10⁴.

C. Genetic feedback loop

Here, we use the LNA to estimate the error in the reduced CME description for genetic feedback loop (41), where the gene concentration is fixed to 1/Ω, i.e., a single gene. We consider the case where the protein P is much more abundant than the gene. Choosing the scalings d₀ = c₀/x and d₁ = c₁/x (where c_i are proportionality constants), it follows by the considerations of Section II B that the steady-state concentration of X₂ (the protein) is proportional to x while that of X₁ (the gene) is a constant; hence, by varying x we have a convenient way to control the ratio of the two concentrations. In particular, in steady-state, the solution of the REs is given by

ϕ_{1} = \frac{2 c_{1} v_{0} + c_{0} Ω v_{1} - \sqrt{(c_{0} Ω v_{1}) (4 c_{1} v_{0} + c_{0} Ω v_{1})}}{2 c_{1} Ω v_{0}},

(84)

ϕ_{2} = {(2 c_{1})}^{- 1} x (\frac{\sqrt{v_{1} (4 c_{1} v_{0} + c_{0} Ω v_{1})}}{\sqrt{c_{0} Ω}} - v_{1}) .

(85)

The LNA relative error in the variance of the fluctuations of the non-abundant gene as given by Eq. (75) is

R = \frac{λ_{0}}{λ_{1} + x λ_{2}},

(86)

where λ = ϕ₂/ϕ₁x, λ₀ = − c₁(−1 + ϕ₁Ω)(c₁(−1 + ϕ₁Ω)(−1 + 2ϕ₁Ω) v₀v₁ − c₀ϕ₁Ω²v₁(v₀ + v₁) + c₀c₁ϕ₁Ω((−1 + ϕ₁Ω) v₀ − ϕ₁v₁) λ), $λ_{1} = (v_{1} + c_{1} ϕ_{1} λ) (c_{0}^{2} ϕ_{1} Ω^{3} v_{1} + c_{0}^{2} c_{1} ϕ_{1} Ω^{2} λ - c_{0}^{2} c_{1} f 1^{2} Ω^{3} λ + c_{1}^{2} (- 1 + ϕ_{1} Ω) ((- 1 + ϕ_{1} Ω) v_{0} - c_{0} ϕ_{1} Ω λ) - c_{1}^{2} ϕ_{1} Ω (- 1 + ϕ_{1} Ω) ((- 1 + ϕ_{1} Ω) v_{0} - c_{0} ϕ_{1} Ω λ))$ and $λ_{2} = (v_{1} + c_{1} ϕ_{1} λ) (c_{1} ϕ_{1} Ω^{2} (1 - ϕ_{1} Ω) v_{0} v_{1} + c_{0} ϕ_{1} Ω^{3} v_{1}^{2} + c_{0} c_{1} ϕ_{1} Ω^{2} v_{1} λ + c_{1}^{2} ϕ_{1} Ω (- 1 + ϕ_{1} Ω) λ ((- 1 + ϕ_{1} Ω) v_{0} - c_{0} ϕ_{1} Ω λ))$ are constants. Note that the ratio of the abundant to the non-abundant species concentrations is proportional to x. Hence, by Eq. (86), it follows that the relative error has a maximum equal to λ₀/λ₁ and decreases monotonically as X₂ becomes more abundant relative to X₁. Note that the form of the LNA estimate of the error in this example is the same as that in the previous example.

In Fig. 10, we plot the true error in the variance of the fluctuations of the non-abundant gene computed using time-averages of SSA (for the full CME) and the analytical solution of the reduced CME in steady-state conditions (a Bernoulli distribution with parameter given by the steady-state solution of Eq. (45)) at Ω = 1. This is compared with LNA’s estimate Eq. (86) which is found to be particularly accurate, as found previously for the enzyme and dimerisation examples. However, unlike the previous examples, for the gene system, in the limit of large Ω, the LNA’s estimate does not become more accurate. The reason is that the LNA is accurate in the deterministic limit in which all species molecule numbers increase to infinity at constant concentration, whereas, in this example, the gene molecule number is fixed to one and only the protein molecule number is taken to infinity.

FIG. 10.

View large Download slide

Relative error in the reduced CME prediction of the variance of gene fluctuations for genetic feedback loop (41) in steady-state conditions. The true relative error at Ω = 1 (red circles) is computed using a time-average of SSA for the full CME and the analytical steady-state solution of the reduced CME and compared with the LNA estimate given by Eq. (86) (blue line). Parameter values are v₀ = 1, v₁ = 0.1, $d_{0} = \frac{1}{x}$ ⁠, $d_{1} = \frac{0.01}{x}$ ⁠, where 1 < x < 10⁴. The relative abundance of protein to gene concentrations is ϕ₂/ϕ₁ (according to the REs).

D. Genetic oscillator

Finally, we use the LNA to estimate the error in the reduced CME description for genetic oscillator (50). We consider the case where the proteins P_i are more abundant than the mRNA, while the gene is restricted to a maximum concentration of $\frac{1}{Ω}$ ⁠, i.e., a single molecule. Choosing the scalings d₁ = c₁/x and k₂ = ⋯ = k_N+1 = c₂/x (where c_i are proportionality constants), it follows by the considerations of Section II B that the steady-state concentrations of P_i (the proteins) is proportional to x while that of D_on (the gene) and M (the mRNA) are constant; hence, by varying x, we have a convenient way to control the ratio of the two concentrations. In particular, in steady-state, the solution of the REs is given by

ϕ_{D_{on}} = \frac{\sqrt{v_{1}^{2} + 4 v_{0} c_{1} k_{1} v_{1} / (d_{0} c_{2} Ω) - v_{1}}}{2 v_{0} c_{1} k_{1} / (d_{0} c_{2})},

(87)

ϕ_{M} = \frac{v_{0}}{d_{0}} ϕ_{D_{on}},

(88)

ϕ_{P_{i}} = \frac{k_{1} x}{c_{2}} ϕ_{M} .

(89)

Given the arbitrarily large number of species, there is no compact analytic expression for the LNA relative error in the variance of the fluctuations of the non-abundant mRNA as given by Eq. (75); however, the error can be calculated by numerical solution of the REs and the Lyapunov equations of the full and reduced systems.

In Fig. 11, we plot the true error in the variance of the fluctuations of the non-abundant mRNA computed using time-averages of SSA (for the full CME) and the solution of the reduced CME in steady-state conditions (computed with the finite-state projection algorithm) at Ω = 1. This is compared with the LNA estimate which is found to be accurate for physically realistic abundances (ratio of protein to mRNA concentrations are commonly larger than a hundred in bacteria; see for example Ref. 23). However, unlike some of the previous examples, for the gene system, in the limit of large Ω, the LNA’s estimate does not become more accurate. The reason is that the LNA is accurate in the deterministic limit in which all species molecule numbers increase to infinity at constant concentration, whereas, p in this example, the gene molecule number is fixed to one and only the protein molecule number is taken to infinity.

FIG. 11.

View large Download slide

Relative error in the reduced CME prediction of the variance of mRNA fluctuations for genetic oscillator (50) in steady-state conditions. The true relative error at Ω = 1 (red circles) is computed using a time-average of the SSA for the full CME and the steady-state solution of the reduced CME computed with the finite-state projection algorithm and compared with the LNA estimate (blue line). Parameter values are N = 10, v₀ = 100, v₁ = 1, d₀ = 1, $d_{1} = \frac{1}{x}$ ⁠, k₁ = 1, $k_{2} = \dots = k_{N + 1} = \frac{1}{x}$ where 1 < x < 10³. The relative abundance of protein to mRNA concentrations is ϕ_{P_i}/ϕ_M (according to the REs) which is the same for all i = 1, …, N.

V. SUMMARY AND DISCUSSION

Summarising, in this paper, we have introduced a novel reduced stochastic description of chemical systems in which some species are abundant. The key intuitive idea is to replace the conditional expectation of the number of molecules of abundant species in the propensities of the exact marginal CME by the solutions of the deterministic REs and hence obtain a reduced CME for the non-abundant species only. Therefore, our method is a hybrid approach. We note that our method is different and simpler than that presented in Refs. 19, 21, and 22 since the latter postulate a more complicated approach than REs to model the abundant species. This relative simplicity indeed leads to three major strengths of our approach over existing approaches: (i) it is easier to implement and computationally more efficient than present approaches; (ii) our reduced CME can be explicitly solved for a number of biochemically relevant examples; (iii) simple rational expressions can be derived which estimate the errors inherent in the hybrid approximation relative to the fully stochastic description. Curiously we also found that the reduced CME at the heart of our hybrid method is exact for some chemical systems, i.e., without requiring the necessity of abundance separation or without restricting the system to purely monomolecular systems (as was found to be the case in Ref. 21 to ensure exactness for the Hellander and Lotstedt model). The major disadvantage of our approach is that its unlikely that it will be able to capture as many features of the fully stochastic model as the more sophisticated approaches mentioned above.

The present work also suggests some new avenues of research. The first is finding exact error bounds for the reduced CME, which could provide useful reassurance for mathematical modellers using this method. A second interesting direction would be to develop a more refined reduction of the CME by replacing the conditional expectation of the number of molecules of abundant species in the propensities of the exact marginal CME by the solutions of effective mesoscopic rate equations (EMREs)²⁶ instead of REs. EMREs have been demonstrated to be more accurate than REs in the sense that the difference between their mean concentration solution and that of the CME is considerably smaller than the difference between the RE solution and that of the CME.²⁷ Hence, a CME reduction based on EMREs is highly likely to be more accurate than the one developed, in this paper, particularly for cases where the compartment volume is small such that even though there is a large ratio of abundant to non-abundant concentrations, the number of molecules of abundant species is small. Finally another interesting area for future work is the relationship between time scale separation and abundance separation. It is clear that latter does not typically imply the type of time scale separation typically used to obtain reduced CMEs (see for example Ref. 16) because it does not lead to a partitioning of reactions into fast and slow types; this is since within our abundance separation method, each limit of a rate constant tending to zero is counterbalanced by the opposite limit of a concentration of an abundant species tending to infinity. Yet it is not difficult to show that our abundance separation limit does lead to a separation of the eigenvalues of the Jacobian and hence point to time scale separation of concentration transients on the deterministic level. Thus, a deeper investigation into the relationship between abundance separation and time scale separation could improve understanding of both types of separation and as well lead to a clearer picture regarding when CMEs can be effectively reduced.

REFERENCES

1.

N.

van Kampen

,

Stochastic Processes in Physics and Chemistry

, 3rd ed. (

Elsevier

,

Amsterdam

,

2007

).

Google Scholar

2.

R.

Grima

,

D.

Schmidt

, and

T.

Newman

,

J. Chem. Phys.

137

,

035104

(

2012

).

https://doi.org/10.1063/1.4736721

Google Scholar

Crossref

PubMed

3.

V.

Shahrezaei

and

P. S.

Swain

,

Proc. Natl. Acad. Sci.

105

,

17256

(

2008

).

https://doi.org/10.1073/pnas.0803850105

Google Scholar

Crossref

4.

S.

Smith

and

V.

Shahrezaei

,

Phys. Rev. E

91

,

062119

(

2015

).

https://doi.org/10.1103/PhysRevE.91.062119

Google Scholar

Crossref

5.

T.

Jahnke

and

W.

Huisinga

,

J. Math. Biol.

54

,

1

(

2007

).

https://doi.org/10.1007/s00285-006-0034-x

Google Scholar

Crossref

PubMed

6.

D. T.

Gillespie

,

J. Phys. Chem.

81

,

2340

(

1977

).

https://doi.org/10.1021/j100540a008

Google Scholar

Crossref

7.

E. L.

Haseltine

and

J. B.

Rawlings

,

J. Chem. Phys.

117

,

6959

(

2002

).

https://doi.org/10.1063/1.1505860

Google Scholar

Crossref

8.

K.

Takahashi

,

H.

Kaizu

,

B.

Hu

, and

M.

Tomita

,

Bioinformatics

20

,

538

(

2004

).

https://doi.org/10.1093/bioinformatics/btg442

Google Scholar

Crossref

PubMed

9.

H.

Salis

and

Y.

Kaznessis

,

J. Chem. Phys.

122

,

054103

(

2005

).

https://doi.org/10.1063/1.1835951

Google Scholar

Crossref

10.

E.

Weinan

,

D.

Liu

, and

E.

Vanden-Eijnden

,

J. Chem. Phys.

123

,

194107

(

2005

).

https://doi.org/10.1063/1.2109987

Google Scholar

Crossref

PubMed

11.

R.

Erban

,

I. G.

Kevrekidis

,

D.

Adalsteinsson

, and

T. C.

Elston

,

J. Chem. Phys.

124

,

084106

(

2006

).

https://doi.org/10.1063/1.2149854

Google Scholar

Crossref

PubMed

12.

S. L.

Cotter

,

K. C.

Zygalakis

,

I. G.

Kevrekidis

, and

R.

Erban

,

J. Chem. Phys.

135

,

094102

(

2011

).

https://doi.org/10.1063/1.3624333

Google Scholar

Crossref

PubMed

13.

S.

Cotter

, “

Constrained approximation of effective generators for multiscale stochastic reaction networks and application to conditioned path sampling

,” e-print arXiv:1506.02446 (unpublished).

14.

S.

Peleš

,

B.

Munsky

, and

M.

Khammash

,

J. Chem. Phys.

125

,

204104

(

2006

).

https://doi.org/10.1063/1.2397685

Google Scholar

Crossref

PubMed

15.

B.

Munsky

and

M.

Khammash

,

J. Chem. Phys.

124

,

044104

(

2006

).

https://doi.org/10.1063/1.2145882

Google Scholar

Crossref

PubMed

16.

Y.

Cao

,

D. T.

Gillespie

, and

L. R.

Petzold

,

J. Chem. Phys.

122

,

014116

(

2005

).

https://doi.org/10.1063/1.1824902

Google Scholar

Crossref

17.

P.

Thomas

,

A. V.

Straube

, and

R.

Grima

,

BMC Syst. Biol.

6

,

39

(

2012

).

https://doi.org/10.1186/1752-0509-6-39

Google Scholar

Crossref

PubMed

18.

P.

Thomas

,

R.

Grima

, and

A. V.

Straube

,

Phys. Rev. E

86

,

041110

(

2012

).

https://doi.org/10.1103/PhysRevE.86.041110

Google Scholar

Crossref

19.

A.

Hellander

and

P.

Lötstedt

,

J. Comput. Phys.

227

,

100

(

2007

).

https://doi.org/10.1016/j.jcp.2007.07.020

Google Scholar

Crossref

20.

A.

Alfonsi

,

E.

Cancès

,

G.

Turinici

,

B.

Di Ventura

, and

W.

Huisinga

, in

ESAIM: Proceedings

(

EDP Sciences

,

2005

), Vol.

14

, pp.

1

–

13

.

21.

T.

Jahnke

,

SIAM Multiscale Model. Simul.

9

,

1646

(

2011

).

https://doi.org/10.1137/110821500

Google Scholar

Crossref

22.

T.

Jahnke

and

M.

Kreim

,

SIAM Multiscale Model. Simul.

10

,

1119

(

2012

).

https://doi.org/10.1137/120871894

Google Scholar

Crossref

23.

Y.

Taniguchi

,

P. J.

Choi

,

G.-W.

Li

,

H.

Chen

,

M.

Babu

,

J.

Hearn

,

A.

Emili

, and

X. S.

Xie

,

Science

329

,

533

(

2010

).

https://doi.org/10.1126/science.1188308

Google Scholar

Crossref

PubMed

24.

B.

Schwanhausser

,

D.

Busse

,

N.

Li

,

G.

Dittmar

,

J.

Schuchhardt

,

J.

Wolf

,

W.

Chen

, and

M.

Selbach

,

Nature

473

,

337

(

2011

).

https://doi.org/10.1038/nature10098

Google Scholar

Crossref

PubMed

25.

C. M.

Bender

and

S. A.

Orszag

,

Advanced Mathematical Methods for Scientists and Engineers: Asymptotic Methods and Perturbation Theory

(

Springer

,

Berlin

,

2010

).

Google Scholar

26.

R.

Grima

,

J. Chem. Phys.

133

,

035101

(

2010

).

https://doi.org/10.1063/1.3454685

Google Scholar

Crossref

PubMed

27.

R.

Ramaswamy

,

N.

Gonzalez-Segredo

,

I. F.

Sbalzarini

, and

R.

Grima

,

Nat. Commun.

3

,

779

(

2012

).

https://doi.org/10.1038/ncomms1775

Google Scholar

Crossref

PubMed

28.

D.

Gillespie

,

J. Chem. Phys.

113

,

297

(

2000

).

https://doi.org/10.1063/1.481811

Google Scholar

Crossref

29.

M.

Voliotis

,

P.

Thomas

,

R.

Grima

, and

C.

Bowsher

, “

Stochastic simulation of biomolecular networks in dynamic environments

,” preprint arXiv:1511.01268v1 (unpublished).

30.

D. F.

Anderson

,

J. Chem. Phys.

127

,

214107

(

2007

).

https://doi.org/10.1063/1.2799998

Google Scholar

Crossref

PubMed

31.

D. P.

Tschudy

and

H. L.

Bonkowsky

,

Mol. Cell. Biochem.

2

,

55

(

1973

).

https://doi.org/10.1007/BF01738678

Google Scholar

Crossref

PubMed

32.

P.

Thomas

,

A. V.

Straube

, and

R.

Grima

,

J. Chem. Phys.

133

,

195101

(

2010

).

https://doi.org/10.1063/1.3505552

Google Scholar

Crossref

PubMed

33.

D.

Gonze

,

J.

Halloy

, and

A.

Goldbeter

,

Proc. Natl. Acad. Sci.

99

,

673

(

2002

).

https://doi.org/10.1073/pnas.022628299

Google Scholar

Crossref

34.

P.

Thomas

,

A. V.

Straube

,

J.

Timmer

,

C.

Fleck

, and

R.

Grima

,

J. Theor. Biol.

335

,

222

(

2013

).

https://doi.org/10.1016/j.jtbi.2013.06.021

Google Scholar

Crossref

PubMed

35.

N.

van Kampen

,

Phys. Lett. A

5

,

333

(

1976

).

https://doi.org/10.1016/0375-9601(76)90398-4

Google Scholar

Crossref

36.

D.

Schnoerr

,

G.

Sanguinetti

, and

R.

Grima

,

J. Chem. Phys.

141

,

024103

(

2014

).

https://doi.org/10.1063/1.4885345

Google Scholar

Crossref

PubMed

37.

J.

Elf

and

M.

Ehrenberg

,

Genome Res.

13

,

2475

(

2003

).

https://doi.org/10.1101/gr.1196503

Google Scholar

Crossref

PubMed

2015

Author(s)

Model reduction for stochastic chemical systems with abundant species

I. INTRODUCTION

II. A REDUCED CHEMICAL MASTER EQUATION

A. A heuristic reduction of the CME

B. The abundant limit

C. Proof of N species abundant convergence

1. Taylor expansion of exact marginal distribution

2. Taylor expansion of the approximate reduced distribution

3. Convergence of full and reduced Taylor series

III. THE ACCURACY OF THE REDUCTION FOR FINITELY LARGE ABUNDANT CONCENTRATIONS AND SPECIAL CASES

A. Stochastic simulations

1. Homodimerisation

2. Genetic feedback loop

3. Metabolic network

4. Genetic oscillator with transcriptional feedback

B. Exact reductions

1. Systems in detailed balance

2. Open Michaelis-Menten reaction with one enzyme molecule

IV. ESTIMATING THE APPROXIMATION ERROR OF THE HYBRID MODEL

A. Open Michaelis-Menten reaction with multiple enzyme molecules

B. Open homodimerisation reaction

C. Genetic feedback loop

D. Genetic oscillator

V. SUMMARY AND DISCUSSION

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

Model reduction for stochastic chemical systems with abundant species

I. INTRODUCTION

II. A REDUCED CHEMICAL MASTER EQUATION

A. A heuristic reduction of the CME

B. The abundant limit

C. Proof of N species abundant convergence

1. Taylor expansion of exact marginal distribution

2. Taylor expansion of the approximate reduced distribution

3. Convergence of full and reduced Taylor series

III. THE ACCURACY OF THE REDUCTION FOR FINITELY LARGE ABUNDANT CONCENTRATIONS AND SPECIAL CASES

A. Stochastic simulations

1. Homodimerisation

2. Genetic feedback loop

3. Metabolic network

4. Genetic oscillator with transcriptional feedback

B. Exact reductions

1. Systems in detailed balance

2. Open Michaelis-Menten reaction with one enzyme molecule

IV. ESTIMATING THE APPROXIMATION ERROR OF THE HYBRID MODEL

A. Open Michaelis-Menten reaction with multiple enzyme molecules

B. Open homodimerisation reaction

C. Genetic feedback loop

D. Genetic oscillator

V. SUMMARY AND DISCUSSION

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only