Rigorous bounds on the stationary distributions of the chemical master equation via mathematical programming
Juan Kuntz
Philipp Thomas
Guy-Bart Stan
Mauricio Barahona
% The stochastic dynamics of biochemical networks is 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. To fill this gap, 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. The bounds obtained provide a computational test for the uniqueness of these distributions. 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 two biochemical networks that exhibit bifurcations to multimodal behaviour: Schl\"ogl's model and a toggle switch model.
main.tex
\documentclass[twocolumn,aip,reprint,jcp %fleqn,
,cha]{revtex4-1}

% Enumitem package to generate concise itemised lists and to substitute arabic numerals by latin letters.
\usepackage{enumitem}

\usepackage{amsthm,mathtools}
%\usepackage[utf8]{inputenc}
\usepackage{amsmath,amssymb}
\usepackage[usenames,dvipsnames,svgnames,table]{xcolor}
\usepackage[colorlinks,linkcolor={blue!50!black},citecolor={blue!50!black}]{hyperref}

\newcommand{\pt}[1]{{\color{magenta}{#1}}}
\newcommand{\ptq}[1]{{\color{blue}{\textbf{[PT: #1]}}}}

% Appendix Lemma label
%\usepackage{chngcntr}
%\usepackage{apptools}
%\AtAppendix{\counterwithin{theorem}{section}}

%
%\setlength{\mathindent}{0.3cm}
% Remove brackets from numbering in List of References
% \makeatletter
% \renewcommand{\@biblabel}[1]{\quad#1.}
% \makeatother

%\numberwithin{equation}{section}


%%% Declaring enviroments %%%
\newtheorem{theorem}{Theorem}
\newtheorem{example}[theorem]{Example}
\newtheorem{assumption}[theorem]{Assumption}
\newtheorem{remark}[theorem]{Remark}
\newtheorem{condition}[theorem]{Condition}
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{lemma}[theorem]{Lemma}
%\newtheorem*{proof}{Proof}

%% Include all macros below
\renewcommand{\cal}[1]{\mathcal{#1}}
\newcommand{\bb}[1]{\mathbb{#1}}
\renewcommand{\r}{\mathbb{R}}
\newcommand{\rp}{\mathbb{R}_{\geq 0}}
\newcommand{\rn}{\mathbb{R}^n}
\newcommand{\cn}{\mathcal{C}^n}

\newcommand{\norm}[1]{\left|\left|{#1}\right|\right|}
\newcommand{\Mag}[1]{\left|{#1}\right|}
\newcommand{\Pb}{\mathbb{P}}
\newcommand{\Eb}{\mathbb{E}}
\newcommand{\Ebb}[1]{\mathbb{E}\left[#1\right]}

\newcommand{\Pbb}[1]{\mathbb{P}\left(#1\right)}
\newcommand{\Pbx}[1]{\mathbb{P}_x\left(#1\right)}
\newcommand{\Ebx}[1]{\mathbb{E}_x\left[#1\right]}
\newcommand{\Pby}[1]{\mathbb{P}_y\left(#1\right)}
\newcommand{\Eby}[1]{\mathbb{E}_y\left[#1\right]}
\newcommand{\Pbl}[1]{\mathbb{P}_\lambda\left(#1\right)} %_\gamma
\newcommand{\Ebl}[1]{\mathbb{E}_\lambda\left[#1\right]} %_\gamma
\newcommand{\Pbp}[1]{\mathbb{P}_\pi\left(#1\right)}
\newcommand{\Ebp}[1]{\mathbb{E}_\pi\left[#1\right]}
\newcommand{\Tr}{\mathrm{Tr}}
\newcommand{\pol}{\r[x]}


%\newcommand{\ms}{\cal{q}}
%\newcommand{\mm}{\cal{M}}
%\newcommand{\m}{\mathbb{M}^+}
\newcommand{\iprod}[2]{\langle #1,#2\rangle}

%% END MACROS SECTION


%% Actually used

\newcommand{\hd}{\#_d}
\newcommand{\hr}{\#_r^w}
\newcommand{\thr}{\tilde{\#}_r^w}

\newcommand{\ri}{\r^{\infty}}
\newcommand{\rs}{\r^S}

\renewcommand{\rhd}{\r^{\hd}}
\newcommand{\rhr}{\r^{\hr}}
\newcommand{\rthr}{\r^{\thr}}



\newcommand{\wmg}[1]{\left|{#1}\right|_w}
\newcommand{\mg}[1]{\left|{#1}\right|}

\newcommand{\s}{\mathcal{S}}

\newcommand{\n}{\mathbb{N}}
\newcommand{\nn}{\mathbb{N}^n} % REMOVE
\newcommand{\ns}{\mathbb{N}^S}
\newcommand{\nts}{\mathbb{N}^{\tilde{S}}}
\newcommand{\nhs}{\mathbb{N}^{\hat{S}}}
\newcommand{\zp}{\mathbb{Z}_+}

\newcommand{\inlineave}[1]{\langle #1 \rangle_p}
\newcommand{\ave}[1]{\left\langle #1 \right\rangle_\pi}
\newcommand{\aver}[1]{\left\langle #1 \right\rangle_{\pi^r}}
\newcommand{\averr}[1]{\left\langle #1 \right\rangle_{\pi_{|r}}}
\newcommand{\avel}[1]{\left\langle #1 \right\rangle_{l^r}}
\newcommand{\aveu}[1]{\left\langle #1 \right\rangle_{u^r}}
\newcommand{\aveq}[1]{\left\langle #1 \right\rangle_q}
\newcommand{\avepr}[1]{\left\langle #1 \right\rangle_{p^r}}
\newcommand{\aveprk}[1]{\left\langle #1 \right\rangle_{p^{r_k}}}
\renewcommand{\a}{\mathcal{A}}
\newcommand{\mmag}[1]{\left| #1 \right|}

\newcommand{\pdc}{\cal{P}^d_c}

\newcommand{\MB}[1]{\textcolor{red}{\textbf{[MB: #1]}}}
\newcommand{\jk}[1]{\textcolor{blue}{\textbf{[JK: #1]}}}
\newcommand{\blue}[1]{\textcolor{blue}{#1}}
\newcommand{\red}[1]{\textcolor{red}{#1}}
%\newcommand{\mbg}[1]{\textcolor{magenta}{#1}}
\newcommand{\mbg}[1]{\textcolor{magenta}{\textbf{[MB2: #1]}}}
\newcommand{\MBB}[1]{\textcolor{TealBlue}{\textbf{[MB2: #1]}}}
\newcommand{\guy}[1]{\textcolor{blue}{\textbf{[GBS: #1]}}}



%%% notation
\renewcommand{\pmb}[1]{#1}

\begin{document}
\title{Rigorous bounds on the stationary distributions of the chemical master equation via mathematical programming}
\author{Juan Kuntz}
\email{juan.kuntz08@imperial.ac.uk}
\affiliation{Department of Mathematics, Imperial College London, London, United Kingdom}
\affiliation{Department of Bioengineering, Imperial College London, London, United Kingdom}


\author{Philipp Thomas}
\email{p.thomas@imperial.ac.uk}
\affiliation{Department of Mathematics, Imperial College London, London, United Kingdom}

\author{Guy-Bart Stan}
\email{Corresponding author: g.stan@imperial.ac.uk}
\affiliation{Department of Bioengineering, Imperial College London, London, United Kingdom}

\author{Mauricio Barahona}
\email{Corresponding author: m.barahona@imperial.ac.uk}
\affiliation{Department of Mathematics, Imperial College London, London, United Kingdom}

\date{\today}

\begin{abstract}
%
The stochastic dynamics of biochemical networks is 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. To fill this gap, 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. The bounds obtained provide a computational test for the uniqueness of these distributions. 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 two biochemical networks that exhibit bifurcations to multimodal behaviour: Schl\"ogl's model and a toggle switch model. 




\end{abstract}


\maketitle


%
\section{Introduction}
%
%
% Significance
Cell-to-cell variability is pervasive in cell biology. A fundamental source of this variability is the fact that biochemical reactions inside living cells often involve species present in very small numbers (i.e., only a few molecules per cell)~\cite{elowitz2002,taniguchi2010,uphoff2016}. Such reactions are key components in gene regulatory and signalling networks involved in cellular adaptation and cell fate decisions \cite{arkin1998,Dandach2010,thomas2014}. 
Mathematically, stochastic reaction networks are modelled using continuous-time Markov chains whose distributions satisfy the chemical master equation (CME). As increasingly accurate single cell measurements become available, it is crucial to develop reliable methods for the analysis of the CME that can facilitate parameter inference~\cite{munsky2015,frohlich2016} to enable the identification of molecular 
mechanisms~\cite{neuert2013} and the design of synthetic cellular circuits~\cite{oyarzun2014,zechner2016}.

% Analytical approaches
Significant efforts have been devoted to investigating the stationary solutions of CMEs, which determine the long time behaviour of the stochastic process~\cite{Meyn1993a}. While exact\cite{hemberg2007,Hemberg2008} and inexact\cite{Asmussen2007} Monte Carlo methods have been developed to sample from stationary solutions of some CMEs, analytical solutions are known only in a few special cases. 
%
%These include birth-death processes\cite{Gardiner2009}, certain classes of unimolecular reaction networks \cite{jahnke2007,lopez2014}, networks obeying the conditions of detailed balance \cite{van1976,danos2013}, deficiency zero networks \cite{Anderson2010a,cappelletti2016}, and small model systems \cite{peccoud1995,bokes2012,Dattani2017}. 
%
More generally, the CME is considered intractable because, aside of systems with finite state space, it consists of infinitely many coupled equations. 

% Lit rev
An approach to circumvent the intractability of the 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, also leading to a system of infinitely many coupled equations that cannot be solved analytically. A variety of \emph{moment closure} schemes, usually requiring assumptions about the unknown solution, are employed to obtain approximations to the moments~\cite{engblom2006,engblom2014,lakatos2015,schnoerr2015,vanKampen1976}. Yet
few of these methods provide quantified approximation errors~\cite{engblom2014}. Independently, mathematical programming techniques have been employed to compute bounds on the moments of Markov processes in various contexts. In those instances, a finite set of moment equations is supplemented by moment inequalities so that the moments can be bounded by solving linear programs (LPs)~\cite{Schwerer1996,Helmes2003} or semidefinite programs (SDPs)~\cite{Hernandez-Lerma2003,Kuntz2016}. Alternatively, in the case that the CME has a unique stationary solution, there are several truncation-based schemes\cite{Tweedie1971,Hart2012,Gupta2017,Dayar2011,Spieler2014,Dayar2011a} that directly approximate the solution (we discuss these further in Sec.~\ref{sec:discussion}). The main drawback of these schemes is that they do not provide estimates of the approximation error and hence we do not know how to verify their accuracy. 

% So, even though, in some cases, some of these schemes are known to produce arbitrarily accurate approximations, there is no known way of verifying their accuracy in practice.

% What is done here?
In this paper, we present two mathematical programming approaches, one that yields upper and lower bounds on the moments of any stationary solution (on the \emph{stationary moments} for short), and another that provides bounds on arbitrary averages of the solutions (or \emph{stationary averages} for short). First, we build on our previous methodology\cite{Kuntz2016} that yields moment bounds of polynomial diffusions using semidefinite programming and adapt it to discrete reaction networks with polynomial or rational propensities. In particular, we consider the first few moment equations, which are generally not closed (that is, underdetermined), and constrain their solution space using semidefinite inequalities satisfied by the moments of any admissible probability distribution. Independently of our work, this same approach was recently proposed in\cite{Sakurai2017,Dowdy2018,Dowdy2017,Ghusinga2017a} for networks with polynomial propensities.

Secondly, we use the approximation schemes introduced in\cite{Kuntz2018a} with a state space truncation guided by a moment bound to obtain upper and lower bounds on arbitrary stationary averages (Theorem~\ref{bounds}). We prove that, in the case of a unique stationary solution, these bounds converge to the corresponding average as the truncation approaches the entire state space. By repeatedly applying this approach, we obtain bounds on the complete stationary solution and its marginals. We show that these bounds converge in total variation to the stationary solution and its marginals as the truncation tends to the entire state-space (Corollaries \ref{compstatcc} and \ref{compstatccm}). Mathematically, the method considers finite but underdetermined systems of stationary equations supplemented by inequalities that involve the moment bounds. Importantly, we show how to compute (resp. bound from above) the total variation distance between the lower (resp. upper) bounds and the exact solution, which quantifies the approximation error. Additionally, we explain how our LP approach doubles up as a test for the uniqueness of the stationary solutions (Corollary \ref{uniqueness-test}), a prerequisite of most other approximation schemes, and we explain how to use our scheme in the non-unique case to obtain converging approximations of the ergodic distributions (however, no computable errors are available for these approximations).

%We illustrate our mathematical programming approaches with applications to stochastic biochemical networks of interest in the literature. Our running example is Schl\"ogl's model, a birth-death process that exhibits a bifurcation and that allows for analytical manipulations. Building on these manipulations, we introduce the full approach applicable to networks that are not analytically tractable. We exemplify its use by computing to high precision the stationary distribution and marginals  of a toggle switch across different parameter regimes.
 
The paper is organised as follows. In Sec.~\ref{sec:stationary}, we introduce the stationary solutions of the CME. We then discuss the problem of bounding the stationary moments in Sec.~\ref{bmom}. We introduce our approach by demonstrating how analytical bounds for the first two moments can be obtained for a simple birth-death process, Sch\"ogl's model, and then develop its generalisation to multi-species networks with rational propensities. In Sec.~\ref{bdist}, we show how to bound entire stationary solutions, their averages and marginals. We first demonstrate the approach in the case of birth-death processes for which distribution bounds can be obtained analytically. Then, we explain how to generalise the approach to arbitrary CMEs. In Sec.~\ref{togglesec}, we exemplify the usefulness of our methods by computing the stationary solution to high verified accuracy of a toggle switch model. We conclude with a discussion of our results in Sec.~\ref{sec:discussion}. 
% For the sake of completeness, we have included Appendix \ref{appendix} an account of the theory pertinent to the stationary distributions of continuous-time chains and their relationship to the stationary solutions of CMEs. Although this theory is not necessary for the methods of this paper, it is crucial to link the bounds we compute to the long-term behaviour of the chains.
%
%
%
\section{Background}
\label{sec:stationary}

Stochastic biochemical kinetics is described by a set of $m$ reactions involving $n$ species $S_1,S_2,...,S_n$:
%
\begin{align}
\label{eq:network}
 R_j:\enskip v_{1j}^-S_1+\dots+v_{nj}^-S_n \xrightarrow{a_j} v_{1j}^+S_1+\dots+v_{nj}^+S_n \\ 
 \quad j=1,\ldots,m,  \nonumber
\end{align}
where $v_{ij}^{\pm}\in\n$ are the stoichiometric coefficients and $a_j:\nn\to[0,\infty)$ is the propensity of reaction $R_j$. 
Under well-mixed conditions, the state of the system is described by the %time-varying 
vector $X(t)=(X_1(t),\dots,X_n(t))$, which compiles the number of molecules of all species at time $t$. 
%
This stochastic process is commonly modelled by a minimal continuous-time Markov chain with rate matrix $Q=(q(x,y))$
%
\begin{equation}
\label{eq:qmatrixsrn}
q(x,y):=\left\{\begin{array}{ll}a_j(x)&\text{if } y=x+v_j,\\
-\sum_{j=1}^ma_j(x)&\text{if } y=x,
\\ 0&\text{otherwise},\end{array}\right.
\end{equation}
%% 
where $v_j:=(v_{1j}^+-v_{1j}^-,\dots,v_{nj}^+-v_{nj}^-)$ denotes the stoichiometric vector consisting of the net changes in molecule numbers produced by reaction $R_j$. 
% If reaction $R_j$ consumes (resp. produces) no molecules, we write $\varnothing$ in the LHS (resp. RHS) of \eqref{eq:network}. 

The state of the system takes values in a (possibly infinite) subset $\s$ of $\n^n$ known as the \emph{state space}. This set must be chosen such that 
%
\begin{align}
\label{eq:qmatrix}
%\begin{array} {c}
& q(x,y) \geq0 \enskip & \forall x\neq y,\\ 
& q(x)  :=  -q(x,x)=\sum_{y \in \s, \, 
y \neq x} q(x,y)<\infty\enskip &\forall x\in\cal{S},
%\end{array}
\end{align}
%
in which case the rate matrix is said to be \emph{totally stable} and \emph{conservative}. 
%In the above, the sum is taken over all states $y$ in $\s$ that are not $x$. 

% In practice, the chain $X$ is often constructed by running the Gillespie Algorithm\cite{Gillespie1976,Gillespie1977,Kendall1950,Doob1945,Doob1942,Feller1940} and it is defined up an \emph{explosion time} at which point the chain leaves the state space, see Appendix \ref{appendix} for details. In the case of a reaction network \eqref{eq:network}, an explosion occurs if and only if the molecule counts of one or more of the species diverges to infinity. If no such explosion occurs so that the chain is defined for all time regardless of its starting location (that is, the initial molecule counts for networks \eqref{eq:network}), then the matrix $Q$ is said to be \emph{regular}. 

If the chain cannot leave the state space in finite time, it is said to be \emph{non-explosive} and the matrix $Q$ is \emph{regular} (see App.~\ref{appendix}). In this case, the collection $p_t:=\{p_t(x)\}_{x\in\s}$ of probabilities $p_t(x)$ of observing the chain in state $x$ at time $t\geq0$ is the only solution of the \emph{chemical master equation} (\emph{CME})
%
\begin{align}\label{eq:CME}
\frac{\mathrm{d}p_t(x)}{\mathrm{d}t} = p_tQ(x),\quad p_0(x)=\lambda(x)\qquad\forall x\in\s,
\end{align}
%
where we use the notation $p_tQ(x):=\sum_{y\in\s}p_t(y)q(y,x)$. Any probability distribution $\pi:=\{\pi(x)\}_{x\in\s}$ %on $\s$ 
that solves the equation
%
\begin{equation}\label{eq:stat}\pi Q(x)=0\qquad\forall x\in\s,\end{equation}
%
is a \emph{stationary solution of the CME}. 
%
The set of all these solutions
%
%
\begin{equation}\label{eq:cmestat}\cal{P}:=\left\{\pi\in\r^{\mmag{\s}} : \begin{array}{ll}
   \pi \, Q(x)=0,& x\in\s,\\
 \pi(\s)= 1, \, \pi(x)\geq0, & \forall x\in\s
 \end{array}\right\};\end{equation}
 %
 %
forms a convex polytope in the $\mmag{\s}$-dimensional space $\r^{\mmag{\s}}$, where $\pi(\s):=\sum_{x\in\s}\pi(x)$, 
%``$\pi\geq0$'' is a shorthand for ``$\pi(x)\geq0$ for all $x\in\s$'', 
and $\mmag{\s}$ denotes the number of states in the state space $\s$ (possibly infinity). For many networks, these solutions determine the long-term behaviour of the chain (see App.~\ref{appendix}).
% We denote the $\pi$-\emph{average} of a real-valued function $f$ on $\s$ as
%
% $$\ave{f}:=\sum_{x\in\s}f(x)\pi(x).$$


% POSITIVE RECURRENT DOES NOT APPEAR IN THE MAIN TEXT.
% If the chain, regardless of its starting location, does not eventually leave the state space in one sense or another, the chain is said to be \emph{positive recurrent}\footnote{Meyn1993a} and the stationary solutions of the CME determine its long-term behaviour
% %in the sense that both its ensemble averages and time averages converge to one of the stationary solutions
% , see Appendix \ref{appendix} for details.  For instance, in the case of the reaction network \eqref{eq:network}, positive recurrence means that the molecule counts do not grow progressively larger and larger approaching infinity, something that is expected of good\footnote{We believe that care should be taken here: models are abstractions of real life and can be far from faithful representations thereof. An unstable model of stable phenomena is not something impossible: it is merely a potential indication of poor modelling choices. This is aggravated by the fact that many fitting algorithms involve automated parameter sweeps, and the algorithms can easily step through parameter values that yield unstable models.} models for most real-life networks. In practice, verifying whether a chain is positive recurrent is done by using a \emph{Foster-Lyapunov criterion} such as that the well-known one that we have included in Appendix \ref{appendix}: Theorem \ref{lyap}. 
%
%
\section{Bounding the stationary moments}
\label{bmom}
%
In this section, we present a semidefinite programming approach that yields bounds of increasing tightness on the stationary moments of reaction networks \eqref{eq:network} with polynomial or rational propensities. We introduce our approach for a model system and obtain analytically bounds on the first two moments in Sec.~\ref{sec:schloegl_moments}. We then present the general semidefinite programming method in Sec.~\ref{sec:mommulti}.
%
%
%
%
\begin{figure}
	\begin{center}
	\includegraphics[width=0.45\textwidth]{1.pdf}
	\vspace{-20pt}
	\end{center}
\caption{\textbf{Outer approximations and bounds for the moments of the stationary distribution of Schl\"ogl's model~\eqref{eq:smdl}. } 
\textbf{(a)} Grey area: The projection on the $y_1$--$y_2$ plane of the outer approximation $\cal{E}^3$ in \eqref{eq:E3_bounds} of the set of stationary moment vectors. The black dots correspond to the upper and lower bounds on the first and second moments obtained in \eqref{eq:e3finboun}.
%of $4$-dimensional vectors that satisfy the first moment equation~\eqref{eq:schmom} involving the first four moments ($z_0$, $z_1$, $z_2$, $z_3$), the positive semidefinitive inequalities~\eqref{eq:schmmm}, and the normalising constraint $z_0=1$. 
\textbf{(b)} By appending further moment equations and inequalities, we tighten the outer approximation~\eqref{eq:spectrahedron} (areas contained within the coloured outlines) of the singleton set of stationary moment vectors (black dot). The corresponding lower and upper bounds on the moments~\eqref{eq:powermoments2} (coloured dots) computed by solving the SDPs \eqref{eq:sdps} close in around the moment vector. Parameter values: $k_1=1$, $k_2=1$, $k_3=4/5$, and $k_4=1$.}
\label{fig:set}
\end{figure}
%
%
%
\subsection{An analytical example: moment bounds for Schl\"ogl's model} 
\label{sec:schloegl_moments}
%
Consider the classic autocatalytic network proposed by Schl\"ogl~\cite{schlogl1972} as a model for a chemical phase transition:
%
\begin{equation}
\label{eq:smdl}
2S\underset{a_2}{\overset{a_1}{\rightleftarrows}}3S,
\qquad \varnothing \underset{a_4}{\overset{a_3}{\rightleftarrows}}S. 
\end{equation}
%
Because the network comprises a single species $S$, the state space $\s$ is $\mathbb{N}$ and we assume that the propensities follow mass-action kinetics:
%
\begin{align}
& a_1(x):=k_1 x (x-1),
\quad a_2(x):= k_2 x (x-1) (x-2) \label{eq:schloegl_props1}\\ 
& a_3(x):= k_3, 
\quad  a_4(x):=k_4x,
\label{eq:schloegl_props2}
\end{align}
%
where $k_1, k_2, k_3, k_4>0$ are the rate constants. The CME has a unique stationary solution $\pi$ and all of its moments are finite, see App.~\ref{appendix}. The first \emph{stationary moment equation} (obtained by plugging $f(x):=x$ into \eqref{eq:steq} of the next section) reads
%
\begin{align}
\label{eq:schmom} 
&b_1z_0-b_2z_1+b_3z_2-b_4z_3 = 0, 
\end{align}
%
where  $b_1:=k_3$, $b_2:=k_1+2k_2+k_4$, $b_3:=k_1+3k_2$, $b_4:=k_2$ are all positive numbers and $z\in\r^4$ is the vector containing the first four moments of $\pi$, $z:=(\ave{1},\ave{x},\ave{x^2},\ave{x^3})$, where $\ave{f}:=\sum_{x\in\s}f(x)\pi(x)$ denotes the $\pi$-average of a real-valued function $f$ on $\s$. Even after noting that  
%
\begin{align}
\label{eq:pos_moments}
\ave{1}=1
\end{align}
%
because $\pi$ is a probability distribution, Eq.~\eqref{eq:schmom} is clearly underdetermined. 

Given any two-dimensional real vector $\mathbf{u}:=(\mathbf{u}_0, \mathbf{u}_1)\in\r^2$, define the affine polynomial  $u(x):=\mathbf{u}_0+\mathbf{u}_1x$ (we say that $\mathbf{u}$ is the \emph{vector of coefficients} of the polynomial $u$). Notice that $u^2$ and $v$ are  non-negative on $[0,\infty)$, where $v(x):=x(u(x))^2$. Because the stationary solution has support contained in $\s=\n\subseteq[0,\infty)$, we have that $\ave{u^2}\geq0$ and $\ave{v}\geq0$. We rewrite these inequalities as
%
\begin{align}\ave{u^2}&=\mathbf{u}_0^2\ave{1}+2\mathbf{u}_0\mathbf{u}_1\pi(x)+\mathbf{u}_1^2\ave{x^2}\nonumber\\
&=\mathbf{u}^TM_3^0(z)\mathbf{u}\geq 0,\label{eq:nfdyuhsanshua1}\\
\ave{v}&=\mathbf{u}_0^2\ave{x}+2\mathbf{u}_0\mathbf{u}_1\ave{x^2}+\mathbf{u}_1^2\ave{x^3}\nonumber\\
&=\mathbf{u}^TM_3^1(z)\mathbf{u}\geq 0,\label{eq:nfdyuhsanshua2}\end{align}
%
where we are employing vector notation and the matrices $M^3_0(y)$ and $M^3_1(y)$ are defined by
%
$$M_3^0(y)  := \begin{bmatrix}y_0 &y_1 \\ y_1 &y_2\end{bmatrix},\qquad 
M_3^1(y) := \begin{bmatrix}y_1 &y_2 \\y_2 &y_3\end{bmatrix},\qquad\forall y\in\r^4.$$
%
Because \eqref{eq:nfdyuhsanshua1}--\eqref{eq:nfdyuhsanshua2} hold for all $\mathbf{u}\in\r^2$, we have that $M^3_{0}(z)$ and $M^3_{1}(z)$ are two positive semidefinite (p.s.d.) matrices or
%
\begin{align}
\label{eq:schmmm}
M_3^0(z)\succeq0,\qquad M_3^1(z)\succeq0,
\end{align}
%
for short. Putting \eqref{eq:schmom}--\eqref{eq:schmmm} together, we have that $z$ belongs to the set
%
\begin{align}
\cal{E}^3=\left\{y\in\r^4:
\begin{array}{l}y_0=1, \\ 
%y_1\geq0, \quad y_2 \geq0,\quad y_3\geq0, \\
b_1y_0-b_2y_1+b_3y_2-b_4y_3 = 0, \\
M_3^0(y) \succeq 0, \quad M_3^1(y) \succeq 0.
\end{array}\right\}.
\label{eq:E3}
\end{align}
%
In other words, the singleton set $\{z\}$ of vectors  whose entries are composed of the first four moments of the stationary solutions of Sch\"ogl's model is contained in $\cal{E}^3$. For this reason, we say that $\cal{E}^3$ is an \emph{outer approximation} of the set of stationary moment vectors of the chain and
%
\begin{align}
l^3_\alpha:=\inf\{y_\alpha:y\in\cal{E}^3\}  & \le \ave{x^\alpha} \le \sup\{y_\alpha:y\in\cal{E}^3\}=:u^3_\alpha
 \label{eq:schloegl_moment_bounds}
 \end{align}
 %
 for all $\alpha=0,1,2,3$. To obtain explicit expressions of these bounds for the first two moments ($\alpha=1,2$), note that Sylvester's criterion implies that the semidefinite inequalities \eqref{eq:schmmm} are equivalent to
%
\begin{align}
y_0\geq0,\quad y_1\geq0, &\quad y_2\geq0,\quad y_3\geq0,\quad \label{eq:schoin1}\\
y_0y_2-y_1^2 &\geq0,  \label{eq:schoin2}  \\
y_1y_3-y_2^2 &\geq0. \label{eq:schoin3} 
\end{align}
%
Inequalities \eqref{eq:schoin1}--\eqref{eq:schoin2} tell us that the moments and the variance must be non-negative, and~\eqref{eq:schoin3} gives us an additional condition involving the first three moments. 

Because of the constraint $y_0=1$ and \eqref{eq:schoin2},  the inequalities $y_0\geq0$ and $y_2\geq0$ in \eqref{eq:schoin1} are redundant, and so we omit them in our description of $\cal{E}^3$. Since $y_0=1$, the moment equation in \eqref{eq:E3} is satisfied if and only if
%
\begin{equation}\label{eq:schoeqr}
b_1-b_2y_1+b_3y_2-b_4y_3 = 0.
\end{equation}
%
If $y_1>0$, then \eqref{eq:schoin3} implies that $y_3\geq0$, and if $y_1=0$ this follows from $y_2\geq0$ and \eqref{eq:schoeqr}. Thus, the inequality $y_3\geq0$ in \eqref{eq:schoin1} is also redundant and we omit it. 

Due to \eqref{eq:schoeqr}, the inequality \eqref{eq:schoin3} is holds if and only if
%
$$b_1y_1+b_3y_1y_2-b_4y_2^2-b_2y_1^2\geq0.$$
%
Suppose that $b_3\geq 2\sqrt{b_2b_4}$ (the other case can be dealt with similarly). Because $y_1\geq0$, the quadratic formula tells us that the above is satisfied if and only if 
% 
$$r_2^-(y_1)\leq y_2\leq r_2^+(y_1)$$
%
with
%
$$r_2^\pm(x):=\frac{b_3x\pm\sqrt{4b_1b_4x+(b_3^2-4b_2b_4)x^2}}{2b_4}.$$
%
In other words, we have reformulated $\cal{E}^3$ in \eqref{eq:E3} as:
%
\begin{align}
\cal{E}^3=\left\{y\in\r^4:
\begin{array}{l}y_0=1,\quad y_1\geq0, \\ 
y_3=(b_1-b_2y_1+b_3y_2)/b_4,\\ 
\max\{y_1^2,r_2^-(y_1)\}\leq y_2\leq r_2^+(y_1).
\end{array}\right\}.
\label{eq:E3_bounds}
\end{align}
%
Fig.~\ref{fig:set}(a) shows the projection of  $\cal{E}^3$ onto the $y_1$-$y_2$ plane for some sample parameter values.

The description \eqref{eq:E3_bounds} leaves clear that $(1,0,0,b_1/b_4) \in \cal{E}^3$. Because $y_1\geq 0$ and $y_2\geq0$ for every vector $y$ in $\cal{E}^3$, the lower bounds in \eqref{eq:schloegl_moment_bounds} for the first two moments are trivial: $l^3_1=l^3_2=0$. The upper bounds however are not: as $r^+(x)> r^-(x)$ for all $x>0$ and both $x\mapsto x^2$ and $r^+$ are monotonically increasing functions, the description \eqref{eq:E3_bounds} implies that the pair $(u^3_1,u^3_2)$ composed of the upper bounds in \eqref{eq:schloegl_moment_bounds} is the northeasternmost point of intersection of the curves $x\mapsto x^2$ and $x\mapsto r^+(x)$. That is, $(u^3_1,u^3_2)=(r_4,r_2^+(r_4))$ where $r_4$ is the rightmost root of the quartic polynomial $x\mapsto x(b_1-b_2x+b_3x^2-b_4x^3)$. It hence follows that
%
\begin{equation}\label{eq:e3finboun}0\leq \ave{x}\leq r_4,\qquad 0\leq \ave{x^2}\leq r^+(r_4).\end{equation}
%
which provides us with lower (albeit uninformative) and upper bounds on the first and second moments. As shown in Fig.~\ref{fig:set}(b), these bounds are rough; we will improve on them in the following section.
%
%
\subsection{Bounding the moments of reaction networks with rational propensities by solving semidefinite programs}
\label{sec:mommulti}
%
\begin{figure*}
	\begin{center}
	\includegraphics[width=\textwidth]{2.pdf}
	\vspace{-20pt}
	\end{center}
\caption{\textbf{The moment bounds for Schl\"ogl's model~\eqref{eq:smdl} converge.} \textbf{(a)} Upper ($u^d_\alpha$, red) and lower ($l^d_\alpha$, blue) bounds on the first three moments ($\alpha = 1, 2, 3$) given in \eqref{eq:powermoments2} computed for increasing number of moment equations $d$ and moment inequalities. The insets show that the difference between the upper and lower bounds decreases to zero as the number of moment equations and inequalities used is increased: both bounds converge to the moment (dashed lines). \textbf{(b)} Using the moment bounds in (a) we compute bounds of three summary statistics: coefficient of variation, variance, and skewness. The insets show the difference of lower and upper bounds on these summary statistics. The parameters for Schl\"ogl's model~\eqref{eq:smdl} used here are  $k_1=6$, $k_2=1/3$, $k_3=50$, $k_4=3$ corresponding to a unimodal stationary solutions, see Fig.~\ref{fig:schodist}.} 
\label{fig:schomom}
\end{figure*}
%
The approach of the previous section can be generalised to any reaction network \eqref{eq:network} with state space $\s\subseteq\nn$ as long as its propensities are rational functions:
%
\begin{equation}\label{eq:numden}a_j(x) := s_j(x)/o(x)\qquad\forall j=1,\dots,m,\end{equation}
%
where $s_1,\dots,s_m,o$ are polynomials on $\rn$ and the common denominator $o$ satisfies $o(x)>0$ for all $x\in\s$. To generalise the approach, we use the \emph{adjoint equation} obtained by multiplying both sides of \eqref{eq:stat} by $f(x)$, summing over $x\in\s$, and applying Fubini's Theorem:
%
\begin{equation}
\label{eq:steq}
\ave{Q f}= \sum_{x\in\s}(Qf)(x)\pi(x):= \sum_{x\in\s} \sum_{y\in\s}q(x,y)f(y)\pi(x)= 0,
\end{equation}
%
as long as $\ave{q\mmag{f}}=\sum_{x\in\s}q(x)\mmag{f(x)}\pi(x)$ is finite, where $\mmag{f(x)}$ denotes the absolute value of $f(x)$, see \cite[Prop. 3]{Glynn2008}.

The key observation is that letting $f$ in \eqref{eq:steq} be any monomial $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\nn$, we obtain a linear equation satisfied by the rational moments
%
\begin{equation}\label{eq:ch7vmom}z_\beta:=\ave{x^\beta/o},\end{equation}
%
of any stationary solution $\pi$ of the CME. We express the power moments as linear combinations of the rational moments: given any $d$ no smaller than the degree $d_o$ of the denominator $o$,
%
$$\ave{x^\alpha}=\sum_{\mmag{\beta}\leq d}\mathbf{o}_\beta \ave{\frac{x^{\alpha+\beta}}{o}}=\sum_{\mmag{\beta}\leq d}\mathbf{o}_\beta z_{\alpha+\beta},$$
%%
where $\mmag{\beta}:=\beta_1+\dots+\beta_n$, the above sums are taken over all $\beta\in\nn$ such that $\mmag{\beta}\leq d$, $z:=(z_\beta)_{\mmag{\beta}\leq d}$ is the vector of the rational moments \eqref{eq:ch7vmom}, and $\mathbf{o}:=(\mathbf{o}_\beta)_{\mmag{\beta}\leq d}$ is the vector of coefficients of the denominator so that $o(x)=\sum_{\mmag{\beta}\leq d}\mathbf{o}_\beta x^\beta$ (using the convention that $\mathbf{o}_\beta=0$ if $\mmag{\beta}>d_o$). As we show in Corollary \ref{boundcor} below, 
%
\begin{align}
 \label{eq:sdps}
 l^d_f : =\inf\{\mathbf{f}^Ty:y\in\cal{E}^d\} \quad \text{and} \quad
 u^d_f : =\sup\{\mathbf{f}^Ty:y\in\cal{E}^d\},
\end{align}
%
bound the \emph{$\alpha$-moment} $\ave{x^\alpha}$, where $\mathbf{f}$ denotes the vector of coefficients of $o(x)x^\alpha$, $\mathbf{f}^Ty:=\sum_{\mmag{\beta}\leq d}\mathbf{f}_\beta y_\beta$ denotes the inner product of the vectors $\mathbf{f}$ and $y$, $\cal{E}^d$ is the set of vectors in $\rhd$ that satisfy both the moment equations and certain positive semidefinite inequalities that we define in the following, and $\#_d:= \binom{n+d}{d}$ is the dimension of the vector of moments $z$. 
%
The set $\cal{E}^d$ is defined by linear equalities and semidefinite inequalities, hence it is a \emph{spectrahedron}\cite{Blekherman2013}. Both the left and right hand side of \eqref{eq:sdps} are SDPs, a class of convex optimisation problems that can be efficiently solved computationally using standard solvers. We flesh out the details in what follows.


\subsection*{The stationary moment equations}
%
Let $f$ be any monomial $x^\alpha$ and define
%
\begin{align*}g(x)&:=o(x)Q f(x)=\sum_{j=1}^m s_j(x)((x+v_j)^\alpha-x^\alpha),\\
h(x)&:=o(x)q(x)f(x)=\sum_{j=1}^m s_j(x)x^\alpha.\end{align*}
%
The polynomials $g$ and $h$ have degrees of at most $\mmag{\alpha}+d_a-1$ and $\mmag{\alpha}+d_a$ respectively, where $d_a$ is the maximum degree of the numerators $s_1,\dots,s_m$ in \eqref{eq:numden}. For this reason, if $\mmag{\alpha}$ is no greater than $d-d_a+1$ and $\pi$ is a stationary solution with finite rational moments of order $d+1$ (meaning that \eqref{eq:ch7vmom} is finite if $\mmag{\beta}\leq d+1$), then $\ave{q\mmag{f}}$ is finite and the adjoint equation \eqref{eq:steq} tells us that
%
\begin{align}\nonumber0=&\ave{Qf}=\ave{\frac{g}{o}}=\ave{\frac{\sum_{\mmag{\beta}\leq d} \mathbf{g}_\beta x^\beta}{o}}\\
&=\sum_{\mmag{\beta}\leq d} \mathbf{g}_\beta \ave{\frac{x^\beta}{o}}=\sum_{\mmag{\beta}\leq d} \mathbf{g}_\beta z_\beta=z^T\mathbf{g}.\label{eq:momeq}\end{align}
%
We call \eqref{eq:momeq} the $\alpha$\emph{-moment equation}. By the \emph{moment equations}, we mean the collection of these equations for all $\alpha$ such that $\mmag{\alpha}\leq d-d_\alpha+1$:
%
\begin{equation}\label{eq:momeqs}z^T\mathbf{oQx^\alpha}=0\qquad\forall \alpha:\mmag{\alpha}\leq d-d_\alpha+1,\end{equation}
%
where we are now writing $oQx^\alpha$ directly instead of $g$ to make the exponent $\alpha$ explicit. For instance, in the case of Schl\"ogl's model \eqref{eq:smdl} with propensities \eqref{eq:schloegl_props1}--\eqref{eq:schloegl_props2}, $o(x)=1$, $d_a=3$, and $\#_d=d+1$, and \eqref{eq:schmom} is \eqref{eq:momeq} with $\alpha=1$.
%
\subsubsection*{The moment constraints}

Given that $\pi$ is a probability distribution, we also have the equation
%
\begin{equation}\label{eq:normalisation_cond}z^T\mathbf{o}=\sum_{\mmag{\beta}\leq d}\mathbf{o}_\beta \ave{\frac{x^\beta}{o}}=\ave{\frac{o}{o}}=\ave{1}=1.\end{equation}
%
%
In general, the system of linear equations formed by \eqref{eq:momeqs}--\eqref{eq:normalisation_cond} is underdetermined. However, just as we did for Schl\"ogl's model in Sec.~\ref{sec:schloegl_moments}, we constrain its solution space by appending to it well-known semidefinite inequalities\cite{Lasserre2009,Blekherman2013,Kuntzthe} satisfied by the rational moments of any unsigned measure (e.g., $\pi$) on a subset of $\n^n$. In our case, we have that
%
\begin{align}
 \label{eq:moment_matrices_posdef}
 %M^0_d(y) \succeq 0  \quad 
 M^i_d(z) \succeq 0\qquad\forall i=0,\dots,n,
 \end{align}
%
where the \emph{localising matrices}\cite{Lasserre2009} $M^0_d(y),\dots,M^n_d(y)$ are defined by
% 
\begin{align*}
 %\label{eq:sdpin}
 [M^0_d(y)]_{\alpha\beta}&:=y_{\alpha+\beta},     & \forall \alpha,\beta:\mmag{\alpha},\mmag{\beta}\leq \lfloor d/2\rfloor,\notag\\
 [M^i_d(y)]_{\alpha\beta}&:=y_{\alpha+\beta+e_i}, & \forall \alpha,\beta:\mmag{\alpha},\mmag{\beta}\leq \lfloor(d-1)/2\rfloor,
\end{align*}
%
with $e_i$ denoting the $i^{th}$ unit vector, $\lfloor b\rfloor$ the largest integer no greater than $b\in\r$, and $y$ any vector in $\rhd$. For Schl\"ogl's model, $M^0_3(y)$ and $M^1_3(y)$ are the matrices in \eqref{eq:schmmm}. The inequalities \eqref{eq:moment_matrices_posdef} follow\cite{Lasserre2009,Blekherman2013,Kuntzthe} from the fact that for any polynomial $g$ of degree $\lfloor d/2\rfloor$ with vector of coefficients $(\mathbf{g}_\beta)_{\mmag{\beta}\leq d}$, it holds that
%
\begin{align}
\mathbf{g}^TM^0_d(z)\mathbf{g}&=\sum_{\mmag{\alpha}\leq\lfloor d/2\rfloor}\sum_{\mmag{\beta}\leq\lfloor d/2\rfloor}\mathbf{g}_\alpha \mathbf{g}_\beta z_{\alpha+\beta}\nonumber\\
&=\ave{\frac{\left(\sum_{\mg{\alpha}\leq\lfloor d/2\rfloor}\mathbf{g}_\alpha x^{\alpha}\right)\left(
\sum_{\mmag{\beta}\leq\lfloor d/2\rfloor}\mathbf{g}_\beta x^{\alpha}\right)}{o}}\nonumber\\
&=\ave{\frac{g^2}{o}}\geq0,\label{eq:semial1}
\end{align}
%
since $o(x)>0$ for all $x\in\s$. Similarly, it can be shown\cite{Lasserre2009,Blekherman2013} that
%
\begin{equation}\label{eq:semial2}\mathbf{g}^TM^i_d(z)\mathbf{g}=\ave{\frac{x_ig^2}{o}}\geq0\quad\forall i=1,\dots,n,\end{equation}
%
because $x_i\geq0$ for all $x\in\s\subseteq\n^n$. As inequalities \eqref{eq:semial1}--\eqref{eq:semial2} hold for all vectors of coefficients $\mathbf{g}$, we have that the corresponding matrices are p.s.d., i.e. \eqref{eq:moment_matrices_posdef}. 
%
\subsubsection*{Bounding the moments via semidefinite programming}
%
The main result of this section is that the spectrahedron $\cal{E}^d$ contains the set of moment vectors (up to order $d$) of the stationary solutions. For this reason, $\cal{E}^d$ is known\cite{Lasserre2009,Kuntz2016} as an \emph{outer approximation} of the set of these stationary moment vectors.
%
\begin{theorem}
\label{inclu}
Let $\cal{P}^{d}$ be the set of stationary solutions of the CME with finite rational moments of order up to $d$: 
%
$$\cal{P}^{d}:=\{\pi\in\cal{P}:z_\beta<\infty\quad\forall\beta:\mmag{\beta}\leq d\},$$
%
in the notation of \eqref{eq:cmestat} and \eqref{eq:ch7vmom}. Assuming that $d$ is no smaller than the degree $d_o$ of the denominator $o$ in \eqref{eq:numden}, for any $\pi \in \cal{P}^{d+1}$, the vector of rational moments $z=(z_\beta)_{\mmag{\beta}\leq d}$ of $\pi$ defined by
%
\begin{equation}\label{eq:momvec}z_\beta:=\ave{x^\beta/o}\qquad \forall\beta:\mmag{\beta}\leq d,\end{equation}
%
belongs to the spectrahedron
%
\begin{align}
 \label{eq:spectrahedron}
 \cal{E}^d:=\left\{y\in \rhd:
 \begin{array}{l} 
  z^T\mathbf{o}=1, \\
 z^T\mathbf{oQx^\alpha}=0\quad\forall \mmag{\alpha}\leq d-d_\alpha+1, \\ 
 M_d^i(y)\succeq0\quad \forall i=0,1,\dots,n.
 \end{array}
 \right\},
\end{align}
%
\end{theorem}
%%
\begin{proof}This follows directly from \eqref{eq:momeqs}--\eqref{eq:moment_matrices_posdef}.
\end{proof}
%

The outer approximation property of $\cal{E}^d$ implies that we obtain bounds on the stationary moments by optimising over $\cal{E}^d$. In particular, by finding the vector in $\cal{E}^d$ with largest (smallest) $\alpha$-entry, we obtain an upper (resp. lower) bound on the $\alpha$-moment.  For example, Fig.~\ref{fig:set}(b) shows the projection onto the $y_1$--$y_2$ plane of the outer approximation $\mathcal{E}^d$ for Schl\"ogl's model with increasing $d$. The northeasternmost (resp. soutwesternmost) vector of these outer approximation yield upper (resp. lower) bounds on the first two moments of Schl\"ogl's model. The moment matrix $M_d^i(y)$ is a principal submatrix of $M_{d+1}^i(y)$ and a matrix is p.s.d. if and only if all of its principal submatrices are p.s.d. For these reasons, these outer approximations \emph{tighten} around the set of stationary moment vectors in the sense that every vector in $\cal{E}^{d+1}$ (appropriately truncated) belongs to $\cal{E}^d$, see again Fig.~\ref{fig:set}(b). Thus, we obtain bounds of increasing quality by raising $d$:
%
\begin{corollary}\label{boundcor}
%
For any polynomial $f$ of degree no greater than $d$ and $\pi\in\cal{P}^{d+1}$ with $\cal{P}^{d+1}$ as in Theorem \ref{inclu},
%
\begin{equation}\label{eq:powermoments}l^d_f \leq\ave{\frac{f}{o}}\leq u^d_f\end{equation}
%
where $l^d_f$ and $u^d_f$ are as in \eqref{eq:sdps}
%

%
Moreover, these bounds form monotonic sequences:
%
\begin{equation}\label{monotone}
 l^d_f \leq l^{d+1}_f\leq \dots  \leq\ave{\frac{f}{o}}\leq \dots\leq u^{d+1}_f  \leq u^d_f.
\end{equation}
%
\end{corollary}
%
\begin{proof}Inequalities \eqref{eq:powermoments} follow from Theorem \ref{inclu} by noting that
%
$$\ave{\frac{f}{o}}=\ave{\frac{\sum_{\mmag{\beta}\leq d}\mathbf{f}_\beta x^\beta}{o}}=\sum_{\mmag{\beta}\leq d}\mathbf{f}_\beta\ave{\frac{x^\beta}{o}}=\mathbf{f}^Tz,$$
%
with $z$ as in \eqref{eq:momvec}. The monotonicity in \eqref{monotone} follows from the definition of $\cal{E}^d$ in \eqref{eq:spectrahedron} and the fact mentioned in the main text that a matrix is p.s.d. if and only if all of its principal submatrices are p.s.d.
\end{proof}

Setting $f(x):=o(x) x^\alpha$ in \eqref{eq:sdps} with $|\alpha|\leq d-d_o$, Corollary \ref{boundcor} tells us that $l_\alpha^{\mmag{\alpha}+d_o},\dots,l_\alpha^{d}$ (resp. $u_\alpha^{\mmag{\alpha}+d_o},\dots,u_\alpha^{d}$) is a monotonically increasing (resp. decreasing) sequence of lower (resp. upper) bounds on the $\alpha$-moment $\ave{x^\alpha}$, where
%
\begin{equation}\label{eq:powermoments2}l_\alpha^d:=l^d_f,\qquad u^d_\alpha:=u^d_f\end{equation}
%
In Sec.~\ref{sec:schloegl_moments}, we derived the \emph{moment bounds} $l^3_1$, $l^3_2$, $u^3_1$, and $u^3_2$, for Schl\"ogl's model through ad hoc analytical manipulations. Instead, we can systematically carry out this bounding procedure numerically:   computing the bounds $l^d_\alpha$ or $u^d_\alpha$ consists of solving the SDPs \eqref{eq:sdps}, see Fig.~\ref{fig:schomom}. As we show in the figure, the moment bounds can be combined to obtain tight bounds on other commonly used statistics such as the variance, coefficient of variation, and skewness.
%
\subsubsection*{Moment bounds may not converge}
%
%
% Computing both lower and upper bounds allows us to quantify the uncertainty remaining in our knowledge of the moment: $\ave{x^\alpha}$ is no more than $u^d_\alpha-l^d_\alpha$ away from $l^d_\alpha$ or $u^d_\alpha$ (or any number in between). If this difference is too large, then \eqref{monotone} tells us we can increase $d$ and obtain sharper (or, at least, no worse) bounds. Suppose that all stationary solutions of the CME have all moments finite. 
An important question is whether we can make the upper and lower bounds arbitrarily close to true value of the moment by increasing the number of moment equations $d$ and size of the corresponding semidefinite constraints (i.e., whether the bounds \emph{converge}). In our experience (e.g., Fig.~\ref{fig:schomom}) the bounds often converge, see also the examples in \cite{Ghusinga2017a,Dowdy2017,Dowdy2018,Sakurai2017}. However, no general guarantee is available for the following reasons.

When the stationary solution is not unique, the lower bounds are limited by the stationary distribution with the smallest moment while the upper bounds are limited by the largest moment. When the stationary solution is unique, the bounds may not converge because the semidefinite conditions we employ are not tailored to distributions with supports on discrete state spaces, but instead support measures defined on the non-negative reals (this can be seen from \eqref{eq:semial1}--\eqref{eq:semial2}). Even in the case of measures with support on $\rp^n$, these bounds are not guaranteed to converge to the moments of such a measure \cite{Lasserre2009,Laurent2009,Landau1987}. An illuminating example demonstrating this phenomenon can be found in\cite[Sec.~IV B]{Dowdy2018}. 

Although more stringent moment conditions tailored to measures on a discrete set $\s$ exist\cite{Lasserre2002a,Lasserre2009}, they are computationally involved (see\cite{Dowdy2018}). However, sharp conditions on the distributions are trivial to construct (e.g., $\pi(x)\geq0$ for all $x$ in $\s$ and $\pi(\s)=1$). The linear programming scheme that is the subject of the remainder of this paper builds on this simple observation.
%
%
\subsubsection*{Computational implementation and numerical considerations}\label{num1}
%
As discussed in\cite{Kuntz2016,Fantuzzi2016a}, the SDPs associated with moment problems tend to be ill-conditioned. The origin of this numerical instability remains an open problem. Our computations indicate that instabilities could be the result of the rapid growth of moments. This leads to moment matrices with elements of vastly different magnitude which can be problematic for the standard double-precision SDP solvers. We mitigated this issue by using the multi-precision solver SDPA-GMP\cite{Nakata2010} as in~\cite{Fantuzzi2016a}. Another approach to deal with this issue is to scale the moments \cite{Kuntz2016,Fantuzzi2016a,Dowdy2018}. Lastly, we should mention that recently developed scalable solvers \cite{Zheng2017} could be adapted to the SDPs considered here.

% as in ~\cite{Fantuzzi2016a} when tackling similar SDPs,  instead of standard double-precision SDP solvers. Another approach to deal with this issue proposed is to scale the decision variables \cite{Kuntz2016,Fantuzzi2016a,Dowdy2018}. Lastly, we should mention that first-order solvers have recently reported success in solving SDPs\cite{Zheng2017} that are dual\cite{Blekherman2013,Laurent2009,Lasserre2009} to ours.


% Alternatively, rewriting a network of polynomial propensities in terms of rational moments (with a judicious choice of an artificial denominator $o$) can also ameliorate the numerical instabilities by increasing the size of the semidefinite constraints without increasing the gap in the orders of magnitude of the moments. 


\vspace{5pt}
\noindent\textbf{Implementation details:} To set up the SDPs \eqref{eq:sdps}, we used MATLAB R2017a and the modelling package YALMIP\cite{Lofberg2004}. We solved the SDPs  with the multi-precision solver SDPA-GMP\cite{Nakata2010} in conjunction with the interface mpYALMIP\cite{Fantuzzi2016}. Alternatively, one may choose to use the modelling package GloptiPoly 3\cite{Henrion2009} and any other SDP solver.
%
\section{Bounding the stationary solutions of the CME}\label{bdist}
%
In this section, we introduce our linear programming scheme that yields bounds on the stationary solutions of CMEs, their marginals, and their averages. Assuming that the CME has a unique stationary solution, we show how iterating our scheme yields converging bounds on the complete solution and its marginals (Corollaries \ref{compstatcc} and \ref{compstatccm}). In the non-unique case, we explain how to adapt our approach to obtain approximations of the ergodic distributions. Additionally, we explain how to use our scheme as a uniqueness test (Corollary \ref{uniqueness-test}). 

The semidefinite programming scheme discussed in the previous section entails constructing outer approximations of the set of stationary moments and optimising over these approximations. Instead, here we construct outer approximations directly of the set \eqref{eq:cmestat} of stationary solutions of the CME and optimise over these other approximations. The approximations involve linear inequalities instead of semidefinite ones like in the previous section. For this reason, the optimisation problems that we must solve here are linear programs instead of semidefinite ones, a simpler subclass of convex optimisation problems for which very mature solvers, e.g.\cite{CPLEX}, are available. We begin in Sec.~\ref{sec:bd} by introducing our approach in the simple setting of the birth-death processes and we generalise it in Sec.~\ref{lpsec} to arbitrary CMEs.
%
\subsection{Semi-analytical bounds for birth-death processes}\label{sec:bd}
\begin{figure*}
	\begin{center}
	\includegraphics[width=0.7\textwidth]{3.pdf}
	\vspace{-20pt}
	\end{center}
\caption{\textbf{Bounding the probability distribution of Schl\"ogl's model. }
%
%
%
%
% 
The stationary distribution of Schl\"ogl's model Eqs.~\eqref{eq:smdl}--\eqref{eq:schloegl_props2} can be unimodal (a,c) or bimodal (b,d).
\textbf{(a)} Shaded areas illustrate the gap between upper and lower bounds on the stationary solution using the truncations $\s_r$ in \eqref{eq:momtrunc} that consist of all states $x$ with $\alpha$-powers $x^\alpha$ smaller than $r$. An upper bound on the first moment ($\alpha=1$) was computed using the approach of Sec.~\ref{sec:mommulti}. The exact unimodal distribution, Eqs.~(\ref{eq:1deg}) and (\ref{eq:schloegl_p0}), is shown for comparison (black line). 
\textbf{(b)} Upper and lower bounds are shown for the bimodal case (red line: scaled by a factor of $1/3$ and computed with $r=52$ and $\alpha=1$; blue: $r=112$, $\alpha=1$; yellow: $r=127$, $\alpha=25$) are compared with the exact solution (black line). 
\textbf{(c)} Approximation error $\varepsilon_r$ of lower bounds as in \eqref{eq:Markov1} is shown for the unimodal case. The approximation error of upper bounds is at most $\varepsilon_r$. The error decreases with truncation size for various values of $\alpha$. Note that our scheme only provides information about the stationary distribution if the error is less than $1$ (dashed line). Moment bounds on the $\alpha^{th}$ moment are computed using the first $25$ moment equations. 
\textbf{(d)} Same as in (c) for the bimodal case.
% 
Parameter values in (a) and (c) are $k_1=6$, $k_2=1/3$, $k_3=50$, $k_4=3$ and the upper bounds on the moments are $u^{25}_1=17.5$ and $u^{25}_{25}=5.79\times 10^{35}$. In (b) and (d), we used $k_1=1/9$, $k_2=1/1215$, $k_3=27/2$, $k_4=59/20$ with upper bounds $u^{25}_1=98.0$ and $u^{25}_{25}=6.37\times 10^{51}$.
%
% 
% given in \eqref{eq:lumeasdef} are com using the truncation \eqref{eq:momtrunc} (with $\alpha=1$ and $r$ values in legend) and 
% \textbf{(b)} Bounds for different parameters yielding a bimodal distribution computed using truncation \eqref{eq:momtrunc} and  $(r,\alpha)=(52,1)$ (red, scaled by a factor of $1/3$), $(r,\alpha)=(112,1)$ (blue), and $(r,\alpha)=(127,25)$ (yellow). Exact distribution (black). %The lower bounds were computed using the upper bound on the $1^{st}$ moment for $r=52$ and $112$ and that on the $25^{th}$ moment for $r=127$.
% %
% \textbf{(b)} Tail bound $\varepsilon_r$ in \eqref{eq:Markov1}, as a function of truncation size, computed using bounds $c=u_1^{25},\dots u^{25}_{25}$ on the first twenty-fifth moments. %The total error computed using the $\ell^1$-distance between the exact solution and the upper bound (yellow) is decreasing with increasing truncation $r$. The total error of the lower bounds obtained using a bound on the $1^{st}$ moment (red) or one on the $25^{th}$ moment (blue line) was computed using \eqref{eq:low}. 
} 
\label{fig:schodist}
\end{figure*}

Birth-death processes are reaction networks \eqref{eq:network} with state space $\s=\n$ whose value changes by $\pm1$ at a time:
%
\begin{equation}
\label{eq:one_step}
\varnothing \xrightarrow{a_+}S\xrightarrow{a_-}\varnothing. 
\end{equation}
The stationary equations $\pi Q=0$ read
% 
\begin{align}\label{eq:introeqs1} &a_-(1)\pi(1)-a_+(0)\pi(0)=0,\\ &a_-(x+1)\pi(x+1)-(a_+(x)+a_-(x))\pi(x)\nonumber\\&+a_+(x-1)\pi(x-1)=0,\label{eq:introeqs2}\end{align}
%
for all $x=1,2\dots$ Assuming that 
%
\begin{equation}\label{eq:drpos}a_-(x)>0\qquad\forall x=1,2,\dots,\end{equation} 
%
these equations have the unique solution \cite{Gardiner2009}:
%
\begin{align}
\label{eq:1deg}
&\pi(x) = \left [ \prod_{z=1}^x\frac{a_+(z-1)}{a_-(z)} \right ] \pi(0) =: \gamma(x)  \, \pi(0),
\end{align}
for $x =1,2,\dots$. The normalising condition allows to solve for $\pi(0)$:
\begin{align}
 \label{eq:normalising}
&\pi(0) = \frac{1}{\sum_{x=0}^\infty \gamma(x)},
\end{align}
%
where $\gamma(0):=1$. Consequently, birth-death processes have at most one stationary solution and it exists if and only if the sum in \eqref{eq:normalising} is finite. 


% $\pi$ belongs to the convex polytope \eqref{eq:cmestat} if and only if
%
\subsubsection*{Upper bounds}
%
For general birth-death processes, no closed-form expression for $\pi(0)$ in \eqref{eq:normalising} is known and, consequently, the solution cannot be computed exactly. However, we can compute an upper bound on $\pi(0)$ by truncating the sum in \eqref{eq:normalising}. In the following, we pick the \emph{state space truncation} (or \emph{truncation} for short) 
%
\begin{equation}\label{eq:momtrunc}\s_r:=\{x\in\n:x^\alpha<r\}.\end{equation}
%
We then have
%
%
\begin{align} 
 & \pi(x)\leq \frac{\gamma(x)}{\sum_{z\in\s_r} \gamma(z)}=:u^r_x, 
  \label{eq:upperbounds}
\end{align}
  %
for all $x\in\s_r$ because $({\sum_{x\in\s_r} \gamma(x)})^{-1}$ bounds $\pi(0)$ from above. 
%
\subsubsection*{Lower bounds}
%
Less obvious, however, is how to obtain lower bounds on $\pi(x)$. To this end, we use a bound on the $\alpha$-moment $\ave{x^\alpha}$ obtained with the SDP approach of Sec.~\ref{bmom}. Markov's inequality then provides us with a bound on the mass $m_r$ outside of the truncation:
%
\begin{align}
 m_r :&= \sum_{x\not\in\s_r} \pi(x) \leq  \frac{\ave{x}^\alpha}{r}\leq \frac{c}{r}=:\varepsilon_r.\label{eq:Markov1}
\end{align}
%
We say that $\varepsilon_r$ is a \emph{tail bound}. A lower bound on $\pi(0)$ then follows from \eqref{eq:1deg}--\eqref{eq:Markov1}:
%
\begin{align}
\pi(0)=& \frac{1-m_r}{\sum_{x\in\s_r}\gamma(x)}\geq \frac{1- \varepsilon_r}{\sum_{x\in\s_r}\gamma(x)},
\end{align}
and hence for all $x\in\s_r$ we have
\begin{align}
& \pi(x)\geq u^r_x (1-\varepsilon_r)=:l^r_x,
\label{eq:lowerbounds}
\end{align}
which is the desired lower bound on $\pi(x)$.
%
%
\subsubsection*{Approximation errors and error bounds}
%
We have shown that
%
\begin{align}
l^r_x  \leq \pi(x) \leq u^r_x, \quad  \forall x \in\s_r.
\label{eq:upper_lower_bounds}
\end{align}
%
Furthermore, it follows from \eqref{eq:1deg}--\eqref{eq:normalising}, that $u^r_x $ converges to $\pi(x)$ as $r$ approaches infinity. Because the tail bound $\varepsilon_r$ tends to zero as $r$ tends to infinity, the lower bound $l^r_x$ also converges to $\pi(x)$. These facts motivate us to approximate $\pi$ with the measures $l^r:=\{l^r(x)\}_{x\in\s}$ and $u^r:=\{u^r(x)\}_{x\in\s}$ obtained by padding these bounds with zeros:
%
\begin{equation}\label{eq:lumeasdef}l^r(x):=\left\{\begin{array}{ll}l^r_{x}&\text{if }x \in\s_r\\0&\text{if }x \not\in\s_r\end{array}\right.,\quad u^r(x):=\left\{\begin{array}{ll}u^r_{x}&\text{if }x \in\s_r\\0&\text{if }x \not\in\s_r\end{array}\right..\end{equation}
%
To quantify the error between $\pi$ and such an approximation $\pi^r$, we use the total variation norm
%
\begin{equation}\label{eq:tvnorm} \norm{\pi-\pi^r}=\sup_{A\subseteq\s}\mmag{\sum_{x\in A}\pi(x)-\sum_{x\in A}\pi^r(x)},\end{equation}
%
so that the mass 
%
\begin{equation}\pi(A):=\sum_{x\in A}\pi(x)\end{equation}
%
that $\pi$ gives to any event $A\subseteq\s$ is the mass $\pi^r(A)$ that the approximation $\pi^r$ awards to this event plus/minus, at most, $\norm{\pi-\pi^r}$. We refer to $\norm{\pi-\pi^r}$ as the \emph{approximation error of} $\pi^r$. Using  \eqref{eq:upperbounds}--\eqref{eq:lumeasdef}, we have that
%
\begin{align}
\norm{\pi-l^r}&=\sup_{A\subseteq\s}\left(\sum_{x\in A}\left(\pi(x)-l^r(x)\right)\right)=\pi(\s)-l^r(\s)\nonumber\\
&=1-\sum_{x\in\s_r}l^r_x=1-(1-\varepsilon_r)\sum_{x\in\s_r}u^r_x=\varepsilon_r,\label{eq:low}\\
\norm{\pi-u^r}&=\max\left\{\sum_{x\in\s_r}(u^r(x)-\pi(x)),\sum_{x\not\in\s_r}\pi(x)\right\} \nonumber\\
&=\max\{1-(1-m_r),m_r\}=m_r\leq\varepsilon_r.\label{eq:up}
\end{align}

Because $\varepsilon_r=c/r$ tends to zero as $r$ tends to infinity, \eqref{eq:low}--\eqref{eq:up} shows that approximations $l^r$ and $u^r$ converge in total variation to $\pi$. We summarise these findings in the following theorem.
%
\begin{theorem}\label{birthdeathth}Consider any birth-death process \eqref{eq:one_step} satisfying \eqref{eq:drpos} and with finite sum \eqref{eq:normalising}. Let $\pi$ be the unique stationary solution of the CME given by \eqref{eq:1deg}--\eqref{eq:normalising}. 
%
\begin{enumerate}[label=(\roman*)]
\item Suppose that $c$ is a constant and $\s_1\subseteq\s_2\subseteq\dots$ is a sequence of increasing truncations such that $m_r\leq c/r=:\varepsilon_r$ holds for all $r=1,2,\dots$, where $m_r$ is as in \eqref{eq:Markov1}. The sequences $\{l^r\}_{r\in\zp}$ and $\{u^r\}_{r\in\zp}$ of approximations defined by \eqref{eq:upperbounds}, \eqref{eq:lowerbounds}, and \eqref{eq:lumeasdef} converge in total variation to $\pi$. Furthermore, the approximations are composed of bounds that vary monotonically with $r$:
%
$$l^r(x)\leq l^{r+1}(x)\leq \dots\leq\pi(x)\leq \dots \leq u^{r+1}(x)\leq u^r(x),$$
%
for all $x\in\s_r$. The approximation error $\norm{\pi-l^r}$ of the lower bounds is $\varepsilon_r$, while that $\norm{\pi-u^r}$ of the upper bounds is $m_r$.  
%
\item If $\pi$ has finite $d+1$-order moment ($\ave{x^{d+1}}<\infty$) and $\alpha\in\{1,\dots,d\}$, then $m_r\leq c/r$ holds for all $r=1,2,\dots$ if $c:=u^d_\alpha$ and $\s_r$ is as in \eqref{eq:momtrunc}.  
\end{enumerate}
\end{theorem}
%
\begin{proof}This follows from \eqref{eq:upperbounds}--\eqref{eq:up} and Corollary \ref{boundcor}.
\end{proof}
%%
\subsubsection*{Application to Schl\"ogl's model}
%
We apply the above semi-analytical theory to compute bounds on the distributions of Schl\"ogl's model \eqref{eq:smdl}. The model presents an ideal test case because it can be solved analytically and its stationary distribution is unique. Defining the birth and death rates $a_+:=a_1+a_3$ and $a_-:=a_2+a_4$ as in \eqref{eq:schloegl_props1}--\eqref{eq:schloegl_props2}, the normalising constant in Eq.~\eqref{eq:normalising} can be computed explicitly 
%
\begin{align}
 \label{eq:schloegl_p0}
 \frac{1}{\pi(0)}= {_2F_2}\left(-\frac{c_1+1}{2},\frac{c_1-1}{2};-\frac{c_2+1}{2},\frac{c_2-1}{2};\frac{k_1}{k_2}\right),
\end{align}
where ${_2F_2}$ denotes the generalised hypergeometric function and $c_1:=\sqrt{1-4 {k_3}/{k_1}}$ and $c_2:= \sqrt{1-4 {k_4}/{k_2}}$. The analytical solution can be unimodal (black line in Fig.~\ref{fig:schodist}a) or bimodal (black line in b) depending on the parameter values. 

We compare our approximations $l^r$ and $u^r$ given in Theorem \ref{birthdeathth} (shades in Fig.~\ref{fig:schodist}a,b) to the analytical solution. As the value of the truncation $r$ is increased, the approximations become tighter. In the unimodal case, the approximation error decreases significantly below $1$ when the truncation cutoff $r$ is larger than the mode (Fig.~\ref{fig:schodist}c). When the truncation is sufficiently large, employing bounds on moments of higher order provides tighter tail bounds and smaller approximation errors (lines in Fig.~\ref{fig:schodist}c). On the other hand, when the truncation is of the same order as the modal value or smaller, 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.~\ref{fig:schodist}d), but the approximation error only decreases when the truncation cutoff $r$ is larger than the second mode. This shows that a large proportion of the distribution may be missed for smaller values of $r$. The choice of the truncation $r$ thus needs to be guided by error bounds.

%
\subsection{Generalisation to arbitrary chains: bounding the stationary distributions via linear programming}\label{lpsec}
%
%
The truncation method of the previous section relies on the special structure of birth-death processes and the fact that the stationary solution can be expressed recursively in terms of $\pi(0)$. For general networks, such a semi-analytic solution cannot be obtained and a fully numerical is needed. We first outline this numerical procedure for the birth-death process and then generalise the concepts introduced.

In the case of a birth-death process, we consider a finite set of the stationary equations that feature only states in the truncated state space $\s_r:=\{0,1,\dots,r-1\}$. These equations are underdetermined since the system of equations \eqref{eq:introeqs1}--\eqref{eq:introeqs2} for $x<r$ involves $r-1$ equations with $r$ unknowns. If $\varepsilon_r$ is a tail bound as in \eqref{eq:Markov1}, the solution satisfies
%
\begin{equation}\label{eq:intromass}
1-\varepsilon_r\leq \pi(\s_r)\leq 1.
\end{equation}
%
For this reason, the \emph{restriction}
%
\begin{equation}\label{eq:restr}\pi_{|r}(x):=\left\{\begin{array}{ll}\pi(x)&\text{if }x \in\s_r\\0&\text{if }x \not\in\s_r\end{array}\right.\qquad\forall x\in\s\end{equation}
%
of $\pi$ to $\s_r$  belongs the convex polytope
%
\begin{equation}
\label{eq:polytopesch}
\left\{\pi^r\in\r^{\mmag{\s}}: \begin{array}{l}\pi^r\text{ satisfies \eqref{eq:introeqs1}, \eqref{eq:introeqs2} with $x<r-1$,}\\\text{\eqref{eq:intromass}, }\pi^r(\s_r^c)=0,\text{ and }\pi^r\geq0\end{array}\right\},
\end{equation}
%
where $\s_r^c$ denotes the set of states $\{r,r+1,\dots\}$ outside of the truncation $\s_r$ and the remainder of the notation is analogous to that in \eqref{eq:cmestat}. For this reason, we have that the stationary solution is bounded as
%
$$\inf\{\pi^r(x):\pi^r\in\text{\eqref{eq:polytopesch}}\}\leq \pi(x)\leq \sup\{\pi^r(x):\pi^r\in\text{\eqref{eq:polytopesch}}\},$$ 
%
for each $x$ in $\s_r$. Because the definition \eqref{eq:polytopesch} only involves linear equations and inequalities, computing each of these bounds numerically consists of solving a linear program. Using  \eqref{eq:1deg}, one can verify that the infimum and supremum are given by the expressions we found using the semi-analytic approach of the previous section, Eqs.~\eqref{eq:lowerbounds} and \eqref{eq:upperbounds}. 
% Note that any $\pi$ in \eqref{eq:polytopesch} satisfies $\pi^r(x)=0$ for all $x\not\in\s_r$. Thus, the feasible set is a  convex polytope in $r$-dimensional space defined by $r-1$ linear equalities and $r+2$ inequalities and computing numerically this infimum and supremum consists of solving two $r$-dimensional LPs. Using \eqref{eq:1deg}, it is straightforward to verify that the infimum (resp. supremum) is $l^r_x$ (resp. $u^r_x$) defined in \eqref{eq:lowerbounds} (resp. \eqref{eq:upperbounds}) with $\alpha$ in \eqref{eq:Markov1} set to $1$.


In the general case, determining the upper and lower bounds above is an optimisation problem that can be solved using linear programming. For this purpose, let $Q$ be a rate matrix satisfying \eqref{eq:qmatrix} and consider the convex polytope 
 %
\begin{equation}\label{eq:brp}\cal{P}_{w,c}:=\left\{\pi\in\r^{\mmag{\s}} : \begin{array}{l}
\pi Q(x)=0, x\in\s,\\
\ave{w}\leq c,\enskip\pi(\s) = 1,\enskip\pi\geq0
\end{array}\right\},\end{equation}
 %
of stationary solutions  that satisfy the \emph{moment bound} 
%
\begin{equation}\label{eq:momboundw}\ave{w}\leq c,\end{equation}
%
where $w$ is a non-negative function and $c$ a constant. We are taking liberties with the name ``moment bound'' as we do not assume in the following that $w$ is a polynomial or even that the state space $\s$ is a subset of $\rn$.  We define our truncations to be the sublevel sets of $w$:
%
\begin{equation}\label{sublevel}\s_r:=\{x\in\s:w(x)<r\}\end{equation}
%
For these truncations to be of use in practice, they must be finite (for all $r\geq0$), in which case we say that $w$ is \emph{norm-like}. 

In the definition of the convex polytope analogous to \eqref{eq:polytopesch}, we include all of the stationary equations that only involve the states in the truncation $\s_r$: $\pi Q(x)=0$ for each $x$ that belongs to 
%
\begin{equation}\label{eq:nr}\cal{N}_r:=\{x\in\s_r:\{z\in\s:q(z,x)\neq 0\}\subseteq\s_r\}.\end{equation}
%
%For instance, in \eqref{eq:upperbounds}--\eqref{eq:up} of the previous section, we used $w(x)=x$, $\s_r=\{0,1,\dots,r-1\}$, and $\cal{N}_r= \{0,1,\dots,r-2\}$. 
For networks \eqref{eq:network} with rate matrix \eqref{eq:qmatrixsrn} we have that $x$ belongs to $\cal{N}_r$ if and only if $x-v_j$ belongs to $\s_r$ for every $j$ such that $x-v_j$ belongs to $\s$ and $a_j(x-v_j)>0$.

A bound on the mass of the tail $m_r$ of any $\pi\in\cal{P}_{w,c}$ follows from Markov's inequality and the moment bound \eqref{eq:momboundw}:
%
\begin{equation}\label{eq:tailboundw}m_r:=\sum_{x\not\in\s_r}\pi(x) \leq \frac{1}{r}\sum_{x\not\in\s_r}w(x)\pi(x)\leq\frac{\ave{w}}{r}\leq \frac{c}{r}.\end{equation}
%
Similarly as with \eqref{eq:polytopesch}, putting the above together we have that, for any $\pi$ in $\cal{P}_{w,c}$, its restriction $\pi_{|r}$ to $\s_r$  in \eqref{eq:restr} belongs to the convex polytope
%
\begin{equation}\label{eq:br}\cal{P}^r_{w,c}:= \left\{\pi^r\in\r^{\mmag{\s}} : \begin{array}{l}
   \pi^r Q(x)=0, \enskip x\in\cal{N}_r,\\
1-\varepsilon_r\leq \pi^r(\s)\leq 1,\\
\aver{w}\leq c,\enskip \pi^r(\s_r^c)=0,\enskip\pi\geq0.
\end{array}\right\}.\end{equation}
%
For this reason, we say that $\cal{P}^r_{w,c}$ is an \emph{outer approximation} of the set $\cal{P}_{w,c}$ of stationary solutions of the CME that satisfy the moment bound \eqref{eq:momboundw}. 
%
\begin{lemma}[Outer approximations of $\cal{P}_{w,c}$]\label{outer} If $\pi$ belongs to $\cal{P}_{w,c}$ in \eqref{eq:brp}, then its restriction $\pi_{|r}$ in \eqref{eq:restr} belongs to $\cal{P}^r_{w,c}$ in \eqref{eq:br}.\end{lemma}
%
\begin{proof} This follows directly from \eqref{eq:tailboundw} and the fact that \eqref{eq:nr} implies that $\pi Q(x)=\pi_{|r} Q(x)$ for each $x$ in $\cal{N}_r$.
\end{proof}


Before discussing how to use the outer approximation to compute bounds on the whole distribution, we first consider the simpler problem of bounding individual stationary averages and probabilities.

%
\subsubsection{Bounding stationary averages and probabilities}\label{sec:piaver}
%
Consider the stationary  average of a real-valued function $f$ on $\s$:
%
\begin{align}
 \ave{f}=\sum_{\s} f(x) \pi(x).
\end{align}
%
If $f$ is the indicator function of a subset $A$ of the state space, the average is the probability $\pi(A):=\sum_{x\in A}\pi(x)$ given to $A$. If $A$ contains only a single state $x$, the average is the stationary probability $\pi(x)$. If $f(x)$ is $x^\alpha$, the average is the $\alpha$- moment. 

Because $\cal{P}^r_{w,c}$ is an outer approximation of $\cal{P}_{w,c}$ (Lemma \ref{outer}), we obtain bounds on an average $\ave{f}$ (with $\pi\in\cal{P}_{w,c}$) by optimising over $\cal{P}^r_{w,c}$. In particular, define
%
\begin{align}\label{eq:lowupthe4}l^r_f:=\inf\{\aver{f}:\pi^r\in\cal{P}^r_{w,c}\},
\notag\\ u^r_f:=\sup\{\aver{f}:\pi^r\in\cal{P}^r_{w,c}\}.\end{align}
%
Since the constraints in \eqref{eq:br} imply that $\pi^r(x)=0$ for all $x$ outside of  $\s_r$ and $\pi^r$ in $\cal{P}^r_{w,c}$, computing $l^r_f$ and $u^r_f$ numerically consists of solving two LPs with $\mmag{\s_r}$ variables, $\mmag{\cal{N}_r}$ equality constraints, and $\mmag{\s_r}+3$ inequality constraints.

Because the restriction $\pi_{|r}$ in \eqref{eq:restr} belongs to $\cal{P}^r_{w,c}$ for any $\pi$ in $\cal{P}_{w,c}$ (Lemma \ref{outer}),  $l^r_f$ and $u^r_f$ bound $\averr{f}$. We use either of these bounds (or any number in between them) as an approximation of the stationary average $\ave{f}$.  This approxmation converges as the truncations approach the entire state space:

\begin{theorem}\label{bounds} Suppose that $\cal{P}_{w,c}$ is non-empty, that $f$ is a real-valued function on $\s$, that $w$ is a norm-like function on $\s$, that 
%
\begin{equation}\label{eq:pqfin}\{z\in\s:q(z,x)\neq 0\}\text{ is finite}\quad\forall x\in\s,\end{equation}
%
and let $l^r_f$ and $u^r_f$ be as in \eqref{eq:lowupthe4}. For any $\pi$ in $\cal{P}_{w,c}$, we have that
%
$$l^r_f\leq\averr{f}\leq u^r_f,\qquad\forall r=1,2,\dots$$
%
Consequently, if $f$ is $\pi$-integrable, we have that
%
\begin{equation}\label{eq:generalaveragecontrollederrorofffunctionsinwspacewithnormlikew}l^r_f-c\left(\sup_{x\not\in\s_r}\frac{\mmag{f(x)}}{w(x)}\right)\leq \ave{f}\leq u^r_f+c\left(\sup_{x\not\in\s_r}\frac{\mmag{f(x)}}{w(x)}\right).\end{equation}
%
Furthermore, if $f$ is non-negative (resp. non-positive) on $\s_r^c$, then 
%
$$l^r_f\leq \ave{f}\qquad\text{ (resp. }\ave{f}\leq u^r_f\text{),}\qquad\forall r=1,2,\dots$$
%
If $f$ eventually grows strictly slower than $w$:
%
\begin{equation}\label{eq:weakstar2}\lim_{r\to\infty}\sup_{x\not\in\s_r}\frac{\mmag{f(x)}}{w(x)}=0,\end{equation}
%
then
%
\begin{align*}&\lim_{r\to\infty}l^r_f=l_f:=\inf\{\ave{f}:\pi\in\cal{P}_{w,c}\},\\
&\lim_{r\to\infty}u^r_f=u_f:=\sup\{\ave{f}:\pi\in\cal{P}_{w,c}\}.\end{align*}
%
%
%
\end{theorem}
%
\begin{proof} The bounds on $\aver{f}$ follow directly from Lemma \ref{outer}. Inequalities \eqref{eq:generalaveragecontrollederrorofffunctionsinwspacewithnormlikew} are then a consequence of the moment bound \eqref{eq:momboundw} and of the following generalisation of Markov's inequality
%
\begin{align*}\sum_{x\not\in\s_r}\rho(x)\mmag{f(x)}&\leq \left(\sup_{x\not\in\s_r}\frac{\mmag{f(x)}}{w(x)}\right)\sum_{x\not\in\s_r}w(x)\rho(x)\\
&\leq\left(\sup_{x\not\in\s_r}\frac{f(x)}{w(x)}\right)\rho(w).\end{align*}
%
The convergence follows from \cite[Theorem 3.2]{Kuntz2018a}.
\end{proof}
%
Condition \eqref{eq:pqfin} means that the chain can reach any given state in a single jump from only finitely many others. For instance, networks \eqref{eq:network} with $Q$-matrices \eqref{eq:qmatrixsrn} satisfy this condition as the chain can only reach state $x$ in a single jump from $x-v_1,\dots,x-v_m$. Note that Eq.~\eqref{eq:generalaveragecontrollederrorofffunctionsinwspacewithnormlikew} allows us to quantify the quality of our approximation of the stationary average $\ave{f}$. 

In the case that the CME has a unique stationary solution $\pi$ satisfying the moment bound \eqref{eq:momboundw} (i.e., $\cal{P}_{w,c}=\{\pi\}$), the feasible points \eqref{eq:br} are good approximations of the stationary solution, meaning that they converge to $\pi$ in \emph{weak*}: 
%
\begin{equation}\label{eq:weakstar}\lim_{r\to\infty}\aver{f}=\ave{f},\end{equation}
%
for each function $f$ satisfying \eqref{eq:weakstar2}. As explained in \cite[Remark 5.10]{Kuntzthe},  weak* convergence implies convergence in total variation if $w$ is norm-like.
%
\begin{corollary}\label{optimalpointsconv}
 Suppose that $w$ is norm-like, that \eqref{eq:pqfin} holds, that $f$ satisfies \eqref{eq:weakstar2} and that there exists a unique $\pi^*\in\cal{P}_{w,c}$ such that
%
$$\langle f\rangle_{\pi^*}=l_f\qquad\text{(resp. }\langle f\rangle_{\pi^*}=u_f\text{)}.  $$
%
For each $r=1,2,\dots$, let $\pi^r$ be an optimal point in $\cal{P}^r_{w,c}$ such that $\aver{f}=l^r_f$ (resp. $\aver{f}=u^r_f$). Then the optimal points $\pi^r$ converge to $\pi^*$ in weak* (and, thus, in total variation) as $r$ tends to infinity. 
\end{corollary}

\begin{proof}Follows from \cite[Corollary 3.3]{Kuntz2018a}.
\end{proof}

The optimal points $\pi^r$ in Corollary \ref{optimalpointsconv} exist because non-emptiness of $\cal{P}_{w,c}$ and the constraints $\pi^r\geq0$, $\pi^r(\s_r^c)=0$, and $\pi^r(\s)\leq 1$ imply that the LPs in \eqref{eq:lowupthe4} consist of optimising a continuous function over a compact non-empty subset of $\r^{\mmag{\s_r}}$. Indeed, LP solvers not only return the optimal value $l^r_f$ (or $u^r_f$), but also an optimal point $\pi^r$ that achieves it. In the case of a unique $\pi$, we have found the optimal points of the LP $\sup\{\aver{1}:\pi^r\in\cal{P}^r_{w,c}\}$ to be particular good approximations of $\pi$. However, we know of no practical way of quantifying the approximation error $\norm{\pi-\pi^r}$. In what follows, we compute upper and lower bounds on $\pi(x)$ for each $x$ in $\s_r$ and combine these to obtain approximations with quantified error bounds. 
%
%
%
\subsubsection{Unique case: bounding the stationary solution}
\label{sec:distributions}
%
Suppose that there is a unique stationary solution $\pi$ of the CME satisfying the moment bound \eqref{eq:momboundw} (i.e, $\cal{P}_{w,c}=\{\pi\}$). By  letting $f$ in \eqref{eq:lowupthe4} be the indicator function $1_x$ ($1_x(y)$ is one if $y=x$ and zero otherwise) of a state $x$ in the truncation $\s_r$  and solving the corresponding LPs 
%
\begin{align}l_x^r&:=\inf\{\pi^r(x):\pi^r\in\cal{P}_{w,c}^r\}\quad\forall x\in\s_r,\label{eq:lowerd}\\
 u^r_x&:=\sup\{\pi^r(x):\pi^r\in\cal{P}_{w,c}^r\}\quad\forall x\in\s_r,\label{eq:upperd}\end{align}
%
we obtain an upper and a lower bound on $\pi(x)$. We use the collections $l^r:=\{l^r(x)\}_{x\in\s}$ and $u^r:=\{u^r(x)\}_{x\in\s}$ of these bounds padded with zeros \eqref{eq:lumeasdef} as approximations of $\pi$. 
%
\begin{corollary}\label{compstatcc}Suppose that $\cal{P}_{w,c}$ is non-empty, that $w$ is a norm-like function on $\s$, that \eqref{eq:pqfin} is satisfied, and that $l^r$ and $u^r$ are as in \eqref{eq:lumeasdef} and \eqref{eq:lowerd}--\eqref{eq:upperd}.
%
\begin{enumerate}[label=(\roman*)]
\item For any $\pi\in\cal{P}_{w,c}$ and $r=1,2\dots$, the approximations $l^r$ and $u^r$ bound $\pi$:
%
$$l^r(x)\leq \pi(x)\quad\forall x\in\s,\quad \pi(x)\leq u^r(x)\quad\forall x\in\s_r,$$
%
and the approximation errors can be quantified as follows:
%
\begin{align}
\norm{l^r-\pi}&=1-l^r(\s_r)=:\varepsilon^l_r,\label{eq:epl}\\
%
\norm{u^r-\pi}&=\max\{u^r(\s_r)-1+m_r,m_r\}\nonumber\\
&\leq \max\{u^r(\s_r)-1+c/r,c/r\}=:\varepsilon^u_r,\label{eq:epu}
\end{align}
%
where $\norm{\cdot}$ denotes the total variation norm in \eqref{eq:tvnorm}.
\item Suppose that $\cal{P}_{w,c}$ is the singleton $\{\pi\}$. The approximation $u^r$ converges pointwise to $\pi$:
%
$$\lim_{r\to\infty}u^r(x)=\pi(x),\qquad\forall x\in\s.$$
%
The approximation $l^r$ converges in weak* (and, thus, in total variation) to $\pi$ as $r$ approaches infinity. 
\end{enumerate}
\end{corollary}
%
\begin{proof}The bounds of $(i)$ follows immediately from Theorem \ref{bounds}. The error expression in \eqref{eq:epl} follows from the bounds and the fact that the total variation norm of an unsigned measure is its mass. Similarly, \eqref{eq:epu} follows from Markov's inequality \eqref{eq:tailboundw} and the fact that the bounds imply that
%
\begin{align*}&\mmag{u^r(A)-\pi(A)}%%=\mmag{u^r(A)-\pi_{|r}(A)-\pi(A\cap\s_r^c)}\\&
\leq \max\{u^r(A)-\pi_{|r}(A),\pi(A\cap\s_r^c)\}\\
&\leq \max\{u^r(\s_r)-1+m_r,m_r\}\\
&=\max\{\mmag{u^r(\s_r)-\pi(\s_r)},\mmag{u^r(\s_r^c)-\pi(\s_r^c)}\},\quad\forall A\subseteq\s.\end{align*}
%
The convergence of the upper bounds in $(ii)$ follows from Theorem \ref{bounds} while that of the lower bounds follows from \cite[Theorem 4.1]{Kuntz2018a}.
\end{proof}
%
In other words, for sufficiently large $r$, $l^r$ and $u^r$ are close to $\pi$. In contrast with the feasible points $\pi^r$ discussed at the of end of the previous section,  we can answer the question ``is $r$ sufficiently large?'' by first computing $l^r$ (or $u^r$) and then evaluating the error using Corollary \ref{compstatcc}~$(i)$. Furthermore, the computable error $\varepsilon^l_r$ of $l^r$ converges to zero as $r$ tends to infinity and so, by increasing $r$, we will always find an approximation $l^r$ that verifiably meets any given error tolerance. We have no proof that the error bound $\varepsilon^u_r$ of the approximation $u^r$ converges to zero (nor that $u^r$ itself converges to $\pi$ in total variation), however we have not yet encountered an example in practice where we did not also observe this convergence for $u^r$ and its error bound.
%
%%
%
%
%
\subsubsection{Unique case: bounding the marginal distributions}
\label{sec:mdist}
%
For high-dimensional state spaces we are often interested in marginal distributions rather than the full stationary solution $\pi$. In particular, let $\{A_i\}_{i\in\cal{I}}$ be a collection of disjoint subsets of $\s$ indexed by $\cal{I}$ such that $\cup_{i\in\cal{I}} A_i=\s$. Consider the \emph{marginal} $\hat{\pi}$ of $\pi$ defined by:
%
\begin{equation}\label{eq:marginal2}\hat{\pi}(i)=\pi(A_i),\qquad \forall i\in\cal{I}.\end{equation}
%
Because the sets $A_1,A_2,\dots$ are disjoint and their union is the entire state space, $\hat{\pi}$ is a probability distribution on $\cal{I}$. For instance, in the case of a  reaction network \eqref{eq:network} with state space $\nn$, if
%
\begin{equation}\label{eq:kthm}A_i:=\n^{k-1}\times\{i\}\times \n^{n-k}\qquad \forall i\in\cal{I}:=\n,
\end{equation}
%
%\begin{array}{ll}A_i:=\{&(x_1,\dots,x_{k-1},i,x_{k+1},\dots,x_n):\\
%&(x_1,\dots,x_{k-1},x_{k+1},\dots,x_n)\in\n^{n-1}\}\end{array}
% for all $i\in\cal{I}:=\n$,
%
then $\hat{\pi}$ is the distribution describing the molecule counts of the $k^{th}$ species. We obtain approximations of the marginals analogous to those  of the entire distribution in \eqref{eq:lumeasdef} and \eqref{eq:lowerd}--\eqref{eq:upperd}:
%
\begin{equation}\label{eq:lumeasdefm}\hat{l}^r(i):=\left\{\begin{array}{ll}\hat{l}^r_{i}&\text{if }i \in\cal{I}_r\\0&\text{if }i \not\in\cal{I}_r\end{array}\right.,\quad \hat{u}^r(i):=\left\{\begin{array}{ll}\hat{u}^r_{i}&\text{if }i \in\cal{I}_r\\0&\text{if }i \not\in\cal{I}_r\end{array}\right..\end{equation}
%
%
where $\cal{I}_r$ is the finite subset of $i\in\cal{I}$ such that $A_i\cap\s_r$ is non-empty and
%
\begin{align}\hat{l}_i^r&:=\inf\{\pi^r(A_i):\pi^r\in\cal{P}_{w,c}^r\}\quad\,\,\forall i\in\cal{I}_r,\label{eq:lowerdm}\\
 \hat{u}^r_i&:=\sup\{\pi^r(A_i):\pi^r\in\cal{P}_{w,c}^r\}\quad\forall i\in\cal{I}_r,\label{eq:upperdm}\end{align}
%
Corollary \ref{compstatcc} holds almost identically for $\hat{l}^r$ and $\hat{u}^r$:
%
\begin{corollary}\label{compstatccm}Suppose that $\cal{P}_{w,c}$ is non-empty, that $w$ is a norm-like function on $\s$, that \eqref{eq:pqfin} is satisfied, that $\{A_i\}_{i\in\cal{I}}$ is a collection of disjoints sets of $\s$ such that $\cup_{i\in\cal{I}} A_i=\s$, and that $\hat{l}^r$ and $\hat{u}^r$ are as in \eqref{eq:lumeasdefm}--\eqref{eq:upperdm}.
%
\begin{enumerate}[label=(\roman*)]
\item For any $\pi\in\cal{P}_{w,c}$ and $r=1,2\dots$, the approximations $\hat{l}^r$ bounds from below $\hat{\pi}$ (defined in \eqref{eq:marginal2}):
%
$$\hat{l}^r(i)\leq \hat{\pi}(i)\quad\forall i\in\cal{I},$$
%
and the approximations satisfy
%
\begin{align*}
||\hat{l}^r-\hat{\pi}||&=1-\hat{l}^r(\cal{I}_r)=:\hat{\varepsilon}^l_r,\\
%
\norm{\hat{u}^r-\hat{\pi}}&\leq \max\{\hat{u}^r(\cal{I}_r)-1+m_r,m_r\}\\
&\leq \max\{\hat{u}^r(\cal{I}_r)-1+c/r,c/r\}=:\hat{\varepsilon}^u_r.\end{align*}
%
\item  Suppose that $\cal{P}_{w,c}$ is the singleton $\{\pi\}$. The approximation $\hat{u}^r$ converges pointwise to $\hat{\pi}$:
%
$$\lim_{r\to\infty}\hat{u}^r(i)=\hat{\pi}(i),\qquad\forall i\in\cal{I}.$$
%
The approximation $\hat{l}^r$ converges in total variation to $\hat{\pi}$ as $r$ approaches infinity.
\end{enumerate}
\end{corollary}
%
\begin{proof}This proof is analogous to that of Corollary \ref{compstatcc}.\end{proof}
%
The difference between Corollary \ref{compstatcc} and \ref{compstatccm} is that $\hat{u}^r$ is not necessarily an upper bound on $\hat{\pi}$, while  $u^r$ is an upper bound on $\pi$. The reason is that in \eqref{eq:marginal2} we are often marginalising over states not included in the truncated space $\s_r$ (e.g., in the case of \eqref{eq:kthm}). The probability mass of these states is bounded by the mass of the tail $m_r$ and hence by $c/r$, see \eqref{eq:tailboundw}. For this reason, if an upper bound is crucial, then we can instead use $\hat{\pi}(i)\leq \hat{u}^r(i)+c/r$.

\subsubsection{A uniqueness test and the non-unique case: approximating the ergodic distributions}\label{nonuniquesec}\label{sec:nonunique}
%
As Theorem \ref{bounds} and Corollaries \ref{compstatcc}--\ref{compstatccm} show, our LP approach yields bounds on the stationary solutions, their averages, and marginals, even if there is more than one of these solutions (i.e., if $\cal{P}_{w,c}$
is not a singleton). In general, however, there will be a non-trivial gap between the lower bounds and the upper bounds as the former will be limited by the stationary solution that achieves the minimum (over $\cal{P}_{w,c}$), while the latter will be limited by the solution that achieves the maximum. This simple observation enables us to turn our LP approach into a test of uniqueness, as follows.

The state space can be decomposed using the \emph{Doeblin Decomposition}
%
\begin{equation}\label{eq:ddc}\s=\left(\bigcup_{j}\cal{C}_j\right)\cup \cal{T},\end{equation}
%
where $\cal{C}_1,\cal{C}_2,\dots$ are the \emph{closed communicating classes} and $\cal{T}$ is the \emph{transient set}. A set $\cal{C}\subseteq\s$ is said to be a closed communicating class if the chain can transit between any pair of states in $\cal{C}$ but cannot leave $\cal{C}$. In terms of the rate matrix, this is the case if and only if given any $x,y\in\cal{C}$ we can find $x_1,x_2,\dots,x_l$ such that
%
\begin{equation}\label{eq:irred}q(x,x_1)q(x_1,x_2)\dots q(x_{l-1},x_l)q(x_l,y)>0,\end{equation}
%
but no such sequence exists if $y$ lies outside of $\cal{C}$\cite{Norris1997}. %If $\cal{C}=\s$, then the chain and its rate matrix are said to be \emph{irreducible}. 
For example, the state space of the network
%
$$ \emptyset\xrightarrow{a_1} 2S_1\xrightarrow{a_2}\emptyset,\qquad S_2\xrightarrow{a_3}\emptyset,$$
%
with mass action-kinetics (e.g., $a_1(x)=1$, $a_2(x)=x_1(x_1-1)$, and $a_3(x)=x_2$) decomposes into the sets
%
\begin{align*}\s=\n^2&=\{(x_1,0):x_1\in\n\text{ is odd}\}&=\cal{C}_1\\
&\cup\{(x_1,0):x_1\in\n\text{ is even}\}&=\cal{C}_2\\
&\cup\{(x_1,x_2):x_1\in\n, x_2\in\zp\}&=\cal{T}.\end{align*}
%The stationary distributions on $\cal{C}_1$ and $\cal{C}_2$ determine the long term behaviour of the chain and are called ergodic distributions.


For the reasons discussed in App.~\ref{appendix}, the Doeblin Decomposition is intimately related to the stationary distributions:
%
\begin{theorem}\label{doeblinc} Let $\cal{C}_1,\cal{C}_2,\dots$ be as in \eqref{eq:ddc} and suppose that $Q$ is regular.
\begin{enumerate}[label=(\roman*)]
\item For each class $\cal{C}_j$, there is at most one stationary solution $\pi_j$ such that $\pi_j(\cal{C}_j)=1$. Whenever it exists, $\pi_j$ is known as the ergodic distribution associated with $\cal{C}_j$.
\item Let $\cal{J}$ be the set of indexes $j$ of the communicating classes $\cal{C}_j$ that have an ergodic distribution $\pi_j$. The set of stationary solutions $\cal{P}$ in \eqref{eq:cmestat} is the set of convex combinations of the ergodic distributions:
%
%\emph{\textrm{conv}}\{\pi_i:i\in\cal{J}\}:=
%
\begin{equation}\label{eq:statconv}\cal{P}=\left\{\sum_{j\in\cal{J}}\theta_j\pi_j: \theta_j\geq0\quad \forall j\in\cal{J}, \quad\sum_{j\in\cal{J}}\theta_j=1\right\}.\end{equation}
\end{enumerate} 
\end{theorem}
%
\begin{proof}Given Theorem \ref{Qstateq} in the App.~\ref{appendix}, $(i)$ and $(ii)$ are shown in Ref~\citep[Theorem~2.44]{Kuntzthe}.\end{proof}
%
\begin{figure*}
	\begin{center}
	\includegraphics[width=\textwidth]{4.pdf} 
	\vspace{-20pt}
	\end{center}
\caption{\textbf{Bounds on the stationary probability of the toggle switch. }  \textbf{(a)}  Lower and upper bounds $l^r$ and $u^r$ exhibit qualitatively similar bimodal distribution shapes. The state truncation parameter was chosen to be $r=46^6$, so that $529$ equations involving $575$ states were included. The approximation error $\varepsilon^l_r$ on the lower bound $l^r$ (see \eqref{eq:epl}) was $0.34$, while the bound on the approximation error $\varepsilon^u_r$ of $u^r$ (see \eqref{eq:epu}) was $0.4$. White areas denote states $x$ outside of our truncation (that is, such that $w(x)\geq 46^6$).
% 
% Errors: $1-1^Tl^{46}=0.3355$, $1^Tu^{46}-1=0.3505$, and $\hat{u}^{10}_6/46^6=0.0537$. 
Computation time was $2.2$ minutes using the solver CPLEX V12.6.3\cite{CPLEX} on a desktop computer with a 3.5GHz processor and 16GB of RAM. 
\textbf{(b)} Computing the difference between these bounds reveals that the largest uncertainties occur near the distribution modes. %Note the difference in the values of the colourmap with those in the colourmap in (a).
% 
\textbf{(c)} The state truncation parameter was increased to $r=76^6$, so that $1444$ equations involving $1520$ states were included. Upper and lower bounds are visually indistinguishable. The error $\varepsilon^l_r$ of $l^r$ was $3.0\times10^{-3}$, while the bound on the approximation error $\varepsilon^u_r$ of $u^r$ was $2.7\times10^{-3}$.  Computation time was $22.6$ minutes on a desktop computer. 
%
\textbf{(d)} Computing the difference between these bounds reveals that these bounds are close to the true distribution with the maximum absolute difference being as small as $10^{-9}$.
% Errors: $1-1^Tl^{76}=2.985E-3$, $1^Tu^{76}-1=1.0332E-4$, and $\hat{u}^{10}_6/76^6=2.641E-3$. 
 Parameters are given in the main text.}
\label{fig:toggle_joint}
\end{figure*}
%
Theorem~\ref{doeblinc}$(ii)$ tells us that the ergodic distributions are the extreme points of the convex polytope $\cal{P}$~\eqref{eq:cmestat} of stationary solutions of the CME. Theorem~\ref{doeblinc}$(i)$ states that these extreme points (vectors in $\mmag{\s}$-dimensional space) are orthogonal to each other (with respect to the usual inner product on $\r^{\mmag{\s}}$). Because these vectors lie in the $\mmag{\s}$-dimensional non-negative orthant, the orthogonality implies that, in the non-unique case, each face of the non-negative orthant contains one of these extreme points. Our uniqueness test follows from this and from the fact that $\cal{P}^r_{w,c}$ is an outer approximation of $\cal{P}_{w,c}$.
%
\begin{corollary}[The uniqueness test]\label{uniqueness-test}Suppose that the premise of Corollary \ref{compstatcc} is satisfied and that $Q$ is regular. If $l^r(x)>0$ for some $x$ in $\s$ and $r=1,2,\dots$, then $\cal{P}_{w,c}$ is a singleton. Conversely, if $\cal{P}_{w,c}$ is the singleton $\{\pi\}$, and $\pi(x)>0$, then $l^r(x)>0$ for sufficiently large $r$.
\end{corollary}
%
\begin{proof} Suppose that $\cal{P}^r_{w,c}$ is not empty or a singleton. Theorem~\ref{doeblinc}$(ii)$ implies that $\cal{P}^r_{w,c}$ contains two ergodic distributions, say $\pi_1$ and $\pi_2$, associated with closed communicating classes $\cal{C}_1$ and $\cal{C}_2$. As these classes are disjoint, $x$ does not belong to one of them, say $\cal{C}_1$. But then Theorem~\ref{doeblinc}$(i)$ implies that $l^r(x)>\pi_1(x)=0$ contradicting the lower bound property of $l^r$ (Corollary~\ref{compstatcc}$(i)$), and hence showing that $\cal{P}^r_{w,c}$ must be empty or a singleton. The remainder follows immediately from the convergence of the bounds in Corollary~\ref{compstatcc}$(ii)$.
\end{proof}

In the non-unique case, $l^r$ or $u^r$ do not provide good approximations of any stationary distribution---Corollary~\ref{uniqueness-test} shows that $l^r$ is trivially zero everywhere, whereas Theorems~\ref{bounds} and~\ref{doeblinc}$(ii)$ show that, for large $r$, the mass of $u^r$ will be no smaller than the number of ergodic distributions. However, Corollary \ref{optimalpointsconv} tells us that any optimal point $\pi^r$ in $\cal{P}^r_{w,c}$ satisfying $\pi^r(x)=u^r(x)$ for any one state $x$ inside the truncation and a closed communicating class $\cal{C}_j$ is a good approximation of the ergodic distribution $\pi_j$ for sufficiently high $r$ (as $\pi^r$ converges to $\pi_j$ in weak$^*$). Importantly, we do not need to know a priori what (if any) class $x$ belongs to. Additionally, by looking at which states $\pi^r$ associates a non-zero probability, we get an idea of what the class $\cal{C}_j$ may be. Once this class is known, and if computable errors are important, we replace $\s$ with $\cal{C}_j$ and proceed as we did in the previous two sections for the unique case.
%
\subsubsection{Choosing the truncation}\label{apriori}
%
%
In this section, we discuss how to choose the truncation $\s_r$ and, the function $w$ and parameter $r$ in its definition \eqref{sublevel}.
% 
For the sake of simplicity, we restrict this section to the case of a unique stationary solution and its  approximation $l^r$ composed of lower bounds (given by \eqref{eq:up} and \eqref{eq:epu}). Analogous arguments hold for the upper bounds. 

The approximation error is $\norm{\pi-l^r}=1-l^r(\s_r)$ (Corollary \ref{compstatcc}) and it can be shown to be at least $c/r$ (c.f. \cite[Lem.~6.12]{Kuntzthe}). Thus, if we wish to guarantee an error smaller than $\epsilon$, we must choose $r$ to be at least $c/\epsilon$. When it comes to $w$, it is best to choose functions such that the size of the truncation \eqref{sublevel} grows slowly with $r$. In other words, functions $w$ that tend to infinity quickly. For instance, for chains that take values in $\n$ and $w(x)=x^\alpha$, higher values of $\alpha$ require smaller truncation sizes $\mmag{\s_r}\approx \sqrt[\alpha]{r}$. Of course, the bounding constants $c$ involved in the moment bound \eqref{eq:momboundw} will be larger for higher moments. However, in our experience, the larger constants only hamper the bounds for small $r$, while for large $r$ they are more than offset by the rapid growth of $w$ (cf. Fig.~\ref{fig:schodist}c,d). 

Generally, one wishes to use functions $w$ defining truncations that cover most of the probability mass.  One could guide the selection by running the scheme a few times using various norm-like functions $w$ to obtain information about the distribution shape. However, an educated guess is often useful, as we illustrate in the example of Sec.~\ref{togglesec}.
%

\subsubsection{Dealing with ill-conditioning by scaling the decision variables}\label{apriori}

For sufficiently large truncations, the orders of magnitude  of the coefficients in the constraints $\pi^r Q(x)=0$ for all $x\in\cal{N}_r$ and $\aver{w}\leq c$ vary greatly. This leads to large round-off errors in the double-precision floating point arithmetic LP solvers typically employ and poor solver performance. A straightforward way to ameliorate this issue is to scale the decision variables so that range of orders of magnitude in the constraint coefficients is smaller. We found that scalings of the form
%
\begin{align*}\tilde{\pi}^r(x)&:=\left\{\begin{array}{ll}-q(x,x)\pi^r(x)&\text{if }q(x,x)\neq0\\\pi^r(x)&\text{if }q(x,x)=0\end{array}\right.\text{ or }\\
 \tilde{\pi}^r(x)&:=\left\{\begin{array}{ll}w(x)\pi^r(x)&\text{if }w(x)\neq0\\\pi^r(x)&\text{if }w(x)=0\end{array}\right.\end{align*}
%
can significantly improve solver performance.
\section{Application: a toggle switch}
\label{togglesec}
%
\begin{figure*}
	\begin{center}
	\includegraphics[width=\textwidth]{5.pdf}
	\vspace{-20pt}
	\end{center}
\caption{\textbf{Bounds on the marginal probability distributions of the toggle switch. }
\textbf{(a)} Lower bounds $\hat{l}^r$ (solid lines) on the marginal probability distributions for increasing state space truncations $r=40^6$ (blue), $45^6$ (red), $50^6$ (yellow) and $120^6$ (purple). Both sets of bounds converge to the stationary distribution of the CME. For comparison, we show simulations performed using the Gillespie Algorithm (black dots). 
% Simulation time: $27$ minutes 
\textbf{(b)} Approximation error $\hat{\varepsilon}^l_r$ of $\hat{l}^r$ (blue) and bound $\hat{\varepsilon}^u_r$ on error of $\hat{u}^r$ (yellow).  
%Lines show the error bounds obtained for the marginal distribution of protein $P_1$, while dots show those the corresponding to the marginal distribution of protein $P_2$. The total error of all estimates decreases with the number of states considered ($r$). 
Parameters are as in Fig.~\ref{fig:toggle_joint}.}
\label{fig:toggle_marginals}
\end{figure*}
%
%
%
%
As an application of the methods discussed in this paper, we consider a network of two mutually repressing genes. Such toggle switches are common motifs in many cell fate decisions genetic circuits\cite{gardner2000,hemberg2007,perez2016}. In particular, we consider the asymmetric case 
%
$$
\varnothing \xrightarrow{a_1} P_1\xrightarrow{a_2}\varnothing, \qquad \varnothing \xrightarrow{a_3} P_2\xrightarrow{a_4}\varnothing,
$$
with propensities
%
\begin{align}
 \label{toggle:prop}
 &a_1(x) = \frac{k_1}{1+(x_2/\theta)^3}, \quad
 &a_2(x) = k_2 x_1, \notag\\
 &a_3(x) = \frac{k_3}{1+x_1}, \quad
 &a_4(x) = k_4 x_2,
\end{align}
%
with $k_1,k_2,k_3,k_4,\theta>0$, which models mutual repression via Hill functions and dilution/degradation via linear unspecific decay, where $x=(x_1,x_2)$ and $x_1$ (resp. $x_2$) denotes the copy number of protein $P_1$ (resp. $P_2$). As shown in App.~\ref{appendix}, the CME has a unique stationary solution $\pi$ and all of its moments are finite. In what follows, we fix $\theta:=1$, $k_1:=30$, $k_2:=k_4:=1$ and $k_3:=10$.
%
%
\subsection{Bounds on the joint distribution}
%
%
%
To obtain a tail bound, we require an upper bound on a moment. Because our example involves rational propensities, we used rational moments as defined in \eqref{eq:momvec} with
%
\begin{align*}
 o(x) := \left({1+x_1}\right) ({1+(x_2/\theta)^3}),
\end{align*}
%
being the product of the denominators of the propensities in \eqref{toggle:prop}. Due to the asymmetric repression of the two proteins, the range of molecule numbers for protein $P_1$ is larger than for protein $P_2$. In an effort to capture the relevant part of the state space in our truncation, we picked
%
$$w(x):=(x_1+2x_2)^6,\qquad\forall x\in\n^2.$$
% 
This choice was guided by quickly obtaining tight bounds on $\ave{x_1}$ and $\ave{x_2}$ using the semidefinite programming scheme of Sec.~\ref{bmom}, which showed that $\ave{x_1}\approx2\ave{x_2}$. We then computed the moment bound 
%%
$$\ave{w}\leq c := 5.0901\times10^8$$
%
by solving the rightmost SDP in \eqref{eq:sdps} with $d:=10$ and $f:=w$. Markov's inequality implies that the mass of the tail $m_r$ is bounded above by $c/r$ depending on the state space truncation parameter $r$.

Using the linear programming approach described in Sec.~\ref{sec:distributions}, we computed upper and lower bounds $u^r$ and $l^r$ on the stationary distribution for small ($r=44$) and large ($r=74$) state space truncations as shown in Fig.~\ref{fig:toggle_joint}. These bounds give an accurate account of the uncertainty in the probabilities using only the first few states of the CME. For small truncations, we find that the maximum absolute discrepancies are found near the distribution modes, as seen in Fig.~\ref{fig:toggle_joint}(a,b). Increasing the truncation, the upper and lower bounds shown in Fig.~\ref{fig:toggle_joint}(c) are nearly indistinguishable: the maximum discrepancy drops under $10^{-9}$ (Fig.~\ref{fig:toggle_joint}(d)). The approximation error (including the tail) is of at most $5.3\times10^{-3}$ (evaluated using Cor.~\ref{compstatcc}$(i)$).
%
\subsection{Bounds on the marginal distributions}

%
\begin{figure*}
	\begin{center}
	\includegraphics[width=0.75\textwidth]{6.pdf}
	\vspace{-20pt}
	\end{center}
\caption{\textbf{Flipping the toggle switch: deterministic vs. stochastic models. } 
Tuning the promoter dissociation constant $\theta$ from small to high values results in a switch-like increase of the production of protein $P_1$ molecules while simultaneously repressing the production of protein $P_2$. In deterministic models such tuning leads to two stable steady states in the transition region (solid red lines) accessible from different initial conditions. In the stochastic model the marginal stationary probabilities (heatmap) displays two coexisting populations. We observe good correspondence between the modes of the distributions and the deterministically stable steady states.
% 
Each marginal is obtained from lower bounds computed using $r = 83$, resulting in $1806$ equations involving $1891$ states. The approximation error $\hat{\varepsilon}^l_r$ was at most $3\times 10^{-3}$ for all parameter values. 
Parameters are as in the main text except for $\theta$.}
\label{fig:toggle_bifurcation}
\end{figure*}
% SIMULATION TIME: 1.6279e3, 0.0032 l1 norm error (double check), needed 2e8 steps
%
The method described in Sec.~\ref{sec:mdist} (with $A_i$ as in \eqref{eq:kthm}) allows us to obtain tight lower bounds $\hat{l}^r$ with computable errors on the marginal distributions. The corresponding bounds for the distributions of protein $P_1$ and $P_2$ are shown in Fig.~\ref{fig:toggle_marginals}(a) and demonstrate the convergence of our approximations. To verify our implementation we compare with histograms obtained by running the Gillespie Algorithm and using \eqref{eq:timeaverages}. In Fig.~\ref{fig:toggle_marginals}(b) we show the corresponding approximation error in terms of the total variation distance between the lower bounds $\hat{l}^r$ and the corresponding marginal. We also show the computable error bound (yellow) on the approximation $\hat{u}^r$ of the marginals

Lastly, we apply the method to gain insights over whole parameter ranges. In particular, increasing the dissociation constant $\theta$ of protein $P_2$ from zero to finite values allows to induce expression of protein $P_1$. To quantify this dependence we use the lower bounds on the marginal distributions shown in Fig.~\ref{fig:toggle_bifurcation}. For intermediate levels of $\theta$, we observe that induction of $P_1$ does not occur gradually but is represented by two coexisting populations corresponding to fully induced and repressed populations, i.e. one with high levels of $P_1$ but low levels of $P_2$ and another one with low levels of $P_1$ but high levels of $P_2$. We find that the modes of the marginal distributions are in good correspondence with the stable solutions of the deterministic steady-state rate equations.

\vspace{5pt}
\noindent\textbf{Implementation details:} To set-up the LPs we used MATLAB R2017a, and to solve them we used the CPLEX V12.6.3\cite{CPLEX}. As in Sec.~\ref{sec:mommulti}, for the SDPs we used YALMIP\cite{Lofberg2004}, SDPA-GMP\cite{Nakata2010}, and mpYALMIP\cite{Fantuzzi2016}.
%

\section{Discussion}
\label{sec:discussion}

We introduced two mathematical programming approaches that yield bounds on the stationary moments and distributions of biochemical reaction networks. These statistical quantities satisfy the stationary moment equations and the stationary CME, respectively, which are typically infinite sets of coupled equations. By considering a finite subset of these equations, the mathematical programming approaches considered here enable us to bound the set of solutions using semidefinite and linear programming. These bounds provide accurate estimates of moments and probabilities with quantifiable errors.

To compute the moment bounds, we append semidefinite inequalities to the moments equations that are satisfied by the moments of the underlying probability distributions. This allow us to formulate an SDP, whose solution is an upper or a lower bound on any given moment. In contrast with moment closure methods \cite{engblom2006,engblom2014,lakatos2015,schnoerr2015}, the proposed method yields moment approximations with quantified errors. We show that repeated applications of our method yield upper (resp. lower) bounds that are monotonically decreasing (resp. increasing) with the number of moment equations and inequalities used (Theorem~\ref{inclu}). As mentioned at the end of Sec.~\ref{sec:mommulti}, there are reaction networks for which these bounds do not converge to the exact moments. In practice, however, they often do converge as we showed using the example in Fig.~\ref{fig:schomom}, see also the examples in the recent studies\cite{Dowdy2017,Dowdy2018,Ghusinga2017a,Sakurai2017}. Sharper bounds certainly apply for restricted state spaces \cite{Lasserre2009,Laurent2009,Blekherman2013,Dowdy2018}, but these refinements are beyond the scope of this paper.

To compute bounds on the stationary distributions (including its marginals or averages), we formulated an LP. In the case where the CME has a unique solution, the method yields converging lower and upper bounds on this solution accompanied by simple error bounds (Corollaries \ref{compstatcc} and \ref{compstatccm}). More generally, the method provides  insights about whether the stationary solution is unique (Corollary \ref{uniqueness-test}). If several solutions to the CME exist, our method provides bounds over the set of possible solutions. As elaborated in Sec.~\ref{sec:nonunique}, the method can be adapted to compute approximations of the ergodic distributions, which are of practical interest, but computable error bounds are not available in this case. 



Several other truncation-based schemes have been proposed to approximate the stationary solutions of the CME. In contrast with ours, these schemes typically assume that the CME has a unique stationary distribution, which has to be verified separately \cite{pauleve2014,Gupta2018}. Perhaps the most extensively studied method is the \emph{truncation--and--augmentation scheme} originally proposed by Seneta\cite{Seneta1967} for discrete time chains. The continuous-time counterpart\cite{Tweedie1971,Hart2012,Gupta2017} of this scheme is known to converge\cite{Hart2012} in total variation for exponentially ergodic chains, monotone chains, and certain generalisations thereof. While hard bounds can be obtained\cite{Hart2012,Gupta2017,Tweedie1998,Meyn1994,Liu2015,Masuyama2017,Masuyama2017a,Liu2018,Liu2018a} on the scheme's approximation error, these bounds can be very conservative and often involve constants that are difficult to compute in practice.

Spieler et al. overcame this issue by iterating the truncation--and--augmentation scheme and applying a tail bound derived from a Foster-Lyapunov criterion to obtain upper and lower bounds on the stationary distribution \cite{Dayar2011,Spieler2014}. Even though the procedure has not been shown to converge, error bounds can be obtained by combining the upper and lower bounds \cite{Dayar2011,Spieler2014} (or as done in this paper). This truncation-based scheme uses tail bounds similar to our method but entails solving only systems of linear equations. This means they are simpler to implement than linear programs and carry smaller computational costs, but they offer no guarantee of convergence. Future improvements in LP and SDP solvers will likely mitigate some of these issues. In fact, the search of moment bounds and that of Lyapunov functions are dual problems in many respects\cite{Glynn2008,Lasserre2009,Laurent2009}.


A distinct feature of our method is that it enables the direct computation of lower and upper bounds on the marginal distributions. In particular, we do not need to compute bounds for each state of joint distribution, but only bounds for each state of the marginal distribution of interest. These distributions are of particular interest for the analysis of high-dimensional networks.

We expect that our approach will be particularly valuable in applications where accurate approximations with quantified errors are needed. These include, for example, estimating distributions of phenotypic switches, determining persister fractions in bacterial populations and inferring model parameters from noisy single cell data \cite{fox2016}.
%
%
%
\section{Acknowledgements}
%
We thank Justine Dattani and Diego Oyarz\'{u}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 by the EPSRC Fellowship EP/M002187/1; and MB by the EPSRC grant EP/N014529/1.
%
%
%
%
\appendix

\section{Minimal continuous-time Markov chains, their long-term behaviour, the stationary distributions, and a Foster-Lyapunov criterion}\label{appendix}

In practice, the chain $X=\{X(t)\}_{t\geq0}$ mentioned in Sec.~\ref{sec:stationary} is often constructed by running the Gillespie Algorithm\cite{Gillespie1976,Kendall1950,Feller1940}: sample  a state $x$ from an \emph{initial distribution} $\lambda:=\{\lambda(x)\}_{x\in\s}$ and start the chain at a state $x$. If $q(x)$ in \eqref{eq:qmatrix} 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 \emph{jumps} from $x$ to $y$ and we call the time at which it jumps the \emph{jump time}). Repeat these steps starting from $y$ instead of $x$. All random variables sampled must be independent of each other, see\cite[Sec.~2.3]{Kuntz2018} for a more detailed construction. 

If $T_n$ denotes the $n^{th}$ jump time, then the limit
%
$$T_\infty:=\lim_{n\to\infty}T_n$$
%
is known as the \emph{explosion time} (or \emph{escape time} or \emph{time of the first infinity}) of the chain because it is the instant by which the chain will have left every finite subset of the state space\cite[Sec.~2]{Kuntz2018}. If no such \emph{explosion} occurs, then $T_\infty$ is infinity. We say that the chain is \emph{non-explosive} if 
%
\begin{equation}\label{eq:nonexp}\Pbl{\{T_\infty=\infty\}}=1,\end{equation}
%
where $\Pb_\lambda$ denotes probability measure underlying the chain (we append the subscript $\lambda$ to emphasise that the chain's starting position was sampled from the distribution $\lambda$). If \eqref{eq:nonexp} holds for every probability distribution $\lambda$ ($\lambda(x)\geq0$ for all $x\in\s$ and $\sum_{x\in\s}\lambda(x)=1$), then the rate matrix $Q$ is said to be \emph{regular}. We collectively refer to the probabilities 
%
$$\{p_t(x)\}_{x\in\s,t\geq0}=\{\Pbl{\{X_t=x,t<T_\infty\}}\}_{x\in\s,t\geq0}$$
%
of observing the process in the state $x:=(x_1,\ldots,x_n)\in\s$ at time $t\geq0$ as the \emph{time-varying law} of the chain. It is\cite[Cor.~2.21]{Kuntzthe} the minimal non-negative solution of the CME \eqref{eq:CME}. 

A probability distribution $\pi:=\{\pi(x)\}_{x\in\s}$ on $\s$ is said to be a \emph{stationary} (or \emph{steady-state} or \emph{equilibrium} or \emph{invariant}) \emph{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:
%
\begin{equation}\label{eq:statdef}\Pbp{\{X_t=x,t<T_\infty\}}=\pi(x),\quad\forall x\in\s,\enskip t\geq0.\end{equation}
%
%
Summing both sides of \eqref{eq:statdef} over $x\in\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:
%
\begin{equation}\label{eq:statnoexp}\Pbp{\{T_\infty=\infty\}}=1.\end{equation}
%
Taking the derivative in time of \eqref{eq:statdef}, we find that stationary distributions are stationary solutions of the CME \eqref{eq:CME} (that is, it belongs to \eqref{eq:cmestat}). The reverse direction is slightly more complicated:
%
\begin{theorem}[Theorem 2.41\cite{Kuntzthe}]\label{Qstateq} Let $X$ be a continuous-time chain with rate matrix $Q$ satisfying \eqref{eq:qmatrix}. A probability distribution $\pi$ on $\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$: \eqref{eq:statnoexp} holds.
\end{theorem}
%
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, \eqref{eq:cmestat} is an \emph{analytical} (as in non-probabilistic) \emph{linear programming characterisation}\cite{Kuntzthe,Kurtz1998} of the set of stationary distributions for regular $Q$. The non-explosivity in Theorem \ref{Qstateq} is crucial: a counterexample is the birth-death process \eqref{eq:one_step} with $a_+(x):=2^{2x}$ and $a_-(x):=2^{2x}/2$. In this case, the sum in \eqref{eq:normalising} is finite showing that the CME has a unique stationary solution $\pi$ given by \eqref{eq:1deg}--\eqref{eq:normalising}. However, \cite[Theorem 11]{Reuter1957} shows that the process is explosive for any initial distribution (including $\pi$) and it follows from \eqref{eq:statnoexp} that no stationary distribution exists.

Stationary distributions are of interest because, if the chain is positive recurrent, then, regardless of the initial distribution $\lambda$, they determine\cite{Meyn1993a} the chain's long term behaviour in the sense that the time-varying law of the chain converges to one of them:
%
\begin{equation}\label{eq:spaceaverages}\lim_{t\to\infty}p_t(x)=\pi(x)\quad\forall x\in\s,\end{equation}
%
and that the fraction of time that the chain spends in any given state tends to the probability that one of the stationary distribution awards to the state:
%
\begin{equation}\label{eq:timeaverages}
\lim_{T\to\infty}\frac{1}{T}\int_0^T1_x(X_t)dt=\pi(x)\quad\forall x\in\s,\qquad\Pb_\lambda\text{-a.s.}
\end{equation}

If the chain starts in one of the closed communicating classes (defined in Sec.~\ref{nonuniquesec}), then it can never escape it. The convergence of the time averages \eqref{eq:timeaverages} then shows that, in the positive recurrent case, there must exist at least one stationary distribution per closed communicating class $\cal{C}_j$ and that the stationary distribution must have support contained in $\cal{C}_j$ ($\pi_j(\cal{C}_j)=1$). This distribution $\pi_j$ is unique and, if the initial distribution $\lambda$ has its mass contained in $\cal{C}_j$, then\cite{Meyn1993a} both the time averages and the space averages converge to $\pi_j$ (in the sense that $\pi$ featuring in both \eqref{eq:spaceaverages} and \eqref{eq:timeaverages} is $\pi_j$). For this reason, $\pi_j$ is known as an \emph{ergodic distribution} of the chain. The transient set $\cal{T}:=\s\backslash\cup_{j}\cal{C}_j$ featuring in the Doeblin Decomposition \eqref{eq:ddc} is called ``transient'' because its definition implies that if the chain leaves a state in $\cal{T}$ it cannot return to this state. For this reason, the chain will eventually leave the set and never return in the sense that its sample paths will enter the closed communicating classes or tend to infinity (in particular $p_t(x)\to0$ as $t\to\infty$ for any $x$ in $\cal{T}$). In the case of a positive recurrent chain, tending to infinity is not an option and so the chain eventually enters one of the closed communicating classes. It then follows from \eqref{eq:statdef} that no stationary distribution $\pi$ such that $\pi(x)>0$ exists for a state $x$ in $\cal{T}$. Bringing this discussion together\cite{Meyn1993a}, we have that stationary distribution $\pi$ in \eqref{eq:timeaverages} is the ergodic distribution $\pi_j$ of the closed communicating class $\cal{C}_j$ that the chain's path eventually enters, while that in \eqref{eq:spaceaverages} is a weighted combination of the ergodic distributions where the weight given to $\pi_j$ is the probability that the chain ever enters $\cal{C}_j$. Theorem \ref{doeblinc} in Sec.~\ref{nonuniquesec} follows from these facts.

In practice, verifying whether a chain is positive recurrent is done by applying a Foster-Lyapunov criterion. For our examples, we will use the following well-known criterion that involves a norm-like function $w$ meaning a non-negative function on $\s$ satisfying \eqref{sublevel}.
%
\begin{theorem}\label{lyap} If there exists constants $K_1\in\r$, $K_2>0$, and a norm-like function $w$ satisfying
%
$$Qw(x):=\sum_{y\in\s}q(x,y)w(y)\leq K_1-K_2w(x)\qquad \forall x\in\s,$$
%
then
%
\begin{enumerate}[label=(\roman*)]
\item The rate matrix $Q$ is regular.
\item There exists at least one stationary distribution.
\item Every stationary distribution $\pi$ satisfies $\ave{w}<\infty.$
\item The stationary distributions determine the long-term behaviour: limits \eqref{eq:spaceaverages}--\eqref{eq:timeaverages} hold for any starting distribution $\lambda$.
\end{enumerate}
\end{theorem}
%
\begin{proof}Part $(i)$ is \cite[Theorem 1.11]{Chen1991} (see also \cite[Theorem 2.1]{Meyn1993b}). Parts $(ii)$--$(iv)$ follow from \cite[Theorems 8.1 and 8.2]{Meyn1993a} and \cite[Theorem 4.6]{Meyn1993b}.\end{proof}

For example, in the case of  Schl\"ogl's model \eqref{eq:smdl}, fixing $w(x):=x^{d-2}$ for some integer $d>2$, we have that
%
$$Qw(x)=g_{d-1}(x)-k_2(d-2)x^{d},$$
%
where $g_{d-1}$ is a polynomial of degree $d-1$. Thus,
%
\begin{align*}Qw(x)&\leq \sup_{x\in\n}\left\{g_{d-1}(x)-\frac{k_2(d-2)}{2}x^{d}\right\}-\frac{k_2(d-2)}{2}x^{d}\\
&\leq \sup_{x\in\n}\left\{g_{d-1}(x)-\frac{k_2(d-2)}{2}x^{d}\right\}-\frac{k_2(d-2)}{2}x^{d-2},\end{align*}
%
and the supremum is finite. Picking $d\geq 3$, Theorem \ref{lyap} then tells us that \eqref{eq:smdl} has a regular rate matrix and at least one stationary distribution, that the moments $\ave{x^1},\dots,\ave{x^{d-2}}$ are finite for any stationary distribution $\pi$, and that the limits \eqref{eq:spaceaverages}--\eqref{eq:timeaverages} hold. Because we can choose ever larger $d$, we have finiteness of all moments. Uniqueness of this distribution follows from \eqref{eq:1deg}--\eqref{eq:normalising}.
 
Similarly, in the case of the toggle switch chain of Sec.~\ref{togglesec}, setting $w(x):=(x_1+x_2)^d$, we have that
%
$$\begin{array}{ll}Qw(x) &=(a_1(x)+a_3(x))((x_1+x_2+1)^d-(x_1+x_2)^d)\\ 
&+(a_2(x)+a_4(x))((x_1+x_2-1)^d-(x_1+x_2)^d) \\ 
&\leq (k_1+k_3)((x_1+x_2+1)^d-(x_1+x_2)^d)\\ 
&+(a_2(x)+a_4(x))((x_1+x_2-1)^d-(x_1+x_2)^d)\\
&\leq g_{d-1}(x_1+x_2)-d(k_2x_1+k_4x_2)(x_1+x_2)^{d-1}\\
&\leq g_{d-1}(x_1+x_2)-d\min\{k_2,k_4\}(x_1+x_2)^d,\end{array}$$
%
where $g_{d-1}$ is a polynomial of degree $d-1$. For this reason, using Theorem \ref{lyap} the same way we did above for Schl\"ogl'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.~\ref{fig:toggle_joint} and Corollary \ref{uniqueness-test} show that it is unique.


\section*{References}
\bibliographystyle{apsrev4-1}
\bibliography{lit}

\end{document}


\grid

main.bbl
%merlin.mbs apsrev4-1.bst 2010-07-25 4.21a (PWD, AO, DPC) hacked
%Control: key (0)
%Control: author (72) initials jnrlst
%Control: editor formatted (1) identically to author
%Control: production of article title (-1) disabled
%Control: page (0) single
%Control: year (1) truncated
%Control: production of eprint (0) enabled
\begin{thebibliography}{73}%
\makeatletter
\providecommand \@ifxundefined [1]{%
 \@ifx{#1\undefined}
}%
\providecommand \@ifnum [1]{%
 \ifnum #1\expandafter \@firstoftwo
 \else \expandafter \@secondoftwo
 \fi
}%
\providecommand \@ifx [1]{%
 \ifx #1\expandafter \@firstoftwo
 \else \expandafter \@secondoftwo
 \fi
}%
\providecommand \natexlab [1]{#1}%
\providecommand \enquote  [1]{``#1''}%
\providecommand \bibnamefont  [1]{#1}%
\providecommand \bibfnamefont [1]{#1}%
\providecommand \citenamefont [1]{#1}%
\providecommand \href@noop [0]{\@secondoftwo}%
\providecommand \href [0]{\begingroup \@sanitize@url \@href}%
\providecommand \@href[1]{\@@startlink{#1}\@@href}%
\providecommand \@@href[1]{\endgroup#1\@@endlink}%
\providecommand \@sanitize@url [0]{\catcode `\\12\catcode `\$12\catcode
  `\&12\catcode `\#12\catcode `\^12\catcode `\_12\catcode `\%12\relax}%
\providecommand \@@startlink[1]{}%
\providecommand \@@endlink[0]{}%
\providecommand \url  [0]{\begingroup\@sanitize@url \@url }%
\providecommand \@url [1]{\endgroup\@href {#1}{\urlprefix }}%
\providecommand \urlprefix  [0]{URL }%
\providecommand \Eprint [0]{\href }%
\providecommand \doibase [0]{http://dx.doi.org/}%
\providecommand \selectlanguage [0]{\@gobble}%
\providecommand \bibinfo  [0]{\@secondoftwo}%
\providecommand \bibfield  [0]{\@secondoftwo}%
\providecommand \translation [1]{[#1]}%
\providecommand \BibitemOpen [0]{}%
\providecommand \bibitemStop [0]{}%
\providecommand \bibitemNoStop [0]{.\EOS\space}%
\providecommand \EOS [0]{\spacefactor3000\relax}%
\providecommand \BibitemShut  [1]{\csname bibitem#1\endcsname}%
\let\auto@bib@innerbib\@empty
%</preamble>
\bibitem [{\citenamefont {Elowitz}\ \emph {et~al.}(2002)\citenamefont
  {Elowitz}, \citenamefont {Levine}, \citenamefont {Siggia},\ and\
  \citenamefont {Swain}}]{elowitz2002}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {M.~B.}\ \bibnamefont
  {Elowitz}}, \bibinfo {author} {\bibfnamefont {A.~J.}\ \bibnamefont {Levine}},
  \bibinfo {author} {\bibfnamefont {E.~D.}\ \bibnamefont {Siggia}}, \ and\
  \bibinfo {author} {\bibfnamefont {P.~S.}\ \bibnamefont {Swain}},\ }\href
  {\doibase 10.1126/science.1070919} {\bibfield  {journal} {\bibinfo  {journal}
  {Science}\ }\textbf {\bibinfo {volume} {297}},\ \bibinfo {pages} {1183}
  (\bibinfo {year} {2002})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Taniguchi}\ \emph {et~al.}(2010)\citenamefont
  {Taniguchi}, \citenamefont {Choi}, \citenamefont {Li}, \citenamefont {Chen},
  \citenamefont {Babu}, \citenamefont {Hearn}, \citenamefont {Emili},\ and\
  \citenamefont {Xie}}]{taniguchi2010}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {Y.}~\bibnamefont
  {Taniguchi}}, \bibinfo {author} {\bibfnamefont {P.~J.}\ \bibnamefont {Choi}},
  \bibinfo {author} {\bibfnamefont {G.-W.}\ \bibnamefont {Li}}, \bibinfo
  {author} {\bibfnamefont {H.}~\bibnamefont {Chen}}, \bibinfo {author}
  {\bibfnamefont {M.}~\bibnamefont {Babu}}, \bibinfo {author} {\bibfnamefont
  {J.}~\bibnamefont {Hearn}}, \bibinfo {author} {\bibfnamefont
  {A.}~\bibnamefont {Emili}}, \ and\ \bibinfo {author} {\bibfnamefont {X.~S.}\
  \bibnamefont {Xie}},\ }\href {\doibase 10.1126/science.1188308} {\bibfield
  {journal} {\bibinfo  {journal} {Science}\ }\textbf {\bibinfo {volume}
  {329}},\ \bibinfo {pages} {533} (\bibinfo {year} {2010})}\BibitemShut
  {NoStop}%
\bibitem [{\citenamefont {Uphoff}\ \emph {et~al.}(2016)\citenamefont {Uphoff},
  \citenamefont {Lord}, \citenamefont {Okumus}, \citenamefont
  {Potvin-Trottier}, \citenamefont {Sherratt},\ and\ \citenamefont
  {Paulsson}}]{uphoff2016}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {S.}~\bibnamefont
  {Uphoff}}, \bibinfo {author} {\bibfnamefont {N.~D.}\ \bibnamefont {Lord}},
  \bibinfo {author} {\bibfnamefont {B.}~\bibnamefont {Okumus}}, \bibinfo
  {author} {\bibfnamefont {L.}~\bibnamefont {Potvin-Trottier}}, \bibinfo
  {author} {\bibfnamefont {D.~J.}\ \bibnamefont {Sherratt}}, \ and\ \bibinfo
  {author} {\bibfnamefont {J.}~\bibnamefont {Paulsson}},\ }\href {\doibase
  10.1126/science.aac9786} {\bibfield  {journal} {\bibinfo  {journal}
  {Science}\ }\textbf {\bibinfo {volume} {351}},\ \bibinfo {pages} {1094}
  (\bibinfo {year} {2016})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Arkin}\ \emph {et~al.}(1998)\citenamefont {Arkin},
  \citenamefont {Ross},\ and\ \citenamefont {McAdams}}]{arkin1998}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {A.}~\bibnamefont
  {Arkin}}, \bibinfo {author} {\bibfnamefont {J.}~\bibnamefont {Ross}}, \ and\
  \bibinfo {author} {\bibfnamefont {H.~H.}\ \bibnamefont {McAdams}},\
  }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal} {Genetics}\
  }\textbf {\bibinfo {volume} {149}},\ \bibinfo {pages} {1633} (\bibinfo {year}
  {1998})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Dandach}\ and\ \citenamefont
  {Khammash}(2010)}]{Dandach2010}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {S.~H.}\ \bibnamefont
  {Dandach}}\ and\ \bibinfo {author} {\bibfnamefont {M.}~\bibnamefont
  {Khammash}},\ }\href {\doibase 10.1371/journal.pcbi.1000985} {\bibfield
  {journal} {\bibinfo  {journal} {PLOS Comput. Biol.}\ }\textbf {\bibinfo
  {volume} {6}},\ \bibinfo {pages} {e1000985} (\bibinfo {year}
  {2010})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Thomas}\ \emph {et~al.}(2014)\citenamefont {Thomas},
  \citenamefont {Popovi{\'{c}}},\ and\ \citenamefont {Grima}}]{thomas2014}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {P.}~\bibnamefont
  {Thomas}}, \bibinfo {author} {\bibfnamefont {N.}~\bibnamefont
  {Popovi{\'{c}}}}, \ and\ \bibinfo {author} {\bibfnamefont {R.}~\bibnamefont
  {Grima}},\ }\href {\doibase 10.1073/pnas.1400049111} {\bibfield  {journal}
  {\bibinfo  {journal} {Proc. Natl. Acad. Sci.}\ }\textbf {\bibinfo {volume}
  {111}},\ \bibinfo {pages} {6994} (\bibinfo {year} {2014})}\BibitemShut
  {NoStop}%
\bibitem [{\citenamefont {Munsky}\ \emph {et~al.}(2015)\citenamefont {Munsky},
  \citenamefont {Fox},\ and\ \citenamefont {Neuert}}]{munsky2015}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {B.}~\bibnamefont
  {Munsky}}, \bibinfo {author} {\bibfnamefont {Z.}~\bibnamefont {Fox}}, \ and\
  \bibinfo {author} {\bibfnamefont {G.}~\bibnamefont {Neuert}},\ }\href
  {\doibase 10.1016/j.ymeth.2015.06.009} {\bibfield  {journal} {\bibinfo
  {journal} {Methods}\ }\textbf {\bibinfo {volume} {85}},\ \bibinfo {pages}
  {12} (\bibinfo {year} {2015})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Fr{\"{o}}hlich}\ \emph {et~al.}(2016)\citenamefont
  {Fr{\"{o}}hlich}, \citenamefont {Thomas}, \citenamefont {Kazeroonian},
  \citenamefont {Theis}, \citenamefont {Grima},\ and\ \citenamefont
  {Hasenauer}}]{frohlich2016}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {F.}~\bibnamefont
  {Fr{\"{o}}hlich}}, \bibinfo {author} {\bibfnamefont {P.}~\bibnamefont
  {Thomas}}, \bibinfo {author} {\bibfnamefont {A.}~\bibnamefont {Kazeroonian}},
  \bibinfo {author} {\bibfnamefont {F.~J.}\ \bibnamefont {Theis}}, \bibinfo
  {author} {\bibfnamefont {R.}~\bibnamefont {Grima}}, \ and\ \bibinfo {author}
  {\bibfnamefont {J.}~\bibnamefont {Hasenauer}},\ }\href {\doibase
  10.1371/journal.pcbi.1005030} {\bibfield  {journal} {\bibinfo  {journal}
  {PLoS Comput. Biol.}\ }\textbf {\bibinfo {volume} {12}},\ \bibinfo {pages}
  {e1005030} (\bibinfo {year} {2016})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Neuert}\ \emph {et~al.}(2013)\citenamefont {Neuert},
  \citenamefont {Munsky}, \citenamefont {Tan}, \citenamefont {Teytelman},
  \citenamefont {Khammash},\ and\ \citenamefont {van
  Oudenaarden}}]{neuert2013}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {G.}~\bibnamefont
  {Neuert}}, \bibinfo {author} {\bibfnamefont {B.}~\bibnamefont {Munsky}},
  \bibinfo {author} {\bibfnamefont {R.~Z.}\ \bibnamefont {Tan}}, \bibinfo
  {author} {\bibfnamefont {L.}~\bibnamefont {Teytelman}}, \bibinfo {author}
  {\bibfnamefont {M.}~\bibnamefont {Khammash}}, \ and\ \bibinfo {author}
  {\bibfnamefont {A.}~\bibnamefont {van Oudenaarden}},\ }\href {\doibase
  10.1126/science.1231456} {\bibfield  {journal} {\bibinfo  {journal}
  {Science}\ }\textbf {\bibinfo {volume} {339}},\ \bibinfo {pages} {584}
  (\bibinfo {year} {2013})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Oyarzún}\ \emph {et~al.}(2014)\citenamefont
  {Oyarzún}, \citenamefont {Lugagne},\ and\ \citenamefont
  {Stan}}]{oyarzun2014}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {D.~A.}\ \bibnamefont
  {Oyarzún}}, \bibinfo {author} {\bibfnamefont {J.-B.}\ \bibnamefont
  {Lugagne}}, \ and\ \bibinfo {author} {\bibfnamefont {G.-B.~V.}\ \bibnamefont
  {Stan}},\ }\href {\doibase 10.1021/sb400126a} {\bibfield  {journal} {\bibinfo
   {journal} {ACS Synth. Biol.}\ }\textbf {\bibinfo {volume} {4}},\ \bibinfo
  {pages} {116} (\bibinfo {year} {2014})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Zechner}\ \emph {et~al.}(2016)\citenamefont
  {Zechner}, \citenamefont {Seelig}, \citenamefont {Rullan},\ and\
  \citenamefont {Khammash}}]{zechner2016}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {C.}~\bibnamefont
  {Zechner}}, \bibinfo {author} {\bibfnamefont {G.}~\bibnamefont {Seelig}},
  \bibinfo {author} {\bibfnamefont {M.}~\bibnamefont {Rullan}}, \ and\ \bibinfo
  {author} {\bibfnamefont {M.}~\bibnamefont {Khammash}},\ }\href {\doibase
  10.1073/pnas.1517109113} {\bibfield  {journal} {\bibinfo  {journal} {Proc.
  Natl. Acad. Sci.}\ }\textbf {\bibinfo {volume} {113}},\ \bibinfo {pages}
  {4729} (\bibinfo {year} {2016})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Meyn}\ and\ \citenamefont
  {Tweedie}(1993{\natexlab{a}})}]{Meyn1993a}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {S.~P.}\ \bibnamefont
  {Meyn}}\ and\ \bibinfo {author} {\bibfnamefont {R.~L.}\ \bibnamefont
  {Tweedie}},\ }\href {\doibase 10.2307/1427521} {\bibfield  {journal}
  {\bibinfo  {journal} {Adv. in Appl. Probab.}\ }\textbf {\bibinfo {volume}
  {25}},\ \bibinfo {pages} {487} (\bibinfo {year}
  {1993}{\natexlab{a}})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Hemberg}\ and\ \citenamefont
  {Barahona}(2007)}]{hemberg2007}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {M.}~\bibnamefont
  {Hemberg}}\ and\ \bibinfo {author} {\bibfnamefont {M.}~\bibnamefont
  {Barahona}},\ }\href {\doibase 10.1529/biophysj.106.099390} {\bibfield
  {journal} {\bibinfo  {journal} {Biophys. J.}\ }\textbf {\bibinfo {volume}
  {93}},\ \bibinfo {pages} {401} (\bibinfo {year} {2007})}\BibitemShut
  {NoStop}%
\bibitem [{\citenamefont {Hemberg}\ and\ \citenamefont
  {Barahona}(2008)}]{Hemberg2008}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {M.}~\bibnamefont
  {Hemberg}}\ and\ \bibinfo {author} {\bibfnamefont {M.}~\bibnamefont
  {Barahona}},\ }\href {\doibase 10.1186/1752-0509-2-42} {\bibfield  {journal}
  {\bibinfo  {journal} {BMC Syst. Biol.}\ }\textbf {\bibinfo {volume} {2}},\
  \bibinfo {pages} {42} (\bibinfo {year} {2008})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Asmussen}\ and\ \citenamefont
  {Glynn}(2007)}]{Asmussen2007}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {S.}~\bibnamefont
  {Asmussen}}\ and\ \bibinfo {author} {\bibfnamefont {W.}~\bibnamefont
  {Glynn}},\ }\href {\doibase 10.1007/978-0-387-69033-9} {\emph {\bibinfo
  {title} {{Stochastic Simulation: Algorithms and Analysis}}}},\ \bibinfo
  {series} {Stochastic Modelling and Applied Probability}, Vol.~\bibinfo
  {volume} {57}\ (\bibinfo  {publisher} {Springer New York},\ \bibinfo
  {address} {New York, NY},\ \bibinfo {year} {2007})\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Engblom}(2006)}]{engblom2006}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {S.}~\bibnamefont
  {Engblom}},\ }\href {\doibase 10.1016/j.amc.2005.12.032} {\bibfield
  {journal} {\bibinfo  {journal} {Appl. Math. Comput.}\ }\textbf {\bibinfo
  {volume} {180}},\ \bibinfo {pages} {498} (\bibinfo {year}
  {2006})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Engblom}\ and\ \citenamefont
  {Pender}(2014)}]{engblom2014}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {S.}~\bibnamefont
  {Engblom}}\ and\ \bibinfo {author} {\bibfnamefont {J.}~\bibnamefont
  {Pender}},\ }\href@noop {} {\  (\bibinfo {year} {2014})},\ \Eprint
  {http://arxiv.org/abs/1406.6164} {arXiv:1406.6164} \BibitemShut {NoStop}%
\bibitem [{\citenamefont {Lakatos}\ \emph {et~al.}(2015)\citenamefont
  {Lakatos}, \citenamefont {Ale}, \citenamefont {Kirk},\ and\ \citenamefont
  {Stumpf}}]{lakatos2015}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {E.}~\bibnamefont
  {Lakatos}}, \bibinfo {author} {\bibfnamefont {A.}~\bibnamefont {Ale}},
  \bibinfo {author} {\bibfnamefont {P.~D.~W.}\ \bibnamefont {Kirk}}, \ and\
  \bibinfo {author} {\bibfnamefont {M.~P.~H.}\ \bibnamefont {Stumpf}},\ }\href
  {\doibase 10.1063/1.4929837} {\bibfield  {journal} {\bibinfo  {journal} {J.
  Chem. Phys.}\ }\textbf {\bibinfo {volume} {143}},\ \bibinfo {pages} {94107}
  (\bibinfo {year} {2015})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Schnoerr}\ \emph {et~al.}(2015)\citenamefont
  {Schnoerr}, \citenamefont {Sanguinetti},\ and\ \citenamefont
  {Grima}}]{schnoerr2015}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {D.}~\bibnamefont
  {Schnoerr}}, \bibinfo {author} {\bibfnamefont {G.}~\bibnamefont
  {Sanguinetti}}, \ and\ \bibinfo {author} {\bibfnamefont {R.}~\bibnamefont
  {Grima}},\ }\href {\doibase 10.1063/1.4934990} {\bibfield  {journal}
  {\bibinfo  {journal} {J. Chem. Phys.}\ }\textbf {\bibinfo {volume} {143}},\
  \bibinfo {pages} {185101} (\bibinfo {year} {2015})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {{Van Kampen}}(1976)}]{vanKampen1976}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {N.~G.}\ \bibnamefont
  {{Van Kampen}}},\ }\href {\doibase 10.1016/0375-9601(76)90398-4} {\bibfield
  {journal} {\bibinfo  {journal} {Phys. Lett. A}\ }\textbf {\bibinfo {volume}
  {59}},\ \bibinfo {pages} {333} (\bibinfo {year} {1976})}\BibitemShut
  {NoStop}%
\bibitem [{\citenamefont {Schwerer}(1996)}]{Schwerer1996}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {E.}~\bibnamefont
  {Schwerer}},\ }\emph {\bibinfo {title} {{A linear programming approach to the
  steady-state analysis of Markov processes}}},\ \href@noop {} {Ph.D. thesis},\
  \bibinfo  {school} {Stanford University} (\bibinfo {year} {1996})\BibitemShut
  {NoStop}%
\bibitem [{\citenamefont {Helmes}\ and\ \citenamefont
  {Stockbridge}(2003)}]{Helmes2003}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {K.}~\bibnamefont
  {Helmes}}\ and\ \bibinfo {author} {\bibfnamefont {R.~H.}\ \bibnamefont
  {Stockbridge}},\ }\href {\doibase 10.1081/STM-120020389} {\bibfield
  {journal} {\bibinfo  {journal} {Stoch. Models}\ }\textbf {\bibinfo {volume}
  {19}},\ \bibinfo {pages} {255} (\bibinfo {year} {2003})}\BibitemShut
  {NoStop}%
\bibitem [{\citenamefont {Hern{\'{a}}ndez-Lerma}\ and\ \citenamefont
  {Lasserre}(2003)}]{Hernandez-Lerma2003}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {O.}~\bibnamefont
  {Hern{\'{a}}ndez-Lerma}}\ and\ \bibinfo {author} {\bibfnamefont {J.~B.}\
  \bibnamefont {Lasserre}},\ }\href {\doibase 10.1007/978-3-0348-8024-4} {\emph
  {\bibinfo {title} {{Markov chains and invariant probabilities}}}}\ (\bibinfo
  {publisher} {Birkh{\"{a}}user},\ \bibinfo {year} {2003})\BibitemShut
  {NoStop}%
\bibitem [{\citenamefont {Kuntz}\ \emph {et~al.}(2016)\citenamefont {Kuntz},
  \citenamefont {Ottobre}, \citenamefont {Stan},\ and\ \citenamefont
  {Barahona}}]{Kuntz2016}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {J.}~\bibnamefont
  {Kuntz}}, \bibinfo {author} {\bibfnamefont {M.}~\bibnamefont {Ottobre}},
  \bibinfo {author} {\bibfnamefont {G.-B.~B.}\ \bibnamefont {Stan}}, \ and\
  \bibinfo {author} {\bibfnamefont {M.}~\bibnamefont {Barahona}},\ }\href
  {\doibase 10.1137/16M107801X} {\bibfield  {journal} {\bibinfo  {journal}
  {SIAM J. Sci. Comput.}\ }\textbf {\bibinfo {volume} {38}},\ \bibinfo {pages}
  {A3891} (\bibinfo {year} {2016})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Tweedie}(1971)}]{Tweedie1971}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {R.~L.}\ \bibnamefont
  {Tweedie}},\ }\href {\doibase 10.2307/3211901} {\bibfield  {journal}
  {\bibinfo  {journal} {J. Appl. Probab.}\ }\textbf {\bibinfo {volume} {8}},\
  \bibinfo {pages} {311} (\bibinfo {year} {1971})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Hart}\ and\ \citenamefont
  {Tweedie}(2012)}]{Hart2012}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {A.~G.}\ \bibnamefont
  {Hart}}\ and\ \bibinfo {author} {\bibfnamefont {R.~L.}\ \bibnamefont
  {Tweedie}},\ }\href {\doibase 10.4236/am.2012.312A301} {\bibfield  {journal}
  {\bibinfo  {journal} {Appl. Math.}\ }\textbf {\bibinfo {volume} {3}},\
  \bibinfo {pages} {2205} (\bibinfo {year} {2012})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Gupta}\ \emph {et~al.}(2017)\citenamefont {Gupta},
  \citenamefont {Mikelson},\ and\ \citenamefont {Khammash}}]{Gupta2017}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {A.}~\bibnamefont
  {Gupta}}, \bibinfo {author} {\bibfnamefont {J.}~\bibnamefont {Mikelson}}, \
  and\ \bibinfo {author} {\bibfnamefont {M.}~\bibnamefont {Khammash}},\ }\href
  {\doibase 10.1063/1.5006484} {\bibfield  {journal} {\bibinfo  {journal} {J.
  Chem. Phys.}\ }\textbf {\bibinfo {volume} {147}},\ \bibinfo {pages} {154101}
  (\bibinfo {year} {2017})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Dayar}\ \emph
  {et~al.}(2011{\natexlab{a}})\citenamefont {Dayar}, \citenamefont {Hermanns},
  \citenamefont {Spieler},\ and\ \citenamefont {Wolf}}]{Dayar2011}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {T.}~\bibnamefont
  {Dayar}}, \bibinfo {author} {\bibfnamefont {H.}~\bibnamefont {Hermanns}},
  \bibinfo {author} {\bibfnamefont {D.}~\bibnamefont {Spieler}}, \ and\
  \bibinfo {author} {\bibfnamefont {V.}~\bibnamefont {Wolf}},\ }\href {\doibase
  10.1002/nla.795} {\bibfield  {journal} {\bibinfo  {journal} {Numer. Linear
  Algebra Appl.}\ }\textbf {\bibinfo {volume} {18}},\ \bibinfo {pages} {931}
  (\bibinfo {year} {2011}{\natexlab{a}})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Spieler}(2014)}]{Spieler2014}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {D.}~\bibnamefont
  {Spieler}},\ }\emph {\bibinfo {title} {{Numerical analysis of long-run
  properties for Markov population models}}},\ \href@noop {} {Ph.D. thesis},\
  \bibinfo  {school} {Saarland University} (\bibinfo {year} {2014})\BibitemShut
  {NoStop}%
\bibitem [{\citenamefont {Dayar}\ \emph
  {et~al.}(2011{\natexlab{b}})\citenamefont {Dayar}, \citenamefont {Sandmann},
  \citenamefont {Spieler},\ and\ \citenamefont {Wolf}}]{Dayar2011a}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {T.}~\bibnamefont
  {Dayar}}, \bibinfo {author} {\bibfnamefont {W.}~\bibnamefont {Sandmann}},
  \bibinfo {author} {\bibfnamefont {D.}~\bibnamefont {Spieler}}, \ and\
  \bibinfo {author} {\bibfnamefont {V.}~\bibnamefont {Wolf}},\ }\href {\doibase
  10.1017/S0001867800005279} {\bibfield  {journal} {\bibinfo  {journal} {Adv.
  in Appl. Probab.}\ }\textbf {\bibinfo {volume} {43}},\ \bibinfo {pages}
  {1005} (\bibinfo {year} {2011}{\natexlab{b}})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Sakurai}\ and\ \citenamefont
  {Hori}(2017)}]{Sakurai2017}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {Y.}~\bibnamefont
  {Sakurai}}\ and\ \bibinfo {author} {\bibfnamefont {Y.}~\bibnamefont {Hori}},\
  }in\ \href {\doibase 10.1109/CDC.2017.8263820} {\emph {\bibinfo {booktitle}
  {IEEE 56th Annual Conference on Decision and Control}}}\ (\bibinfo {year}
  {2017})\ pp.\ \bibinfo {pages} {1206--1211}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Dowdy}\ and\ \citenamefont
  {Barton}(2018)}]{Dowdy2018}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {G.~R.}\ \bibnamefont
  {Dowdy}}\ and\ \bibinfo {author} {\bibfnamefont {P.~I.}\ \bibnamefont
  {Barton}},\ }\href {\doibase 10.1063/1.5009950} {\bibfield  {journal}
  {\bibinfo  {journal} {J. Chem. Phys.}\ }\textbf {\bibinfo {volume} {148}},\
  \bibinfo {pages} {084106} (\bibinfo {year} {2018})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Dowdy}\ and\ \citenamefont
  {Barton}(2017)}]{Dowdy2017}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {G.~R.}\ \bibnamefont
  {Dowdy}}\ and\ \bibinfo {author} {\bibfnamefont {P.~I.}\ \bibnamefont
  {Barton}},\ }\href {\doibase 10.1016/B978-0-444-63965-3.50375-5} {\bibfield
  {journal} {\bibinfo  {journal} {Computer Aided Chemical Engineering}\
  }\textbf {\bibinfo {volume} {40}},\ \bibinfo {pages} {2239} (\bibinfo {year}
  {2017})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Ghusinga}\ \emph {et~al.}(2017)\citenamefont
  {Ghusinga}, \citenamefont {Vargas-Garcia}, \citenamefont {Lamperski},\ and\
  \citenamefont {Singh}}]{Ghusinga2017a}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {K.~R.}\ \bibnamefont
  {Ghusinga}}, \bibinfo {author} {\bibfnamefont {C.~A.}\ \bibnamefont
  {Vargas-Garcia}}, \bibinfo {author} {\bibfnamefont {A.}~\bibnamefont
  {Lamperski}}, \ and\ \bibinfo {author} {\bibfnamefont {A.}~\bibnamefont
  {Singh}},\ }\href {\doibase 10.1088/1478-3975/aa75c6} {\bibfield  {journal}
  {\bibinfo  {journal} {Phys. Biol.}\ }\textbf {\bibinfo {volume} {14}},\
  \bibinfo {pages} {04LT01} (\bibinfo {year} {2017})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Kuntz}\ \emph
  {et~al.}(2018{\natexlab{a}})\citenamefont {Kuntz}, \citenamefont {Thomas},
  \citenamefont {Stan},\ and\ \citenamefont {Barahona}}]{Kuntz2018a}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {J.}~\bibnamefont
  {Kuntz}}, \bibinfo {author} {\bibfnamefont {P.}~\bibnamefont {Thomas}},
  \bibinfo {author} {\bibfnamefont {G.~B.}\ \bibnamefont {Stan}}, \ and\
  \bibinfo {author} {\bibfnamefont {M.}~\bibnamefont {Barahona}},\ }\href@noop
  {} {\  (\bibinfo {year} {2018}{\natexlab{a}})},\ \Eprint
  {http://arxiv.org/abs/1810.03658} {arXiv:1810.03658} \BibitemShut {NoStop}%
\bibitem [{\citenamefont {Schl{\"{o}}gl}(1972)}]{schlogl1972}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {F.}~\bibnamefont
  {Schl{\"{o}}gl}},\ }\href {\doibase 10.1007/BF01379769} {\bibfield  {journal}
  {\bibinfo  {journal} {Z. Phys. A}\ }\textbf {\bibinfo {volume} {253}},\
  \bibinfo {pages} {147} (\bibinfo {year} {1972})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Glynn}\ and\ \citenamefont
  {Zeevi}(2008)}]{Glynn2008}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {P.~W.}\ \bibnamefont
  {Glynn}}\ and\ \bibinfo {author} {\bibfnamefont {A.}~\bibnamefont {Zeevi}},\
  }in\ \href {\doibase 10.1214/074921708000000381} {\emph {\bibinfo {booktitle}
  {Markov Processes and Related Topics: A Festschrift for Thomas G. Kurtz}}}\
  (\bibinfo  {publisher} {Institute of Mathematical Statistics},\ \bibinfo
  {year} {2008})\ pp.\ \bibinfo {pages} {195--214}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Blekherman}\ \emph {et~al.}(2013)\citenamefont
  {Blekherman}, \citenamefont {Parrilo},\ and\ \citenamefont
  {Thomas}}]{Blekherman2013}%
  \BibitemOpen
  \bibinfo {editor} {\bibfnamefont {G.}~\bibnamefont {Blekherman}}, \bibinfo
  {editor} {\bibfnamefont {P.~A.}\ \bibnamefont {Parrilo}}, \ and\ \bibinfo
  {editor} {\bibfnamefont {R.~R.}\ \bibnamefont {Thomas}},\ eds.,\ \href
  {\doibase 10.1137/1.9781611972290.fm} {\emph {\bibinfo {title} {{Semidefinite
  optimization and convex algebraic geometry}}}}\ (\bibinfo  {publisher}
  {SIAM},\ \bibinfo {year} {2013})\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Lasserre}(2009)}]{Lasserre2009}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {J.~B.}\ \bibnamefont
  {Lasserre}},\ }\href {\doibase 10.1142/p665} {\emph {\bibinfo {title}
  {{Moments, positive polynomials and their applications}}}}\ (\bibinfo
  {publisher} {Imperial College Press},\ \bibinfo {year} {2009})\BibitemShut
  {NoStop}%
\bibitem [{\citenamefont {Kuntz}(2017)}]{Kuntzthe}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {J.}~\bibnamefont
  {Kuntz}},\ }\emph {\bibinfo {title} {{Deterministic approximation schemes
  with computable errors for the distributions of Markov chains}}},\ \href@noop
  {} {Ph.D. thesis},\ \bibinfo  {school} {Imperial College London} (\bibinfo
  {year} {2017})\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Laurent}(2009)}]{Laurent2009}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {M.}~\bibnamefont
  {Laurent}},\ }in\ \href {\doibase 10.1007/978-0-387-09686-5_7} {\emph
  {\bibinfo {booktitle} {Emerging applications of algebraic geometry}}}\
  (\bibinfo  {publisher} {Springer-Verlag},\ \bibinfo {year} {2009})\ pp.\
  \bibinfo {pages} {157--270}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Landau}(1987)}]{Landau1987}%
  \BibitemOpen
  \bibinfo {editor} {\bibfnamefont {H.}~\bibnamefont {Landau}},\ ed.,\ \href
  {\doibase 10.1090/psapm/037} {\emph {\bibinfo {title} {{Moments in
  mathematics}}}},\ Proceedings of Symposia in Applied Mathematics\ (\bibinfo
  {publisher} {American Mathematical Society},\ \bibinfo {year}
  {1987})\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Lasserre}(2002)}]{Lasserre2002a}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {J.~B.}\ \bibnamefont
  {Lasserre}},\ }\href {\doibase 10.1090/S0002-9947-01-02898-7} {\bibfield
  {journal} {\bibinfo  {journal} {Trans. Amer. Math. Soc.}\ }\textbf {\bibinfo
  {volume} {354}},\ \bibinfo {pages} {631} (\bibinfo {year}
  {2002})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Fantuzzi}\ \emph {et~al.}(2016)\citenamefont
  {Fantuzzi}, \citenamefont {Goluskin}, \citenamefont {Huang},\ and\
  \citenamefont {Chernyshenko}}]{Fantuzzi2016a}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {G.}~\bibnamefont
  {Fantuzzi}}, \bibinfo {author} {\bibfnamefont {D.}~\bibnamefont {Goluskin}},
  \bibinfo {author} {\bibfnamefont {D.}~\bibnamefont {Huang}}, \ and\ \bibinfo
  {author} {\bibfnamefont {S.~I.}\ \bibnamefont {Chernyshenko}},\ }\href
  {\doibase 10.1137/15M1053347} {\bibfield  {journal} {\bibinfo  {journal}
  {SIAM J. Appl. Dyn. Syst.}\ }\textbf {\bibinfo {volume} {15}},\ \bibinfo
  {pages} {1962} (\bibinfo {year} {2016})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Nakata}(2010)}]{Nakata2010}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {M.}~\bibnamefont
  {Nakata}},\ }in\ \href {\doibase 10.1109/CACSD.2010.5612693} {\emph {\bibinfo
  {booktitle} {EEE International Symposium on Computer-Aided Control System
  Design}}}\ (\bibinfo {year} {2010})\ pp.\ \bibinfo {pages}
  {29--34}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Zheng}\ \emph {et~al.}(2017)\citenamefont {Zheng},
  \citenamefont {Fantuzzi},\ and\ \citenamefont
  {Papachristodoulou}}]{Zheng2017}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {Y.}~\bibnamefont
  {Zheng}}, \bibinfo {author} {\bibfnamefont {G.}~\bibnamefont {Fantuzzi}}, \
  and\ \bibinfo {author} {\bibfnamefont {A.}~\bibnamefont
  {Papachristodoulou}},\ }\href@noop {} {\  (\bibinfo {year} {2017})},\ \Eprint
  {http://arxiv.org/abs/1708.04174} {arXiv:1708.04174} \BibitemShut {NoStop}%
\bibitem [{\citenamefont {L{\"{o}}fberg}(2004)}]{Lofberg2004}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {J.}~\bibnamefont
  {L{\"{o}}fberg}},\ }in\ \href {\doibase 10.1109/CACSD.2004.1393890} {\emph
  {\bibinfo {booktitle} {IEEE International Symposium on Computer Aided Control
  System Design}}}\ (\bibinfo {year} {2004})\ pp.\ \bibinfo {pages}
  {284--289}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Fantuzzi}\ and\ \citenamefont
  {Fuentes}(2016)}]{Fantuzzi2016}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {G.}~\bibnamefont
  {Fantuzzi}}\ and\ \bibinfo {author} {\bibfnamefont {F.}~\bibnamefont
  {Fuentes}},\ }\href@noop {} {\enquote {\bibinfo {title} {{mpYALMIP (Version
  1.1.2), https://github.com/giofantuzzi/mpYALMIP}},}\ } (\bibinfo {year}
  {2016})\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Henrion}\ \emph {et~al.}(2009)\citenamefont
  {Henrion}, \citenamefont {Lasserre},\ and\ \citenamefont
  {L{\"{o}}fberg}}]{Henrion2009}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {D.}~\bibnamefont
  {Henrion}}, \bibinfo {author} {\bibfnamefont {J.~B.}\ \bibnamefont
  {Lasserre}}, \ and\ \bibinfo {author} {\bibfnamefont {J.}~\bibnamefont
  {L{\"{o}}fberg}},\ }\href {\doibase 10.1080/10556780802699201} {\bibfield
  {journal} {\bibinfo  {journal} {Optim. Methods Softw.}\ }\textbf {\bibinfo
  {volume} {24}},\ \bibinfo {pages} {761} (\bibinfo {year} {2009})}\BibitemShut
  {NoStop}%
\bibitem [{\citenamefont {CPLEX}(2009)}]{CPLEX}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {I.~I.}\ \bibnamefont
  {CPLEX}},\ }\href@noop {} {\enquote {\bibinfo {title} {{V12.1: User's Manual
  for CPLEX}},}\ } (\bibinfo {year} {2009})\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Gardiner}(2009)}]{Gardiner2009}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {C.}~\bibnamefont
  {Gardiner}},\ }\href@noop {} {\emph {\bibinfo {title} {{Stochastic
  Methods}}}},\ \bibinfo {edition} {4th}\ ed.\ (\bibinfo  {publisher}
  {Springer-Verlag},\ \bibinfo {year} {2009})\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Norris}(1997)}]{Norris1997}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {J.~R.}\ \bibnamefont
  {Norris}},\ }\href {\doibase 10.1017/CBO9780511810633} {\emph {\bibinfo
  {title} {{Markov chains}}}}\ (\bibinfo  {publisher} {Cambridge University
  Press},\ \bibinfo {year} {1997})\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Gardner}\ \emph {et~al.}(2000)\citenamefont
  {Gardner}, \citenamefont {Cantor},\ and\ \citenamefont
  {Collins}}]{gardner2000}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {T.~S.}\ \bibnamefont
  {Gardner}}, \bibinfo {author} {\bibfnamefont {C.~R.}\ \bibnamefont {Cantor}},
  \ and\ \bibinfo {author} {\bibfnamefont {J.~J.}\ \bibnamefont {Collins}},\
  }\href {\doibase 10.1038/35002131} {\bibfield  {journal} {\bibinfo  {journal}
  {Nature}\ }\textbf {\bibinfo {volume} {403}},\ \bibinfo {pages} {339}
  (\bibinfo {year} {2000})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Perez-Carrasco}\ \emph {et~al.}(2016)\citenamefont
  {Perez-Carrasco}, \citenamefont {Guerrero}, \citenamefont {Briscoe},\ and\
  \citenamefont {Page}}]{perez2016}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {R.}~\bibnamefont
  {Perez-Carrasco}}, \bibinfo {author} {\bibfnamefont {P.}~\bibnamefont
  {Guerrero}}, \bibinfo {author} {\bibfnamefont {J.}~\bibnamefont {Briscoe}}, \
  and\ \bibinfo {author} {\bibfnamefont {K.~M.}\ \bibnamefont {Page}},\ }\href
  {\doibase 10.1371/journal.pcbi.1005154} {\bibfield  {journal} {\bibinfo
  {journal} {PLoS Comput. Biol.}\ }\textbf {\bibinfo {volume} {12}},\ \bibinfo
  {pages} {e1005154} (\bibinfo {year} {2016})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Paulev{\'{e}}}\ \emph {et~al.}(2014)\citenamefont
  {Paulev{\'{e}}}, \citenamefont {Craciun},\ and\ \citenamefont
  {Koeppl}}]{pauleve2014}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {L.}~\bibnamefont
  {Paulev{\'{e}}}}, \bibinfo {author} {\bibfnamefont {G.}~\bibnamefont
  {Craciun}}, \ and\ \bibinfo {author} {\bibfnamefont {H.}~\bibnamefont
  {Koeppl}},\ }\href {\doibase 10.1007/s00285-013-0686-2} {\bibfield  {journal}
  {\bibinfo  {journal} {J. Math. Biol.}\ }\textbf {\bibinfo {volume} {69}},\
  \bibinfo {pages} {55} (\bibinfo {year} {2014})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Gupta}\ and\ \citenamefont
  {Khammash}(2018)}]{Gupta2018}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {A.}~\bibnamefont
  {Gupta}}\ and\ \bibinfo {author} {\bibfnamefont {M.}~\bibnamefont
  {Khammash}},\ }\href {\doibase 10.1137/17M1134299} {\bibfield  {journal}
  {\bibinfo  {journal} {SIAM J. Appl. Dyn. Syst}\ }\textbf {\bibinfo {volume}
  {17}},\ \bibinfo {pages} {1213} (\bibinfo {year} {2018})}\BibitemShut
  {NoStop}%
