A finite state projection algorithm for the stationary solution of the chemical master equation

TL;DR

This paper introduces the stationary Finite State Projection (sFSP) algorithm, enabling efficient approximation of stationary solutions of the chemical master equation in systems biology, even for very large state-spaces.

Contribution

The paper develops the sFSP method, extending FSP to accurately estimate stationary distributions of CMEs using finite linear systems and tensor train techniques.

Findings

sFSP provides accurate stationary distribution approximations

Error bounds can be minimized by expanding the truncated state-space

Efficiently solves problems with over 100 million states

Abstract

The chemical master equation (CME) is frequently used in systems biology to quantify the effects of stochastic fluctuations that arise due to biomolecular species with low copy numbers. The CME is a system of ordinary differential equations that describes the evolution of probability density for each population vector in the state-space of the stochastic reaction dynamics. For many examples of interest, this state-space is infinite, making it difficult to obtain exact solutions of the CME. To deal with this problem, the Finite State Projection (FSP) algorithm was developed by Munsky and Khammash (Jour. Chem. Phys. 2006), to provide approximate solutions to the CME by truncating the state-space. The FSP works well for finite time-periods but it cannot be used for estimating the stationary solutions of CMEs, which are often of interest in systems biology. The aim of this paper is to develop a version of FSP which we refer to as the stationary FSP (sFSP) that allows one to obtain accurate approximations of the stationary solutions of a CME by solving a finite linear-algebraic system that yields the stationary distribution of a continuous-time Markov chain over the truncated state-space. We derive bounds for the approximation error incurred by sFSP and we establish that under certain stability conditions, these errors can be made arbitrarily small by appropriately expanding the truncated state-space. We provide several examples to illustrate our sFSP method and demonstrate its efficiency in estimating the stationary distributions. In particular, we show that using a quantised tensor train (QTT) implementation of our sFSP method, problems admitting more than 100 million states can be efficiently solved.

Tables

Table 1. Table 1: Application of sFSP on the gene-expression network. The transition rate matrix Q ¯ i subscript ¯ 𝑄 𝑖 \overline{Q}_{i} is constructed in C++ while its stationary distribution is found in Matlab.

Iteration	Cut-offs		State-space size	Convergence factor	CPU Time (seconds)
$i$	$C_{l,i}$	$C_{r,i}$	$n_i$	${\gamma }_i=r_{\text{out}}^{\left(i\right)}{C}_{r,i}$	Constructing ${\overline{Q}}_i$	Finding ${\overline{\pi }}_i$
$1$	$1860$	$2340$	$1,008,240$	$2.541\times {10}^3$	13.6	18.7
$2$	$1620$	$2580$	$2,016,480$	$0.278$	27.5	46.4
$3$	$1380$	$2820$	$3,024,720$	$7.473\times {10}^{-5}$	40.4	75.6
$4$	$1140$	$3060$	$4,032,960$	$2.292\times {10}^{-9}$	53.6	110.7
$5$	$900$	$3300$	$5,041,200$	$7.336\times {10}^{-15}$	68.3	159.6

Table 2. Table 2: Application of sFSP on the Toggle-Switch network. The transition rate matrix Q ¯ i subscript ¯ 𝑄 𝑖 \overline{Q}_{i} is constructed in C++ while its stationary distribution is found in Matlab.

Iteration	Cut-offs		State-space size	Convergence factor	CPU Time (seconds)
$i$	$C_{l,i}$	$C_{r,i}$	$n_i$	${\gamma }_i=r_{\text{out}}^{\left(i\right)}{C}_{r,i}$	Constructing ${\overline{Q}}_i$	Finding ${\overline{\pi }}_i$
$1$	$860$	$1360$	$555,250$	$1.03\times {10}^3$	7.6	13.1
$2$	$610$	$1610$	$1,110,500$	$3.87\times {10}^2$	14.8	24.9
$3$	$360$	$1860$	$1,665,750$	$4.14\times {10}^1$	23.3	37.2
$4$	$110$	$2110$	$2,221,000$	$1.37\times {10}^{-1}$	30.63	55.9
$5$	$0$	$2360$	$2,785,980$	$1.51\times {10}^{-53}$	37.66	77.1

Table 3. Table 3: Reactions for the Pap-Switch. Here [LRP] = 100 [LRP] 100 \textnormal{[LRP]}=100 denotes the total number of LRP molecues and x = ( x 1 , … , x 5 ) 𝑥 subscript 𝑥 1 … subscript 𝑥 5 x=(x_{1},\dots,x_{5}) denotes the copy-numbers of the five species ordered as G 1 , G 2 , G 3 , G 4 subscript 𝐺 1 subscript 𝐺 2 subscript 𝐺 3 subscript 𝐺 4 G_{1},G_{2},G_{3},G_{4} and P ​ a ​ p ​ I 𝑃 𝑎 𝑝 𝐼 PapI . Propensities of reactions 2 , 4 , 6 2 4 6 2,4,6 and 8 8 8 contain a term for the repression of LRP unbinding by P ​ a ​ p ​ I 𝑃 𝑎 𝑝 𝐼 PapI molecules.

Table 4. Table 4: Application of sFSP on the Pap-Switch network. The transition rate matrix Q ¯ i subscript ¯ 𝑄 𝑖 \overline{Q}_{i} is constructed in C++ while its stationary distribution is found in Matlab.

Iteration	Cut-offs		State-space size	Convergence factor	CPU Time (seconds)
$i$	$C_{l,i}$	$C_{r,i}$	$n_i$	${\gamma }_i=r_{\text{out}}^{\left(i\right)}{C}_{r,i}$	Constructing ${\overline{Q}}_i$	Finding ${\overline{\pi }}_i$
$1$	$0$	$10$	$44$	$7.72\times {10}^{-1}$	0.00062	0.0344
$2$	$0$	$16$	$68$	$9.64\times {10}^{-2}$	0.00117	0.0566
$3$	$0$	$22$	$92$	$1.75\times {10}^{-2}$	0.001719	0.0517
$4$	$0$	$28$	$116$	$6.26\times {10}^{-6}$	0.002906	0.0519
$5$	$0$	$34$	$140$	$6.31\times {10}^{-9}$	0.003498	0.0464
$6$	$0$	$40$	$164$	$2.23\times {10}^{-12}$	0.00255	0.0486

Table 5. Table 5: Reactions for the triple-repressor model. Here x = ( x 1 , … , x 9 ) 𝑥 subscript 𝑥 1 … subscript 𝑥 9 x=(x_{1},...,x_{9}) denotes the copy-numbers of the 9 network species ordered as G A 1 subscript superscript 𝐺 1 𝐴 G^{1}_{A} , G B 1 subscript superscript 𝐺 1 𝐵 G^{1}_{B} , G C 1 subscript superscript 𝐺 1 𝐶 G^{1}_{C} , M A subscript 𝑀 𝐴 M_{A} , M B subscript 𝑀 𝐵 M_{B} , M C subscript 𝑀 𝐶 M_{C} , P A subscript 𝑃 𝐴 P_{A} , P B subscript 𝑃 𝐵 P_{B} and P C subscript 𝑃 𝐶 P_{C} . Note that G A 0 subscript superscript 𝐺 0 𝐴 G^{0}_{A} is the species denoting that Gene A is in the OFF state and hence its copy-number is simply ( 1 − x 1 ) 1 subscript 𝑥 1 (1-x_{1}) . The interpretation for species G B 0 subscript superscript 𝐺 0 𝐵 G^{0}_{B} and G C 0 subscript superscript 𝐺 0 𝐶 G^{0}_{C} is similar.

No.	Reaction	Propensity
1	$G_A^0⟶G_A^1$	${\lambda }_1\left(x\right)=\left(10+1.5x_7\right)\left(1-x_1\right)$
2	$G_A^1⟶G_A^0$	${\lambda }_2\left(x\right)=\left(7+2x_9\right)x_1$
3	$G_B^0⟶G_B^1$	${\lambda }_3\left(x\right)=\left(9+4x_8\right)\left(1-x_2\right)$
4	$G_B^1⟶G_B^0$	${\lambda }_4\left(x\right)=\left(10+4x_7\right)x_2$
5	$G_C^0⟶G_C^1$	${\lambda }_5\left(x\right)=\left(11+1.5x_9\right)\left(1-x_3\right)$
6	$G_C^1⟶G_C^0$	${\lambda }_6\left(x\right)=\left(9+2x_8\right)x_3$
7	$G_A^1⟶G_A^1+M_A$	${\lambda }_7\left(x\right)=1.5x_1$
8	$G_B^1⟶G_B^1+M_B$	${\lambda }_8\left(x\right)=1x_2$
9	$G_C^1⟶G_C^1+M_C$	${\lambda }_9\left(x\right)=1.1x_3$
10	$M_A⟶\varnothing$	${\lambda }_{10}\left(x\right)=0.5x_4$
11	$M_B⟶\varnothing$	${\lambda }_{11}\left(x\right)=0.3x_5$
12	$M_C⟶\varnothing$	${\lambda }_{12}\left(x\right)=0.425x_6$
13	$M_A⟶M_A+P_A$	${\lambda }_{13}\left(x\right)=9.5x_4$
14	$M_B⟶M_B+P_B$	${\lambda }_{14}\left(x\right)=11x_5$
15	$M_C⟶M_C+P_C$	${\lambda }_{15}\left(x\right)=10x_6$
16	$P_A⟶\varnothing$	${\lambda }_{16}\left(x\right)=14.5x_7$
17	$P_B⟶\varnothing$	${\lambda }_{17}\left(x\right)=15x_8$
18	$P_C⟶\varnothing$	${\lambda }_{18}\left(x\right)=11x_9$

Table 6. Table 6: Application of sFSP on the triple-repressor model.

Iteration	Upper bounds		State-space size	Convergence factor	CPU Time (minutes)
$i$	mRNAs $U_{m,i}$	Proteins $U_{p,i}$	$\left|{ℰ}_i\right|$	${\gamma }_i=r_{\text{out}}^{\left(i\right)}{U}_{p,i}$	t
$1$	$4$	$4$	$32,768$	$13.5607$	3.66
$2$	$8$	$4$	$262,144$	$47.3662$	6.77
$3$	$8$	$8$	$2,097,152$	$2.4899$	27.67
$4$	$16$	$8$	$16,777,216$	$5.2869$	60.09
$5$	$16$	$16$	$134,217,728$	$0.0036$	69.66

Table 7. Table 7: Comparison of the sFSP estimated stationary distribution π ¯ ¯ 𝜋 \overline{\pi} and the SSA estimated stationary distribution π ^ ^ 𝜋 \widehat{\pi} for the triple-repressor model. Computed ℓ 1 subscript ℓ 1 \ell_{1} distance ‖ π ¯ − π ^ ‖ ℓ 1 subscript norm ¯ 𝜋 ^ 𝜋 subscript ℓ 1 \|\overline{\pi}-\widehat{\pi}\|_{\ell_{1}} and CPU times to generate SSA samples are shown for three sample sizes 10 5 , 10 6 superscript 10 5 superscript 10 6 10^{5},10^{6} and 10 7 superscript 10 7 10^{7} .

No. of SSA samples	${\Vert \overline{\pi }-\hat{\pi }\Vert }_{{\ell }_1}$	CPU Time
${10}^5$	0.5969	12 minutes
${10}^6$	0.2461	117 minutes
${10}^7$	0.091	1076 minutes

Equations

$$Q^T\pi =0,$$

$$\sum _{i=1}^d{\nu }_{ik}{\mathbf{X}}_i⟶\sum _{i=1}^d{\nu }_{ik}^{\prime }{\mathbf{X}}_i.$$

$${\lambda }_k(x_1,\dots ,x_d)={\theta }_k\prod _{i=1}^d\frac{x_i(x_i-1)\dots (x_i-{\nu }_{ik}+1)}{{\nu }_{ik}!},$$

$$\mathbb{Q}f(x)=\sum _{k=1}^K{\lambda }_k(x)\left(f(x+{\zeta }_k)-f(x)\right),$$

$$\text{ for each }x\in \mathcal{E}\text{ and }k=1,\dots ,K,\text{if}{\lambda }_k(x)>0\text{ then }(x+{\zeta }_k)\in \mathcal{E}.$$

