Rigorous bounds on the stationary distributions of the chemical master equation via mathematical programming

TL;DR

This paper introduces mathematical programming methods to approximate stationary distributions of the chemical master equation with computable error bounds, enabling verification of accuracy and handling non-uniqueness in biochemical network models.

Contribution

It develops semidefinite and linear programming techniques to bound stationary distributions and moments, providing a rigorous, verifiable approximation framework for complex biochemical systems.

Findings

Provides tighter bounds on moments of stationary distributions.

Yields converging approximations of stationary distributions and ergodic distributions.

Demonstrates methodology on various biochemical models.

Abstract

The stochastic dynamics of biochemical networks are usually modelled with the chemical master equation (CME). The stationary distributions of CMEs are seldom solvable analytically, and numerical methods typically produce estimates with uncontrolled errors. Here, we introduce mathematical programming approaches that yield approximations of these distributions with computable error bounds which enable the verification of their accuracy. First, we use semidefinite programming to compute increasingly tighter upper and lower bounds on the moments of the stationary distributions for networks with rational propensities. Second, we use these moment bounds to formulate linear programs that yield convergent upper and lower bounds on the stationary distributions themselves, their marginals and stationary averages. The bounds obtained also provide a computational test for the uniqueness of the distribution. In the unique case, the bounds form an approximation of the stationary distribution with a computable bound on its error. In the non-unique case, our approach yields converging approximations of the ergodic distributions. We illustrate our methodology through several biochemical examples taken from the literature: Schl\"ogl's model for a chemical bifurcation, a two-dimensional toggle switch, a model for bursty gene expression, and a dimerisation model with multiple stationary distributions.

Figures

Equations

$$R_j:v_{1j}^{-}{S}_1+\cdots +v_{nj}^{-}{S}_n\stackrel{a_j}{\to }{v}_{1j}^{+}{S}_1+\cdots +v_{nj}^{+}{S}_n$$

$$j=1,\dots ,m,$$

$$q(x,y):=\sum _{j=1}^m{a}_j(x)(1_{x+v_j}(y)-1_x(y)),$$

$$\displaystyle 1_{y}(x):=\begin{cases*}1&if $x=y$\\ 0&otherwise.\end{cases*}$$

$$q(x,y)$$

$$-q(x,x)$$

$$\frac{\mathrm{d}{p}_t(x)}{\mathrm{d}t}=p_{t}Q(x),p_0(x)=\lambda (x)\forall x\in \mathcal{S},$$

$$p_{t}Q(x):=\sum _{y\in \mathcal{S}}{p}_t(y)q(y,x),$$

$$\pi Q(x)=0\forall x\in \mathcal{S},$$

$${\ell }^1:=\left\{(\rho (x){)}_{x\in \mathcal{S}}:\sum _{x\in \mathcal{S}}\left|\rho (x)\right|<\infty \right\}$$

$$\displaystyle\mathcal{P}:=\left\{\pi\in\ell^{1}:\begin{array}[]{l}\pi\,Q(x)=0,\,\,\forall x\in\mathcal{S},\\ \pi(\mathcal{S}):=\sum_{x\in\mathcal{S}}\pi(x)=1,\\ \pi\geq 0\end{array}\right\},$$

$${⟨f⟩}_{\pi }:=\sum _{x\in \mathcal{S}}f(x)\pi (x).$$

$$2S\underset{a_2}{\stackrel{a_1}{\rightleftarrows }}3S,\varnothing \underset{a_4}{\stackrel{a_3}{\rightleftarrows }}S,$$

$$a_1(x):=k_{1}x(x-1),a_2(x):=k_{2}x(x-1)(x-2)$$

$$a_3(x):=k_3,a_4(x):=k_{4}x,$$

$${⟨Q{x}^k⟩}_{\pi }:=\sum _{x^{\prime }\in \mathcal{S}}\sum _{x\in \mathcal{S}}q(x^{\prime },x)x^k\pi (x^{\prime })=0,k\in \mathbb{N}.$$

$${⟨Qx⟩}_{\pi }=b_1{⟨1⟩}_{\pi }-b_2{⟨x⟩}_{\pi }+b_3{⟨x^2⟩}_{\pi }-b_4{⟨x^3⟩}_{\pi }=0,$$

$${⟨1⟩}_{\pi }=1,$$

$${⟨x^2⟩}_{\pi }-{⟨x⟩}_{\pi }^2\ge 0.$$

$${⟨f^2⟩}_{\pi }$$

$${⟨x{f}^2⟩}_{\pi }$$

$$z:={\left(\begin{array}{c}{⟨1⟩}_{\pi },{⟨x⟩}_{\pi },{⟨x^2⟩}_{\pi },{⟨x^3⟩}_{\pi }\end{array}\right)}^T\in {\mathbb{R}}^4.$$

$${⟨f^2⟩}_{\pi }$$

$${⟨x{f}^2⟩}_{\pi }$$

$$M_3^0(y):=\left[\begin{array}{cc}{y}_0& y_1\\ y_1& y_2\end{array}\right],M_3^1(y):=\left[\begin{array}{cc}{y}_1& y_2\\ y_2& y_3\end{array}\right],$$

$$M_3^0(z)⪰0,M_3^1(z)⪰0.$$

$$z_0\ge 0,z_1\ge 0,$$

$$z_0{z}_2-z_1^2$$

$$z_1{z}_3-z_2^2$$

$$\displaystyle\mathcal{E}^{3}=\left\{y\in\mathbb{R}^{4}:\begin{array}[]{l}y_{0}=1\\ b_{1}y_{0}-b_{2}y_{1}+b_{3}y_{2}-b_{4}y_{3}=0\\ M_{3}^{0}(y)\succeq 0\\ M_{3}^{1}(y)\succeq 0\end{array}\right\}.$$

$$L_{\alpha }^3:=\inf \{y_{\alpha }:y\in {\mathcal{E}}^3\}$$

$$\alpha =0,1,2,3.$$

$$b_1{y}_1+b_3{y}_1{y}_2-b_4{y}_2^2-b_2{y}_1^2\ge 0.$$

$$r_2^{-}(y_1)\le y_2\le r_2^{+}(y_1)$$

Full text

I Introduction

Cell-to-cell variability is pervasive in cell biology. A fundamental source of this variability is the fact that biochemical reactions inside cells often involve only a few molecules per cell Elowitz et al. (2002); Taniguchi et al. (2010); Uphoff et al. (2016). Such reactions are key components in gene regulatory and signalling networks involved in cellular adaptation and cell fate decisions Arkin et al. (1998); Dandach and Khammash (2010); Thomas et al. (2014). Mathematically, stochastic reaction networks are modelled using continuous-time Markov chains whose distributions satisfy the chemical master equation (CME). As the availability of accurate single cell measurements widens, it is crucial to develop reliable methods for the analysis of the CME that facilitate parameter inference Munsky et al. (2015); Fröhlich et al. (2016), both for the identification of molecular mechanisms Neuert et al. (2013) and for the design of synthetic cellular circuits Strelkowa and Barahona (2010, 2011); Oyarzún et al. (2014); Zechner et al. (2016); Tomazou et al. (2018).

Significant effort has been devoted to investigating the stationary solutions of CMEs, which determine the long time behaviour of the stochastic process Meyn and Tweedie (1993a). While exactHemberg and Barahona (2007, 2008) Monte Carlo methods have been developed to sample from stationary solutions of some CMEs, analytical solutions are known only in a few special cases. In general, the CME is considered intractableSchnoerr et al. (2017) because, aside of systems with finite state space, it consists of an infinite set of coupled equations.

An approach to circumvent the intractability of the full CME is to compute moments of its stationary solutions. However, moment computations are only exact for networks of unimolecular reactions; in all other cases, the equations of lower moments involve higher moments, leading to an infinite system of coupled equations that cannot be solved analytically. Moment closure schemes, usually requiring assumptions about the unknown solution, are thus employed to approximate the moments Engblom (2006); Engblom and Pender (2014); Lakatos et al. (2015); Schnoerr et al. (2015); Van Kampen (1976). Yet few of these methods provide quantified approximation errors Engblom and Pender (2014). Another approach is provided by mathematical programming techniques, which have been employed to compute bounds on the moments of Markov processes in various contexts. In such schemes, the finite set of moment equations is supplemented by moment inequalities and the moments are bounded by solving linear programs (LPs) Schwerer (1996); Helmes and Stockbridge (2003) or semidefinite programs (SDPs) Hernández-Lerma and Lasserre (2003); Kuntz et al. (2016). Alternatively, when the CME has a unique stationary solution, several truncation-based schemesHart and Tweedie (2012); Gupta et al. (2017); Dayar et al. (2011a); Spieler (2014); Dayar et al. (2011b) have been proposed to approximate the solution, although in most cases they do not provide estimates of the error they introduce.

Here, we present two different mathematical programming approaches that yield bounds on, and approximations of, the stationary solutions of the CME. Our first approach builds on our previous work Kuntz et al. (2016, 2017); Kuntz (2017), and uses semidefinite programming to obtain upper and lower bounds on the moments of stationary solutions of networks with polynomial and rational propensities. The scheme constrains the possible solutions of a truncated, underdetemined set of moment equations by appending semidefinite inequalities that are satisfied by all probability distributions on the state space. Independently of this work, similar approaches have been recently proposed for networks with polynomial propensities Kuntz et al. (2017); Sakurai and Hori (2017); Dowdy and Barton (2017); Ghusinga et al. (2017). Here we extend the mathematical framework to rational networks of interest in biochemistry, and we state rigorous, checkable mathematical conditions for the validity of the approach.

The second approach employs linear programming to obtain lower and upper bounds on stationary averages of the CME, using state space truncations guided by moment bounds computed with our first (SDP) approach. In the case of a unique stationary solution, we prove that the LP bounds converge to the true average as the truncation approaches the entire state space. Because stationary averages can be tailored to bound distributions, the LP bounds provide approximations that converge in total variation to the stationary solution and its marginals with a computable approximation error (see also Ref. Kuntz et al. (2018a)). Additionally, the LP bounds provide a computational test for the uniqueness of the stationary solution, a prerequisite for most other approximation schemes. In the non-unique case, the scheme provides converging approximations of the ergodic distributions.

The paper is organised as follows. In Section II, we introduce definitions regarding stationary solutions of the CME. Section III presents the SDP method to bound moments: first, the conceptual framework is introduced analytically for a simple birth-death process that displays a chemical bifurcation, followed by the general computational approach for multi-species networks with rational propensities using semidefinite programming. Section IV presents the LP approach to bound and approximate entire stationary solutions, their averages and marginals: first, the mathematical framework is introduced through semi-analytical expressions for birth-death processes, followed by the general computational approach for rational networks using linear programming. In Section V, we apply the methods to three additional examples: a toggle switch, a model of bursty gene expression with negative feedback, and a model with multiple stationary solutions. We conclude with a discussion in Section VI. For completeness, Appendix A presents theoretical results linking stationary distributions of continuous-time chains and the stationary solutions of CMEs, and Appendix B presents a Foster-Lyapunov criterion that guarantees existence and finiteness of moments of the stationary distribution for the examples in the paper.

II Notation and definitions

Stochastic biochemical kinetics under well-mixed conditions are usually described by a set of $m$ reactions $R_{j}$ involving $n$ species $S_{1},S_{2},...,S_{n}$:

[TABLE]

where $v_{ij}^{\pm}\in\mathbb{N}$ denote the stoichiometric coefficients, and $a_{j}:\mathbb{N}^{n}\to[0,\infty)$ is the propensity of reaction $R_{j}$.

Formally, the state of the system is described by the random variable $X(t)=(X_{1}(t),\dots,X_{n}(t))\in\mathbb{N}^{n}$, a vector with components representing the number of molecules of each species at time $t$. The dynamical process is modelled with a minimal continuous-time Markov chainNorris (1997) with rate matrix $Q=(q(x,y))$ defined by:

[TABLE]

where $v_{j}:=(v_{1j}^{+}-v_{1j}^{-},\dots,v_{nj}^{+}-v_{nj}^{-})$ denotes the stoichiometric vector containing the net changes in molecule numbers produced by reaction $R_{j}$, and $1_{y}$ denotes the indicator function of state $y$:

[TABLE]

The state of the system takes values in a subset $\mathcal{S}\subseteq\mathbb{N}^{n}$, known as the state space, with (possibly infinite) cardinality $\left|\mathcal{S}\right|$. The set $\mathcal{S}$ must be chosen such that

[TABLE]

in which case $Q$ is said to be totally stable and conservative.