\bibitem [{\citenamefont {Seneta}(1967)}]{Seneta1967}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {E.}~\bibnamefont
  {Seneta}},\ }\href {\doibase 10.1017/S0305004100042006} {\bibfield  {journal}
  {\bibinfo  {journal} {Proc. Camb. Phil. Soc.}\ }\textbf {\bibinfo {volume}
  {63}},\ \bibinfo {pages} {983} (\bibinfo {year} {1967})}\BibitemShut
  {NoStop}%
\bibitem [{\citenamefont {Tweedie}(1998)}]{Tweedie1998}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {R.~L.}\ \bibnamefont
  {Tweedie}},\ }\href {\doibase 10.1017/S0021900200016181} {\bibfield
  {journal} {\bibinfo  {journal} {J. Appl. Probab.}\ }\textbf {\bibinfo
  {volume} {35}},\ \bibinfo {pages} {517} (\bibinfo {year} {1998})}\BibitemShut
  {NoStop}%
\bibitem [{\citenamefont {Meyn}\ and\ \citenamefont
  {Tweedie}(1994)}]{Meyn1994}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {S.~P.}\ \bibnamefont
  {Meyn}}\ and\ \bibinfo {author} {\bibfnamefont {R.~L.}\ \bibnamefont
  {Tweedie}},\ }\href {\doibase 10.1214/aoap/1177004900} {\bibfield  {journal}
  {\bibinfo  {journal} {Ann. Appl. Probab.}\ }\textbf {\bibinfo {volume} {4}},\
  \bibinfo {pages} {981} (\bibinfo {year} {1994})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Liu}(2015)}]{Liu2015}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {Y.}~\bibnamefont
  {Liu}},\ }\href {\doibase 10.1007/s11425-015-5019-z} {\bibfield  {journal}
  {\bibinfo  {journal} {Sci. China Math.}\ }\textbf {\bibinfo {volume} {58}},\
  \bibinfo {pages} {2633} (\bibinfo {year} {2015})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Masuyama}(2017{\natexlab{a}})}]{Masuyama2017}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {H.}~\bibnamefont
  {Masuyama}},\ }\href {\doibase 10.15807/jorsj.60.271} {\bibfield  {journal}
  {\bibinfo  {journal} {J. Oper. Res. Soc. Japan}\ }\textbf {\bibinfo {volume}
  {60}},\ \bibinfo {pages} {271} (\bibinfo {year}
  {2017}{\natexlab{a}})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Masuyama}(2017{\natexlab{b}})}]{Masuyama2017a}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {H.}~\bibnamefont
  {Masuyama}},\ }\href {\doibase 10.1016/J.LAA.2016.10.014} {\bibfield
  {journal} {\bibinfo  {journal} {Linear Algebra Appl.}\ }\textbf {\bibinfo
  {volume} {514}},\ \bibinfo {pages} {105} (\bibinfo {year}
  {2017}{\natexlab{b}})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Liu}\ \emph {et~al.}(2018)\citenamefont {Liu},
  \citenamefont {Li},\ and\ \citenamefont {Masuyama}}]{Liu2018}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {Y.}~\bibnamefont
  {Liu}}, \bibinfo {author} {\bibfnamefont {W.}~\bibnamefont {Li}}, \ and\
  \bibinfo {author} {\bibfnamefont {H.}~\bibnamefont {Masuyama}},\ }\href
  {\doibase 10.1016/J.ORL.2018.05.001} {\bibfield  {journal} {\bibinfo
  {journal} {Oper. Res. Lett.}\ }\textbf {\bibinfo {volume} {46}},\ \bibinfo
  {pages} {409} (\bibinfo {year} {2018})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Liu}\ and\ \citenamefont {Li}(2018)}]{Liu2018a}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {Y.}~\bibnamefont
  {Liu}}\ and\ \bibinfo {author} {\bibfnamefont {W.}~\bibnamefont {Li}},\
  }\href {\doibase 10.1017/apr.2018.28} {\bibfield  {journal} {\bibinfo
  {journal} {Adv. Appl. Probab.}\ }\textbf {\bibinfo {volume} {50}},\ \bibinfo
  {pages} {645} (\bibinfo {year} {2018})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Fox}\ \emph {et~al.}(2016)\citenamefont {Fox},
  \citenamefont {Neuert},\ and\ \citenamefont {Munsky}}]{fox2016}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {Z.}~\bibnamefont
  {Fox}}, \bibinfo {author} {\bibfnamefont {G.}~\bibnamefont {Neuert}}, \ and\
  \bibinfo {author} {\bibfnamefont {B.}~\bibnamefont {Munsky}},\ }\href
  {\doibase 10.1063/1.4960505} {\bibfield  {journal} {\bibinfo  {journal} {J.
  Chem. Phys.}\ }\textbf {\bibinfo {volume} {145}},\ \bibinfo {pages} {74101}
  (\bibinfo {year} {2016})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Gillespie}(1976)}]{Gillespie1976}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {D.~T.}\ \bibnamefont
  {Gillespie}},\ }\href {\doibase 10.1016/0021-9991(76)90041-3} {\bibfield
  {journal} {\bibinfo  {journal} {J. Comput. Phys.}\ }\textbf {\bibinfo
  {volume} {22}},\ \bibinfo {pages} {403} (\bibinfo {year} {1976})}\BibitemShut
  {NoStop}%