$$\displaystyle Q_{ij}=\left\{\begin{array}[]{cc}-\sum_{k=1}^{K}\lambda_{k}(x_{i})&\textnormal{ if }i=j\\ \lambda_{k}(x_{i})&\textnormal{ if }x_{j}=x_{i}+\zeta_{k}\textnormal{ for some }k\\ 0&\textnormal{ otherwise}.\end{array}\right.$$

$$p(t,x)=\mathbb{P}\left(X(t)=x\right)$$

$$\frac{dp(t,x)}{dt}=$$

$$p(t,A)=\sum _{x\in A}p(t,x)$$

$$p(t)=(p(t,x_0),p(t,x_1),p(t,x_2),\dots )$$

$$\frac{dp}{dt}=Q^{T}p(t).$$

$$p(t)=\exp (Q^{T}t)p(0)\text{for any}t\ge 0,$$

$$\underset{t\to \infty }{\lim }\Vert p(t)-\pi {\Vert }_{{\ell }_1}=0,$$

$$\Vert p(t)-\pi {\Vert }_{{\ell }_1}\le C{e}^{-\rho t}.$$

$$\frac{d{p}_n}{dt}=Q_n^T{p}_n(t).$$

$${ϵ}_n(t):=1-1^T{p}_n(t)=1-1^T\exp (Q_n^{T}t)p_n(0)\ge 0.$$

$$\displaystyle\widetilde{Q}_{n}=\left[\begin{array}[]{cc}Q_{n}&c_{n}\\ {\bf 0}&0\end{array}\right],$$

$$c_{n,i}=\sum _{k=1,(x_{j_i}+{\zeta }_k)\notin {\mathcal{E}}_n}^K{\lambda }_k(x_{j_i}).$$

$${\tilde{p}}_n(t)=(p_n(t),{ϵ}_n(t)),$$

$${\overline{Q}}_n=Q_n+c_n{b}_l,$$

$$r_{\text{out}}^{(n)}=c_n^T{\overline{\pi }}_n.$$

$$\mathcal{E}={\mathcal{E}}_b\times {\mathbb{N}}_0^{d_f},$$

$$\mathbb{Q}V(x)\le C_1-C_{2}V(x),$$

$$V(x)=1+⟨v,x⟩,$$

$$\sum _{k=1}^K{\lambda }_k(x)⟨v,{\zeta }_k⟩\le C_1-C_2(1+⟨v,x⟩)\text{ for all }x\in \mathcal{E}.$$

$$\sum _{k=1}^K{\lambda }_k(x)⟨v,{\zeta }_k{⟩}^2\le C_3+C_4(1+⟨v,x⟩)\text{ for all }x\in \mathcal{E}.$$

$$\Vert \mu {\Vert }_V=\sum _{x\in \mathcal{E}}|\mu (x)|V(x).$$

$$\mathcal{B}({\mathcal{E}}_n)=\{x\in {\mathcal{E}}_n:{\lambda }_k(x)>0\text{ and }(x+{\zeta }_k)\notin {\mathcal{E}}_n\text{ for some }k=1,\dots ,K\}.$$

$${\gamma }_V^{(n)}=r_{\text{out}}^{(n)}\Vert {\mathcal{E}}_n{\Vert }_V,$$

$$\Vert {\mathcal{E}}_n{\Vert }_V=V(x_{\ell })+\underset{x\in \mathcal{B}({\mathcal{E}}_n)}{\max }V(x)$$

Full text

Keywords: stochastic reaction networks; the Chemical Master Equation; Finite State Projection; stationary distribution; ergodicity; irreducibility; tensor trains

Mathematical Subject Classification (2010): 60J22; 60J27; 60H35; 65C40; 92E20

1 Introduction

Many intracellular reaction networks consist of biomolecular species that are typically present in low copy-numbers. The reactions involving these species fire intermittently at random times, rather than continuously. Hence deterministic descriptions of the reaction dynamics, based on Ordinary Differential Equations (ODEs), become highly inaccurate [1]. It is now well-known that macroscopic properties of the system can be heavily influenced by the intrinsic noise or randomness that arises due to the random timing of reactions [2]. Consequently stochastic formulations of the reaction dynamics, based on continuous-time Markov chains (CTMCs), has become a popular approach for studying the effects of intrinsic noise [3]. In this paper we provide a tool for the analysis of such models.

In the CTMC model of a reaction network, the state at any time is the vector of copy-number counts of all the species. When the number of network species is $d$, the dynamics evolves on a discrete state-space $\mathcal{E}$ which is a subset of the $d$-dimensional nonnegative integer lattice $\mathbb{N}^{d}_{0}$ and this subset must be large enough to include all the states that are accessible by the random dynamics. The effects of intrinsic noise on the reaction network are generally studied using the probability distribution $p(t)$ of the random state-vector $X(t)$ at time $t$. It is known that the time-evolution of this probability distribution is given by a system of coupled ODEs, known as the Chemical Master Equation (CME) in the literature (see (2.7)). For each state in $\mathcal{E}$ there is an ODE in the CME that captures the inflow and outflow of probability at that state. If this state-space $\mathcal{E}$ is finite, then the CME is a finite system of linear ODEs which can in principle be solved to yield the probability distribution $p(t)$. However in many examples of biological interest, the state-space $\mathcal{E}$ is infinite, making the CME impossible to solve. A common approach in such cases is to estimate the CME solution by computing the empirical distribution of the samples obtained by simulating the CTMC using Monte Carlo methods such as Gillespie’s Stochastic Simulation Algorithm (SSA) [4]. This simulation-based approach can be very time-consuming and the estimates suffer from statistical errors due to finitely many samples being used. In particular the low-probability events are sparsely sampled by Monte Carlo simulations, which can lead to incorrect representations of the CME solution. Such problems can be avoided by using the Finite State Projection (FSP) method developed by Munsky and Khammash [5], that directly solves the CME by truncating the state-space $\mathcal{E}$ to a manageable size (see Section 2.2). The solution obtained is approximate but FSP provides an iterative way to ensure that the approximation error is within some pre-specified tolerance level.

The truncated state-space needed by FSP to solve the CME accurately is still exorbitantly large for many problems of interest. For example, consider a simple gene-expression network where ten protein species are interacting with each other. Typically each protein in a cell has copy-numbers of the order of several thousands. So even if we have a conservative upper-bound of 1000 on the copy-number of each protein, the size of the state-space required for FSP is of the order of $10^{30}$, which is beyond the computational and storage capacity of modern day computers. This combinatorial explosion in the state-space size is often called the “curse of dimensionality” and it presents a major challenge in making the CME practically solvable. Several advanced numerical techniques have been developed that address this challenge by adapting the FSP. These techniques include Krylov Subspace approximations [6], Tensor-Train representations [7], and using sparse grids and aggregation methods [8]. Unlike these methods which attempt to solve the exact version of CME, there also exist a body of methods that aim to solve simplified versions of CME, which are derived by approximating the CTMC dynamics by a Stochastic Differential Equation (SDE) or a Piecewise-deterministic Markov Process (PDMP) (see [9, 10, 11]). Such dynamical approximations only hold for finite time-periods, and the assumptions on species copy-numbers and reaction propensities they require, are not always satisfied by networks encountered in systems biology.

For many biological applications, one is interested in the steady-state behavior which is captured by the stationary probability distribution $\pi$ to which the solution $p(t)$ of the CME converges to as $t\to\infty$. For CTMCs whose state-space $\mathcal{E}$ is finite and not too large, estimation of the stationary distribution $\pi$ is a simple linear-algebraic problem (see (1.1)). However in situations where the state-space is very large or infinite, this linear-algebraic problem cannot be practically solved, and we need to estimate $\pi$ by other means. The methods mentioned above for estimating $p(t)$ only work over finite time-intervals and they would generally fail to provide an accurate estimate of the stationary distribution $\pi$. The reason for this failure depends on the method being used. The dynamical approximations based on PDMPs or SDEs introduce an error that can become unbounded in the limit $t\to\infty$, and the Monte Carlo simulation based approach for estimating $\pi$ is highly undesirable due to statistical errors and the computational costs associated with these simulations over large time-intervals. The FSP algorithm also cannot be used for estimating $\pi$ because this method introduces an absorbing state to catch all the transitions that leave the truncated state-space (see Figure 1B). However in the limit $t\to\infty$, all the probability mass flows into this absorbing state, and so the obtained probability distributions are unable to capture the true stationary distribution $\pi$. We revisit this point later in this section and also explain it in detail in Section 2.2.

The aim of this paper is to present a FSP-like method for accurately estimating the stationary distribution $\pi$. This method also involves truncating the state-space but rather than solving a linear system of ODEs for probabilities over the truncated state-space (as in FSP), our method estimates the true stationary distribution $\pi$ by computing the stationary distribution of a suitably defined CTMC over the truncated state-space. As the latter step can be accomplished by solving a linear-algebraic system, rather than a system of ODEs, the computational complexity of our method is much lower than that of FSP. Consequently it can be successfully applied on a larger class of networks. We call our method the stationary Finite State Projection (or sFSP) algorithm and we provide theoretical arguments to establish its accuracy under certain stability conditions which are usually satisfied by networks in systems and synthetic biology. Even though sFSP can be applied on larger systems than FSP, the combinatorial explosion of state-space sizes still limits the range of applicability of sFSP severely. As was the case with FSP, this issue can be somewhat resolved by adapting sFSP to work with quantized tensor-train (QTT) representations [7], sparse grids and aggregation methods [8]. We illustrate this point with a computational example where sFSP is applied to the QTT representation of the CME (see Section 5). We remark here that QTT representations have already been successfully employed for obtaining approximations of the stationary distribution for reaction networks satisfying certain graph-theoretic criteria [12, 13]. However these criteria are highly-restrictive and it will become evident that the sFSP based approach is more versatile.

We now describe the problem of estimating stationary distributions in more detail. Henceforth let $|A|$ denote the size of any set $A$, and let ${\bf 0}$ and ${\bf 1}$ denote the vector of all zeros and all ones respectively111The dimension of these vectors will be clear from the context.. The stochastic model of a reaction network (see Section 2.1) represents the dynamics as a CTMC over a discrete state-space $\mathcal{E}\subset\mathbb{N}^{d}_{0}$. Such a CTMC can be described by its $|\mathcal{E}|\times|\mathcal{E}|$ transition rate matrix $Q$ (see [14]), whose diagonal entries are non-positive, off-diagonal entries are non-negative and all the rows sum to zero. The stationary distribution for this CTMC can be described by a non-negative vector222Throughout the paper we assume that vector and matrix indices start from [math] rather than $1$. $\pi=(\pi_{0},\pi_{1},\dots)$, which is in the left null-space of transition rate matrix $Q$, i.e.

[TABLE]

and its components sum to $1$ (i.e. ${\bf 1}^{T}\pi=\sum_{i}\pi_{i}=1$). Such a stationary distribution may not be unique and if $|\mathcal{E}|=\infty$ then it may not even exist (see [14]). In our recent work, we have dealt extensively with the issue of computationally verifying the existence and uniqueness of the stationary distribution corresponding to the CTMC models for a large class of biomolecular reaction networks (see [15] and [16]). Assuming that the existence and uniqueness of the stationary distribution $\pi$ has been ascertained for the network, our aim here is to estimate $\pi$ numerically. We are primarily interested in situations where $\mathcal{E}$ is infinite, and so the direct computation of $\pi$ using (1.1) is computationally impossible.

It is natural to try to estimate $\pi$ by solving a finite, truncated version of the linear-algebraic system (1.1). This truncated version can be obtained by first identifying a truncated state-space and then projecting the CTMC dynamics on this truncated state-space. Thereafter the stationary distribution for the projected CTMC, found by solving the corresponding linear-algebraic system of the form (1.1), serves as an estimate of the true stationary distribution $\pi$. An important issue that arises here is how to handle the outgoing transitions from the truncated state-space, so that the obtained estimate of $\pi$ is accurate. In the FSP approach [5], these outgoing transitions are preserved but their target states are collapsed into a single absorbing state (see Figure 1B). This leads to the “probability leakage” problem which can be managed over finite time-intervals but not in the asymptotic $t\to\infty$ regime. This problem manifests itself in the fact that the only stationary distribution for the projected CTMC would be the one that puts all the probability-mass at the absorbing state. Obviously this does not capture the true stationary distribution and hence modifications to the FSP approach are necessary to circumvent the probability-leakage problem. One such modification that has been tried is motivated by the use of “reflected” boundary conditions in the study of Fokker-Planck equations [12]. In this approach all the outgoing transitions from the truncated state-space are simply eliminated by setting their propensities to zero. It has been shown that this reflected version of FSP yields accurate estimates of the stationary distribution for some reaction network examples [17, 12]. However there is no theoretical guarantee that this approach will work well in general.

The method sFSP that we present in this paper modifies the FSP in another way. It preserves the outgoing transitions from the truncated state-space, but rather than channeling them to an absorbing state (as in FSP), it redirects them to a designated state within the truncated state-space (see Figure 1C). This modification is simple to implement and its appealing feature is that for a wide range of biomolecular reaction networks, we can theoretically guarantee that the stationary distribution of the projected CTMC converges to the actual stationary distribution $\pi$ as the truncated state-space expands to the full state-space $\mathcal{E}$. Moreover we derive bounds for the approximation error incurred by this approach, in terms of the outflow rate of all the outgoing transitions evaluated at the estimated stationary distribution (see Theorem 3.1). These results provide the theoretical basis for our method which expands the truncated state-space iteratively to recover a “good” approximation of $\pi$. Note that our approach for estimating the stationary distribution is very different from the stochastic complementation approach proposed in [18]. This complementation approach is generally difficult to implement for infinite state-space CTMCs and it only yields the conditional stationary distribution which can then be used to derive upper and lower bounds for the true stationary probabilities. However such bounds are not guaranteed to be close to each other. In contrast, our method allows one to estimate the true stationary probabilities directly.

For our method sFSP to work we require that the original CTMC representing the reaction network satisfies a couple of stability conditions. The first condition is that the state-space $\mathcal{E}$ needs to be irreducible i.e. all the states in $\mathcal{E}$ must be accessible from each other via a sequence of positive-propensity reactions. The second condition is a Foster-Lyapunov criterion (see [19]) which ensures that the original CTMC is exponentially ergodic i.e. the solution $p(t)$ of the CME converges to the stationary distribution $\pi$ exponentially fast. We elaborate these stability conditions in Section 3.1 and there we also explain how these conditions can be easily checked for a wide range of networks arising in systems and synthetic biology, using the computational procedures developed in our recent papers [15] and [16]. This makes the proposed sFSP method broadly applicable and of interest to the growing community of researchers working with stochastic models of biomolecular reaction networks.

This paper is organized as follows. In Section 2 we describe the stochastic model and the original FSP method [5]. In Section 3 we present and mathematically analyze our stationary Finite State Projection (or sFSP) algorithm. A simple implementation of sFSP is presented in Section 4 while its QTT implementation is presented in Section 5. These sections also include the computational examples that illustrate the respective implementations. Finally in Section 6 we conclude and discuss directions for future research.

2 Preliminaries

2.1 The stochastic model

We now formally describe the CTMC model of a reaction network. Suppose this network has $d$ species, called $\mathbf{X}_{1},\dots,\mathbf{X}_{d}$, and $K$ reactions of the form

[TABLE]

Here $\nu_{ik}$ and $\nu^{\prime}_{ik}$ are nonnegative integers denoting the number of molecules of species $\mathbf{X}_{i}$ that are consumed and produced by the $k$-th reaction. The state of the system at any time is the vector $x=(x_{1},\dots,x_{d})\in\mathbb{N}^{d}_{0}$ of molecular counts of all the $d$ species. When the $k$-th reaction fires, the state is displaced by the stoichiometric vector $\zeta_{k}\in\mathbb{Z}^{d}$ whose $i$-th component is $\zeta_{ik}=(\nu^{\prime}_{ik}-\nu_{ik})$. At any state $x$, the rate of the $k$-th reaction is $\lambda_{k}(x)$, where $\lambda_{k}:\mathbb{N}^{d}_{0}\to[0,\infty)$ is the propensity function for this reaction. Commonly mass action kinetics (see [3]) is assumed, where each $\lambda_{k}$ is given by

[TABLE]

with the positive parameter $\theta_{k}$ being the associated rate constant. We model the reaction dynamics as a CTMC which jumps from state $x$ after a random waiting time which is exponentially distributed with rate $\lambda_{0}(x):=\sum_{k=1}^{K}\lambda_{k}(x)$, and this jump is in direction $\zeta_{k}$ with probability $\lambda_{k}(x)/\lambda_{0}(x)$. Formally this CTMC can be specified by its generator333The generator of a Markov process is an operator which specifies the rate of change of the probability distribution of the process (see Chapter 4 in [20] for details). $\mathbb{Q}$ defined as

[TABLE]

where $f$ is any bounded real-valued function on $\mathbb{N}^{d}_{0}$.

From now on we suppose that there is a nonempty state-space $\mathcal{E}\subset\mathbb{N}^{d}_{0}$ on which the CTMC evolves i.e.

[TABLE]

In other words, if at state $x\in\mathcal{E}$, reaction $k$ has a positive probability of firing then the resulting state $(x+\zeta_{k})$ must also be in $\mathcal{E}$. As $\mathcal{E}$ is at most countable, it can be enumerated. This means that we can find a one-to-one and onto map $\phi$ from $\mathcal{E}$ to the set $\{0,1,\dots,|\mathcal{E}|-1\}$. Once such an enumeration is fixed, the set $\mathcal{E}$ can be expressed as $\mathcal{E}=\{x_{0},x_{1},\dots\}$, where $x_{i}=\phi^{-1}(i)$. Then the CTMC generator $\mathbb{Q}$ can be expressed as the transition rate matrix $Q=[Q_{ij}]$ given by444Here we assume for convenience that all stoichiometry vectors ($\zeta_{k}$-s) are distinct.

[TABLE]

Let $(X(t))_{t\geq 0}$ be the CTMC with this transition rate matrix and some initial state $X(0)\in\mathcal{E}$. For any state $x\in\mathcal{E}$, let

[TABLE]

be the probability that the CTMC is in state $x$ at time $t$. These probabilities evolve in time according to the Chemical Master Equation (CME) given by

[TABLE]

for each $x\in\mathcal{E}$. Note that this system has as many ODEs as the number of elements in the state-space $\mathcal{E}$, which is generally infinite or very large.

Let $p(t)$ be the probability distribution defined by

[TABLE]

for any $A\subset\mathcal{E}$. The vectorized form of $p(t)$ w.r.t. enumeration $\phi$ is simply given by

[TABLE]

and using this form we can express the CME as

[TABLE]

If the number of states in $\mathcal{E}$ is finite, then this first-order system can in principle be solved by exponentiating the matrix $Q^{T}$, i.e. the solution is given by

[TABLE]

where $p(0)$ is the vectorized form of the probability distribution of the initial state $X(0)$. However this approach is infeasible for large state-spaces and in such cases, the Finite State Projection (FSP) method [5] can be used to approximately solve the CME (see Section 2.2).

In many biological applications, rather than the finite-time behavior, one is interested in the properties of the system after it has settled down, or in other words, the CTMC $(X(t))_{t\geq 0}$ has reached a steady-state which is characterized by a stationary distribution $\pi$ satisfying (1.1), that is essentially a fixed-point for the CME (2.9). We say that the CTMC $(X(t))_{t\geq 0}$ is ergodic if this fixed-point is unique and globally attracting in the sense that for any initial probability distribution $p(0)$, the solution $p(t)$ of (2.9) satisfies

[TABLE]

where $\|p(t)-\pi\|_{\ell_{1}}=\sum_{x\in\mathcal{E}}|p(t,x)-\pi(x)|$ denotes the $\ell_{1}$-distance555Generally ergodicity is defined using the total variation distance between probability distributions. However for a discrete state-space $\mathcal{E}$ the total variation distance among probability distributions is exactly half of the distance computed using the $\ell_{1}$ norm. So we work with the $\ell_{1}$ norm in this paper. between probability measures $p(t)$ and $\pi$. Furthermore the CTMC is called exponentially ergodic if the convergence in (2.10) is exponentially fast, i.e. there exist positive constants $C$ and $\rho$ such that for any $t>0$

[TABLE]

Here the constant $C$ may depend on the initial distribution $p(0)$ but the constant $\rho$ does not (see [19] for example).

2.2 The Finite State Projection Algorithm

In the FSP method, approximate solutions of the CME (2.9) are obtained by restricting it to a truncated state-space. Suppose this truncated subset is given by a finite set $\mathcal{E}_{n}\subset\mathcal{E}$ of size $n=|\mathcal{E}_{n}|$. Using the same enumeration $\phi$ as in Section 2.1, we can express the set $\mathcal{E}_{n}$ as $\mathcal{E}_{n}=\{x_{j_{1}},x_{j_{2}},\dots,x_{j_{n}}\}$. Letting $Q_{n}$ to be the matrix formed by the rows and columns of matrix $Q$ in the set $J_{n}:=\{j_{1},\dots,j_{n}\}$, we approximate (2.9) by the $n$-dimensional linear system

[TABLE]

The solution of this system is simply $p_{n}(t)=\exp(Q_{n}^{T}t)p_{n}(0)$, where $p_{n}(0)$ is the $n\times 1$ containing the components of vector $p(0)$ in the set $J_{n}$.

Let ${\bf 1}$ be the $n$-dimensional vector of all ones. We assume that the initial state $X(0)$ can only take values in $\mathcal{E}_{n}$ and so ${\bf 1}^{T}p_{n}(0)=1$. It is easy to check that all the rows of matrix $Q_{n}$ have a non-positive sum, which implies that for any $t\geq 0$

[TABLE]

Results in [5] show that $\epsilon_{n}(t)$ quantifies the “error” between the actual solution of CME $p(t)$ and its approximation $p_{n}(t)$. For any fixed $t$, this error $\epsilon_{n}(t)$ decreases monotonically with increasing values of $n$. Moreover as $n\to\infty$ and the truncated state-space $\mathcal{E}_{n}$ approaches the full state-space $\mathcal{E}$, we have $\epsilon_{n}(t)\to 0$. In the FSP algorithm of [5], the final time $t_{f}$ is fixed and the system (2.11) is solved in the time-interval $[0,t_{f}]$ with some truncated state-space $\mathcal{E}_{n}$. Thereafter the error $\epsilon_{n}(t_{f})$ is evaluated and if this value is below some pre-specified tolerance level $\epsilon$, then the algorithm is terminated. Otherwise the truncated state-space is expanded to include more states and the same is process is repeated. After finitely many such iterations, the truncated state-space becomes large enough to ensure that the tolerance criterion is met.

Another way to formulate the FSP method is to consider the projected CTMC over the state-space $\widetilde{\mathcal{E}}_{n}=\mathcal{E}_{n}\cup\{x_{A}\}$, whose transitions among the states in $\mathcal{E}_{n}$ are same as the original CTMC, but any outgoing transitions from the set $\mathcal{E}_{n}$ are absorbed in the state $x_{A}$ (see Figure 1B), which serves as a proxy for all states in the set $\mathcal{E}^{c}_{n}=\{x\in\mathcal{E}:x\notin\mathcal{E}_{n}\}$. Enumerating the elements of $\widetilde{\mathcal{E}}_{n}$ as $\widetilde{\mathcal{E}}_{n}=\{x_{j_{1}},\dots,x_{j_{n}},x_{A}\}$, the $(n+1)\times(n+1)$ transition rate matrix $\widetilde{Q}_{n}$ for this projected CTMC is given by

[TABLE]

where $c_{n}$ is the $n$-dimensional column vector whose $i$-th component is

[TABLE]

This choice of $c_{n}$ ensures that all the rows of matrix $\widetilde{Q}_{n}$ sum to [math] and hence $\widetilde{Q}_{n}$ is a valid transition rate matrix. Let $\widetilde{p}_{n}(t)$ be the solution of the CME (2.9) corresponding to rate matrix $\widetilde{Q}_{n}$ and with initial value $\widetilde{p}_{n}(0)=(p_{n}(0),0)$. Then it can be shown that for any $t\geq 0$ we can express $\widetilde{p}_{n}(t)$ as

[TABLE]

which proves that the FSP approximation error $\epsilon_{n}(t)$ at time $t$ is exactly the amount of probability-mass that has been absorbed by the state $x_{A}$ in the time-interval $[0,t]$.

One can show that typically for any fixed truncated state-space $\mathcal{E}_{n}$, we would have $\epsilon_{n}(t)\to 1$ as $t\to\infty$, which says that eventually all the probability mass gets absorbed by the state $x_{A}$. Therefore even if the original CTMC is ergodic and the solution $p(t)$ of the CME (2.7) converges to $\pi$ as $t\to\infty$, the approximate solution $p_{n}(t)$ obtained by solving the FSP system (2.11), will not be close to $\pi$ for large times and in fact $p_{n}(t)$ converges to a vector of all zeros at $t\to\infty$. This is also evident from the stationary distribution $\widetilde{\pi}_{n}$ that can be computed by finding a non-zero solution to the linear-algebraic system (1.1) corresponding to the matrix $\widetilde{Q}_{n}$ (see (2.14)). This would yield a stationary distribution of the form $\widetilde{\pi}_{n}=({\bf 0},1)$, which assigns all the mass to the absorbing state $x_{A}$ and hence $\widetilde{\pi}_{n}$ cannot be close to $\pi$. This shows that the FSP approach is not conducive for the estimation of stationary distribution $\pi$.

3 The stationary FSP method

In this section we present our method sFSP for estimating the stationary distribution $\pi$ for the CTMC model of a reaction network. This is accomplished by constructing a projected CTMC over the truncated state-space and computing the stationary distribution of this new CTMC. Keeping the same notation as in Section 2.2, this projected CTMC over the truncated state-space $\mathcal{E}_{n}=\{x_{j_{1}},\dots,x_{j_{n}}\}\subset\mathcal{E}$ is constructed by redirecting the transitions that leave $\mathcal{E}_{n}$ to some designated state $x\in\mathcal{E}_{n}$ (see Figure 1C). Let the $n\times n$ matrix $Q_{n}$ and the $n\times 1$ vector $c_{n}$ be as in Section 2.2. Then the $n\times n$ transition rate matrix $\overline{Q}_{n}$ for this CTMC is simply given by

[TABLE]

where $l$ corresponds to the address of the designated state (i.e. $x_{j_{l}}=x$) and $b_{l}$ is the $1\times n$ vector whose $l$-th component is $1$ and the rest are all zeros. Essentially, $\overline{Q}_{n}$ is formed by adding the non-negative vector $c_{n}$ to the $l$-th column of matrix $Q_{n}$. All the rows of matrix $\overline{Q}_{n}$ sum to [math] and hence $\overline{Q}_{n}$ is a valid transition rate matrix and so our projected CTMC is well-defined. Our method sFSP estimates the stationary distribution $\pi$ by computing the finite-dimensional stationary distribution $\overline{\pi}_{n}$ for the projected CTMC with transition rate matrix $\overline{Q}_{n}$. Using $\overline{\pi}_{n}$, we can also compute the overall outflow rate at the estimated stationary distribution by

[TABLE]

This quantity will play a key role in bounding the sFSP approximation error whose direct computation is impossible.

3.1 Analysis of sFSP

The aim of this section is to demonstrate that under certain conditions, that are commonly satisfied by biological reaction networks, the sFSP approximation error can be made arbitrarily small by picking a truncated state-space $\mathcal{E}_{n}$, that is large enough. Moreover it is possible to check if $\mathcal{E}_{n}$ is large enough by computing a convergence factor which is defined by suitably scaling the outflow rate $r^{(n)}_{\textnormal{out}}$. The main results of this section are collected in Theorem 3.1 and they provide the theoretical basis for our sFSP method.

Before we present our result we need to discuss some preliminary concepts. The state-space $\mathcal{E}$ of the original CTMC $(X(t))_{t\geq 0}$ is called irreducible if this CTMC has a positive probability of reaching any state in $\mathcal{E}$ from any other state in $\mathcal{E}$, in a finite time. More formally, the state-space $\mathcal{E}$ is irreducible, if for any $x,y\in\mathcal{E}$ we have $\mathbb{P}(X(t)=y|X(0)=x)>0$ for some $t>0$. In our setting of reaction networks, this is equivalent to saying that between any two states $x,y\in\mathcal{E}$ there exists a sequence of positive-propensity reactions $k_{1},\dots,k_{n}$ that takes the dynamics from $x$ to $y$. For this to hold we must have $y=x+\sum_{i=1}^{n}\zeta_{k_{i}}$ and at each intermediate state $z_{j}=(x+\sum_{i=1}^{j-1}\zeta_{k_{i}})$ the next reaction in the sequence ($k_{j}$) has a positive propensity of firing ($\lambda_{k_{j}}(z_{j})>0$). When only finitely many states are accessible by the reaction dynamics, irreducible state-spaces can be easily found by manipulating the transition rate matrix $Q$ (see [21]). However when infinitely many states are accessible, finding irreducible state-spaces within the infinite lattice becomes a complicated task. In a recent work [15] we address this challenge and develop a computational procedure that can find all the irreducible state-spaces for a large class of biological reaction networks. In particular, for most networks of interest each irreducible state-space has the form666To obtain this form relabeling of species may be required.

[TABLE]

where $\mathcal{E}_{b}$ is a finite set in $\mathbb{N}^{d_{b}}$, and $d_{b},d_{f}$ are non-negative integers summing up to the total number of species $d$. Here $\mathcal{E}_{b}$ contains the dynamics of $d_{b}$ bounded species whose copy-numbers are required to satisfy a positive mass-conversation relation. A typical example is a gene-expression network where the gene of interest has many (say $d_{b}$) activity modes. To represent the dynamics we need to represent each such mode by a different network species, but all these species will be bounded and their copy-numbers will evolve in a finite set $\mathcal{E}_{b}$, because the gene of interest has a fixed copy-number (see the Pap-Switch example in Section 4.3 for instance). The species that are not bounded are free777Apart from free and bounded species, there may also exist another type of species, called restricted species, whose dynamics essentially mimics the dynamics of free species according to some affine map. However these restricted species can be easily eliminated to obtain a dynamically equivalent network and hence we ignore such species here (see [15] for more details). to have any copy-number and hence the state-space for their dynamics is taken to be the full non-negative integer orthant $\mathbb{N}^{d_{f}}_{0}$.

Note that the property of ergodicity (see Section 2.1) will obviously fail if there do not exist any stationary distributions or there exist more than one stationary distributions for the CTMC $(X(t))_{t\geq 0}$. If the state-space $\mathcal{E}$ is finite, then its irreducibility is sufficient to guarantee that the stationary distribution exists uniquely and the CTMC is exponentially ergodic (see [21]). However when $\mathcal{E}$ is infinite, its irreducibility can only guarantee the uniqueness of a stationary distribution but the existence of this distribution must be checked by other means, for example, using the results in [22] and [19]. In particular Theorem 7.1 in [19] guarantees the existence of a stationary distribution along with exponential ergodicity, if we can construct a function $V:\mathcal{E}\to[1,\infty)$ which is norm-like (i.e. $V(x)\to\infty$ as $\|x\|\to\infty$) and for some $C_{1},C_{2}>0$, the following holds for all $x\in\mathcal{E}$:

[TABLE]

where $\mathbb{Q}$ is the CTMC generator given by (2.4). This condition is called the Foster-Lyapunov criterion in the literature and it describes the tendency of the CTMC to experience a drift towards some finite set in the state-space with a force that is proportional to the distance from this finite set, measured according to $V$. In [16] it is shown that for many biomolecular reaction networks, a linear Foster-Lyapunov function

[TABLE]

satisfying (3.19) can be constructed. Here $v\in\mathbb{R}^{d}$ is a positive vector which is chosen using simple Linear Programming and $\langle\cdot,\cdot\rangle$ denotes the standard inner product in $\mathbb{R}^{d}$. Observe that for the linear function $V(x)$ (3.20), the drift condition (3.19) is simply

[TABLE]

As demonstrated in [16], often for biological reaction networks the vector $v$ can be chosen in such a way that along with this drift condition, the following diffusivity condition is also satisfied - for some $C_{3},C_{4}>0$

[TABLE]

When (3.21) and (3.22) hold simultaneously, then in addition to exponential ergodicity, one can also guarantee other desirable properties like finiteness of all statistical moments of the stationary distribution $\pi$ and convergence of all the moments of the CTMC to their steady-state values as time approaches infinity (see Theorem 5 in [16]).

To study the sFSP approximation error we need to work with the norm prescribed by the Foster-Lyapunov function $V$. For any signed measure $\mu$ on $\mathcal{E}$, this norm is given by

[TABLE]

Note that this norm is tighter than the $\ell_{1}$ norm because $\|\mu\|_{V}\geq\|\mu\|_{\ell_{1}}$ as $V\geq 1$. Let $\mathcal{B}(\mathcal{E}_{n})$ denote the boundary of the truncated state-space $\mathcal{E}_{n}$, which includes all those states in $\mathcal{E}_{n}$ for which there exists a positive-propensity reaction that takes the dynamics outside $\mathcal{E}_{n}$, i.e.

[TABLE]

Based on the outflow rate $r^{(n)}_{\textnormal{out}}$ given by (3.17), we define the convergence factor as

[TABLE]

where

[TABLE]

and $x_{\ell}$ is the designated state. Our next result will show that the convergence factor $\gamma^{(n)}_{V}$ is a useful diagnostic tool to assess the approximation error $\|\pi-\overline{\pi}_{n}\|_{V}$ of sFSP. Note that unlike the approximation error, $\gamma^{(n)}_{V}$ can be explicitly computed from the sFSP output $\overline{\pi}_{n}$ if the Foster-Lyapunov function $V$ is known. In situations where $V$ is unknown, the definition of $\gamma^{(n)}_{V}$ can often be suitably modified to preserve its diagnostic purpose (see Remark 3.3).

We now come to the main result of our paper.

Theorem 3.1

Suppose that state-space $\mathcal{E}$ is irreducible for the original CTMC with transition rate matrix $Q$, and there exists a Foster-Lyapunov function $V:\mathcal{E}\to[1,\infty)$ satisfying (3.19). Also assume that $\{\mathcal{E}_{n}:n=1,2,\dots\}$ is a sequence of finite sets that is increasing (i.e. $\mathcal{E}_{n_{1}}\subset\mathcal{E}_{n_{2}}$ if $n_{1}<n_{2}$) and that covers the full state-space $\mathcal{E}$ in the limit $n\to\infty$. Fix a designated state $x_{\ell}\in\mathcal{E}_{1}$ and let $\overline{Q}_{n}$ be the transition rate matrix of our projected CTMC with state-space $\mathcal{E}_{n}$, defined according to (3.16). Then we have the following:

• (A)

The stationary distribution $\overline{\pi}_{n}$ for the projected CTMC exists uniquely.

• (B)

As $n\to\infty$, $\overline{\pi}_{n}$ converges to the stationary distribution $\pi$ for the original CTMC, in the $\ell_{1}$ metric, i.e.

[TABLE]

• (C)

There exists a positive constant $M$ such that for any $n$

[TABLE]

where $\gamma^{(n)}_{V}$ is the convergence factor defined by (3.24).

• (D)

Suppose that the Foster-Lyapunov function $V$ has the linear form (3.20) and the positive vector $v$ is such that both (3.21) and (3.22) are satisfied. Furthermore assume that the sequence of sets $\{\mathcal{E}_{n}\}$ grows uniformly w.r.t. function $V$ which means that for some constant $\theta\in(0,1)$ we have

[TABLE]

where $\mathcal{B}(\mathcal{E}_{n})$ is the boundary of $\mathcal{E}_{n}$ defined by (3.23). Then there exists a constant $M^{\prime}>0$ for which the converse of (3.27) also holds, i.e. for each $n$

[TABLE]

Furthermore, the convergence factor $\gamma^{(n)}_{V}$ converges to [math] as $n\to\infty$.

Remark 3.2

It will become evident from the proof that if the Foster-Lyapunov function $V$ and constants $C_{1},C_{2}$ in (3.19) are known, then a constant $M$ satisfying part (C) can be explicitly computed using the results in Meyn and Tweedie [23]. Hence part (C) provides a computable upper-bound for the approximation error $\|\pi-\overline{\pi}_{n}\|_{V}$. Similarly the constant $M^{\prime}$ satisfying part (D) may be explicitly computed from constants $C_{1},\dots,C_{4}$ in (3.21) and (3.22), and the constant $\theta$ that appears in (3.28). The tightness of the error bounds obtained from these explicitly computable constants remains to be investigated. Nevertheless parts (C) and (D) are useful in demonstrating that up to a constant, the magnitude of the uncomputable approximation error $\|\pi-\overline{\pi}_{n}\|_{V}$ can be assessed by computing the convergence factor $\gamma^{(n)}_{V}$. In other words, if $\gamma^{(n)}_{V}\leq\epsilon$ then $\|\pi-\overline{\pi}_{n}\|_{V}\leq M\epsilon$, and similarly if $\gamma^{(n)}_{V}\geq\epsilon$ then $\|\pi-\overline{\pi}_{n}\|_{V}\geq M^{\prime}\epsilon$, where $M$ and $M^{\prime}$ are the optimal constants for which parts (C) and (D) hold.

Remark 3.3

Note that computation of the convergence factor $\gamma^{(n)}_{V}$ (3.24) requires knowledge of the Foster-Lyapunov function $V$ which is undesirable from the point of view of applications. However it is possible to circumvent this problem, if one has information about the form of $V$ and the shape of finite sets $\{\mathcal{E}_{n}\}$. For this one needs to pick a sequence $\{\beta_{n}\}$ such that for some constants $\alpha,\alpha^{\prime}>0$

[TABLE]

holds for each $n$, with $\|\mathcal{E}_{n}\|_{V}$ defined by (3.25). Then one can define the convergence factor as

[TABLE]

with the outflow rate $r^{(n)}_{\textnormal{out}}$ given by (3.17), and parts (C) and (D) will hold with the substitutions, $\gamma^{(n)}_{V}\to\gamma_{n}$, $M\to M\alpha$ and $M^{\prime}\to M^{\prime}\alpha^{\prime}$. For example, if $V$ has the linear form (3.20), then one can define $\beta_{n}$ in the same way as $\|\mathcal{E}_{n}\|_{V}$ but with $V(x)$ replaced by any norm $\|x\|$ on $\mathbb{R}^{d}$.

Proof. We start by proving part (A). The stationary distribution $\overline{\pi}_{n}$ for the projected CTMC certainly exists because the transition rate matrix $\overline{Q}_{n}$ is finite (see [21]). This stationary distribution can be found by solving the linear-algebraic system (1.1) with transition-rate matrix $\overline{Q}_{n}$. We now prove by contradiction the uniqueness of this stationary distribution. Suppose that this uniqueness does not hold. Then there would exist at least two disjoint non-empty irreducible state-spaces (say $A$ and $B$) for the projected CTMC within the state-space $\mathcal{E}_{n}$. This implies that if the projected CTMC starts in set $A$ then it remains in this set for all times, and there is a positive probability for this CTMC to reach any state in $A$ from any other state in $A$ in a finite time. The same holds true for set $B$. Certainly one of these sets, say $A$, will not contain the designated state $x_{\ell}$ but this leads to a contradiction due to the following reasons. Since the state-space $\mathcal{E}$ is irreducible for the original CTMC, there exists a sequence of reactions $k_{1},\dots,k_{m}$ that takes the original CTMC from any state $x\in A$ to the designated state $x_{\ell}$ with a positive probability. If all the intermediate states that arise in this reaction path (recall $z_{j}$-s from above) lie within the set $\mathcal{E}_{n}$, then the same sequence of reactions will also take the projected CTMC from state $x\in A$ to state $x_{j_{l}}$, which is a contradiction because $A$ is an irreducible state-space not containing $x_{j_{l}}$. On the other hand if one of the intermediate states lies outside $\mathcal{E}_{n}$, then the last reaction, say $k_{q}$, in the sequence that takes the dynamics outside $\mathcal{E}_{n}$ will be redirected to the designated state $x_{\ell}$ in the projected CTMC and hence again we have a contradiction because $k_{1},\dots,k_{q}$ is a positive-probability sequence of reactions that takes the projected CTMC from state $x\in A$ to state $x_{j_{l}}\notin A$. Therefore the stationary distribution $\overline{\pi}_{n}$ for the projected CTMC is unique, and this completes the proof of part (A).

We now prove part (B). Clearly the assertion of part (B) is trivial when the full state-space $\mathcal{E}$ is finite and so we assume that $\mathcal{E}$ is infinite from now on. Let $\{\mathcal{E}_{n}\}$ be a sequence of sets as stated in the proposition and let $\phi:\mathcal{E}\to\mathbb{N}_{0}$ be an enumeration of $\mathcal{E}$ satisfying

[TABLE]

Such an enumeration exists because $\{\mathcal{E}_{n}\}$ is an increasing sequence of sets that cover the set $\mathcal{E}$ in the limit $n\to\infty$. Note that as each $\mathcal{E}_{n}$ is a finite set, condition (3.31) ensures that

[TABLE]

Now consider the $\mathbb{N}_{0}$-valued, one-dimensional process $(Y(t))_{t\geq 0}$ given by $Y(t)=\phi(X(t))$ for each $t\geq 0$, where $(X(t))_{t\geq 0}$ is the original CTMC with transition rate matrix $Q$ and generator $\mathbb{Q}$ (see (2.4)). As $\phi$ is a one-to-one and onto map, the process $(Y(t))_{t\geq 0}$ is also a CTMC and its generator is given by

[TABLE]

where $g$ is a bounded real-valued function on $\mathbb{N}_{0}$ and $f$ is the bounded real-valued function on $\mathbb{N}^{d}_{0}$ defined by $f(x)=g(\phi(x))$.

Irreducibility of state-space $\mathcal{E}$ for $(X(t))_{t\geq 0}$ implies the irreducibility of state-space $\mathbb{N}_{0}$ for $(Y(t))_{t\geq 0}$. Let $V:\mathbb{N}^{d}_{0}\to[0,\infty)$ be the norm-like Foster-Lyapunov function satisfying (3.19) and define the function $\widehat{V}:\mathbb{N}_{0}\to[0,\infty)$ by $\widehat{V}(i)=V(\phi^{-1}(i))$. Then using (3.32) and (3.19) we can deduce that $\widehat{V}$ is a norm-like function satisfying

[TABLE]

Therefore $\widehat{V}$ is a Foster-Lyapunov function for CTMC $(Y(t))_{t\geq 0}$ with generator $\widehat{\mathbb{Q}}$ and hence this CTMC is exponentially ergodic due to Theorem 7.1 in [19]. Let $\widehat{\pi}$ and $\widehat{\pi}_{n}$ be the probability distributions on $\mathbb{N}_{0}$ and $\{0,1,\dots,|\mathcal{E}_{n}|-1\}$ defined by

[TABLE]

Then $\widehat{\pi}$ is the stationary distribution for the CTMC $(Y(t))_{t\geq 0}$ and $\widehat{\pi}_{n}$ is the stationary distribution of this CTMC projected onto the finite state-space $\{0,1,\dots,|\mathcal{E}_{n}|-1\}$ by redirecting all the outgoing transitions to the designated state $\phi(x_{\ell})$. Theorem 3.3 in [24] proves

[TABLE]

using resolvent forms (see (3.40)). This limit is equivalent to (3.26) and this proves part (B).

We will now prove part (C). Without loss of generality we can assume that $\mathcal{E}_{n}=\{0,1,\dots,n-1\}$. Define an infinite vector

[TABLE]

whose first $n$ elements are $\vartheta_{1}=Q_{n}^{T}\overline{\pi}_{n}$, where $Q_{n}$ denotes the $n\times n$ northwest sub-matrix of $Q$. Recall that matrix $\overline{Q}_{n}$ is given by (3.16) and the outflow rate $r^{(n)}_{\textnormal{out}}$ is defined by (3.17). As $\overline{Q}^{T}_{n}\overline{\pi}_{n}={\bf 0}$ we can write $\vartheta_{1}$ as

[TABLE]

which shows that the $n\times 1$ vector $\vartheta_{1}$ has only one non-zero entry which is equal to $-r^{(n)}_{\textnormal{out}}$ and it is at the position corresponding to the designated state $x_{\ell}$. Since $Q{\bf 1}={\bf 0}$ we have ${\bf 1}^{T}\vartheta_{n}={\bf 0}$ which implies that

[TABLE]

One can check that all entries of the infinite vector $\vartheta_{2}$ are non-negative and only those entries are non-zero that correspond to the states in the boundary set $\mathcal{B}(\mathcal{E}_{n})$ (see (3.23)) of $\mathcal{E}_{n}$. Therefore, viewing $\vartheta_{n}$ as a signed measure over $\mathcal{E}$, we can express it as

[TABLE]

where $\mu_{1}$ and $\mu_{2}$ are probability measures on $\mathcal{E}$, supported on $\{x_{\ell}\}$ and $\mathcal{B}(\mathcal{E}_{n})$ respectively. With a slight abuse of notation, we will denote the vector-version of $\mu_{i}$ also as $\mu_{i}$.

Define the sFSP approximation error in vector form as

[TABLE]

and since $Q^{T}\pi={\bf 0}$ we get the following equation from (3.37)

[TABLE]

One can verify that $\epsilon_{n}$ is the unique solution of this linear system with the constraint $\langle{\bf 1},\epsilon_{n}\rangle=0$. For any $\beta>0$, let $R_{\beta}$ denote the $\beta$-resolvent matrix corresponding to the transition rate matrix $Q$. It is defined by

[TABLE]

where ${\bf I}$ denotes the identity matrix. It is known (see [24]) that $R_{\beta}$ is a positive matrix satisfying $R_{\beta}{\bf 1}={\bf 1}$, $\pi^{T}R_{\beta}=\pi^{T}$ and

[TABLE]

One can regard $R_{\beta}$ as the transition matrix of a discrete-time Markov chain over $\mathcal{E}=\{x_{0},x_{1},\dots\}$ whose unique stationary distribution is $\pi$.

Expressing the Foster-Lyapunov function $V$ as the vector $V=(V(x_{0}),V(x_{1}),\dots)$ we can write the drift condition (3.19) as

[TABLE]

This relation along with (3.41) and the positivity of $R_{\beta}$ implies

[TABLE]

Letting $\lambda=(1+C_{2}/\beta)^{-1}$ and $C=C_{1}/(\lambda\beta)$ we obtain

[TABLE]

Note that $\lambda\in(0,1)$. Theorem 6.1 in [23] shows that we can explicitly compute constants $C^{\prime}>0$ and $\rho\in(0,1)$, such that for any probability distribution $\mu$ over $\mathcal{E}$ we have

[TABLE]

where $R^{m}_{\beta}$ denotes the $m$-th power of the matrix $R_{\beta}$. Transposing (3.39), multiplying both sides by $R_{\beta}$ and using (3.41) and (3.38) we get

[TABLE]

One can write $\epsilon_{n}$ as

[TABLE]

where $\epsilon_{j}$ is the solution to

[TABLE]

for $j=1,2$. This solution can be expressed as

[TABLE]

and using (3.42) we get

[TABLE]

Therefore

[TABLE]

where $M=C^{\prime}\beta^{-1}\rho(1-\rho)^{-1}$. As $\mu_{1}$ and $\mu_{2}$ are probability distributions supported on $\{x_{\ell}\}$ and $\mathcal{B}(\mathcal{E}_{n})$, we have $(\|\mu_{1}\|_{V}+\|\mu_{2}\|_{V})\leq\|\mathcal{E}_{n}\|_{V}$ (see (3.25)). This proves part (C) of the theorem.

We now prove part (D). Here we assume that the Foster-Lyapunov function $V$ has the linear form (3.20) and both (3.21) and (3.22) are satisfied. Note that by rescaling the positive vector $v$ in (3.20) if necessary, we can assume that

[TABLE]

As $\mathbb{Q}V(x)=\sum_{k=1}^{K}\lambda_{k}(x)\langle v,\zeta_{k}\rangle$ from condition (3.21) we obtain

[TABLE]

for each $x\in\mathcal{E}$. Transposing (3.39), multiplying both sides by vector $V$ and taking absolute values we get

[TABLE]

Using (3.44) we can upper-bound the l.h.s. as

[TABLE]

Since $\vartheta_{n}$ is given by (3.38), with $\mu_{1}$ and $\mu_{2}$ being probability distributions supported on $\{x_{\ell}\}$ and $\mathcal{B}(\mathcal{E}_{n})$ respectively, we can lower-bound the r.h.s. of (3.45) as

[TABLE]

The uniform growth condition (3.28), along with the fact that $V(x_{\ell})$ does not depend on $n$, ensures that there exists a positive constant $\theta^{\prime}$ such that

[TABLE]

for each $n$, and hence obtain the lower-bound

[TABLE]

This relation along with (3.45) and (3.46) yield

[TABLE]

which is sufficient to prove (3.29) as $\|\epsilon_{n}\|_{\ell_{1}}\leq\|\epsilon_{n}\|_{V}$.

We now prove the second assertion of part (D), i.e. $\gamma^{(n)}_{V}\to 0$ as $n\to\infty$. For this we first demonstrate that the square of the linear Foster-Lyapunov $V$ will also satisfy the drift condition (3.19). To see this note that for any $x\in\mathcal{E}$

[TABLE]

Using (3.21) and (3.22) we obtain

[TABLE]

As $V(x)$ is a semi-norm, the quadratic term will dominate the linear term for all $x$ outside some compact set and hence the drift condition (3.19) will be satisfied by function $V^{2}(x)$ for some constants $\widehat{C}_{1},\widehat{C}_{2}>0$. This drift condition also ensures that (see [19]) there exists a constant $L$ such that

[TABLE]

for each $n$. Now using Cauchy-Schwarz inequality we get

[TABLE]

As $n\to\infty$, part (B) shows that $\|\epsilon_{n}\|_{\ell_{1}}\to 0$ and hence $\|\epsilon_{n}\|_{V}\to 0$ as well. Now (3.29) proves that $\gamma^{(n)}_{V}\to 0$ and this concludes the proof of the theorem.

$\Box$

3.2 The sFSP Algorithm

Theorem 3.1 proves that under certain conditions, the sFSP approximation error, measured in a certain norm, converges to [math] as $n\to\infty$ and the truncated state-space $\mathcal{E}_{n}$ expands to the fully state-space $\mathcal{E}$. Moreover for any $\mathcal{E}_{n}$ the magnitude of the approximation error can be judged by computing the convergence factor $\gamma_{n}$ defined according to (3.30) with the sequence $\{\beta_{n}\}$ chosen as in Remark 3.3. These results form the basis of our stationary Finite State Projection (sFSP) algorithm, that is presented as Algorithm 1. This algorithm takes as input a $d$-species reaction network $\mathcal{R}$, specified as a set of $K$ reactions with propensity functions $\lambda_{1},\dots,\lambda_{K}$ and stoichiometric vectors $\zeta_{1},\dots,\zeta_{K}$. It is required that the CTMC describing the reaction kinetics admits a Foster-Lyapunov function satisfying (3.19) and its state-space $\mathcal{E}$ is irreducible. These conditions can be checked using the results in [15] and [16] as discussed before.

Algorithm 1 starts by picking an increasing sequence of finite state-space truncations $\{\mathcal{E}_{i}:i=1,2,\dots\}$ as in Theorem 3.1, a sequence $\{\beta_{i}:i=1,2,\dots\}$ as in Remark 3.3, and a designated state $x_{\ell}\in\mathcal{E}_{1}$. Thereafter for each iteration cycle $i$, the transition rate matrix $\overline{Q}_{i}$ for the projected CTMC over the truncated state-space $\mathcal{E}_{i}$ is constructed and its stationary distribution $\overline{\pi}_{i}$ is found by solving the linear-algebraic system (1.1) for matrix $\overline{Q}_{i}$. Next the outflow rate $r^{(i)}_{\textnormal{out}}$ and the convergence factor $\gamma_{i}=r^{(i)}_{\textnormal{out}}\beta_{i}$ are computed. If this convergence factor is below an acceptable threshold level $\epsilon$ (chosen in step 4 of Algorithm 1), then sFSP terminates after returning $\overline{\pi}_{i}$ as the estimate of the true stationary distribution $\pi$. Otherwise if $\gamma_{i}\geq\epsilon$, then the algorithm goes into the new iteration cycle with the expanded truncated state-space $\mathcal{E}_{i+1}$.

4 sFSP Algorithm: Simple Implementation

In this section we present the simple implementation of sFSP akin to to the classical FSP [5], where the multi-dimensional state-space is explicitly enumerated, and accordingly the transition rate matrix for the projected CTMC is constructed and its stationary distribution vector is computed. The performance of sFSP depends crucially on the choice of finite state-space truncations $\mathcal{E}_{i}$-s and their enumerating functions $\phi_{i}$-s. We now discuss these choices for our implementation of sFSP.

4.1 State-space enumeration and truncation

The basic ingredient of our state-space enumeration strategy is the Cantor Pairing function (see [25]) which is the bijective map between $\mathbb{N}^{2}_{0}$ to $\mathbb{N}_{0}$ defined by

[TABLE]

Under this bijection, the elements in $\mathbb{N}^{2}_{0}$ are mapped to $\mathbb{N}_{0}$ by moving along the anti-diagonals, which are the straight lines given by $x_{1}+x_{2}=k$ (see Figure 2A). This map is easy to invert and for any $z\in\mathbb{N}_{0}$, $(x_{1},x_{2})=\Phi^{-1}_{2}(z)$ can be computed as $x_{1}=v-x_{2}$ and $x_{2}=z-v(v+1)/2$, where

[TABLE]

Henceforth we define $\Phi_{1}$ as the identity map on $\mathbb{N}_{0}$. By composition, one can extend the Cantor function to obtain a bijection from $\mathbb{N}^{n}_{0}$ to $\mathbb{N}_{0}$ for any positive integer $n$. Such a bijective map $\Phi_{n}$ can be defined recursively as

[TABLE]

Similarly the inverse $\Phi^{-1}_{n}:\mathbb{N}_{0}\to\mathbb{N}^{n}_{0}$ of this map can also be defined recursively as

[TABLE]

where $(z_{1},z_{2})=\Phi^{-1}_{2}(z)$.

Consider the situation where the irreducible state-space $\mathcal{E}$ has the form (3.18) with $d_{b}=0$ and $d=d_{f}$. In this case, $\mathcal{E}$ is just the $d$-dimensional non-negative integer orthant $\mathbb{N}^{d}_{0}$ and we enumerate it using the Cantor function $\Phi_{d}$. An explicit formula for $\Phi_{d}$ can be obtained (see [25]) as

[TABLE]

where ${n\choose k}=\frac{n!}{k!(n-k)!}$ denotes the binomial coefficient. This formula shows that for any $C_{l},C_{r}\in\mathbb{N}_{0}$ with $C_{l}\leq C_{r}$, the following set

[TABLE]

is non-empty, and we call this set a trapezoidal truncation of $\mathbb{N}^{d}_{0}$ with left cut-off point $C_{l}$ and right cut-off point $C_{r}$. For $d=2$, we plot such a set in Figure 2B and it simply consists of all the states $(x_{1},x_{2})$ whose component-sum $x_{1}+x_{2}$ is between $C_{l}$ and $C_{r}$. This may not be exactly true is higher-dimensions ($d>2$) but still one can think of a trapezoidal truncation as the set of states whose component-sum is within certain bounds. Note that in our setting of reaction networks, the component-sum of a state represents the total molecular count of all the species. In many biomolecular reaction networks this total molecular count is within certain tight bounds even though each species can individually have high copy-number variation. This is mainly because the species are often in competition with each other, through mechanisms such as mutual repression or interconversion, which ensures that the total molecular count is tightly regulated. This property makes trapezoidal truncations very appealing for our purpose of estimating stationary distributions. This point is nicely illustrated by the Toggle-Switch example considered in Section 4.3.2.

We now consider the situation where the irreducible state-space $\mathcal{E}$ has the form (3.18) for some finite non-empty set $\mathcal{E}_{b}\subset\mathbb{N}^{d_{b}}_{0}$. Let $N_{b}=|\mathcal{E}_{b}|$ and we fix an enumeration of this set as $\mathcal{E}_{b}=\{e_{0},\dots,e_{N_{b}-1}\}$. This enables us to define an enumeration over the full state-space $\mathcal{E}=\mathcal{E}_{b}\times\mathbb{N}^{d_{f}}_{0}$ by

[TABLE]

where $e=e_{j}\in\mathcal{E}_{b}$ and $x\in\mathbb{N}^{d_{f}}_{0}$. One can easily see that this map is a bijection between $\mathcal{E}_{b}\times\mathbb{N}^{d_{f}}_{0}$ and $\mathbb{N}_{0}$, and its inverse is given by

[TABLE]

where $j$ is the remainder in the division of $z$ by $N_{b}$ and $q$ is the corresponding quotient. For the state-space $\mathcal{E}_{b}\times\mathbb{N}^{d_{f}}_{0}$ we define the trapezoidal truncation as

[TABLE]

where $C_{l}$ and $C_{r}$ are non-negative integers satisfying $C_{l}\leq C_{r}$ as before.

We now come to the definitions of finite state-space truncations $\mathcal{E}_{i}$-s and their enumerating functions $\phi_{i}$-s. Let $\{C_{l,i}:i=1,2,\dots\}$ and $\{C_{r,i}:i=1,2,\dots\}$ be monotonic sequences of non-negative integers that satisfy $C_{l,i}\leq C_{r,i}$ for each $i$ along with the limits

[TABLE]

For each $i=1,2,\dots$ we define the finite state-space truncation $\mathcal{E}_{i}$ as $\mathcal{T}(C_{l,i},C_{r,i})$. Note that monotonicity of the left and right cut-off sequences $\{C_{l,i}\}$ and $\{C_{r,i}\}$ along with (4.52) ensures that $\{\mathcal{E}_{i}:i=1,2,\dots\}$ is an increasing sequence of finite sets that covers the full state-space $\mathcal{E}$ in the limit $i\to\infty$, as demanded by the sFSP Algorithm 1. Assuming that the Foster-Lyapunov function $V$ has the linear form (3.20), we can choose the sequence $\{\beta_{i}:i=1,2,\dots\}$ (see Remark 3.3) in step 2 of Algorithm 1 as $\beta_{i}=C_{r,i}$.

In the case where the full state-space $\mathcal{E}$ is the non-negative integer orthant $\mathbb{N}^{d}_{0}$, the size of the truncated state-space $\mathcal{E}_{i}$ is

[TABLE]

and we enumerate the set $\mathcal{E}_{i}$ using the map $\phi_{i}:\mathcal{E}_{i}\to\{0,1,\dots,n_{i}-1\}$ given by

[TABLE]

Based on this enumeration the transition rate matrix $\overline{Q}_{i}$ for the projected CTMC over the truncated state-space $\mathcal{E}_{i}$ (see step 5 in Algorithm 1) can be constructed with Algorithm 2. In the other situation where the irreducible state-space $\mathcal{E}$ has the form (3.18) for some finite non-empty set $\mathcal{E}_{b}=\{e_{0},\dots,e_{N_{b}-1}\}$ with $N_{b}=|\mathcal{E}_{b}|$ elements, the size of the truncated state-space $\mathcal{E}_{i}$ is

[TABLE]

and we enumerate the set $\mathcal{E}_{i}$ using the map $\phi_{i}:\mathcal{E}_{i}\to\{0,1,\dots,n_{i}-1\}$ given by

[TABLE]

where $\Psi$ is the map defined by (4.49). The transition rate matrix $\overline{Q}_{i}$ for the projected CTMC over the truncated state-space $\mathcal{E}_{i}$ can be constructed using Algorithm 2 with some minor changes.

4.2 Implementation Details

We now provide some details on our computer implementation of sFSP Algorithms 1 and 2, and discuss the related issues. Note that the size $n_{i}$ of the truncated state-space $\mathcal{E}_{i}$ can be very large, causing problems in storing the $n_{i}\times n_{i}$ transition rate matrix $\overline{Q}_{i}$, and also in solving the linear-algebraic system (1.1) to obtain $\overline{\pi}_{i}$. Note however that out of $n^{2}_{i}$ entries in matrix $\overline{Q}_{i}$, at most $n_{i}(K+1)$ entries can be non-zero, where $K$ is the number of reactions which is typically much smaller than $n_{i}$. Hence $\overline{Q}_{i}$ is an extremely sparse matrix and this sparsity can be exploited for storing matrix $\overline{Q}_{i}$ and for finding the vector $\overline{\pi}_{i}$.

Another issue that commonly arises is that for states with large components, the propensity functions take very high values which causes the matrix $\overline{Q}_{i}$ to have very large entries. This creates numerical issues while solving the linear-algebraic system (1.1) for computing $\overline{\pi}_{i}$. A simple way to circumvent this problem is to scale the matrix $\overline{Q}_{i}$ by its diagonal entries and apply the same scaling to the solution of the linear-algebraic system to recover $\overline{\pi}_{i}$. In other words, matrix $\overline{Q}_{i}$ is constructed by modifying Algorithm 2 by setting $Q_{mm}$ to $-1$ in step 7 and by replacing $\lambda_{k}(y_{m})$ with $\lambda_{k}(y_{m})/\lambda_{0}(y_{m})$ in steps 12 and 14. Such a scaling is allowed because the state-space $\mathcal{E}$ is irreducible for the original CTMC and hence any $y_{m}\in\mathcal{E}$ cannot be an absorbing state and so $\lambda_{0}(y_{m})=\sum_{k=1}^{K}\lambda_{k}(y_{m})$ is nonzero. While constructing matrix $\overline{Q}_{i}$ we must also store the values $\lambda_{0}(y_{m})$ for $m=0,1,\dots,n_{i}$. These values help in recovering $\overline{\pi}_{i}$ from the solution $\widehat{\pi}_{i}$ of the linear-algebraic system solved in step 6 of Algorithm 1

[TABLE]

Of course $\overline{\pi}_{i}$ is then normalized (step 7 of Algorithm 1) to ensure that its component-sum is $1$ and it represents a valid stationary distribution.

In our setup we implement the main sFSP method (Algorithm 1) in Matlab but we delegate the construction of the transition rate matrix $\overline{Q}_{i}$ to a C++ program that implements Algorithm 2. Once constructed, this matrix is imported into the sFSP Matlab program as a sparse matrix. The linear-algebraic system (1.1) for this matrix is solved by computing the eigenvector corresponding to the smallest-magnitude eigenvalue (i.e. [math]) using the eigs function in Matlab. This function performs an Arnoldi iterative procedure [26] to efficiently compute a subset of eigenvalues and eigenvectors for large sparse matrices. It also allows us to pass a starting vector for the Arnoldi procedure. In our implementation we use the stationary distribution vector $\overline{\pi}_{i-1}$ obtained in iteration $(i-1)$ as the starting vector in iteration $i$888In the first iteration $i=1$, the starting vector is chosen to correspond to the uniform stationary distribution over the first state-space truncation $\overline{\mathcal{E}}_{1}$. For the sFSP implementation considered in this section, we use the scaled version of matrix $\overline{Q}_{i}$ as described above.

4.3 Computational Examples

In this section we illustrate our simple implementation of sFSP using examples from systems biology. In all the considered examples, sFSP is applicable because with results in [15] and [16] we can verify that the theoretical conditions required by sFSP (see Theorem 3.1) are satisfied. Moreover for all the examples, we fix the acceptable threshold level $\epsilon$ (see step 4 of Algorithm 1) to be $10^{-10}$, and we specify the increasing family of trapezoidal state-space truncations $\{\mathcal{E}_{i}=\mathcal{T}(C_{l,i},C_{r,i})\}$ via a pair of monotonic cut-off sequences $\{C_{l,i}\}$ and $\{C_{r,i}\}$ that satisfy $C_{l,i}\leq C_{r,i}$ for each $i$ along with the limits (4.52). The choice of these sequences can have a big influence on the overall performance of sFSP and especially the number of iterations it needs to terminate. Recall that $C_{l,i}$ and $C_{r,i}$ can be interpreted as bounds on the component-sum of states in the trapezoidal truncation $\mathcal{E}_{i}$ (see Section 4.1). Therefore we can use crudely estimated values of the mean and standard deviation of the state component-sum at stationary, as a guidance for selecting these cut-off sequences. These crude estimates can be obtained with a few sample trajectories of the original CTMC generated with Gillespie’s SSA [4].

Since the two main steps of sFSP, viz. constructing the rate matrix $\overline{Q}_{i}$ and solving the linear-algebraic for $\overline{\pi}_{i}$, are performed on two separate computing platforms (C++ and Matlab), we will report the CPU times999All the computations for this simple implementation of sFSP were performed on an Apple machine with 2.9 GHz Intel Core i5 processor. for both these steps individually for each iteration $i$. The total CPU time needed for an iteration is approximately the sum of these two times, and we will plot it along with the convergence factor $\gamma_{i}$, as a function of the iteration counter $i$, to show how they change as the truncated state-space $\mathcal{E}_{i}$ expands in size. For the computation of convergence factors we choose $\beta_{i}=C_{r,i}$ in step 2 of Algorithm 1 (see Section 4.1).

4.3.1 Gene-expression network

Our first example is the gene-expression network given in [27], where molecules of the messenger RNA or mRNA (denoted by $M$) are created by a gene, and these mRNA molecules catalytically produce molecules of some protein (denoted by $P$). Molecules of both these species can degrade spontaneously. This two-species network has the following four reactions:

[TABLE]

The propensity functions are given by mass-action kinetics (2.3) and $\theta_{i}$-s denote the associated rate constants. We assume that the values of these rate constants are given by $\theta_{1}=50$, $\theta_{2}=4$, $\theta_{3}=0.5$ and $\theta_{4}=0.2$.

For the CTMC model of this network, the state-space $\mathcal{E}=\mathbb{N}^{2}_{0}$ is irreducible and so it can be enumerated with the Cantor Pairing function $\Phi_{2}$ (see Section 4.1). We apply sFSP to this network to obtain an estimate of the stationary probability distribution $\pi$. The cut-off sequences $\{C_{l,i}\}$ and $\{C_{r,i}\}$ that define the trapezoidal truncation $\mathcal{E}_{i}=\mathcal{T}(C_{l,i},C_{r,i})$ for iteration $i$ are chosen as

[TABLE]

where $\widehat{\mu}=2100$ and $\widehat{\sigma}=120$, are crudely estimated values of the mean and standard deviation of the state component-sum at stationarity, and these are obtained with a few SSA-generated trajectories of the CTMC. The designated state we select for sFSP is $(0,\widehat{\mu})$, which corresponds to [math] mRNA molecules and $\widehat{\mu}=2100$ protein molecules.

The performance of sFSP on the gene-expression network is summarized in Table 1, where for each iteration $i$, the cut-off values ($C_{l,i}$ and $C_{r,i}$), the truncated state-space size ($n_{i}=|\mathcal{E}_{i}|$), the convergence factor $\gamma_{i}$ and the CPU times for the two main sFSP steps are provided. One can see that sFSP terminated in $5$ iterations and overall it required 615 seconds of CPU time. To assess the accuracy of sFSP, we also estimate $\pi$ using $10^{6}$ CTMC trajectories simulated with SSA in the time-interval $[0,100]$. This SSA-based estimation was implemented in C++ and it needed 7246 seconds of CPU time which is much higher than the 615 seconds needed for sFSP.

Note that the size of the truncated state-space $n_{i}$ is increasing linearly with $i$ and hence the size of the $n_{i}\times n_{i}$ rate matrix $\overline{Q}_{i}$ is increasing quadratically with $i$. So we would expect the CPU time for constructing $\overline{Q}_{i}$ and solving the linear-algebraic system for $\overline{\pi}_{i}$, to also increase quadratically with $i$. However this is not the case and the two CPU times only increase linearly (see Table 1). This is because matrix $\overline{Q}_{i}$ is extremely sparse with only $n_{i}K$ non-zero entries, where $K=4$ is the number of reactions. This sparsity is exploited in our implementation of sFSP for both constructing the matrix and solving the linear-algebraic system.

The linear increase in the required CPU time can be seen from Figure 3A. Here the convergence factor $\gamma_{i}$ is also plotted in log-scale and the almost linear decay shows that the convergence factor drops exponentially to zero as the truncated state-space expands iteratively. Such an exponential decay is perhaps due to the fact that the joint stationary distribution is unimodal, as indicated by the contour plot in Figure 3B. This unimodality is also visible from the marginal distribution plots in Figure 3C. These sFSP-estimated marginal distribution plots are compared with the SSA-estimated distributions in Figure 3C and they show a close match.

4.3.2 Toggle-Switch network

We now consider the example of the genetic toggle-switch network proposed by Gardner et. al. [28]. This network has two species ${\bf X}_{1}$ and ${\bf X}_{2}$ that are competing by repressing each other’s production. This repression is modeled through propensities given by nonlinear Hill functions [29]. The network has four simple reactions

[TABLE]

where the propensity functions $\lambda_{i}$-s are given by

[TABLE]

Here $x_{1}$ and $x_{2}$ denote the copy-numbers of ${\bf X_{1}}$ and ${\bf X_{2}}$ respectively. For our computations we set $\alpha_{1}=500$, $\alpha_{2}=0.3$, $\alpha_{3}=200$, $\alpha_{4}=0.4$, $\beta=1.5$ and $\gamma=1$.

For the CTMC model of this network, the state-space $\mathcal{E}=\mathbb{N}^{2}_{0}$ is irreducible, and we apply sFSP with trapezoidal truncations using the cut-off sequences $\{C_{l,i}\}$ and $\{C_{r,i}\}$

[TABLE]

at iteration $i$, where $\widehat{\mu}=1110$ and $\widehat{\sigma}=500$, are crude SSA-based estimates of the mean and standard deviation of the state component-sum at stationarity. The designated state we select for sFSP is $(0,\widehat{\mu})$, which corresponds to [math] molecules of ${\bf X_{1}}$ and $\widehat{\mu}=1110$ molecules of ${\bf X_{2}}$.

The performance of sFSP on the Toggle-Switch network is summarized in Table 2, where for each iteration $i$, the cut-off values, the truncated state-space size, the convergence factor and the CPU times for the two main sFSP steps are provided. In this example, sFSP terminated in $5$ iterations and overall it required 322 seconds of CPU time. In comparison, the SSA-based estimation of $\pi$, implemented in C++, using $6\times 10^{6}$ CTMC trajectories simulated in the time-interval $[0,100]$, needed 42192 seconds of CPU time.

As in the previous example, the required CPU time increases almost linearly with iteration $i$ and the convergence factor $\gamma_{i}$ decreases slowly for the first four iterations and then plummets to nearly [math] in the fifth iteration (see Figure 4A). The contour plot for the joint stationary distribution estimated by sFSP is shown in Figure 4B and it indicates that this distribution is bimodal with each mode corresponding to one of the species being dominant. In Figure 4C we plot the sFSP-estimated marginal stationary distributions for the copy-numbers of the two species and compare them with the SSA-estimated marginal stationary distributions. One can clearly see that unlike sFSP, SSA fails to adequately capture the stationary distribution in the low-probability regions of the state-space even though a large sample of size $6$ million is used. These statistical errors and other numerical issues associated with computing very low probabilities, may explain the slight discrepancy in the sFSP and SSA estimated marginal distribution for species ${\bf X_{2}}$ (see Figure 4C).

4.3.3 Pap-Switch network

We now consider the Pap epigenetic switch whose finite-time CME was solved in [5] with the FSP method. This stochastic switch is responsible for deciding whether or not E. coli will develop hairlike structures called pili. The Pap-switch network is illustrated in Figure 5A and it consists of a single pap operon $G$ that can exist in four states $G_{1},G_{2},G_{3}$ and $G_{4}$ determined by the binding sites occupied by the leucine-responsive regulatory protein (LRP) molecules. When the operon is in state $G_{2}$, it can produce a local regulatory protein called PapI which represses the unbinding of the LRP molecules from the operon binding sites. This $PapI$ protein is allowed to degrade spontaneously at a certain rate. As in [5] we assume that the number of LRP molecules is fixed at $100$. The dynamics of the copy-numbers of the five species $G_{1},G_{2},G_{3},G_{4}$ and $PapI$ in the Pap-Switch network can be modeled with $10$ reactions described in Table 3.

For the CTMC model of this network, the state-space $\mathcal{E}=\mathcal{E}_{b}\times\mathbb{N}_{0}$ is irreducible, where

[TABLE]

is the finite set which contains the dynamics of the copy-numbers $(x_{1},x_{2},x_{3},x_{4})$ of the four operon states $G_{1},G_{2},G_{3}$ and $G_{4}$. The copy-numbers of $PapI$ can take values in the whole set of non-negative integers $\mathbb{N}_{0}$. The state-space of the form $\mathcal{E}=\mathcal{E}_{b}\times\mathbb{N}_{0}$ can be enumerated using the function $\Psi$ (see (4.49)) with $N_{b}=4$. Similarly the trapezoidal truncations $\mathcal{E}_{i}$-s can be defined as (4.51). In our application of sFSP for this network we construct these truncations using the cut-off sequences $\{C_{l,i}\}$ and $\{C_{r,i}\}$ specified by

[TABLE]

at iteration $i$, where $\widehat{\mu}=4$ and $\widehat{\sigma}=3$, are coarse SSA-based approximations of the mean and standard deviation of the $PapI$ copy-numbers. Note that due to the low copy-numbers involved we fix the left cut-off point $C_{l,i}$ to be zero for all the iterations. Also the designated state we select for sFSP is $(1,0,0,0,0)$, which corresponds to the operon being in state $G_{1}$ and $PapI$ having [math] molecules.

The performance of sFSP on the Pap-Switch network is summarized in Table 4, where for each iteration $i$, the cut-off values, the truncated state-space size, the convergence factor and the CPU times for the two main sFSP steps are provided. For this network, sFSP took $6$ iterations to terminate and overall it required only $0.304$ seconds of CPU time. By contrast, the SSA-based estimation of the stationary distribution, implemented in C++, with $10^{6}$ CTMC trajectories generated in the time-period $[0,100]$, required 134 seconds of CPU time.

In this example, the sizes of the truncated state-spaces are very small and so sFSP executes very quickly, causing the CPU times to vary non-monotonically with iteration $i$ while the convergence factor decreases almost exponentially (see Figure 5B). In Figure 5C we plot the sFSP-estimated stationary distributions for $PapI$ copy-numbers at each operon state $G_{1},G_{2},G_{3}$ and $G_{4}$. These are compared with the corresponding SSA-estimated stationary distributions and it can be seen from Figure 5C that the match is almost perfect.

4.3.4 Self-activated gene expression

We end this section with a simple but instructive example borrowed from [30]. Consider a gene whose protein output ${\bf X}$ can activate its own expression through a nonlinear feedback loop. A simple reaction network model for this would be

[TABLE]

where the propensity function for the degradation reaction is linear $\lambda_{2}(x)=\gamma x$ while the propensity function for the production reaction is given by a Hill-type function

[TABLE]

Here $x$ denotes the copy-number of protein ${\bf X}$. For our computations we set $k_{1}=20$, $k_{2}=125$, $\alpha=5$, $m=70$ and $\gamma=1$.

As there is only one species, the trapezoidal truncation $\mathcal{E}_{i}=\mathcal{T}(C_{l,i},C_{r,i})$ is simply the set $\mathcal{E}_{i}=\{C_{l,i},C_{l,i}+1,\dots,C_{r,i}\}$. We choose $C_{l,i}=0$ and $C_{r,i}=(5+i)$ at iteration $i$, and apply sFSP on this example with designated state [math]. The results are shown in Figure 6. Note that the stationary distribution is bimodal, with a small peak around $20$ and a larger peak around $145$. Most of the stationary probabilities are concentrated in two disjoint regions $R_{1}=\{0,\dots,55\}$ and $R_{2}=\{80,\dots,210\}$ around the two peaks. Observe that the end-points $x_{1}=55$ and $x_{2}=210$ of these two regions are inflection or turning points for the behavior of the convergence factor $\gamma_{i}$ with increasing iteration counter $i$ or expanding truncated state-space $\mathcal{E}_{i}$. The convergence factor $\gamma_{i}$ decays exponentially before $x_{1}$ and after $x_{2}$, but in the intermediate region $I=\{x_{1}+1,\dots,x_{2}-1\}$ it shows a gradual increase. Further computations reveal that for iterations corresponding to this intermediate region, the outflow rate $r^{(i)}_{\textnormal{out}}$ remains approximately constant, and so the convergence factor $\gamma_{i}$ increases slowly due to scaling by the cut-off value $C_{r,i}$. This relationship between bimodality of the stationary distribution and non-monotonicity of the convergence factor $\gamma_{i}$ is very interesting and should be investigated in a greater detail elsewhere.

5 sFSP Algorithm: QTT Implementation

The second implementation is motivated by the recently developed Quantized Tensor-Train (QTT) version of FSP [7], which works with QTT representations of the transition rate matrix and its stationary distribution vector. The use of such representations expands the range of applicability of sFSP and we demonstrate this by applying sFSP on a network which is much larger than the networks considered in Section 4.3.

5.1 The CME in QTT form

A tensor is essentially a multi-dimensional generalization of a two-dimensional matrix or a one-dimensional vector. A $d$-dimensional tensor $T$ of size ${\bm{n}}=n_{1}\times\dots\times n_{d}$, represents a structured collection of real numbers given by

[TABLE]

Each dimension of this tensor $T$ is also called its mode, and $n_{1},\ldots,n_{d}$ denote the mode sizes. The tensor $T$ can also be viewed as a real-valued function over the $d$-dimensional hyper-rectangle

[TABLE]

which is a subset of the non-negative integer orthant $\mathbb{N}^{d}_{0}$.

Tensors are particularly well suited to express the CME since the system already has a physical interpretation as tensors, where each species corresponds to one tensor mode and for any mode $k$, its size $n_{k}$ serves as the strict upper-bound for the allowable copy-numbers for species $\mathbf{X}_{k}$. As in FSP [5], consider a CME over the truncated state-space $\mathcal{E}_{\bm{n}}$ (see (2.11) for example). The probability distribution $p_{\bm{n}}(t)$ of the random state-vector at time $t$ can be represented as a $d$-dimensional tensor of size ${\bf n}$ and the matrix $Q^{T}_{\bm{n}}$ that captures its rate of change can be represented as a $2d$-dimensional tensor of size ${\bf n}\times{\bf n}$.

The tensor train (TT) representation of a $d$-dimensional tensor $T$ with size ${\bm{n}}=n_{1}\times\dots\times n_{d}$ is given by

[TABLE]

where $r_{0}=r_{d}=1$ and for each $j=1,\dots,d$, $U_{j}$ is a three-dimensional tensor with size $r_{j-1}\times n_{j}\times r_{j}$. The tensors $U_{1}$ to $U_{d}$ are called the core tensors and $r_{1},\ldots r_{d-1}$ are referred to as tensor ranks. The TT-representation can potentially provide a high compression of the tensor, especially if the ranks are low. Most basic matrix-vector operations (like matrix-vector product, dot product, outer product etc.) can be applied directly on the compressed TT format (for details see [31]). The complexity of these basic operations as well as the storage cost can be bound by $n_{\textnormal{max}}r_{\textnormal{max}}^{2}d$ where $n_{\textnormal{max}}=\max\{n_{1},\ldots n_{d}\}$ and $r_{\textnormal{max}}=\max\{r_{1},\ldots,r_{d-1}\}$. Any tensor can be decomposed into the TT format by the TT-SVD algorithm [31], which is based on the Singular Value Decomposition (SVD) for matrices. The TT format can be extended to the quantized tensor train (QTT) format which provides another layer of compression by dividing each mode of the tensor into several virtual modes that are then further compressed using tensor trains (see [32] and [33]).

In [7] the authors show how the matrix $Q^{T}_{\bm{n}}$ for the CME (2.11) over the truncated state-space $\mathcal{E}_{\bm{n}}$ (5.54) can be directly constructed in the QTT format and thereafter used for efficiently solving the FSP and obtaining the transient CME solution $p_{\bm{n}}(t)$. The main observation underlying the QTT construction of $Q^{T}_{\bm{n}}$ is that one can think of this matrix in terms of the spatial shift operator $\bm{S}_{\zeta_{k}}$, shifting a probability density tensor $p$ by the stoichiometry vector $\zeta_{k}$ for reaction $k$, and a multiplication operator $\bm{M}_{\lambda_{k}}$, multiplying a probability density tensor $p$ by the propensity function $\lambda_{k}$ for reaction $k$, i.e.

[TABLE]

for any $x\in\mathcal{E}_{\bm{n}}$. Using these operators along with the identity operator $\mathbb{I}$, the matrix $Q^{T}_{\bm{n}}$ can be expressed as

[TABLE]

and this form can be exploited for efficiently constructing the QTT representation of $Q^{T}_{\bm{n}}$. As explained in [7], for mass action kinetics, the operator $\bm{M}_{\lambda}$ can be constructed by taking the outer products of state-vectors in $\mathcal{E}_{\bm{n}}$ and the appropriate vector of ones $\mathbf{1}$, while the operator $\bm{S}_{\zeta_{k}}$ can be constructed as a matrix of zeros with a shifted diagonal of ones.

5.2 Implementation Details

In our QTT implementation of sFSP, we use a similar expression as (5.55) to construct the QTT representation of the transpose $\overline{Q}^{T}_{\bm{n}}$ of the transition rate matrix $\overline{Q}_{\bm{n}}$ (see (3.16)) for our projected CTMC over the truncated state-space $\mathcal{E}_{\bm{n}}$, where all the outgoing transitions are redirected to the designated state ${\bf 0}$ of all zeros. Using the QTT representation of $\overline{Q}^{T}_{\bm{n}}$, the corresponding linear-algebraic system (1.1) is directly solved in QTT format to yield the stationary probability distribution $\overline{\pi}_{\bm{n}}$ in QTT format.

For solving the linear-algebraic system, we use the inverse iteration approach (see [34]) which is known to have very good convergence properties and work well with tensor algebra [35]. In this approach a linear system of the form $Ax={\bf 0}$, for a singular matrix $A$, is solved by iteratively solving the linear systems

[TABLE]

starting with some initial guess $x_{0}$. The solution $x_{j}$ is suitably normalized before commencing iteration $(j+1)$. Generally this procedure requires very few iterations (like 2 or 3) to converge, and this convergence can be judged by checking that the distance between subsequent solutions $\|x_{j}-x_{j-1}\|$ is below some threshold level $\delta$.

In our setup we implement the QTT version of sFSP method (Algorithm 1) in Matlab, using Version 2.2 of the qtt-toolbox developed by I. Oseledets, S. Dolgov, V. Kazeev, O. Lebedeva, and T. Mach [36]. In particular the linear systems that arise in the inverse iteration procedure are solved using the function dmrg_solve3.m from this toolbox. For each sFSP iteration $i$, the initial guess for the inverse iteration procedure is chosen based on the estimate obtained in iteration $(i-1)$, as mentioned in Section 4.2. For the computational example we consider next, we found that only two inverse iterations were always sufficient to yield a convergent solution of the linear-algebraic system (1.1) for the threshold level $\delta=10^{-4}$.

5.3 A Computational Example

We now illustrate our QTT implementation on a toy example with features similar to the Repressilator network given by Elowitz and Leibler [37], which has three gene-expression modules (say A, B and C) that interact by mutual inhibition of each other in a cyclic fashion i.e. $A$ represses $B$, $B$ represses $C$ and $C$ represses $A$ (see Figure 7A). This inhibition is carried out by the corresponding proteins ($P_{A}$, $P_{B}$ and $P_{C}$) and it is achieved by enhancing the rate at which the inhibited gene becomes inactive (OFF) from an active (ON) state. Each protein also activates its own production by increasing the rate at which its gene switches ON from the OFF state. The mRNAs ($M_{A}$, $M_{B}$ and $M_{C}$) associated with the genes are only transcribed when the corresponding gene is in the ON state. Overall this network consists of $9$ species and $18$ reactions described in Table 5. These $9$ species include the indicators for the three genes being in the ON state ($G_{A}^{1}$, $G_{B}^{1}$ and $G_{C}^{1}$), the three mRNAs ($M_{A}$, $M_{B}$ and $M_{C}$) and finally the three proteins ($P_{A}$, $P_{B}$ and $P_{C}$).

For the CTMC model of this network, the state-space $\mathcal{E}=\mathcal{E}_{b}\times\mathbb{N}^{6}_{0}$ is irreducible, where

[TABLE]

is the finite set which contains the dynamics of the copy-numbers $(x_{1},x_{2},x_{3})$ of the three genes being in the ON state. The copy-numbers of all the mRNAs and proteins can take values in the whole set of non-negative integers $\mathbb{N}_{0}$. We apply sFSP on the 3-gene network with the finite truncated state-space $\mathcal{E}_{i}$ for sFSP iteration $i$ chosen as $\mathcal{E}_{i}=\mathcal{E}_{\bm{n}_{i}}$ (see (5.54)) with

[TABLE]

Here $U_{m,i}$ and $U_{p,i}$ denote the strict upper-bounds for the copy-numbers of all the mRNAs and proteins respectively. The convergence factor $\gamma_{i}$ is computed for this example using $\beta_{i}=\max\{U_{m,i},U_{p,i}\}$ in step 2 of Algorithm 1. Due to limitations posed by the qtt-toolbox and our computational hardware, we fix the acceptable threshold level $\epsilon$ (see step 4 of Algorithm 1) to be $10^{-2}$ instead of $10^{-10}$ used previously.

The performance of sFSP on this triple-repressor network is summarized in Table 6, where for each iteration $i$, the upper-bounds ($U_{m,i}$ and $U_{p,i}$), the truncated state-space size ($|\mathcal{E}_{i}|$), the convergence factor $\gamma_{i}$ and the CPU times are provided. One can see that sFSP terminated in $5$ iterations and overall it required around 168 minutes of CPU time101010All the computations for this QTT implementation of sFSP were performed on a Lenovo T440 machine with 1.6 GHz Intel i5-4200U processor with 8GB of RAM. Note that the copy-numbers of all the species are relatively small in this example. However due to the large number of species, the size of the final truncated state-space $\mathcal{E}_{5}$ is several times larger than the truncated state-spaces encountered in the examples considered before, for the simple implementation of sFSP. To assess the accuracy of sFSP, we also estimate $\pi$ using $10^{6}$ CTMC trajectories simulated with SSA in the time-interval $[0,200]$. As in the previous examples, we plot the CPU times and the convergence factors at all the sFSP iterations in Figure 7B, the contour plots for the various joint stationary distributions estimated by sFSP in Figure 7C, and the estimated marginal stationary distributions for the all the $9$ species in Figure 8. These marginal stationary distributions are also compared with the corresponding SSA-estimated marginal stationary distributions and one can see that the match is quite good.

The SSA-based estimation with $10^{6}$ trajectories needed around 117 minutes of CPU time, based on a C++ implementation, which is slightly faster than sFSP (168 minutes). However we must note that even though this SSA-based estimation captures the marginal distributions very well (see Figure 8), it is unable to capture the full stationary distribution because the state-space is high-dimensional and the size of the final truncated state-space $\mathcal{E}_{5}$ for sFSP suggests that the support of the true stationary distribution is much larger ($>$130 million) than the number of SSA samples ($1$ million) being used for the estimation. To illustrate this point, we compute the $\ell_{1}$ distance between the sFSP estimated stationary distribution $\overline{\pi}$ and the stationary distribution $\widehat{\pi}$ estimated with $10^{5},10^{6}$ and $10^{7}$ SSA samples. The results are shown in Table 7 along with the associated CPU times for generating the SSA samples. Notice that as the number of SSA samples increases, the $\ell_{1}$ distance $\|\overline{\pi}-\widehat{\pi}\|_{\ell_{1}}$ decreases sharply, which strongly suggests that sFSP is an accurate approximation of the true stationary distribution $\pi$. However this $\ell_{1}$ distance is significant when $\widehat{\pi}$ is estimated with $1$ million SSA samples, which implies that $\widehat{\pi}$ is quite inaccurate. If we use $10^{7}$ SSA samples to estimate $\widehat{\pi}$ then the accuracy improves but the total CPU time required is approximately 18 hours, that is $6.4$ times larger than the time needed for sFSP.

6 Conclusion

In this paper we presented a new method for estimating the stationary probability distributions of continuous-time Markov chain (CTMC) models of reaction networks based on suitable truncations of the CME. The method which we call the stationary Finite State Projection (sFSP) algorithm is similar to the Finite State Projection (FSP) algorithm[5], with the crucial difference being that instead of introducing an absorbing state, we redirect all the outgoing transitions from the truncated state-space to a designated state within the truncated state-space (see Figure 1C). This simple modification creates a projected CTMC over the truncated state-space, whose stationary distribution can be obtained by solving a finite linear-algebraic system. We provided theoretical arguments to establish that this stationary distribution estimated from the projected CTMC is unique, converges to the true stationary distribution as the truncated state-space expands to the full state-space and for any truncated state-space the error between the estimated stationary distribution and the true stationary distribution can be assessed by computing the overall rate of outgoing transitions at the estimated stationary distribution (see Theorem 3.1). These results form the basis of our sFSP method. We illustrated the efficiency and accuracy of this method using several examples. These examples indicated that sFSP can easily outperform the stochastic simulation-based approach for estimating the stationary distribution, both in terms of computational speed as well as accuracy. This is not unexpected, as stochastic simulations are expensive to perform over large time-intervals, and the stationary distribution they estimate suffers from statistical errors that can be significant in regions of the state-space where the probabilities are extremely low. These issues do not arise in sFSP and this makes it an appealing method for estimating stationary distributions of CTMCs representing reaction networks.

There are several ways to improve and extend sFSP. Like FSP, this method is iterative in nature and the number of iterations it requires to converge depends on the specifics of the implementation of sFSP. In this paper we discussed two such implementations. In the first implementation the state-space was explicitly enumerated with Cantor pairing functions and then truncated in trapezoidal shapes (see Section 4), while the second implementation was based on the recently developed quantized tensor train (QTT) version of CME where each state-space truncation is a hyper-rectangle (see Section 5). Both these implementations will benefit from better state-space truncation schemes that adapt to the problem at hand. One way to do this would be to use Lyapunov function theory or use stationary moment bounds to construct optimal state-space truncations (see [18], [38] and [16]). Observe that unlike FSP which solves a linear system of ODEs, sFSP only requires solving a linear-algebraic system which is computationally much easier. Hence sFSP can handle a wider range of networks in comparison to FSP. Indeed with the QTT implementation, finding stationary distributions for problems with state truncations exceeding 100 million states was shown to be feasible. A possible approach for enhancing the feasibility of sFSP to even larger problems would be to integrate it with sparse grids and aggregation methods [8]. Note that at the core of sFSP, is the problem of finding vectors in the one-dimensional null-spaces of large, but extremely sparse matrices (see Section 4.2). This sparsity and the structure of the matrices that arise, make this problem quite amenable to parallel-computing approaches [39].

Acknowledgments

The authors would like to thank Prof. Sean Meyn (University of Florida) and Prof. Brian Munsky (Colorado State University) for their helpful comments and suggestions.