If the Markov chain cannot leave the state space in finite time, the matrix $Q$ is said to be regular (see Appendix A). In this case, the collection of probabilities $p_{t}(x)$ of observing the chain in state $x$ at time $t\geq 0$ is the only solution of the chemical master equation (CME)

[TABLE]

where we define the vector $p_{t}:=(p_{t}(x))_{x\in\mathcal{S}}$ and

[TABLE]

Following standard convention, probability distributions and measures are defined throughout as row vectors.

Any probability distribution $\pi:=(\pi(x))_{x\in\mathcal{S}}$ that solves the equation

[TABLE]

is called a stationary solution of the CME. By definition, each stationary solution belongs to the space

[TABLE]

of absolutely summable real sequences indexed by states in $\mathcal{S}$. The set of all stationary solutions forms a convex polytope in $\ell^{1}$:

[TABLE]

where $\pi\geq 0$ is shorthand for $\pi(x)\geq 0,\forall x\in\mathcal{S}$. For most networks of interest the stationary solutions determine the long-term behaviour of the chain (see Appendix A). In many cases, we will be interested in obtaining the $\pi$-average of a real-valued function $f$ on $\mathcal{S}$:

[TABLE]

For example, the $k$-th stationary moment is $\left\langle x^{k}\right\rangle_{\pi}$.

III Bounding the stationary moments of the CME

As a first use of optimisation techniques, we present a systematic approach that yields bounds of increasing tightness on the stationary moments of reaction networks with polynomial or rational propensities, and we give rigorous sufficient conditions for its validity. The moment bounds obtained in this section will be used in conjunction with linear programming to bound the full stationary distributions of the CME in Sec. IV.

To motivate our optimisation approach, we first present the mathematical formulation through a simple example, for which explicit analytical expressions can be obtained (Sec. III.1). For more complex systems, the approach can be implemented computationally in a systematic manner through a general semidefinite programming method (Sec. III.2). Readers interested in the computational approach (and not the theory behind it) should skip Sec. III.1 and go directly to Sec. III.2.

III.1 A simple analytical example: moment bounds for Schlögl’s model

To illustrate the mathematical framework, consider the classic autocatalytic network with a single species $S$ proposed by Schlögl Schlögl (1972) as a model for a chemical bifurcation:

[TABLE]

with mass-action propensities:

[TABLE]

where $k_{1},k_{2},k_{3},k_{4}>0$ are rate constants. Here, $n=1$ and the state space is $\mathcal{S}=\mathbb{N}$.

The CME of this network has a unique stationary solution $\pi$ and all of its moments are finite (see Appendix B). The reaction network, as encoded in the rate matrix $Q$, imposes certain relationships between the stationary moments $\left\langle x^{k}\right\rangle_{\pi}$. Such relations form an infinite system of coupled stationary moment equations:

[TABLE]

Except in particular instances, it is not possible to solve this coupled system exactly, and closure approximations are usually adopted by neglecting higher order moments.

An alternative approach is that of mathematical programming, which includes additional constraints in the form of inequalities that must be fulfilled by the moments of distributions. Including such inequalities allows us to obtain feasible regions for the solutions of the system, and hence rigorous bounds for the moments. Increasing the number of inequalities considered, restricts the feasible region further and makes the bounds tighter.

Let us consider the first moment equation for (13):

[TABLE]

where $b_{1}:=k_{3}$, $b_{2}:=k_{1}+2k_{2}+k_{4}$, $b_{3}:=k_{1}+3k_{2}$, $b_{4}:=k_{2}$ are positive numbers. Even after noting that

[TABLE]

as $\pi$ is a probability distribution, Eq. (17) is underdetermined. However, further relationships between moments can be added to constrain the system. For example, the non-negativity of the variance implies the inequality:

[TABLE]

Further inequalities involving higher moments can be built systematically from polynomial functions as follows.

Consider the polynomial $f(x):=f_{0}+f_{1}x$, where $\mathbf{f}:=(f_{0},f_{1})^{T}\in\mathbb{R}^{2}$ is a (column) vector of polynomial coefficients. Clearly, $f^{2}(x)$ and $x\,f^{2}(x)$ are non-negative on $x\in[0,\infty)$. Hence it follows that

[TABLE]

Let us define the vector of moments

[TABLE]

Then the inequalities  (20)–(21) are written compactly as

[TABLE]

where the matrices $M^{3}_{0}(y)$ and $M^{3}_{1}(y)$ are defined by

[TABLE]

for any four-dimensional vector $y=(y_{0},y_{1},y_{2},y_{3})^{T}\in\mathbb{R}^{4}$. Since (22)–(23) hold for all $\mathbf{f}\in\mathbb{R}^{2}$, we have that $M^{3}_{0}(z)$ and $M^{3}_{1}(z)$ are positive semidefinite (p.s.d.):

[TABLE]

From Sylvester’s criterion, (24) is equivalent to

[TABLE]

Hence (26) recovers (19), whereas (27) gives an additional condition involving the first three moments.

Putting (17)–(24) together, we conclude that the singleton set $\{z\}$ of vectors whose entries are composed of the first four moments of the stationary solutions of Schlögl’s model belongs to the set

[TABLE]

We say that $\mathcal{E}^{3}$ is an outer approximation of the set of stationary moment vectors. Hence the stationary moments have the following lower and upper bounds:

[TABLE]

Such bounds are usually handled computationally, but it is illustrative to obtain explicit expressions in this simple case. Combine (27) with the equalities in (32) to get

[TABLE]

Assuming $b_{3}\geq 2\sqrt{b_{2}b_{4}}$ (the other case is analogous), we use the quadratic formula to obtain the bound

[TABLE]

Hence $\mathcal{E}^{3}$ can be rewritten equivalently as:

[TABLE]

Figure 1(a) shows the projection of $\mathcal{E}^{3}$ onto the $y_{1}$-$y_{2}$ plane. From (40), it is clear that $(1,0,0,b_{1}/b_{4})\in\mathcal{E}^{3}$, and the lower bounds in (33) for the first two moments are trivial: $(L^{3}_{1},L^{3}_{2})=(0,0)$. The upper bounds, however, are not. Since $r^{+}_{2}(x)>r^{-}_{2}(x)$ for all $x>0$, the upper bounds in (33) are obtained by the northeasternmost intersection of $y_{1}^{2}$ and $r^{+}_{2}(y_{1})$: $(U^{3}_{1},U^{3}_{2})=(r_{4},r_{2}^{+}(r_{4}))$, where $r_{4}$ is the rightmost root of $x(b_{1}-b_{2}x+b_{3}x^{2}-b_{4}x^{3})=0$. In summary, we get:

[TABLE]

As seen in Fig. 1(a), the analytical bounds based on $\mathcal{E}^{3}$ are rough. However, we show in the following section how to obtain tighter bounds systematically by appending further moment equations and inequalities and solving the associated optimisations over higher order sets $\mathcal{E}^{d}$.

III.2 The general approach: Bounding the moments of rational CMEs by solving semidefinite programs

The approach in the previous section can be applied to any reaction network (1) with $n$ species and state space $\mathcal{S}$, as long as the propensities of its $m$ reactions are rational (or polynomial) functions, i.e., they can be rewritten as

[TABLE]

where $b_{1},\dots,b_{m}$ and $s$ are polynomials on $\mathbb{R}^{n}$ and the common denominator $s$ satisfies $s(x)>0,\forall x\in\mathcal{S}$.

To deal with multiple species, we use standard multi-index notation: $x^{\alpha}:=x_{1}^{\alpha_{1}}x^{\alpha_{2}}\dots x^{\alpha_{n}}_{n}$, where $\alpha$ is the multi-index $(\alpha_{1},\alpha_{2},\dots,\alpha_{n})\in\mathbb{N}^{n}$ and $\left|\alpha\right|:=\alpha_{1}+\alpha_{2}+\dots+\alpha_{n}$ is the degree of the monomial. Polynomial functions are expressed in terms of such monomials. For instance, the denominator of (42) is:

[TABLE]

where $d_{s}$ is the degree of $s$ and we define the (column) vector of coefficients

[TABLE]

We also define the stationary rational moments

[TABLE]

which are directly related to the raw moments:

[TABLE]

The following checkable assumption is a sufficient condition for our general SDP approach to apply to generic reaction networks with rational propensities.

Assumption 1** (Order of the approximation and finiteness of moments).**

Recall that $d_{s}$ is the degree of the denominator $s$ in (42). Let us denote the order of the approximation by an integer $d\geq d_{s}$, and let us compile the stationary rational moments (45) up to order $d$ into the (column) vector

[TABLE]

of dimension $\#_{d}:=\binom{n+d}{n}$.

We assume that all stationary solutions of the CME have finite rational moments up to order $d+1$, i.e.,

[TABLE]

This requirement can be verified using a Foster-Lyapunov criterion as detailed in Appendix B.

Remark 2**.**

From (46) and Assumption 1 it follows that raw moments $\left\langle x^{\alpha}\right\rangle_{\pi}$ with $\left|\alpha\right|\leq d+1-d_{s}$ are also finite.

In order to write down each $\alpha$-moment, it is helpful to define the associated polynomial function $g_{\alpha}(x)$:

[TABLE]

with degree $d_{g_{\alpha}}=\left|\alpha\right|+d_{b}-1$, where $d_{b}:=\max\{d_{b_{i}}\}$ is the maximum degree of the numerators in (42). We also define the (column) vector of polynomial coefficients

[TABLE]

where ${(g_{\alpha})}_{\beta}=0$ if $\left|\beta\right|>d_{g_{\alpha}}$.

The finiteness of moments guarantees that a subset of moment equations will hold, as stated in the following lemma.

Lemma 3** (The moment equations).**

If Assumption 1 is satisfied and $\pi\in\mathcal{P}$, then the $\alpha$-moment equation

[TABLE]

holds for every $\alpha\in\mathbb{N}^{n}$ such that $\left|\alpha\right|\leq d-d_{b}+1$.

Proof.

Consider the adjoint of (8):

[TABLE]

which, by Fubini’s theorem, is valid for any $f$ as long as $\sum_{x\in\mathcal{S}}\left|q(x,x)f(x)\right|\pi(x)$ is finite Glynn and Zeevi (2008). Because

[TABLE]

then Assumption 1 implies

[TABLE]

for all $\left|\alpha\right|\leq d-d_{b}+1$. Setting $f(x):=x^{\alpha}$ in (51), we get:

[TABLE]

∎

In addition to the moment equations (50), the moments $z$ satisfy additional constraints. Firstly, since $\pi$ is a probability distribution, we have:

[TABLE]

Furthermore, the rational moments satisfy well-known semidefinite inequalitiesLasserre (2009); Blekherman et al. (2013); Kuntz (2017). Specifically, the localising matrices are positive semidefinite:

[TABLE]

where the $M^{i}_{d}(y)$ are defined by

[TABLE]