\bibitem [{\citenamefont {Kendall}(1950)}]{Kendall1950}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {D.~G.}\ \bibnamefont
  {Kendall}},\ }\href@noop {} {\bibfield  {journal} {\bibinfo  {journal} {J. R.
  Stat. Soc. Ser. B Stat. Methodol.}\ }\textbf {\bibinfo {volume} {12}},\
  \bibinfo {pages} {116} (\bibinfo {year} {1950})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Feller}(1940)}]{Feller1940}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {W.}~\bibnamefont
  {Feller}},\ }\href {\doibase 10.1090/S0002-9947-1940-0002697-3} {\bibfield
  {journal} {\bibinfo  {journal} {Trans. Amer. Math. Soc.}\ }\textbf {\bibinfo
  {volume} {48}},\ \bibinfo {pages} {488} (\bibinfo {year} {1940})}\BibitemShut
  {NoStop}%
\bibitem [{\citenamefont {Kuntz}\ \emph
  {et~al.}(2018{\natexlab{b}})\citenamefont {Kuntz}, \citenamefont {Thomas},
  \citenamefont {Stan},\ and\ \citenamefont {Barahona}}]{Kuntz2018}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {J.}~\bibnamefont
  {Kuntz}}, \bibinfo {author} {\bibfnamefont {P.}~\bibnamefont {Thomas}},
  \bibinfo {author} {\bibfnamefont {G.~B.}\ \bibnamefont {Stan}}, \ and\
  \bibinfo {author} {\bibfnamefont {M.}~\bibnamefont {Barahona}},\ }\href@noop
  {} {\  (\bibinfo {year} {2018}{\natexlab{b}})},\ \Eprint
  {http://arxiv.org/abs/1801.09507} {arXiv:1801.09507} \BibitemShut {NoStop}%
\bibitem [{\citenamefont {Kurtz}\ and\ \citenamefont
  {Stockbridge}(1998)}]{Kurtz1998}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {T.~G.}\ \bibnamefont
  {Kurtz}}\ and\ \bibinfo {author} {\bibfnamefont {R.~H.}\ \bibnamefont
  {Stockbridge}},\ }\href {\doibase 10.1137/S0363012995295516} {\bibfield
  {journal} {\bibinfo  {journal} {SIAM J. Control Optim.}\ }\textbf {\bibinfo
  {volume} {36}},\ \bibinfo {pages} {609} (\bibinfo {year} {1998})}\BibitemShut
  {NoStop}%