with $e_{i}$ denoting the $i^{th}$ unit vector and $y\in\mathbb{R}^{\#_{d}}$.

The inequalities (53) follow from the fact that for any polynomial function $f(x)$ of degree $\lfloor d/2\rfloor$ with (column) vector of coefficients $\mathbf{f}=(f_{\beta})_{\left|\beta\right|\leq\lfloor d/2\rfloor}$, we have

[TABLE]

Similarly, it can be shownLasserre (2009); Blekherman et al. (2013) that

[TABLE]

Since (54)–(55) hold for any vector $\mathbf{f}$, the matrices $M^{i}_{d}(z)$ are positive semidefinite.

Remark 4**.**

In the case of Schlögl’s model (13), we had $s(x)=1$, $d_{b}=3$ and $\#_{d}=d+1$. Eq. (17) is the moment equation (50) with $\alpha=1$, and the matrices in (24) with $d=3$ are $M^{0}_{3}(y)$ and $M^{1}_{3}(y)$.

III.2.1

Bounding the moments

We can then establish the following lemma regarding outer approximations of the set of rational moments.

Lemma 5** (Outer approximations of the set of rational moments).**

If Assumption 1 is satisfied and $\pi\in\mathcal{P}$, the vector of rational moments $z$ belongs to the spectrahedron

[TABLE]

where $d_{b}$ is the maximum degree of the numerators in (42), and $\mathbf{s}$, $\mathbf{g_{\alpha}}$ and $M_{d}^{i}(y)$ are defined in (43),(49) and (53), respectively.

Proof.

The proof follows from Lemma 3, (52)–(53). ∎

In summary, the vectors of stationary moments of order $d$ are contained in a feasible set $\mathcal{E}^{d}$, defined by linear equalities and inequalities, which constitutes an outer approximation to the set of moment vectors. There are several implications of this lemma.

Firstly, the outer approximation property implies that extremal points of $\mathcal{E}^{d}$ provide bounds on the stationary moments. Specifically, the vector in $\mathcal{E}^{d}$ with largest (resp. smallest) $\alpha$-entry provides an upper (resp. lower) bound on the $\alpha$-moment. For example, Fig. 1(b) shows the projection of $\mathcal{E}^{d}$ with increasing $d$ onto the $y_{1}$–$y_{2}$ plane for Schlögl’s model. The northeasternmost (resp. soutwesternmost) vector of these outer approximations yield upper (resp. lower) bounds on the first two moments of Schlögl’s model.

Secondly, note that the moment matrix $M_{d}^{i}(y)$ is a principal submatrix of $M_{d+1}^{i}(y)$. Since a matrix is p.s.d. if and only if all of its principal submatrices are p.s.d. and $\mathcal{E}^{d+1}$ includes all moment equations in $\mathcal{E}^{d}$, then it follows that every vector in $\mathcal{E}^{d+1}$ (appropriately truncated) belongs to $\mathcal{E}^{d}$. As a result, the outer approximations tighten around the set of stationary moment vectors with bounds of increasing quality as the order of the approximation $d$ is increased, as seen inFig. 1(b). These two observations are summarised in the following corollary.

Corollary 6** (Monotonic moment bounds).**

*Suppose that Assumption 1 is satisfied and $\pi\in\mathcal{P}$. If $f$ is a polynomial of degree $d_{f}\leq d$, then *

[TABLE]

[TABLE]

Proof.

Note that $z\in\mathcal{E}^{d}$ because

[TABLE]

Hence (60) follows from (61)–(62). As explained in the main text, the monotonicity of the bounds follows from the definition of $\mathcal{E}^{d}$ and the fact that a matrix is p.s.d. if and only if all of its principal submatrices are p.s.d and $\mathcal{E}^{d+1}$ includes all moment equations in $\mathcal{E}^{d}$. ∎

Applying these results to $\left\langle x^{\alpha}\right\rangle_{\pi}$, the $\alpha$-moment of the CME, is straightforward. Let $f(x):=s(x)x^{\alpha}$ and choose $d\geq|\alpha|+d_{s}$ to obtain the bounds

[TABLE]

Corollary 6 establishes that outer approximations $\mathcal{E}^{d}$ of increasing order can be used to compute a monotonically increasing (resp. decreasing) sequence of lower (resp. upper) bounds for $\left\langle x^{\alpha}\right\rangle_{\pi}$:

[TABLE]

Remark 7** (The sequence of moment bounds is monotonic but may not converge).**

The monotonicity of the bounds does not imply that the gap between the bounds $(U_{\alpha}^{d}-L_{\alpha}^{d})$ will converge to zero as $d\to\infty$. Although in our experience the bounds often converge numerically, there is no general guarantee for several reasons. Firstly, the stationary solution may not be unique and, in that case, the lower bounds are limited by the stationary solution with the smallest moment while the upper bounds are limited by that with the largest moment. Even if the solution is unique, the bounds may not converge because the semidefinite conditions are tailored to distributions with support on the non-negative real orthant but not to distributions with support on discrete state spaces Dowdy and Barton (2018), for which more stringent conditions can be produced at a higher computational cost Lasserre (2002, 2009); Dowdy and Barton (2018).

III.2.2 Computing moment bounds via semidefinite programming

Given a reaction network with rational propensities and a polynomial $f$ of degree $d_{f}$, the moment bounds (61)–(62) are the extreme points of the linear functional $y\mapsto\mathbf{f}^{T}y$ over the set $\mathcal{E}^{d}$, which is defined by linear equalities and semidefinite inequalities. Hence, computing the bounds amounts to solving a semidefinite program (SDP), a convex optimisation problem for which there exist efficient computational tools. Therefore, instead of ad hoc analytical manipulations, like those leading to (41), a general procedure by constructing and solving the SDPs systematically is implemented as follows:

Rewrite the reaction propensities in the form (42) removing all common factors and setting $s$ to be the lowest common denominator. 2. 2.

Verify the existence of stationary solutions $\pi$ and choose the order of the approximation $d$, an integer $d\geq d_{f}$ for which the $d+1$ stationary moments are finite (Assumption 1) using a Foster-Lyapunov criterion (Theorem 31 in Appendix B). 3. 3.

Compute the bounds $L_{f}^{d}$ and $U_{f}^{d}$ by solving the two SDPs (61)–(62). We set up the SDPs using the modelling package YALMIPLöfberg (2004), and solve them using the multi-precision solver SDPA-GMPNakata (2010) with the interface mpYALMIPFantuzzi and Fuentes (2016). Examples of computation times are given in the figure captions.

SDPs involving high order moments can be numerically ill-conditioned Kuntz et al. (2016); Fantuzzi et al. (2016); Dowdy and Barton (2018). Although the origin of this numerical instability remains an open problem, our computations suggest that it could be the result of the rapid growth of moments, which leads to ill-conditioned moment matrices. Such disparity is problematic for standard double-precision SDP solvers but we have mitigated it with the multi-precision solver SDPA-GMPNakata (2010) as in Ref. Fantuzzi et al. (2016). Alternatively, one can scale the moments Kuntz et al. (2016); Fantuzzi et al. (2016); Dowdy and Barton (2018), or adapt recently developed specialised solvers Papp and Yildiz (2017); Zheng et al. (2018). 4. 4.

Evaluate the error of the bounds by computing the gap $U^{d}_{f}-L^{d}_{f}$. If the gap is unsatisfactorily large, increase the order of the approximation and return to Step 2 to compute new bounds.

Corollary 6 guarantees that the bounds will not loosen as we increase $d$, yet the bounds may stagnate (Remark 7). In this case, we recommend breaking the impasse by employing the LP approach of the next section (see Fig. 7(b)).

As an example of this procedure, Fig. 2(a) shows how the computed upper and lower bounds $(L^{d}_{\alpha},U^{d}_{\alpha})$ for the first three stationary moments of Schlögl’s model become tighter as we increase $d$, the order of the approximation. Fig. 2(b) combines the moment bounds to obtain bounds on commonly used statistics, e.g., variance, coefficient of variation, and skewness.

IV Bounding and approximating the stationary solutions of the CME

The semidefinite programming scheme in the previous section allows us to compute bounds of the stationary moments of the CME by constructing and optimising over outer approximations of the set of stationary moments. In this section, we go further and introduce a linear programming scheme that yields bounds on the full stationary solutions of the CME, their marginals and their averages by constructing and optimising over outer approximations of the set of stationary solutions of the CME. To do so, we use a moment bound obtained in the previous section. Note that the approximations in this section only involve linear inequalities (instead of semidefinite ones). Hence the optimisations to be solved are linear programs (instead of SDPs), a simpler subclass of convex optimisation problems for which mature, industrial solvers CPLEX (2009) are available.

As for the SDP scheme above, we introduce the mathematical approach through a simple semi-analytic example (i.e., birth-death processes) in Sec. IV.1, and then present the computational framework for general CMEs in Sec. IV.2. Readers interested in the computational implementation (and not the mathematical background) should skip Sec. IV.1 and go directly to Sec. IV.2.

IV.1 A simple example: bounding the stationary solution of birth-death processes

Birth-death processes are one-species reaction networks ($n=1$) with state space $\mathcal{S}=\mathbb{N}$ whose value $x\in\mathbb{N}$ changes by $\pm 1$ in each reaction:

[TABLE]

The specific birth-death process is defined by the functional form of the given propensities $a_{+}(x)$ and $a_{-}(x)$.

The stationary equations $\pi Q=0$ then read

[TABLE]

Assuming non-vanishing death rates

[TABLE]

it is well known that the unique stationary solution is:Gardiner (2009)

[TABLE]

with $\pi(0)$ given by the normalisation condition:

[TABLE]

where we have introduced the notation for sums over sets:

[TABLE]

Hence birth-death processes have at most one stationary solution, which exists if and only if $\gamma(\mathcal{S})$ is finite.

IV.1.1 Semi-analytical approach for bounds and approximation

For most birth-death processes, no closed-form expression for $\gamma(\mathcal{S})$ is known and, consequently, the stationary solution cannot be computed exactly. However, we can obtain upper and lower bounds for the distribution, as follows.

Let us consider a state space truncation

[TABLE]

with size $\left|\mathcal{S}_{r}\right|=\lceil r^{1/\alpha}\rceil$ controlled by the parameters $\alpha\in\mathbb{Z}_{+}$ and $r>0$. Let $\mathcal{S}_{r}^{c}$ denote its complement, i.e., the set of states outside of the truncation $\mathcal{S}_{r}$.

Let us assume that we have available an upper bound on the stationary $\alpha$-moment:

[TABLE]

Note that for rational propensities, such a bound can be computed with the SDP scheme in Sec. III.

Upper bound:

An easy upper bound on $\pi(0)$ is obtained by truncating the sum in (69) to get:

[TABLE]

whence it follows that

[TABLE]

Lower bound:

Using (72), we obtain a bound on the probability mass $m_{r}$ outside of the truncation:

[TABLE]

We say that $\varepsilon_{r}$ is a tail bound.

A lower bound on $\pi(0)$ then follows from (68)–(74):

[TABLE]

whence we obtain a lower bound for $\pi(x)$:

[TABLE]

Convergent bounds:

We have thus shown that

[TABLE]

and it is easy to see that both bounds converge to the stationary solution as the size of the truncation grows: as $r\to\infty$, both $u^{r}_{x}\to\pi(x)$ and $l^{r}_{x}\to\pi(x)$. This follows from (68)–(69) and $\varepsilon_{r}\to 0$.

Approximating the distribution and the approximation error:

Motivated by these facts, we define the two following measures (lower and upper bounds padded with zeros),

[TABLE]

and introduce them as approximations for $\pi(x)$:

[TABLE]

We quantify the approximation error of $\widetilde{\pi}$ with the total variation norm:

[TABLE]

where $\pi(A)$ and $\widetilde{\pi}(A)$ are sums over sets, defined in (70).

For $u^{r}$ and $l^{r}$, the approximation error can be chararacterised further. Using (73)–(83), we have:

[TABLE]

Since $\varepsilon_{r}=c/r\to 0$ as $r\to\infty$, it thus follows that both $l^{r}$ and $u^{r}$ converge in total variation to $\pi$.

We summarise these findings in the following theorem:

Theorem 8** (Bounds and approximations of the stationary solution of birth-death processes).**

Consider any birth-death process (64) with non-vanishing decay rates (67) and finite sum $\gamma(\mathcal{S})$ (69), such that it has a unique stationary solution $\pi$ (68). Suppose that $\pi$ satisfies the moment bound (72) and let $\mathcal{S}_{r}\subseteq\mathcal{S}$ be the truncation (71) of the state space controlled by the parameters $r,\alpha\in\mathbb{Z}_{+}$, with tail bound  $m_{r}=\pi(\mathcal{S}_{r}^{c})\leq c/r=\varepsilon_{r}$. Then the following hold:

1. (i)

The values of the distribution over the truncation are bounded above and below:

[TABLE]

where $u^{r}_{x}=\gamma(x)/\gamma(\mathcal{S}_{r})$ and $l^{r}_{x}=u^{r}_{x}(1-\varepsilon_{r})$. 2. (ii)

The measures $l^{r}=(l^{r}(x))_{x\in\mathcal{S}}$ and $u^{r}=(u^{r}(x))_{x\in\mathcal{S}}$ defined in (80)–(83) approximate the solution with approximation errors:

[TABLE] 3. (iii)

The bounds vary monotonically with $r$:

[TABLE]

and the sequences of approximations converge in total variation to $\pi$:

[TABLE]

Proof.

This follows from (73)–(86) and Corollary 6. ∎

Remark 9**.**

If Assumption 1 holds, then $\pi$ satisfies the moment bound (72) with $c=U^{d}_{\alpha}$, where $\alpha\in\{1,\dots,d-d_{s}\}$ and $U^{d}_{\alpha}$ is defined in (62).

**An application of Theorem 8:

Schlögl’s model**

To illustrate our results, we apply Theorem 8 to compute bounds on the unique stationary solution of Schlögl’s model (13)–(15). This model is a birth-death process for which an explicit analytical stationary solution can be obtained, thus allowing us to test the results directly without any simulations.

Through some analytical manipulations, the solution of Schlögl’s model can be obtained explicitly in terms of:

[TABLE]

where ${{}_{2}F_{2}}$ denotes the generalised hypergeometric function; $c_{1}:=\sqrt{1-4{k_{3}}/{k_{1}}}$; and $c_{2}:=\sqrt{1-4{k_{4}}/{k_{2}}}$. The stationary solution goes from being unimodal (black line, Fig. 3(a)) to bimodal (black line, Fig. 3(b)) depending on the parameter values, analogously to a bifurcation.

In Fig. 3(a)–(b), we compare the approximations $l^{r}$ and $u^{r}$ given in Theorem 8 (colour shades) to the analytical solution (black lines) of the unimodal and bimodal cases. As the size of the truncation (controlled by the parameter $r$) is increased, the bounds tighten around the analytical solution. In Fig. 3(c)–(d) we show that the approximation error tends to zero as the size of the truncation $\left|\mathcal{S}_{r}\right|=\lceil r^{1/\alpha}\rceil$ is increased. In the unimodal case (Fig. 3(c)), the approximation error decreases rapidly when $\left|\mathcal{S}_{r}\right|$ is larger than the mode. Furthermore, when the truncation is sufficiently large, employing bounds on higher order moments (larger $\alpha$) provides tighter tail bounds and smaller approximation errors. On the other hand, if the size of the truncation is smaller than the mode, using higher order moments does not necessarily improve the approximation error. A similar dependence of the approximation error is observed in the bimodal case (Fig. 3(d)), but the approximation error only decreases when the truncation size is larger than the second (larger) mode. This example shows how the ability to compute error bounds can reveal the presence of modes outside of the truncation.

IV.1.2 Reformulation of the bounds as optimisations

The truncation method leading to Theorem 8 relies on the detailed balanced structure of birth-death processes. However, this semi-analytic method is not generalisable to arbitrary reaction networks. Instead, the bounds can be reformulated as an equivalent (and generalisable) optimisation problem, as follows.

Consider a truncation $\mathcal{S}_{r}$ (71) controlled by the parameters $\alpha,r\in\mathbb{Z}_{+}$ with tail bound $m_{r}=\pi(\mathcal{S}_{r}^{c})\leq c/r=\varepsilon_{r}$.

Definition 10** (Restriction of $\pi$ to $\mathcal{S}_{r}$).**

The restriction of $\pi$ to $\mathcal{S}_{r}$ is:

[TABLE]

The restriction $\pi_{|r}$ belongs to the convex polytope:

[TABLE]

where $\mathcal{N}_{r}:=\{0,1,\dots,\lceil r^{1/\alpha}\rceil-1\}$, which follows directly from the definition of the restriction (88); the tail bound (74); and the fact that the stationary equations (65)–(66) with $x<\lceil r^{1/\alpha}\rceil-1$ only involve states inside of the truncation.

From the definition of the polytope (89), we can show that the bounds of the stationary solution in Theorem 8 (i) are recovered by optimising over $\mathcal{P}^{r}$, as stated in the following theorem.

Theorem 11** (Bounds and LP formulation).**

The bounds $l^{r}$ and $u^{r}$ in (80)–(83) are obtained by optimising over the polytope $\mathcal{P}^{r}$:

[TABLE]

Proof.

We only present the argument for the upper bounds—the proof for the lower bounds is analogous. If $x\not\in\mathcal{S}_{r}$, the result is trivial. A distribution $\pi^{r}$ satisfies (89) if and only if $\pi^{r}(x)=\gamma(x)\pi^{r}(0),\forall x\in\mathcal{S}_{r}$, where $\gamma(x)$ is given in (68). Therefore, we have

[TABLE]

which follows from $\pi^{r}(0)=\pi^{r}(\mathcal{S}_{r})/\gamma(\mathcal{S}_{r})\leq 1/\gamma(\mathcal{S}_{r}).$ However, $u^{r}$ clearly satisfies all constraints in (89), including the moment constraint:

[TABLE]

which follows from $\pi(x)=\gamma(x)/\gamma(\mathcal{S})$ and the inequality

[TABLE]

Hence $u^{r}\in\mathcal{P}^{r}$ and together with (90), this completes the proof. ∎

Importantly, the definition (89) only involves linear equations and inequalities. Therefore optimising over the polytope consists of solving a linear program, a class of optimisations for which there exist powerful computational platforms and algorithms. Furthermore, this optimisation reformulation can be extended seamlessly to arbitrary networks, as expanded in the next section.

IV.2 Generalisation to arbitrary networks via linear programming

We now generalise the optimisation approach to obtain bounds and approximations with controlled errors of the stationary solutions of arbitrary reaction networks.

Let us consider a reaction network (1) with state space $\mathcal{S}$, rate matrix $Q$ satisfying (4)–(5), and stationary solutions $\pi\in\ell^{1}$ that form the polytope $\mathcal{P}$ (12).

To characterise the solutions of the CME, we choose a norm-like function $w$, which plays the same role as the moment function ($x^{\alpha}$) in Section IV.1.2.

Definition 12** (Norm-like function).**

A function $w:\mathcal{S}\to\mathbb{R}$ is norm-like if it is non-negative

[TABLE]

and has finite sublevel sets:

[TABLE]

Furthermore, we require that the growth of $w$ be dominated by the stationary solution $\pi$, so that its expectation with respect to $\pi$ is finite. We summarise these requirements in the following checkable assumption.

Assumption 13** (Existence of CME solutions and moment bound).**

We assume that the CME has at least one stationary solution $\pi$, and that we have available a norm-like function $w$ with sublevel sets $\mathcal{S}_{r}$ such that every stationary solution $\pi$ satisfies

[TABLE]

where $c$ is a known constant. This inequality can be thought of as a generalisation of (72), hence we refer to it as a moment bound.

The existence of the stationary solutions can be verified on a case by case basis using a Foster-Lyapunov criterion (e.g., Theorem 31 in App. B).

Regarding the moment bound, in the case of networks with rational propensities satisfying Assumption 1, $w$ can be chosen to be any norm-like rational function with numerator of degree $d$ and the bounding constant $c$ can then be computed using the SDP approach of Sec. III. For general reaction networks, the moment bound can be obtained using Foster-Lyapunov criteria Glynn and Zeevi (2008) (see Appendix B).

In analogy with Lyapunov theory, the sublevel sets (92) of $w$ play an important role in characterising the stationary solutions of the CME. Specifically, we use the sublevel sets $\mathcal{S}_{r}$ as our state space truncations, noting that (93) allows us to establish a bound on the mass of the tail of the distribution outside of $\mathcal{S}_{r}$:

[TABLE]

Just as in the previous section, this choice yields a sequence of increasing truncations that approach the entire state space:

[TABLE]

For each truncation, let $\mathcal{N}_{r}$ denote the set of states $x\in\mathcal{S}_{r}$ that cannot be reached in a single jump from outside of the truncation:

[TABLE]

and the associated convex polytope:

[TABLE]

which, analogously to (89), includes all the stationary equations that only involve states in $\mathcal{S}_{r}$.

We then have the following lemma:

Lemma 14** (Outer approximations of $\mathcal{P}$).**

Suppose that Assumption 13 holds and let $\pi_{|r}$ be the restriction of $\pi$ to $\mathcal{S}_{r}$, as defined in (88). If $\pi\in\mathcal{P}$, then $\pi_{|r}\in\mathcal{P}^{r}$.

Proof.

This follows directly from (94) and the fact that $\pi Q(x)=\pi_{|r}\,Q(x),\,\forall x\in\mathcal{N}_{r}$. ∎

The outer approximation property means that optimising over $\mathcal{P}^{r}$ provides convergent bounds on the averages of functions $f$ on the state space, as summarised in the following theorem.

Theorem 15** (Convergent bounds of stationary averages).**

Consider a reaction network (1) with state space $\mathcal{S}$, rate matrix $Q$ satisfying (2)–(5), and stationary solutions $\pi$ forming the set $\mathcal{P}$ (12) and suppose that Assumption 13 holds.

If $\pi\in\mathcal{P}$ and $f:\mathcal{S}\mapsto\mathbb{R}$ is any real-valued function, then we can bound its averages over the restrictions:

[TABLE]

where $\pi_{|r}$ is the restriction of $\pi$ to $\mathcal{S}_{r}$ defined in (88) and the bounds are given by:

[TABLE]

If we have additional information on $f$, we have the following bounds on the full $\pi$-averages:

1. (i)

If $f(x)\geq 0,\,\forall x\not\in\mathcal{S}_{r}$, then $l^{r}_{f}\leq\left\langle f\right\rangle_{\pi},\enskip\forall r\in\mathbb{Z}_{+}$. 2. (ii)

If $f(x)\leq 0,\,\forall x\not\in\mathcal{S}_{r}$, then $\left\langle f\right\rangle_{\pi}\leq u^{r}_{f},\enskip\forall r\in\mathbb{Z}_{+}$. 3. (iii)

If $\left\langle\left|f\right|\right\rangle_{\pi}<\infty$ (i.e., $f$ is $\pi$-integrable), then

[TABLE] 4. (iv)

If the growth of $f$ is stricly dominated by $w$ as the size $r$ of the truncations $\mathcal{S}_{r}$ increases, i.e.,

[TABLE]

then $f$ is $\pi$-integrable for all $\pi\in\mathcal{P}$ and the sequences of bounds from below $(l^{r}_{f})_{r\in\mathbb{Z}_{+}}$ and above $(u^{r}_{f})_{r\in\mathbb{Z}_{+}}$ converge:

[TABLE]

Proof.

Eq. (97) follows directly from Lemma 14.

(i) and (ii) follow from (97) and $\left\langle f\right\rangle_{\pi}=\left\langle f\right\rangle_{\pi_{|r}}+\sum_{x\not\in\mathcal{S}_{r}}f(x)\pi(x)$.

(iii) is a consequence of (97), the moment bound (93), and the following inequality:

[TABLE]

(iv) has two parts. First, the $\pi$-integrability of $f$ follows from

[TABLE]

Second, the convergence of the bounds follows from the fact that every subsequence of $(l^{r}_{f})_{r\in\mathbb{Z}_{+}}$ has a converging subsequence $(l^{r_{k}}_{f})_{k\in\mathbb{Z}_{+}}$ with limit $\pi(f)$, where $\pi\in\mathcal{P}$ (see Remark 16 and Theorem 3.5 in Ref. Kuntz et al. (2018a)). By definition, $l_{f}\leq\pi(f)$, hence taking limits in (99) shows that $l_{f}\geq\pi(f)$ and the result follows. The case of the upper bounds is identical. For details, see Corollary 3.6$(i)$ in Ref. Kuntz et al. (2018a). ∎

Remark 16** (LP computation).**

The bounds $l^{r}_{f}$ and $u^{r}_{f}$ in (15) are obtained by solving two linear programmes (LPs) with $\left|\mathcal{S}_{r}\right|$ variables, $\left|\mathcal{N}_{r}\right|$ equality constraints, and $\left|\mathcal{S}_{r}\right|+3$ inequality constraints. LP solvers return the optimal value $l^{r}_{f}$ (or $u^{r}_{f}$) and an optimal point $\pi^{*,r}$, such that $\langle f\rangle_{\pi^{*,r}}=l^{r}_{f}$ (or $\left\langle f\right\rangle_{\pi^{r}}=u^{r}_{f}$). The optimal points exist because the LPs are optimisations of a continuous function over a compact non-empty subset of $\mathbb{R}^{\left|\mathcal{S}_{r}\right|}$.

Theorem 15 provides a general framework to obtain bounds (15) that can be used as approximations of stationary averages $\left\langle f\right\rangle_{\pi}$ with a quantifiable error given by (i)–(iii); furthermore, under the conditions in (iv), the approximations converge to $\left\langle f\right\rangle_{\pi}$ as the truncations approach the entire state space $\mathcal{S}$ of the reaction network.

IV.2.1 The case of a unique distribution: bounds and approximations

Throughout this section, we assume that the CME has a unique stationary solution $\pi$, i.e.,

[TABLE]

In this case, the results of Theorem 15 can be strengthened. and In particular, the feasible points $\pi^{r}\in\mathcal{P}^{r}$ are good approximations of the stationary solution in the sense that they converge to $\pi$ in weak∗, as detailed in the following corollary.

Corollary 17** (Convergence of bounds and feasible points for a unique solution).**

Let us assume that the conditions of Theorem 15(iv) hold and, in addition, that $\mathcal{P}=\{\pi\}$ consists of a single stationary solution $\pi$. Then the upper and lower bounds (15) converge to the average:

[TABLE]

and any sequence of feasible points $(\pi^{r})_{r\in\mathbb{Z}_{+}}$ belonging to the outer approximations $(\mathcal{P}^{r})_{r\in\mathbb{Z}_{+}}$ (i.e., $\pi^{r}\in\mathcal{P}_{r}$, $\forall r\in\mathbb{Z}_{+}$) converges to $\pi$ in weak∗:

[TABLE]

for any function $g$ that satisfies (100).

Proof.

This proof is similar to that of Theorem 15(iv). For full details, see Corollary 3.6(iii) in Ref. Kuntz et al. (2018a). ∎

Remark 18**.**

When $w$ is norm-like, convergence in weak implies convergence in total variation—see Appendix B of Ref. Kuntz et al. (2018a) for a proof.*

Remark 19**.**

Corollary 17 shows that, given a sufficiently large truncation, the feasible points $\pi^{r}$ provide arbitrarily accurate approximations to $\pi$, yet with no quantifiable bound on the approximation error $||\pi-\pi^{r}||$.

Bounding and approximating a unique stationary solution with quantifiable error:

Corollary 17 can still be used to obtain approximations with quantified errors of the distribution $\pi$ itself. Specifically, by applying (103) repeatedly with specific $f$s, the indicator function at each point of the truncation. This allows us to compute bounds on each value of the distribution $\pi(x)$ in $\mathcal{S}_{r}$.

Given a truncation $\mathcal{S}_{r}$, let us define the set of indicator functions $\{1_{x}\}_{x\in\mathcal{S}_{r}}$, one for every state in the truncation, where each $1_{x}$ is defined as (3). For each function $1_{x}$ in the set (i.e., for each state in the truncation), we compute the bounds (97) by solving the LPs (15) with $f=1_{x}$. Hence we obtain lower and upper bounds on $\left\langle 1_{x}\right\rangle_{\pi_{|r}}=\pi(x)$, for all states in the truncation:

[TABLE]

As in (80)–(83), we then collect these bounds, pad them with zeros, and define two approximations for $\pi$:

[TABLE]

These approximations of $\pi$ have controlled errors, as summarised in the following corollary.

Corollary 20** (Upper and lower bounding approximations of the unique solution).**

If Assumption 13 holds and $\mathcal{P}=\{\pi\}$, then the approximations $l^{r}$ (110) and $u^{r}$ (113) fulfill the following:

1. (i)

The approximations $l^{r}$ and $u^{r}$ bound $\pi$ from below and above, respectively,

[TABLE]

with quantified approximation errors $\varepsilon^{l}_{r}$ and $\varepsilon^{u}_{r}$:

[TABLE]

where $\left|\left|{\cdot}\right|\right|$ denotes the total variation norm (84). 2. (ii)

As the truncation size $r$ approaches infinity (and $\mathcal{S}_{r}$ approaches $\mathcal{S}$), the approximation $u^{r}$ converges pointwise to the unique $\pi$,

[TABLE]

and the approximation $l^{r}$ converges to $\pi$ in weak∗* (104) to $\pi$.*

Proof.

The bounds (114)–(115) follow directly from Theorem 15. The error (116) follows from (114)–(115) and the fact that the total variation norm of an unsigned measure is its mass. Similarly, (117) follows from (94) and

[TABLE]

Theorem 15(iv) shows that $l^{r}$ and $u^{r}$ converge pointwise to $\pi$. For any $f$ satisfying (100) and $r,r^{\prime}\in\mathbb{Z}_{+}$, we have that

[TABLE]

Using the pointwise convergence of $l^{r}$ and (102), we can pick $r^{\prime}$ such that the second sum is arbitrarily small and, subsequently, an $r$ such that the first sum is arbitrarily small. Hence the weak∗ convergence of $l^{r}$ follows. See Theorem 4.1 in Ref. Kuntz et al. (2018a) for details. ∎

Corollary 20 states that, for sufficiently large $r$, $l^{r}$ and $u^{r}$ are close to $\pi$. In contrast with the feasible points $\pi^{r}$, we can answer the question “is $r$ sufficiently large?” by evaluating the errors $\varepsilon^{l}_{r}$ (116) and $\varepsilon^{u}_{r}$ (117). Since

[TABLE]

we will always find an approximation $l^{r}$ that verifiably meets any given error tolerance by increasing $r$.

Remark 21**.**

Although we have no proof that $\varepsilon^{u}_{r}$ converges to zero (nor that $u^{r}$ itself converges to $\pi$ in total variation), all the examples we have encountered in practice exhibit convergence of $u^{r}$ and $\varepsilon^{u}_{r}$.

To instead answer the question ‘when is $r$ too small?’ (i.e., to establish how large the partition $\mathcal{S}_{r}$ must be to guarantee a given approximation error), the following proposition is of use.

Proposition 22** (Achievable approximation errors).**

Under the same conditions as in Corollary 20, the errors of the approximations $l^{r}$ (110) and $u^{r}$ (113) cannot be made smaller than the tail bound or the mass of the tail, respectively, i.e.,

[TABLE]

*where $m_{r}=\pi(\mathcal{S}^{c}_{r})\leq\varepsilon_{r}=c/r$. *

Proof.

The inequality (120) follows directly from (117). For (119), recall that there exists at least one optimal point $\pi^{*,r}$ such that $\pi^{*,r}(x)=l^{r}(x)$ for any $x\in\mathcal{S}_{r}$ (Remark 16). It is straightforward to verify that

[TABLE]

which implies $\pi^{*,r}(\mathcal{S}_{r})=1-\varepsilon_{r}$ (due to the minimality of $\pi^{*,r}(x)$). Since $l^{r}$ bounds from below all feasible points of $\mathcal{P}^{r}$, it follows that $l^{r}(\mathcal{S}_{r})\leq\pi^{*,r}(\mathcal{S}_{r})=1-\varepsilon_{r}$. Combined with (116), this gives (119). ∎

In other words, the approximation error of the lower bounds is no smaller than the tail bound, whereas that of the upper bounds is no smaller than the tail mass.

Remark 23**.**

The inequalities (119)–(120) are sharp for birth-death processes with $w(x)=x^{\alpha}$ (see (85)–(86)).

Approximating marginal distributions:

For high-dimensional state spaces, we are often interested in marginal distributions rather than the full multivariate solution $\pi$ defined on $\mathcal{S}$. A marginalisation is associated with a partition of the state space into a collection of disjoint subsets:

[TABLE]

and the marginal distribution is defined with respect to each subset:

[TABLE]

Because $\{A_{i}\}_{i\in\mathcal{I}}$ is a partition of $\mathcal{S}$, $\hat{\pi}$ is a probability distribution on $\mathcal{I}$.

Typically, we are interested in marginalising the distribution of a reaction network with $n$ species (and state space $\mathcal{S}=\mathbb{N}^{n}$) over a subset of species. For instance, if we are interested in the molecule counts of species $k$, we consider the following (infinite) set of subsets:

[TABLE]

whose union trivially recovers the entire state space. Associated with this set $\{A_{i}\}_{i\in\mathbb{N}}$ we then have the marginal distribution $\hat{\pi}$

[TABLE]

which, in this case, corresponds to the (univariate) distribution describing the molecule counts of the $k^{th}$ species.

The marginal distribution $\hat{\pi}$ can also be bounded and approximated following a similar procedure to the one described above for the full distribution. Using the indicator functions $1_{A_{i}}$ as the functions $f$, we solve the analogous LPs:

[TABLE]

for all the subsets $A_{i}$ that intersect with the truncation, i.e., $\mathcal{I}_{r}=\{i\in\mathcal{I}:A_{i}\cap\mathcal{S}_{r}\neq\varnothing\}$.

As before, we construct two approximations by padding (123)–(124) with zeros:

[TABLE]

which are the analogues for the marginal distribution of the approximations to the entire distribution (110)–(113), and have similar (but not identical) properties, as summarised in the following two corollaries.

Corollary 24** (Lower bounding approximation of the marginal distribution).**

Let us assume that the conditions of Corollary 20 hold. If $\{A_{i}\}_{i\in\mathcal{I}}$ is a partition of $\mathcal{S}$, then the associated marginal distribution $\hat{\pi}$ (121) is lower bounded by the approximation $\hat{l}^{r}$ defined in (127):

[TABLE]

with approximation error

[TABLE]

Furthermore, as the truncation size is increased ($r\to\infty$ and $\mathcal{S}_{r}$ approaches $\mathcal{S}$), $\hat{l}^{r}$ converges to $\hat{\pi}$ in total variation.

Proof.

This proof is analogous to that of Corollary 20 except that one needs to use Corollary 4.3 in Ref. Kuntz et al. (2018a) (with $g(x,y)=1_{A_{y}}(x)$ instead of Theorem 4.1 in Ref. Kuntz et al. (2018a). ∎

Corollary 25** (Convergent approximation of the marginal distribution).**

Under the same conditions as in Corollary 24, $\hat{u}^{r}$ defined in (130) approximates the marginal distribution $\hat{\pi}$ with error bounded by

[TABLE]

Furthermore, as the truncation $\mathcal{S}_{r}$ approaches $\mathcal{S}$, $\hat{u}^{r}$ converges pointwise to $\hat{\pi}$:

[TABLE]

Proof.

The proof is analogous to that of Corollary 20. ∎

Note that $\hat{u}^{r}$ does provide a controlled approximation of the marginal distribution, as it is a pointwise convergent approximation to $\hat{\pi}$ with a guaranteed, computable error bound $\hat{\varepsilon}^{u}_{r}$ (132).

Remark 26** (Upper bounds for the marginal distribution).**

The approximation $\hat{u}^{r}(i)$ bounds the marginal $\hat{\pi}(i)$ if and only if $A_{i}\cap\mathcal{S}_{r}^{c}=\varnothing$, i.e., when the set $A_{i}$ is fully contained inside the truncation $\mathcal{S}_{r}$. Hence $\hat{u}^{r}(i)$ does not provide an upper bound if the truncation does not include all the space of the marginalised variables.

However, using the fact that the probability mass of $A_{i}\cap\mathcal{S}_{r}^{c}$ is bounded by the mass of the tail $m_{r}$ (94), we have the following easy (but loose) upper bounds:

[TABLE]

IV.2.2 Non-uniqueness, ergodic distributions and a uniqueness test

Theorem 15 shows that our LP optimisation over the polytopes $\mathcal{P}^{r}$ yields bounds on the stationary averages, even if there are multiple stationary solutions. In the non-unique case, however, the gap between the lower bounds and the upper bounds will reflect the fact that the extreme points of $\pi\mapsto\left\langle f\right\rangle_{\pi}$ over $\mathcal{P}$ can be achieved by different solutions in the polytope. Yet it is possible to characterise further the set of solutions and the extreme points in terms of the ergodic distributions of the CME, and use this description to turn our LP approach into a test of uniqueness, as we show below.

To see how multiple stationary solutions of the CME can arise, consider the simple reaction network

[TABLE]

with mass action kinetics. It is clear that its state space $\mathcal{S}=\mathbb{N}^{2}$ decomposes into three disjoint sets:

[TABLE]

where $\mathcal{C}_{1}$ and $\mathcal{C}_{2}$ are closed communicating classes and $\mathcal{T}$ contains the remaining states. A set $\mathcal{C}\subseteq\mathcal{S}$ is a closed communicating class Norris (1997) if the chain can transit between any pair of states in $\mathcal{C}$ but cannot leave $\mathcal{C}$. Another common source of multiple closed communicating classes are conservation laws in reaction networksVallabhajosyula et al. (2005). For instance, the reactions

[TABLE]

conserve the quantity $n=x_{1}+2x_{2}$. Hence there exists a different closed communicating class for every $n\in\mathbb{N}$.

The closed communicating classes are intimately related to the stationary solutions, as summarised in the following theorem that compiles some facts that are broadly known in the literature.

Theorem 27** (Ergodic distributions and communicating classes Meyn and Tweedie (1993a)).**

Consider a reaction network (1) with rate matrix $Q$ satisfying (2)–(5), assume $Q$ is regular, and decompose the state space as

[TABLE]

where $\mathcal{C}_{j}$ are closed communicating classes and $\mathcal{T}$ contains the remaining states.

1. (i)

For each $\mathcal{C}_{j}$, there is at most one stationary solution $\pi_{j}$; hence $\pi_{j}(\mathcal{C}_{j})=1$. Whenever it exists, $\pi_{j}$ is known as the ergodic distribution associated with $\mathcal{C}_{j}$. 2. (ii)

Let $\mathcal{J}$ be the set of indexes $j$ of the ergodic distributions $\pi_{j}$. The set of stationary solutions $\mathcal{P}$ (12) is the set of convex combinations of the ergodic distributions:

[TABLE]

Proof.

If $Q$ is regular, the stationary solutions of the CME are the stationary distributions of the chain (Theorem 30 in Appendix A). Using this fact, $(i)$ can be found in most books on continuous-time chains (e.g., Ref. (Norris, 1997, Th. 3.5.2)), and $(ii)$ is given in Ref.(Meyn and Tweedie, 1993a, Th. 3.4). ∎

Theorem 27 states that the ergodic distributions $\pi_{j}$ are orthogonal to each other and that they are the extreme points of the convex polytope $\mathcal{P}$ (12) of stationary solutions of the CME. Because $\mathcal{P}$ is contained in the non-negative orthant of $\ell^{1}$, it follows that each face of the non-negative orthant contains at most one of the ergodic distributions.

Using this fact, we obtain a computational test of the uniqueness of stationary solutions, as summarised in the following corollary.

Corollary 28** (A uniqueness test).**

If Assumption 13 holds and $Q$ is regular, then

[TABLE]

Proof.

The proof is by contradiction. Suppose that $\mathcal{P}$ is not a singleton. Theorem 27$(ii)$ implies that $\mathcal{P}$ contains two or more ergodic distributions, $\pi_{j}$, each associated with a different closed communicating class, $\mathcal{C}_{j}$. Let us consider, e.g., a state $x\in\mathcal{C}_{1}$. Since the classes $\mathcal{C}_{j}$ are disjoint, then $\pi_{j}(x)=0,\forall j\neq 1$, which contradicts the lower bound property of $l^{r}$ (Theorem 15$(i)$). Hence, $\mathcal{P}$ must be a singleton. The converse follows from the convergence of the bounds in Corollary 20 (ii). ∎

In other words, if $\pi$ is unique, then the lower bound of any (and all) states in the support of $\pi$ is non-zero for sufficiently large $r$. Conversely, finding a single non-zero lower bound for any $x\in\mathcal{S}$ provides a proof of uniqueness of the distribution. Hence, if there is more than one ergodic distribution, all the lower bounds are zero for all states in the state space $\mathcal{S}$.

When $\pi$ is unique, $l^{r}$ or $u^{r}$ are good approximations of the stationary solution. However, this is not so in the non-unique case. Indeed, it is easy to show that, in the non-unique case, the lower and upper bounds are always loose: Corollary 28 shows that the lower bounds are trivially zero everywhere ($l^{r}=0$), whereas Theorems 15$(ii)$ and 27 imply that for large $r$, the mass of $u^{r}$ will be no smaller than the number of ergodic distributions:

[TABLE]

hence the upper bound is not tight.

However, our LP framework can still be used to obtain approximations of the ergodic distributions by using sequences of feasible points $\pi^{r}$. To this end, we require the following generalisation of Corollary 17.

Corollary 29** (Convergent approximations of ergodic distributions).**

Let $f$ be any function that satisfies (100) (i.e., $f$ is dominated by the norm-like function $w$ as the size $r$ of the truncations increases), and let $(\pi^{*,r})_{r\in\mathbb{Z}_{+}}$ be a sequence of optimal points such that

[TABLE]

for all $r\in\mathbb{Z}_{+}$.

1. (i)

If there exists a unique point $\pi^{*}\in\mathcal{P}$ such that

[TABLE]

then the sequence of optimal points $\pi^{*,r}$ converges to $\pi^{*}$ in weak∗* as $r\to\infty$.* 2. (ii)

If $f$ is the indicator function of a state (or of a subset) that is contained in a closed communicating class $\mathcal{C}_{j}$ with associated ergodic distribution $\pi_{j}$, then the sequence of optimal points $(\pi^{*,r})_{r\in\mathbb{Z}_{+}}$ converges to $\pi_{j}$ in weak as $r\to\infty$.*

Proof.

The proof of $(i)$ is similar to that of Theorem 15(iv). See Corollary 3.6$(ii)$ and Remark 3.7 in Ref. Kuntz et al. (2018a) for details. $(ii)$ follows from $(i)$, Theorem 27(ii), and the fact that indicator functions satisfy (100). ∎

Corollary 29 provides a rationale for how to use our computational framework to obtain approximations of the ergodic distributions $\pi_{j}$ in the non-unique case. Importantly, we do not need to know a priori what the closed communicating classes are. Using the indicator function for a chosen state $x$, we obtain the sequence of optimal points $\pi^{r}$ satisfying $\pi^{r}(x)=u^{r}(x)$. Should $x$ belong to a closed communicating class $\mathcal{C}_{j}$ with ergodic distribution $\pi_{j}$, Corollary 29$(ii)$ shows that $\pi^{r}$ will converge to $\pi_{j}$ as $r$ tends to infinity. Indeed, by looking at the states for which $\pi^{r}(x)>0$, we can in principle deduce which communicating class $\mathcal{C}_{j}$ (if any) the state belongs. Once the class is known, we replace $\mathcal{S}$ with $\mathcal{C}_{j}$ and proceed as for the unique case to obtain bounds on $\pi_{j}$. See Section V.3 for an example of the application of this procedure.

IV.2.3 Computational implementation and numerical considerations

Let us consider a given reaction network with rational propensities. In order to obtain approximations of its stationary solutions with controlled error smaller than a tolerance $\epsilon$, we proceed as follows:

Verify the existence of stationary solutions $\pi$ and the finiteness of their moments (Assumption 1) using a Foster-Lyapunov criterion (Theorem 31 in App. B). 2. 2.

Choose a norm-like rational function $w$ and define the truncations $\mathcal{S}_{r}$ as the sublevel sets (92) controlled by $r$.

We have found it best to choose functions $w$ that define truncations that cover most of the probability mass and that tend quickly to infinity, so that the size of the truncation grows slowly with $r$. For example, if we take $w(x)=x^{\alpha}$, then higher values of $\alpha$ induce smaller truncation sizes $\left|\mathcal{S}_{r}\right|\approx\sqrt[\alpha]{r}$. To guide the selection of $w$, one can run the scheme with various $w$ to gain information about the shape of the distribution. 3. 3.

Use the SDP approach of Sec. III.2 to find a moment bound (93) with bounding constant $c$ satisfied by all stationary solutions (Assumption 13). In particular, we employ YALMIPLöfberg (2004), SDPA-GMPNakata (2010), and mpYALMIPFantuzzi and Fuentes (2016) to formulate and solve the SDP (62) with $f:=w$. See details in Sec. III.2.2. 4. 4.

Choose an initial truncation size $r$ based on the achievable errors established in Proposition 22. To guarantee an error smaller than our tolerance $\epsilon$, we must choose an initial $r>c/\epsilon$. 5. 5.

Solve the LPs (106) to obtain the upper and lower approximations $l_{r}$ (110) and $u_{r}$ (113). Here we use the dual simplex algorithm of CPLEX V12.6.3CPLEX (2009) to solve the LPs. The tool YALMIPLöfberg (2004) is convenient to formulate the LPs but, if speed is a priority, the model can be fed straight to the solver to avoid computational overheads.

For large truncations, the coefficients in the constraints of the LPs span many orders of magnitude, leading to round-off errors in double-precision arithmetic and poor solver performance. One way to ameliorate this issue is to scale the decision variables; in particular, scaling $\pi^{r}(x)$ by $-q(x,x)$ or $w(x)$ often significantly improves solver performance. 6. 6.

Evaluate the error of the approximations $l_{r}$ and $u_{r}$ using (116) and (117), respectively. If the error is larger than our tolerance $\epsilon$, we increase the truncation size $r$ and return to the previous step. 7. 7.

In addition to approximating the full distribution, we can apply the above steps to compute other measures of interest by changing the LPs and associated errors in Steps 5–6:

• •

If we want to approximate a marginal distribution, we solve the LPs (123)–(124) and quantify the error using (131)–(132).

• •

If we are interested in a particular stationary average $\left\langle f\right\rangle_{\pi}$, we instead solve the LPs (15) and control the error using the bounds in Theorem 15$(i)$–$(iii)$. 8. 8.

As the particular reaction network could have several stationary solutions, we check in Step 5 for non-trivial lower bounds, $l_{r}(x)>0$. If we find one such bound, the solution is unique (Corollary 28). Otherwise, we investigate further the uniqueness question by increasing $r$ and recomputing the lower bounds to examine the presence of communicating classes as discussed in Corollary 29.

Our computations were carried out on a desktop computer with a 3.5GHz processor and 16GB of RAM.

V Application to biological examples

We now present the application of the methodology to three examples. First, we showcase how to obtain tight bounds on the stationary solution (and marginals) of a two-dimensional toggle switch. Second, we consider a model of bursty gene expression with negative feedback, through which we explore the capabilities of our method to deal with promoter switching noise. Third, we demonstrate the application of our methods to the non-unique case with a dimerisation network. The code used to compute the approximations and bounds for this last example is available at Kun .

V.1 A toggle switch

Toggle switches are common motifs in many cell-fate decision genetic circuitsGardner et al. (2000); Hemberg and Barahona (2007); Perez-Carrasco et al. (2016). A simple such circuit consists of two mutually repressing genesGardner et al. (2000). In particular, we consider the asymmetric case with mutual repression modelled via Hill functions and dilution/degradation modelled via linear decay:

[TABLE]

The state space of the CME is $x\in\mathcal{S}=\mathbb{N}^{2}$ with $x=(x_{1},x_{2})$, where $x_{1}$ and $x_{2}$ denote the number of protein $P_{1}$ and $P_{2}$, respectively, and the propensities of the reactions are:

[TABLE]

where the $k_{i}>0$ are kinetic constants and $\theta>0$ is the dissociation constant of $P_{1}$.

We follow the steps detailed in Section IV.2.3 to obtain bounds and approximations for this reaction network. First, we show that a stationary solution $\pi$ exists and that all of the moments of every solution are finite using a Foster-Lyapunov criterion (App. B).

We pick the norm-like function

[TABLE]

and compute the moment bound

[TABLE]

by solving the SDP (62) with $d:=10$ and $f:=w$ (solver time = $3.6$ minutes).

We then solve the LPs (106)–(107) and compute the bounding approximations $u^{r}$ and $l^{r}$ of the stationary solutions. The fact that the lower bounds $l^{r}$ are non-zero provides us with a proof of uniqueness of the stationary solution. Figure 4 shows the bounds for small ($r=45$) and large ($r=75$) state space truncations. The maximum absolute discrepancies are found near the modes, and by increasing the size of the truncation, the upper and lower bounds become nearly indistinguishable—the maximum discrepancy drops under $10^{-4}$ (Fig. 4(d)). Overall, the total approximation error including the tail is less than $2.6\times 10^{-3}$, as given by (116)–(117).

We have also used our method to obtain approximations on the marginal distributions of the number of proteins $P_{1}$ ($x_{1}$) and $P_{2}$ ($x_{2}$) (Sec. IV.2.1). The results in Fig. 5(a) show that the bounds get tighter for truncations of increasing $r$ (although, as discussed in Remark 26, the upper bound (133) remains loose). Note, however, that the approximation $\hat{u}^{r}$ in Fig. 5(b) rapidly approaches the Gillespie numerical simulations. Fig. 5(c) shows that the errors of both $\hat{l}^{r}$ and $\hat{u^{r}}$ can be made arbitrarily small by increasing the truncation size (Corollaries 24–25).

Finally, we apply the method to chart the change of the stationary solution as a function of a parameter. In particular, the dissociation constant of protein $P_{2}$ ($\theta$) can be thought of as a bifurcation parameter: increasing $\theta$ allows for higher expression of protein $P_{1}$. Fig. 6 presents the lower bounds $\hat{l}$ on the marginals of both proteins. At small values of $\theta$, we observe a single population with high numbers of $P_{2}$ repressing $P_{1}$. At large values of $\theta$, the opposite happens: the population we observe has high numbers of $P_{1}$ repressing $P_{2}$. For intermediate $\theta$, we observe coexistence of both populations. Indeed, we find that the modes of the marginal distributions are in good correspondence with the stable solutions of the deterministic steady-state rate equations.

V.2 Bursty gene expression with negative feedback

As a second example, consider a model of bursty production of a protein that regulates (negatively) its own expression. The model Kumar et al. (2014) involves a promoter that switches between active ($G_{\text{on}}$) and inactive ($G_{\text{off}}$) states, and the protein $P$ it encodes. When the promoter is on, the protein is expressed in bursts of size $b$, a geometrically distributed random variable Shahrezaei and Swain (2008) with mean $\langle b\rangle$. The protein represses its own production by switching off the promoter:

[TABLE]

The state space of the CME is $x=(x_{1},x_{2})\in\mathcal{S}=\{0,1\}\times\mathbb{N}$, where $x_{1}=\{0,1\}$ is a binary variable describing the off/on state of the promoter and $x_{2}\in\mathbb{N}$ represents the protein count. The propensities are

[TABLE]

where the $k_{i}>0$ are reaction rate constants. In App. B we show that the network has a unique stationary solution $\pi$ and that all of its moments are finite.

This example provides an interesting test case for SDP methods since the protein noise is particularly large: $CV(x_{2})$, the coefficient of variation of $x_{2}$, grows Kumar et al. (2014) with the burst size $\langle b\rangle$. Therefore, we expect that getting tight bounds for the CV will entail the use of a large number of moment equations. To investigate the effect of such large noise on the efficacy of our SDP method, we compute the following bounds:

[TABLE]

where we use (61)–(62) and we append the following equalities to our SDP (59):

[TABLE]

Figure 7(a) shows how the bounds (145) get tighter as we increase the number of moment equations in our SDP calculations. As expected, for small mean burst sizes ($\langle b\rangle=1$), the bounds become tight with $10$ moment equations, but tightening the bounds becomes difficult when the burst size is larger ($\langle b\rangle=10,100$).

Computing tight bounds for large $\langle b\rangle$ with the naive SDP approach would thus require a prohibitive number of moment equations. However, we can apply the LP method of Sec. IV to overcome this limitation. To do this, use SDP to compute a (cheaper) loose upper bound on the sixth moment $L^{13}_{x_{2}^{6}}=2.3\times 10^{13}\leq\left\langle x_{2}^{6}\right\rangle_{\pi}\leq U^{13}_{x_{2}^{6}}=4.5\times 10^{13}$ (time=1 minute per bound). We then set $w(x):=x_{2}^{6}$ and $c:=U^{13}_{x_{2}^{6}}$ in Theorem 15 to obtain:

[TABLE]

which we combine as in (145) to obtain much tighter bounds on CV($x_{2}$). Figure 7(b) shows the convergence of these tight bounds for CV($x_{2}$) with $\langle b\rangle=100$ as the truncation size is increased. These results exemplify the fact that it is enough to obtain loose SDP bounds on a higher order moment in order to obtain arbitrarily tight LP bounds on lower order moments (Theorem 15).

Finally, we exemplify in Figure 7(c) another use of our capability to bound marginal distributions following the steps in Section IV.2.3. In this case, we marginalise over the on/off promoter variable ($x_{1}$) and we compute upper bounds $\hat{u}^{r}$ and lower bounds $\hat{l}^{r}$ (visually indistinguishable in Fig. 7(c)) on $\hat{\pi}(x_{2})$, the distribution of protein counts, for three different burst sizes, $\langle b\rangle=1,10,100$. As expected, the protein distribution widens considerably (yet still with light tails) as the burst size increases. To compute these bounds, we set

[TABLE]

and solve the LPs (123)–(124) to obtain (127)–(130). Note that this marginalisation is over the complete, untruncated domain of the marginalised variable $x_{1}=\{0,1\}$. As a result, the $\hat{u}^{r}$ do provide upper bonds in this case (Remark 26).

V.3 Dimerisation network with multiple stationary distributions

To illustrate the use of our method on CMEs with multiple stationary solutions (see Section IV.2.2), we consider the reversible dimerisation network

[TABLE]

with state space $\mathcal{S}=\mathbb{N}$ and mass action kinetics: $a_{1}(x):=k_{1}$ and $a_{2}(x):=k_{2}\,x_{1}(x_{1}-1)$. The reactions preserve the parity of the number of molecules; hence the even numbers $\mathcal{C}_{0}$ and the odd numbers $\mathcal{C}_{1}$ are closed communicating classes. Furthermore, because the network obeys detailed balanceVan Kampen (1976), it is straightforward to obtain analytical expressions for the two ergodic distributions:

[TABLE]

where $\mu:=\sqrt{k_{1}/k_{2}}$ and the normalising constants are given by $Z_{0}:=\cosh(\mu)$ and $Z_{1}:=\sinh(\mu)$.

To illustrate the use of our tools in such a non-unique case, suppose we were not aware of the above facts and instead apply our computational procedure. First, we use a Foster-Lyapunov criterion and find that the rate matrix is regular; at least one stationary solution exists; all of the moments of each stationary solution are finite; and, for each closed communicating class, there exists an ergodic distribution with support in that class (App. B).

Next, we obtain bounds on the stationary solutions. Using the norm-like function

[TABLE]

we solve the SDP (62) with $d:=7$ and $f:=w$ (solver time = $3$ seconds) to obtain the moment bound

[TABLE]

We then solve the LPs (106)–(107) to compute upper and lower bounds, $u^{r}$ and $l^{r}$ (Fig. 8(a)). Note that the lower bounds remain trapped at zero, indicating that the stationary solution is non-unique (as given by Theorem 28). Furthermore, we observe that the mass of the upper bounds approaches two as the truncation grows (i.e., $u^{r}(\mathcal{S}_{r})\to 2$ as $r\to\infty$), hinting that there exist two ergodic distributions corresponding to two closed communicating classes (c.f. (136)).

To identify the communicating classes, we note that the optimal points $\pi^{*,r}_{0}$ and $\pi^{*,r}_{1}$ of the linear programs

[TABLE]

approach distributions with support on the even and odd numbers, respectively (Fig. 8(b)). These observations indicate that the communicating classes are the even numbers $\mathcal{C}_{0}$ and the odd numbers $\mathcal{C}_{1}$ (Corollary 29).

To verify this claim, we compute two different sets of upper and lower bounds over $\mathcal{C}_{0}$ and $\mathcal{C}_{1}$ separately. When computed over the subsets $\mathcal{C}_{0}$ and $\mathcal{C}_{1}$ separately, the lower bounds are non-zero, thus attesting to the uniqueness of each stationary distribution over its communicating class. Using $r=25^{7}$, the total variation error of the lower bounds (116) is approximately $8\times 10^{-3}$, while that of the upper bounds (117) is bounded above by $8\times 10^{-3}$. In fact, the actual errors of the upper bounds computed from the analytical expressions (147) are $<10^{-4}$, i.e., substantially smaller than the guaranteed bound.

VI Discussion

We have introduced two mathematical programming approaches that yield bounds on: (i) the stationary moments, and (ii) the stationary distributions of biochemical reaction networks. These statistical quantities typically satisfy infinite sets of coupled equations: (i) the stationary moment equations and (ii) the stationary CME. Both our approaches consider a subset of these equations and employ: (i) semidefinite programming and (ii) linear programming to bound the set of solutions. The bounds we obtain provide converging estimates of moments and probabilities with quantifiable errors.

Regarding our first method, which provides bounds for stationary moments, recently, and independently of our work, SDP-based procedures have been proposed by several authors Kuntz et al. (2017); Dowdy and Barton (2017); Ghusinga et al. (2017); Sakurai and Hori (2017). Our work differs from those works in two ways: firstly, our results apply to networks with both polynomial and rational propensities, a wider class of networks of interest in biochemistry, beyond the mass action models considered in Refs. Dowdy and Barton (2017); Ghusinga et al. (2017); Sakurai and Hori (2017); secondly, we give mathematically precise conditions for the validity of the method (Assumption 1) and we explain how these conditions can be verified in practice. To the best of our knowledge, our second approach, the LP bounding and approximation procedure for probability distributions, has not appeared in the CME literature (see Ref. Kuntz et al. (2018a) for a discussion of related methodologies in the optimisation literature). Importantly, both methods are tightly interlinked: our second method uses SDP moment bounds to formulate the LPs, in order to obtain controlled approximations of stationary solutions and marginals.

For some CMEs, the SDP approach might need to include a large number of moments to obtain accurate estimates of lower order moments (see Figs. 2 and 7). Such large SDPs pose a computational challenge for larger networks, as the number of moments $\#_{d}$ is $\binom{n+d}{d}$, where $n$ denotes the number of species in the network and $d$ is the maximum moment order, and thus explodes combinatorially with the number of species. Similar costs are encountered when using moment-closure methods Lakatos et al. (2015); Schnoerr et al. (2015). In contrast with moment closure methods Engblom (2006); Engblom and Pender (2014); Lakatos et al. (2015); Schnoerr et al. (2015), however, the proposed SDP method to bound moments yields approximations with quantified errors. Furthermore, we show that repeated applications of our SDP method yield upper (resp. lower) bounds that are monotonically decreasing (resp. increasing) as the number of moment equations and inequalities is increased (Theorem 5). Although, as mentioned in Sec. III.2, there are reaction networks for which the bounds do not converge to the exact moments, they often converge in practice (Fig. 2 and other examples in Refs. Dowdy and Barton (2017, 2018); Ghusinga et al. (2017); Sakurai and Hori (2017, 2018)). In addition, when tight SDP bounds prove computationally too expensive, our LP approach can be used to tighten bounds on moments of interest employing a loose SDP bound on a higher moment (Fig. 7). Lastly, it is possible to obtain sharper SDP bounds for restricted state spaces Lasserre (2009); Laurent (2009); Blekherman et al. (2013); Dowdy and Barton (2018), but these refinements are beyond the scope of this paper.

As stated above, the LP approach produces convergent bounds on the stationary solutions (including their marginals or averages). To do so, it uses a moment bound obtained with the SDP method. It is worth remarking that, while we have limited ourselves to rational networks where the moment bound can be obtained using SDPs, the LP approach can be extended beyond rational propensities by using Foster-Lyapunov criteria Glynn and Zeevi (2008). If the CME has a unique solution, our LP method yields converging lower and upper bounds on this solution and easy-to-evaluate error bounds (Corollaries 20, 24 and 25). If the CME has multiple solutions, our method provides bounds over the set of possible solutions. Furthermore, the procedure can be adapted to infer the closed communicating classes and to compute converging approximations of (and bounds on) the corresponding ergodic distributions. Additionally, if we are unsure whether the stationary solution is unique, our method provides a uniqueness test (Corollary 28) that settles the question.

Although LP solvers are highly mature and scalable, the applicability of our LP approach can present computational challenges. Firstly, as discussed in Sec. IV.2.3 (Step 5), the LPs can become ill-conditioned if the truncation is large, although this issue is mitigated by scaling the variables and by ongoing improvements in LP solvers. Secondly, although the computational cost of solving an LP depends on the algorithm, the cost per bound is at least $\mathcal{O}(\left|\mathcal{S}_{r}\right|)$, where $\left|\mathcal{S}_{r}\right|$ is the size of the truncation. For the purpose of computing the entire distribution, we need $\mathcal{O}(\left|\mathcal{S}_{r}\right|)$ such bounds; hence the cost is at least $\mathcal{O}(\left|\mathcal{S}_{r}\right|^{2})$. If computing a marginal distribution where $k$ species remain, we need $\mathcal{O}(\left|\mathcal{S}_{r}\right|^{k/n})$ bounds with an cost of at least $\mathcal{O}(\left|\mathcal{S}_{r}\right|^{1+k/n})$. If only a stationary average is of interest, the cost is at least $\mathcal{O}(\left|\mathcal{S}_{r}\right|)$, since only two bounds per average need to be computed. Note also that the truncation size typically grows combinatorially in the number of unbounded species, e.g., the number of states for a simplex truncation $\{x\in\mathbb{N}^{n}:x_{1}+\dots+x_{n}\leq M\}$ is $\binom{n+M}{n}$, where $M$ is an upper cut-off for the species count. Hence the cost of LPs suffers a combinatorial explosion in the number of species, as for all truncation-based methods Tweedie (1971); Hart and Tweedie (2012); Gupta et al. (2017); Dayar et al. (2011a); Spieler (2014); Dayar et al. (2011b).

Several other truncation-based schemes have been proposed to approximate the stationary solutions of the CMEHart and Tweedie (2012); Gupta et al. (2017); Dayar et al. (2011a, b). In contrast with ours, those schemes typically assume that the CME has a unique stationary distribution, which has to be verified separately Paulevé et al. (2014); Gupta and Khammash (2018). Perhaps most extensively studied is the truncation-and-augmentation (TA) scheme, originally proposed by SenetaSeneta (1967) for discrete-time chains. Its continuous-time counterpartTweedie (1971); Hart and Tweedie (2012); Gupta et al. (2017) converges in total variation for exponentially ergodic chains, monotone chains, and certain generalisations Hart and Tweedie (2012). However, bounds on the TA approximation error can be conservative and often involve constants that are difficult to compute in practice Hart and Tweedie (2012); Gupta et al. (2017); Tweedie (1998); Meyn and Tweedie (1994); Liu (2015); Masuyama (2017a, b); Liu et al. (2018); Liu and Li (2018). Spieler et al.Dayar et al. (2011a); Spieler (2014) overcame this issue by iterating the TA scheme and applying a tail bound derived from a Foster-Lyapunov criterion to bound the stationary distribution. Spieler’s truncation-based scheme is thus closest to ours. However, their scheme entails solving only systems of linear equations which, although cheaper to compute and simpler to implement than LPs, offer no guarantee of convergence and are only applicable in the unique case.

Another distinct feature of our method is that it enables the direct computation of bounds on the marginal distributions, without the need to compute bounds for each state of the joint distribution. Marginal distributions are of particular interest for the analysis of high-dimensional networks and for inference of model parameters from single cell data. Since our approach yields upper and lower bounds on the marginals, it can be used to bound the likelihood or likelihood ratios from experimental observations. This would be useful to extend the work in Ref. Fox et al. (2016) avoiding error redistribution using water filling methods, and aiding by selecting the size of truncations that are sufficient for parameter identifiability. Similar bounds on acceptance ratios could be used in Metropolis-Hastings algorithms to extend the applicability of our method to Bayesian inference. We therefore expect that our approach will be valuable not only for estimating distributions, but also for estimating model parameters from noisy single cell data where accurate approximations with quantified errors are needed.

VII Acknowledgements

We thank Justine Dattani and Diego Oyarzún for stimulating discussions, and Michela Ottobre and Jure Vogrinc for insights on the stability and long-term behaviour of continuous-time chains. JK gratefully acknowledges support through a BBSRC PhD studentship (BB/F017510/1); PT through a Fellowship of the Royal Commission for the Exhibition of 1851; GBS through the EPSRC Fellowship EP/M002187/1; and MB through the EPSRC grant EP/N014529/1 funding the EPSRC Centre for Mathematics of Precision Healthcare.

Appendix A Minimal continuous-time Markov chains, their long-term behaviour and stationary distributions

In practice, the chain $X=(X(t))_{t\geq 0}$ in Sec. II is often constructed by running the Gillespie AlgorithmGillespie (1976); Kendall (1950); Feller (1940), i.e., one starts the chain from a state $x$ sampled from an initial distribution $\lambda:=(\lambda(x))_{x\in\mathcal{S}}$. If $q(x):=-q(x,x)$ in (4)–(5) is zero, leave the chain at the $x$ for all time. Otherwise, wait an exponentially distributed amount of time with mean $1/q(x)$, sample $y\neq x$ from the probability distribution $(q(x,y)/q(x))_{y\neq x}$, and update the chain’s state to $y$ (we say that the chain jumps from $x$ to $y$ and we call the time at which it jumps the jump time). Repeat these steps starting from $y$ instead of $x$. All random variables sampled must be independent of each other.

If $T_{n}$ denotes the $n^{th}$ jump time, then the limit

[TABLE]

is known as the explosion time of the chain, i.e., the first instant by which the chain has left every finite subset of the state space(Kuntz et al., 2018b, Sec. 2.3). If no such explosion occurs, then $T_{\infty}=\infty$, and we say that the chain is non-explosive:

[TABLE]

where $\mathbb{P}_{\lambda}$ denotes the probability measure underlying the chain (the subscript $\lambda$ emphasises the fact that the starting state was sampled from the distribution $\lambda$). If (149) holds for every probability distribution $\lambda$ ($\lambda(x)\geq 0,\forall x\in\mathcal{S},\lambda(\mathcal{S})=1$), then the rate matrix $Q$ is regular.

We collectively refer to the probabilities

[TABLE]

of observing the process in the state $x:=(x_{1},\ldots,x_{n})\in\mathcal{S}$ at time $t\geq 0$ as the time-varying law of the chain. The time-veraying law is the minimal non-negative solution of the CME (6) (see Ref.(Kuntz, 2017, Cor. 2.21))

A probability distribution $\pi:=\{\pi(x)\}_{x\in\mathcal{S}}$ on $\mathcal{S}$ is said to be a stationary (or steady-state or invariant) distribution of the chain if sampling the chain’s starting position from $\pi$ ensures that it will be distributed according to $\pi$ for all time:

[TABLE]

Summing both sides of (150) over $x\in\mathcal{S}$ and taking the limit $t\to\infty$, we find that the chain is non-explosive when its starting location is sampled from a stationary distribution:

[TABLE]

Taking the derivative in time of (150), we find that stationary distributions are stationary solutions of the CME (6) (that is, it belongs to (12)). The reverse direction is slightly more complicated:

Theorem 30** (Theorem 2.41Kuntz (2017)).**

Let $X$ be a continuous-time chain with rate matrix $Q$ satisfying (4)–(5). A probability distribution $\pi$ on $\mathcal{S}$ is a stationary distribution of $X$ if and only if it is a stationary solution of the CME and the chain is non-explosive when initialised with law $\pi$ (i.e., (151) holds).

In particular, assuming that $Q$ is regular, $\pi$ is a stationary distribution if and only if it is a stationary distribution of the chain. In other words, (12) is an analytical (as in non-probabilistic) linear programming characterisationKuntz (2017); Kurtz and Stockbridge (1998) of the set of stationary distributions for regular $Q$. The non-explosivity in Theorem 30 is crucial: a counterexample is the birth-death process (64) with $a_{+}(x):=2^{2x}$ and $a_{-}(x):=2^{2x}/2$. In this case, the sum in (69) is finite showing that the CME has a unique stationary solution $\pi$ given by (68)–(69). However, (Reuter, 1957, Theorem 11) shows that the process is explosive for any initial distribution (including $\pi$) and it follows from (151) that no stationary distribution exists.

Stationary distributions are of interest because, if the chain is stable, then, regardless of the initial distribution $\lambda$, they determineMeyn and Tweedie (1993a) the chain’s long term behaviour. i.e., the time-varying law of the chain converges in total variation to a stationary distribution, $\pi$:

[TABLE]

Furthermore, the empirical distribution $\rho_{T}$, which denotes the fraction of the time interval $[0,T]$ that the chain will spend in a state $x$

[TABLE]

also converges to $\pi$

[TABLE]

In general, the stationary distributions featuring in (152) and (153) depend on the initial distribution $\lambda$ and on the starting location $X_{0}$, respectively.

If the chain starts in one of the closed communicating classes (defined in Sec. IV.2.2), then it can never escape the class. The convergence of the empirical distribution in (153) then implies that, in the stable case, there must exist at least one stationary distribution per closed communicating class $\mathcal{C}_{j}$ and that the stationary distribution must have support contained in $\mathcal{C}_{j}$ ($\pi_{j}(\mathcal{C}_{j})=1$). This distribution $\pi_{j}$ is unique and, if the initial distribution $\lambda$ has its mass contained in $\mathcal{C}_{j}$, thenMeyn and Tweedie (1993a) both the time-varying law and the empirical distribution converge to $\pi_{j}$ (in the sense that the $\pi$ featuring in both (152) and (153) is $\pi_{j}$). For this reason, $\pi_{j}$ is known as an ergodic distribution of the chain. The definition of the set $\mathcal{T}:=\mathcal{S}\backslash\cup_{j}\mathcal{C}_{j}$ featuring in the decomposition (134) implies that the chain visits any given state $x$ in $\mathcal{T}$ at most finitely many times (in particular $1_{x}(X_{t})\to 0$ as $t\to\infty$ $\mathbb{P}_{\lambda}$-almost surely). It follows that the chain’s paths must eventually either leave $\mathcal{T}$ or diverge to infinity. In the case of a stable chain, tending to infinity is not an option and so the chain eventually enters one of the closed communicating classes. It then follows from (150) that no stationary distribution $\pi$ such that $\pi(x)>0$ exists for a state $x$ in $\mathcal{T}$. Bringing this discussion togetherMeyn and Tweedie (1993a), we have that stationary distribution $\pi$ in (153) is the ergodic distribution $\pi_{j}$ of the closed communicating class $\mathcal{C}_{j}$ that the chain’s path eventually enters, while that in (152) is a weighted combination of the ergodic distributions where the weight given to $\pi_{j}$ is the probability that the chain ever enters $\mathcal{C}_{j}$. Theorem 27 in Sec. IV.2.2 follows from these facts.

Appendix B A Foster-Lyapunov criterion

In practice, verifying whether a chain is stable is done by applying a Foster-Lyapunov criterion. For our examples, we will use the following well-known criterion:

Theorem 31** (Foster-Lyapunov criterion Chen (1991); Meyn and Tweedie (1993a, b)).**

If there exist constants $K_{1}\in\mathbb{R}$, $K_{2}>0$, and a norm-like function $w$ (Definition 12) that satisfy

[TABLE]

then the following hold:

1. (i)

The rate matrix $Q$ is regular. 2. (ii)

There exists at least one stationary distribution. 3. (iii)

For each closed communicating class, there exists an ergodic distribution with support in that class. 4. (iv)

Every stationary distribution $\pi$ satisfies $\left\langle w\right\rangle_{\pi}\leq K_{1}/K-2<\infty.$ 5. (v)

The stationary distributions determine the long-term behaviour of the chain:

1. a)

for every deterministic starting distribution there exists a stationary distribution $\pi$ satisfying (152); 2. b)

for any starting distribution $\lambda$ and for $\mathbb{P}_{\lambda}$-almost every path, there exists a stationary distribution $\pi$ satisfying (153).

Proof.

Part $(i)$ is Ref. (Chen, 1991, Theorem 1.11) (see also Ref. (Meyn and Tweedie, 1993b, Theorem 2.1)). Parts $(ii)$–$(v)$ follow from Ref. (Meyn and Tweedie, 1993a, Theorems 8.1 and 8.2) and Ref. (Meyn and Tweedie, 1993b, Theorem 4.6).∎

We have used this criterion for our examples in the main text, as follows.

Schlögl’s model:

In the case of Schlögl’s model (13), fixing $w(x):=x^{d-2}$ for an integer $d>2$, we have

[TABLE]

where $g_{d-1}$ is a polynomial of degree $d-1$. Thus,

[TABLE]

and the supremum is finite. Taking $d\geq 3$, Theorem 31 tells us that (13) has a regular rate matrix and at least one stationary distribution; that the moments $\left\langle x^{1}\right\rangle_{\pi},\dots,\left\langle x^{d-2}\right\rangle_{\pi}$ are finite for any stationary distribution $\pi$; and that the limits (152)–(153) hold. Since we can choose ever larger $d$, we have finiteness of all moments. Uniqueness of the distribution follows from (68)–(69).