\bibitem [{\citenamefont {Reuter}(1957)}]{Reuter1957}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {G.~E.~H.}\
  \bibnamefont {Reuter}},\ }\href {\doibase 10.1007/BF02392391} {\bibfield
  {journal} {\bibinfo  {journal} {Acta Math.}\ }\textbf {\bibinfo {volume}
  {97}},\ \bibinfo {pages} {1} (\bibinfo {year} {1957})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Chen}(1991)}]{Chen1991}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {M.}~\bibnamefont
  {Chen}},\ }\href {\doibase 10.2307/3214868} {\bibfield  {journal} {\bibinfo
  {journal} {J. Appl. Probab.}\ }\textbf {\bibinfo {volume} {28}},\ \bibinfo
  {pages} {305} (\bibinfo {year} {1991})}\BibitemShut {NoStop}%
\bibitem [{\citenamefont {Meyn}\ and\ \citenamefont
  {Tweedie}(1993{\natexlab{b}})}]{Meyn1993b}%
  \BibitemOpen
  \bibfield  {author} {\bibinfo {author} {\bibfnamefont {S.~P.}\ \bibnamefont
  {Meyn}}\ and\ \bibinfo {author} {\bibfnamefont {R.~L.}\ \bibnamefont
  {Tweedie}},\ }\href {\doibase 10.2307/1427522} {\bibfield  {journal}
  {\bibinfo  {journal} {Adv. Appl. Probab.}\ }\textbf {\bibinfo {volume}
  {25}},\ \bibinfo {pages} {518} (\bibinfo {year}
  {1993}{\natexlab{b}})}\BibitemShut {NoStop}%
\end{thebibliography}%