Toggle switch:

In the case of the toggle switch chain of Sec. V.1, setting $w(x):=(x_{1}+x_{2})^{d}$, we have

[TABLE]

where $g_{d-1}$ is a polynomial of degree $d-1$. For this reason, proceeding as we did above for Schlögl’s model, we have that $Q$ is regular, that the chain does have a stationary distribution, and that all of the moments are finite of each of the stationary distributions are finite. The non-trivial lower bounds in Fig. 4 and Corollary 28 show that it is unique.

Bursty gene model:

For the bursty gene expression model of Sec. V.2, let $w(x):=x_{2}^{d}$ and $(p(k))_{k\in\mathbb{N}}$ denotes the geometric distribution, $p(k)=(1-p(0))^{k}p(0)$ with $p(0)=1/(1+\langle b\rangle)$), and $\langle b^{l}\rangle$ denotes its $l^{th}$ moment. We then have

[TABLE]

where $g_{d-1}$ is a polynomial of degree $d-1$ (note that all moments of a geometric random variable are finite). Because the state space is $\{0,1\}\times\mathbb{N}$, $w$ is norm-like. Thus, proceeding as we did above for Schlögl’s model, we have that $Q$ is regular, that there exists at least one stationary distribution, and that each stationary distribution has all moments finite. For each set of parameter values, we solved LP (110) with $x=(0,0)$ and obtained a non-trivial lower on $\pi((0,0))$. Corollary 28 then showed that the stationary distribution is unique.

Dimerisation network with multiple stationary distributions:

For the network (146), choosing $w(x):=x^{d-1}$ with integer $d>1$, we have

[TABLE]

where $g_{d-1}$ is a polynomial of degree $d-1$. Proceeding as for Schlögl’s model above, we have that $Q$ is regular, that there exists at least one stationary distribution, and that each stationary distribution has all moments finite.