A finite state projection method for steady-state sensitivity analysis of stochastic reaction networks

TL;DR

This paper introduces a new computational method for estimating steady-state sensitivities in stochastic reaction networks, avoiding time-consuming simulations by leveraging the stationary Finite State Projection algorithm and solving a Poisson equation.

Contribution

It presents a novel approach combining the stationary FSP and Poisson equation solutions to efficiently compute steady-state sensitivities without simulations.

Findings

Accurately estimates steady-state sensitivities

Reduces computational time compared to simulation-based methods

Validated on multiple example networks

Abstract

Consider the standard stochastic reaction network model where the dynamics is given by a continuous-time Markov chain over a discrete lattice. For such models, estimation of parameter sensitivities is an important problem, but the existing computational approaches to solve this problem usually require time-consuming Monte Carlo simulations of the reaction dynamics. Therefore these simulation-based approaches can only be expected to work over finite time-intervals, while it is often of interest in applications to examine the sensitivity values at the steady-state after the Markov chain has relaxed to its stationary distribution. The aim of this paper is to present a computational method for the estimation of steady-state parameter sensitivities, which instead of using simulations, relies on the recently developed stationary Finite State Projection (sFSP) algorithm [J. Chem. Phys. 147, 154101 (2017)] that provides an accurate estimate of the stationary distribution at a fixed set of parameters. We show that sensitivity values at these parameters can be estimated from the solution of a Poisson equation associated with the infinitesimal generator of the Markov chain. We develop an approach to numerically solve the Poisson equation and this yields an efficient estimator for steady-state parameter sensitivities. We illustrate this method using several examples.

Tables

Table 1. Table A.1: Reactions in the gene-expression network. The propensities are parametrized by θ = ( θ 1 , θ 2 , θ 3 , θ 4 ) 𝜃 subscript 𝜃 1 subscript 𝜃 2 subscript 𝜃 3 subscript 𝜃 4 \theta=(\theta_{1},\theta_{2},\theta_{3},\theta_{4}) .

No.	Reaction	Propensity
1	$\varnothing \to {𝒮}_1$	${\lambda }_1(x,\theta )={\theta }_1$
2	${𝒮}_1\to {𝒮}_1+{𝒮}_2$	${\lambda }_2(x,\theta )={\theta }_2{x}_1$
3	${𝒮}_1\to \varnothing$	${\lambda }_3(x,\theta )={\theta }_3{x}_1$
4	${𝒮}_2\to \varnothing$	${\lambda }_4(x,\theta )={\theta }_4{x}_2$

Table 2. Table A.2: Truncated state-spaces and designated states used in the gene-expression network. The cut-off values ( C l subscript 𝐶 𝑙 C_{l} and C r subscript 𝐶 𝑟 C_{r} ) for the trapezoidal truncations ( 4.52 ) are provided.

No.	Cut-offs		State-space size	Designated state
$i$	$C_{l,i}$	$C_{r,i}$	$n_i$	$x_{\ell }$
$1$	$3280$	$4280$	$3,784,784$	$(10,3270)$
$2$	$2780$	$4780$	$7,565,784$	$(10,2770)$
$3$	$2280$	$5280$	$11,346,784$	$(10,2270)$
$4$	$1780$	$5780$	$15,127,784$	$(10,1770)$
$5$	$1280$	$6280$	$18,908,784$	$(10,1270)$

Table 3. Table A.3: Sensitivity analysis for the gene-expression network. This table lists the PE-estimated sensitivity values w.r.t. all four parameters θ 1 , … , θ 4 subscript 𝜃 1 … subscript 𝜃 4 \theta_{1},\dots,\theta_{4} for two objective functions f ​ ( x ) = x 1 𝑓 𝑥 subscript 𝑥 1 f(x)=x_{1} and f ​ ( x ) = x 2 𝑓 𝑥 subscript 𝑥 2 f(x)=x_{2} with expected steady-state values 𝔼 θ ​ ( x 1 ) = 180 subscript 𝔼 𝜃 subscript 𝑥 1 180 \mathbb{E}_{\theta}(x_{1})=180 and 𝔼 θ ​ ( x 2 ) = 3600 subscript 𝔼 𝜃 subscript 𝑥 2 3600 \mathbb{E}_{\theta}(x_{2})=3600 respectively. The corresponding percentage relation error RE % percent \% ( 4.51 ) for the sensitivity estimates is also provided. To perform all these computations the total processing time was 1142 1142 1142 seconds on four processors and the total memory requirement was 166’278 MB.

Objective	Sensitive	Sensitivity	RE%	Sensitive	Sensitivity	RE%
Function	Parameter	Values (PE)		Parameter	Values (PE)
$x_1$	${\theta }_1$	$2.0000$	$1.9609\cdot {10}^{-8}$	${\theta }_3$	$-360.0000$	$1.8802\cdot {10}^{-8}$
	${\theta }_2$	$1.2186\cdot {10}^{-8}$	$1.2186\cdot {10}^{-8}$	${\theta }_4$	$-2.5096\cdot {10}^{-7}$	$2.5096\cdot {10}^{-7}$
$x_2$	${\theta }_1$	$40.000$	$7.6378\cdot {10}^{-8}$	${\theta }_3$	$-7200.0$	$7.5121\cdot {10}^{-8}$
	${\theta }_2$	$900.00$	$9.8708\cdot {10}^{-8}$	${\theta }_4$	$-18000$	$9.9957\cdot {10}^{-10}$

Table 4. Table A.4: Reactions in the Toggle-Switch network. The propensities are parametrized by θ = ( θ 1 , θ 2 , θ 3 , θ 4 , θ 5 , θ 6 ) 𝜃 subscript 𝜃 1 subscript 𝜃 2 subscript 𝜃 3 subscript 𝜃 4 subscript 𝜃 5 subscript 𝜃 6 \theta=(\theta_{1},\theta_{2},\theta_{3},\theta_{4},\theta_{5},\theta_{6}) .

No.	Reaction	Propensity
1	$\varnothing \to {𝒮}_1$	${\lambda }_1(x,\theta )=\frac{{\theta }_1}{1+x_2^{{\theta }_5}}$
2	${𝒮}_1\to \varnothing$	${\lambda }_2(x,\theta )={\theta }_2{x}_1$
3	$\varnothing \to {𝒮}_2$	${\lambda }_3(x,\theta )=\frac{{\theta }_3}{1+x_1^{{\theta }_6}}$
4	${𝒮}_2\to \varnothing$	${\lambda }_4(x,\theta )={\theta }_4{x}_2$

Table 5. Table A.5: Truncated state-spaces and designated states used in the toggle-switch network. The cut-off values ( C l subscript 𝐶 𝑙 C_{l} and C r subscript 𝐶 𝑟 C_{r} ) for the trapezoidal truncations ( 4.52 ) are provided.

No.	Cut-offs		State-space size	Designated state
$i$	$C_{l,i}$	$C_{r,i}$	$n_i$	$x_{\ell }$
$1$	$0$	$360$	$65341$	$(235,115)$
$2$	$0$	$860$	$371,091$	$(235,115)$
$3$	$0$	$1360$	$926,841$	$(235,115)$
$4$	$0$	$2360$	$2,788,341$	$(235,115)$
$5$	$0$	$3360$	$5,649,841$	$(235,115)$
$6$	$0$	$4360$	$9,511,341$	$(235,115)$

Table 6. Table A.6: Sensitivity analysis for the toggle-switch network. This table lists the PE-estimated sensitivity values w.r.t. all six parameters θ 1 , … , θ 6 subscript 𝜃 1 … subscript 𝜃 6 \theta_{1},\dots,\theta_{6} for two objective functions f ​ ( x ) = x 1 𝑓 𝑥 subscript 𝑥 1 f(x)=x_{1} and f ​ ( x ) = x 2 𝑓 𝑥 subscript 𝑥 2 f(x)=x_{2} with expected steady-state values 𝔼 θ ​ ( x 1 ) ≈ 0.015148 subscript 𝔼 𝜃 subscript 𝑥 1 0.015148 \mathbb{E}_{\theta}(x_{1})\approx 0.015148 and 𝔼 θ ​ ( x 2 ) = 496.23 subscript 𝔼 𝜃 subscript 𝑥 2 496.23 \mathbb{E}_{\theta}(x_{2})=496.23 respectively. The corresponding percentage relation error RE % percent \% ( 4.51 ) for the sensitivity estimates is also provided. To perform all these computations the total processing time was 782 782 782 seconds on a single processor and the total memory requirement was 45’050 MB.

Objective	Sensitive	Sensitivity	Sensitive	Sensitivity
Function	Parameter	Values (PE)	Parameter	Values (PE)
$x_1$	${\theta }_1$	$3.0677\cdot {10}^{-5}$	${\theta }_4$	$0.057690$
	${\theta }_2$	$-5.1177\cdot {10}^{-3}$	${\theta }_5$	$-0.095117$
	${\theta }_3$	$-1.1531\cdot {10}^{-4}$	${\theta }_6$	$6.1345\cdot {10}^{-7}$
$x_2$	${\theta }_1$	$-7.5857\cdot {10}^{-3}$	${\theta }_4$	$-1254.8$
	${\theta }_2$	$1.2645$	${\theta }_5$	$23.522$
	${\theta }_3$	$2.5095$	${\theta }_6$	$-9.8566\cdot {10}^{-3}$

Table 7. Table A.7: Reactions in the deficiency zero network. The propensities are parametrized by the vector θ = ( θ 1 , θ 2 , θ 3 , θ 4 , θ 5 , θ 6 ) 𝜃 subscript 𝜃 1 subscript 𝜃 2 subscript 𝜃 3 subscript 𝜃 4 subscript 𝜃 5 subscript 𝜃 6 \theta=(\theta_{1},\theta_{2},\theta_{3},\theta_{4},\theta_{5},\theta_{6}) .

No.	Reaction	Propensity
1	$\varnothing \to {𝒮}_1$	${\lambda }_1(x,\theta )={\theta }_1$
2	${𝒮}_1\to \varnothing$	${\lambda }_2(x,\theta )={\theta }_2{x}_1$
3	$\varnothing \to {𝒮}_2$	${\lambda }_3(x,\theta )={\theta }_3$
4	${𝒮}_2\to \varnothing$	${\lambda }_4(x,\theta )={\theta }_4{x}_2$
5	$2{𝒮}_2\to {𝒮}_4$	${\lambda }_5(x,\theta )={\theta }_5{x}_2\left(x_2-1\right)$
6	${𝒮}_4\to 2{𝒮}_2$	${\lambda }_6(x,\theta )={\theta }_6{x}_4$
7	${𝒮}_1+{𝒮}_2\to {𝒮}_3$	${\lambda }_7(x,\theta )={\theta }_7{x}_1{x}_2$
8	${𝒮}_3\to {𝒮}_1+{𝒮}_2$	${\lambda }_8(x,\theta )={\theta }_8{x}_3$

Table 8. Table A.8: Truncated state-spaces and designated states used in the deficiency zero network. The cut-off values ( C l subscript 𝐶 𝑙 C_{l} and C r subscript 𝐶 𝑟 C_{r} ) for the trapezoidal truncations ( 4.52 ) are provided.

No.	Cut-offs		State-space size	Designated state
$i$	$C_{l,i}$	$C_{r,i}$	$n_i$	$x_{\ell }$
$1$	$0$	$15$	$3,876$	$(10,0,0,0)$
$2$	$0$	$20$	$10,626$	$(10,0,0,0)$
$3$	$0$	$25$	$23,751$	$(10,0,0,0)$
$4$	$0$	$30$	$46,376$	$(10,0,0,0)$
$5$	$0$	$35$	$82,251$	$(10,0,0,0)$
$6$	$0$	$40$	$135,751$	$(10,0,0,0)$
$7$	$0$	$45$	$211,876$	$(10,0,0,0)$
$8$	$0$	$50$	$316,251$	$(10,0,0,0)$
$9$	$0$	$55$	$455,126$	$(10,0,0,0)$

Table 9. Table A.9: Sensitivity analysis for the deficiency zero network. This table lists the PE-estimated sensitivity values w.r.t. all eight parameters θ 1 , … , θ 8 subscript 𝜃 1 … subscript 𝜃 8 \theta_{1},\dots,\theta_{8} for four objective functions f ​ ( x ) = x 1 𝑓 𝑥 subscript 𝑥 1 f(x)=x_{1} , f ​ ( x ) = x 2 𝑓 𝑥 subscript 𝑥 2 f(x)=x_{2} , f ​ ( x ) = x 3 𝑓 𝑥 subscript 𝑥 3 f(x)=x_{3} and f ​ ( x ) = x 4 𝑓 𝑥 subscript 𝑥 4 f(x)=x_{4} with expected steady-state values 𝔼 θ ​ ( x 1 ) = 5.625 , 𝔼 θ ​ ( x 2 ) = 5 , 𝔼 θ ​ ( x 3 ) = 1.0547 formulae-sequence subscript 𝔼 𝜃 subscript 𝑥 1 5.625 formulae-sequence subscript 𝔼 𝜃 subscript 𝑥 2 5 subscript 𝔼 𝜃 subscript 𝑥 3 1.0547 \mathbb{E}_{\theta}(x_{1})=5.625,\mathbb{E}_{\theta}(x_{2})=5,\mathbb{E}_{\theta}(x_{3})=1.0547 and 𝔼 θ ​ ( x 4 ) = 1.3636 subscript 𝔼 𝜃 subscript 𝑥 4 1.3636 \mathbb{E}_{\theta}(x_{4})=1.3636 respectively. The corresponding percentage relation error RE % percent \% ( 4.51 ) for the sensitivity estimates is also provided. To perform all these computations the total processing time was 4386 4386 4386 seconds on a single processor and the total memory requirement was 26’458 MB.

Objective	Sensitive	Sensitivity	RE%	Sensitive	Sensitivity	RE%
Function	Parameter	Values (PE)		Parameter	Values (PE)
$x_1$	${\theta }_1$	$1.2450$	$8.0956\cdot {10}^{-8}$	${\theta }_5$	$-3.2680\cdot {10}^{-9}$	$3.2680\cdot {10}^{-9}$
	${\theta }_2$	$-7.0312$	$8.2053\cdot {10}^{-8}$	${\theta }_6$	$1.7825\cdot {10}^{-10}$	$1.7825\cdot {10}^{-10}$
	${\theta }_3$	$-1.5332\cdot {10}^{-9}$	$1.5332\cdot {10}^{-9}$	${\theta }_7$	$-5.3773\cdot {10}^{-10}$	$5.3773\cdot {10}^{-10}$
	${\theta }_4$	$7.6662\cdot {10}^{-9}$	$7.6662\cdot {10}^{-9}$	${\theta }_8$	$2.0159\cdot {10}^{-11}$	$2.0159\cdot {10}^{-11}$
$x_2$	${\theta }_1$	$-1.4302\cdot {10}^{-9}$	$1.4302\cdot {10}^{-9}$	${\theta }_5$	$-7.7320\cdot {10}^{-9}$	$7.7320\cdot {10}^{-9}$
	${\theta }_2$	$8.0446\cdot {10}^{-9}$	$8.0446\cdot {10}^{-9}$	${\theta }_6$	$4.2174\cdot {10}^{-10}$	$4.2174\cdot {10}^{-10}$
	${\theta }_3$	$1.0000$	$3.5816\cdot {10}^{-7}$	${\theta }_7$	$-9.3908\cdot {10}^{-10}$	$9.3908\cdot {10}^{-10}$
	${\theta }_4$	$-5.0000$	$3.5429\cdot {10}^{-7}$	${\theta }_8$	$3.5205\cdot {10}^{-11}$	$3.5205\cdot {10}^{-10}$
$x_3$	${\theta }_1$	$0.23438$	$1.0319\cdot {10}^{-7}$	${\theta }_5$	$-7.6189\cdot {10}^{-9}$	$7.6189\cdot {10}^{-9}$
	${\theta }_2$	$-1.3184$	$1.0260\cdot {10}^{-7}$	${\theta }_6$	$4.1558\cdot {10}^{-10}$	$4.1558\cdot {10}^{-10}$
	${\theta }_3$	$0.21094$	$1.8204\cdot {10}^{-8}$	${\theta }_7$	$0.35156$	$6.7294\cdot {10}^{-9}$
	${\theta }_4$	$-1.0547$	$1.3920\cdot {10}^{-8}$	${\theta }_8$	$-0.013184$	$6.6909\cdot {10}^{-9}$
$x_4$	${\theta }_1$	$4.1305\cdot {10}^{-10}$	$4.1305\cdot {10}^{-10}$	${\theta }_5$	$2.2727$	$4.8438\cdot {10}^{-7}$
	${\theta }_2$	$-2.3228\cdot {10}^{-9}$	$2.3228\cdot {10}^{-9}$	${\theta }_6$	$-0.12397$	$4.8415\cdot {10}^{-7}$
	${\theta }_3$	$0.54545$	$9.6426\cdot {10}^{-7}$	${\theta }_7$	$-5.3900\cdot {10}^{-10}$	$5.3900\cdot {10}^{-10}$
	${\theta }_4$	$-2.7273$	$9.6111\cdot {10}^{-7}$	${\theta }_8$	$2.0204\cdot {10}^{-11}$	$2.0204\cdot {10}^{-11}$

Table 10. Table A.10: Reactions in the simple linear network. The propensities are parametrized by the vector θ = ( θ 1 , θ 2 , θ 3 , θ 4 ) 𝜃 subscript 𝜃 1 subscript 𝜃 2 subscript 𝜃 3 subscript 𝜃 4 \theta=(\theta_{1},\theta_{2},\theta_{3},\theta_{4}) .

No.	Reaction	Propensity
1	${𝒮}_1\to {𝒮}_2$	${\lambda }_1(x,\theta )={\theta }_1{x}_1$
2	${𝒮}_2\to {𝒮}_1$	${\lambda }_2(x,\theta )={\theta }_2{x}_2$
3	${𝒮}_2\to {𝒮}_3$	${\lambda }_3(x,\theta )={\theta }_3{x}_2$
4	${𝒮}_3\to {𝒮}_2$	${\lambda }_4(x,\theta )={\theta }_4{x}_3$

Table 11. Table A.11: Sensitivity analysis for the simple linear network. This table lists the PE-estimated sensitivity values w.r.t. all four parameters θ 1 , … , θ 4 subscript 𝜃 1 … subscript 𝜃 4 \theta_{1},\dots,\theta_{4} for three objective functions f ​ ( x ) = x 1 𝑓 𝑥 subscript 𝑥 1 f(x)=x_{1} , f ​ ( x ) = x 2 𝑓 𝑥 subscript 𝑥 2 f(x)=x_{2} and f ​ ( x ) = x 3 𝑓 𝑥 subscript 𝑥 3 f(x)=x_{3} with expected steady-state values 𝔼 θ ​ ( x 1 ) = 4.4444 , 𝔼 θ ​ ( x 2 ) = 2.2222 formulae-sequence subscript 𝔼 𝜃 subscript 𝑥 1 4.4444 subscript 𝔼 𝜃 subscript 𝑥 2 2.2222 \mathbb{E}_{\theta}(x_{1})=4.4444,\mathbb{E}_{\theta}(x_{2})=2.2222 and 𝔼 θ ​ ( x 3 ) = 3.3333 subscript 𝔼 𝜃 subscript 𝑥 3 3.3333 \mathbb{E}_{\theta}(x_{3})=3.3333 respectively. The corresponding percentage relation error RE % percent \% ( 4.51 ) for the sensitivity estimates is also provided. To perform all these computations the total processing time was only 0.2 0.2 0.2 seconds on a single processor and the total memory requirement was 105 MB.

Objective	Sensitive	Sensitivity	RE%	Sensitive	Sensitivity	RE%
Function	Parameter	Values (PE)		Parameter	Values (PE)
$x_1$	${\theta }_1$	$-0.24691$	$1.0137\cdot {10}^{-5}$	${\theta }_3$	$-49.383$	$1.8538\cdot {10}^{-5}$
	${\theta }_2$	$0.12346$	$1.0242\cdot {10}^{-5}$	${\theta }_4$	$74.074$	$9.3173\cdot {10}^{-6}$
$x_2$	${\theta }_1$	$0.098765$	$8.0763\cdot {10}^{-6}$	${\theta }_3$	$-24.691$	$1.8642\cdot {10}^{-5}$
	${\theta }_2$	$-0.049383$	$8.1807\cdot {10}^{-6}$	${\theta }_4$	$37.037$	$9.2132\cdot {10}^{-6}$
$x_3$	${\theta }_1$	$0.14815$	$1.1511\cdot {10}^{-5}$	${\theta }_3$	$74.074$	$1.8572\cdot {10}^{-5}$
	${\theta }_2$	$-0.074074$	$1.1616\cdot {10}^{-5}$	${\theta }_4$	$-111.11$	$9.2826\cdot {10}^{-6}$

Equations

$$S_{\theta }(f,t):=\frac{\partial }{\partial \theta }\mathbb{E}(f(X_{\theta }(t))).$$

$$S_{\theta }(f):=\underset{t\to \infty }{\lim }\frac{\partial }{\partial \theta }\mathbb{E}(f(X_{\theta }(t))).$$

$$S_{\theta }(f)=\frac{\partial }{\partial \theta }⟨f,{\pi }_{\theta }⟩,$$

$$⟨f,{\pi }_{\theta }⟩=\sum _{x\in \mathcal{E}}f(x){\pi }_{\theta }(x).$$

$$\sum _{i=1}^M{\nu }_{ik}{\mathcal{S}}_i⟶\sum _{i=1}^M{\nu }_{ik}^{\prime }{\mathcal{S}}_i,$$

$${\lambda }_k(x_1,\dots ,x_M)={\theta }_k\prod _{i=1}^M{1\mathrm{l}}_{\{x_i\ge {\nu }_{ik}\}}\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)(f(x+{\zeta }_k)-f(x)),$$

$$\displaystyle Q_{ij}=\left\{\begin{array}[]{cc}-\sum_{k=1}^{K}\lambda_{k}(x_{i})&\textnormal{ if }i=j\\ \sum_{k\in\mathcal{K}_{ij}}\lambda_{k}(x_{i})&\textnormal{ if }i\neq j\end{array}\right.$$

$$p_i(t)=\mathbb{P}(X(t)=x_i)$$

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

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

$$\underset{t\to \infty }{\lim }\Vert p(t)-\pi {\Vert }_{{\ell }_1}:=\sum _{i=0}^{|\mathcal{E}|}|p_i(t)-{\pi }_i|=0.$$

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

$$\mathbb{Q}V(x)\le C_1-C_{2}V(x)\text{for all}x\in \mathcal{E},$$

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

$$\sum _{k=1}^K{\lambda }_k(x)⟨v,{\zeta }_k⟩\le (C_1-C_2)-C_2⟨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⟨v,x⟩\text{for all}x\in \mathcal{E},$$

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

$${\overline{Q}}_n=Q_n+c_n{e}_{\ell }^T,$$

$$\underset{n\to \infty }{\lim }\Vert \pi -{\pi }_n{\Vert }_{{\ell }_1}=0.$$

$$-{\mathbb{Q}}_{\theta }g=f-⟨f,{\pi }_{\theta }⟩,$$

$$|f(x)|\le C(1+\Vert x{\Vert }^r)\text{for each}x\in \mathcal{E},$$

$$S_{\theta }(f)=\sum _{k=1}^K\sum _{x\in \mathcal{E}}\frac{\partial {\lambda }_k(x,\theta )}{\partial \theta }(g(x+{\zeta }_k)-g(x)){\pi }_{\theta }(x)$$

$$\Vert h{\Vert }_V=\underset{x\in \mathcal{E}}{\sup }\frac{|h(x)|}{V(x)}.$$

$$V_m(x)=1+⟨v,x{⟩}^m,$$

$$\Vert {\pi }_{\theta }{\Vert }_{V_m}:=\sum _{x\in \mathcal{E}}{V}_m(x){\pi }_{\theta }(x)<\infty .$$

$$\Psi (y,h,t)=\mathbb{E}\left(h(Y_{\theta }(t))|Y_{\theta }(0)=y\right).$$

$$\underset{t\to \infty }{\lim }\Psi (y,h,t)=⟨h,{\pi }_{\theta }⟩$$

$$|\Psi (y,h,t)-⟨h,{\pi }_{\theta }⟩|\le \beta V_m(y)e^{-ct}\text{for all}t\ge 0\text{and}y\in \mathcal{E},$$

$$S_{\theta }(f,t)=\sum _{k=1}^K\mathbb{E}\left[{\int }_0^t\frac{\partial {\lambda }_k(X_{\theta }(s),\theta )}{\partial \theta }\left(\Psi (X_{\theta }(s)+{\zeta }_k,f,t-s)-\Psi (X_{\theta }(s),f,t-s)\right)ds\right].$$

Full text

Keywords: stochastic reaction networks; the Chemical Master Equation; Finite State Projection; stationary distribution; sensitivity analysis; parameter sensitivity; Poisson equation

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

1 Introduction

Stochastic models of reaction networks are commonly used in systems and synthetic biology to model reaction dynamics within cells, where some biomolecular reactants are typically present in low copy-numbers [1]. Such models capture the random timing of reactions and allows one to investigate the role of this randomness in causing cell-to-cell variability and shaping the macroscopic properties of clonal cell-populations [2].

In a stochastic model the reaction dynamics is represented as a continuous-time Markov chain (CTMC), which keeps track of the molecular or copy-number counts of all the reacting species [3]. Therefore the CTMC evolves on a discrete state-space $\mathcal{E}\subset\mathbb{N}^{d}_{0}$, where $\mathbb{N}_{0}$ is the set of nonnegative integers and $d$ is the number of reacting species. It is known that the dynamics of the probability distribution of the CTMC is given by the Chemical Master Equation (CME) that consists of an ODE describing the inflow and outflow of probabilities at each state in the state-space $\mathcal{E}$. Therefore if this state-space is infinite in size, as is the case in many examples of interest, then the CME is nearly impossible to solve exactly. However approximate solutions to the CME over finite time-periods can be obtained using the Finite-State Projection (FSP) that projects the probability dynamics on a finite truncated state-space and solves the resulting system of ODEs [4]. Consider the situation when the CTMC is ergodic and so the CME solution converges to a unique stationary distribution as time tends to infinity [5]. In such a setting, the classical FSP is inappropriate for estimating the stationary distribution but recently a modification of this method has been proposed that is able to estimate the stationary distribution accurately under certain conditions [6]. This method is called the stationary Finite State Projection (sFSP) and it is described in detail in Section 2.3.

In many applications it is of interest to quantitatively determine the influence of some parameter $\theta$ (like temperature, extracellular ligand concentration, cell-volume etc.) on an output of the form $\mathbb{E}(f(X_{\theta}(t)))$, where $\mathbb{E}$ is the expectation operator, $(X_{\theta}(t))_{t\geq 0}$ is the $\theta$-dependent CTMC that describes the reaction dynamics, and $f$ is some real-valued function on the CTMC state-space $\mathcal{E}$. This parameter influence is often measured by estimating the infinitesimal sensitivity value

[TABLE]

Computing such sensitivities w.r.t. various parameters is useful for many applications, like investigating robustness properties of networks [7], finding critical reactions, parameter inference [8] and controlling a system’s output [9]. Suppose that the CTMC $(X_{\theta}(t))_{t\geq 0}$ is ergodic and its unique $\theta$-dependent stationary distribution is $\pi_{\theta}$. The goal of this paper is to develop a method to numerically estimate the steady-state sensitivity defined by

[TABLE]

Under certain mild conditions, this limit can be shown to exist and it can be represented as [10]

[TABLE]

where r.h.s. is the expectation of function $f$ under the stationary distribution $\pi_{\theta}$, i.e.

[TABLE]

Many methods [11, 12, 13, 14, 15, 16, 17] exist for estimating finite-time sensitivities of the form (1.1), but they rely on simulations of the dynamics $(X_{\theta}(t))_{t\geq 0}$ obtained with Gillespie’s stochastic simulation algorithm (SSA) [18] or its counterparts [19, 20]. As simulations can only be performed over finite time-periods, these methods cannot be naturally extended to estimate steady-state sensitivities (1.2). It is of course possible to approximate the steady-state sensitivity $S_{\theta}(f)$ with the finite-time sensitivity $S_{\theta}(f,T)$ for a very large time-value $T$. However this approximation is only accurate if the distribution of the random state $X_{\theta}(T)$ of the dynamics is sufficiently close to the stationary distribution $\pi_{\theta}$, and it is very difficult to determine how large $T$ needs to be for this to hold. Moreover when $T$ is large, simulations of the dynamics over time-period $[0,T]$ become computationally very expensive, and a large number of such simulations will be required to obtain a statistically useful estimate of $S_{\theta}(f,T)$. This is because the variance of finite-time sensitivity estimators generally blows-up as $T\to\infty$. For the Likelihood Ratio (LR) sensitivity estimation method [12], this variance blow-up problem can be circumvented by appropriately centering the estimator [21], but such strategies have not been found for other sensitivity estimation methods.

Motivated by the problems faced by simulation-based methods in estimating steady-state sensitivity $S_{\theta}(f)$, our aim in this paper is to present a simulation-free approach for estimating this sensitivity value. Our method relies on approximate computation of the stationary distribution $\pi_{\theta}$ with sFSP [6], and a novel mathematical result which shows how steady-state sensitivity can be evaluated from the solution of a Poisson Equation associated with the generator (see Section 3.1) of the CTMC $(X_{\theta}(t))_{t\geq 0}$. Solving the Poisson Equation is numerically challenging but we develop a Basis Function Method (BFM) that is able to efficiently obtain projection of the solution on the linear space spanned by user-specified basis functions (see Section 3.3). We demonstrate that if the collection of basis function is sufficiently large, then BFM can recover the solution of the Poisson Equation almost exactly, and this yields a very accurate estimate of the steady-state sensitivity $S_{\theta}(f)$. We refer to our method as the Poisson Estimator (PE) and we illustrate it with a number of examples. We also compare it with the simulation-based integrated centered Likelihood Ratio (IntCLR) [21] method that was recently proposed for steady-state sensitivity estimation (see Section 4).

We provide a well-documented open source C++ implementation of our method (see Section 3.4), which includes parallelized solvers for the stationary distribution (with sFSP) and for the Poisson Equation (with BFM). We must point out that apart from estimating sensitivities, solving the Poisson Equation for Markov chains has many other applications, for example in computing the optimal policy function for Markov Decision Processes [22, 23].

2 Preliminaries

2.1 The Stochastic Reaction Network model

Consider a reaction network involving $M$ biochemical species, denoted by $\mathcal{S}_{1},\dots,\mathcal{S}_{M}$. These species interact via $K$ reactions, and each reaction $k$ has the form

[TABLE]

where $\nu_{ik}$ $(\nu^{\prime}_{ik})$ is a nonnegative integer specifying the number of $\mathcal{S}_{i}$ molecules consumed (produced) by reaction $k$. Hence the change in the number of $\mathcal{S}_{i}$ molecules, caused by the firing of reaction $k$ is given by the integer $\zeta_{ik}=(\nu^{\prime}_{ik}-\nu_{ik})$ and $\zeta_{k}=(\zeta_{1k},\dots,\zeta_{Mk})$ is called the stoichiometry vector for reaction $k$. To each reaction $k$ we associate a propensity function $\lambda_{k}:\mathbb{N}^{M}_{0}\to[0,\infty)$, that specifies the firing rate of this reaction as $\lambda_{k}(x_{1},\dots,x_{M})$, when the molecular count or copy-number of species $\mathcal{S}_{i}$ is $x_{i}$. Commonly this rate function is given by mass-action kinetics [3]

[TABLE]

where $\theta_{k}$ is a positive rate constant.

In the classical continuous-time Markov chain (CTMC) model of a reaction network [3], the state at time $t$ is simply the vector of species copy-numbers $X(t)=(X_{1}(t),\dots,X_{M}(t))\in\mathbb{N}^{M}_{0}$ at time $t$. When the state is $X(t)=x$, the rate of firing of the $k$-th reaction is $\lambda_{k}(x)$ and if it fires before any other reaction, the state moves to $(x+\zeta_{k})$. In other words, the generator111A Markov process can be uniquely characterized by its generator, which is an operator specifying the rate of change of the distribution of the Markov process. of the CTMC $(X(t))_{t\geq 0}$ is given by

[TABLE]

where $f$ is any bounded real-valued function on the state-space $\mathcal{E}\subset\mathbb{N}^{M}_{0}$ of the CTMC. We shall assume that this state-space $\mathcal{E}$ is non-empty and closed under the CTMC dynamics, i.e. if $x$ is a state in $\mathcal{E}$ and $k$ is some reaction with $\lambda_{k}(x)>0$ then $(x+\zeta_{k})$ is also a state in $\mathcal{E}$.

As the state-space $\mathcal{E}$ is countable, we can construct a bijection $\phi:\mathcal{E}\longrightarrow\{0,1\dots,|\mathcal{E}|\}$, where $|\mathcal{E}|$ denotes the number of elements in $\mathcal{E}$. Letting $x_{i}=\phi^{-1}(i)$ for $i=0,1,\dots$, we can represent the state-space as $\mathcal{E}=\{x_{0},x_{1},\dots\}$, and express the operator $\mathbb{Q}$ (2.7) as the transition rate matrix $Q=[Q_{ij}]$

[TABLE]

where $\mathcal{K}_{ij}=\{k=1\dots,K:x_{j}=x_{i}+\zeta_{k}\}$ is the set of reactions that can take the state from $x_{i}$ to $x_{j}$ in a single firing. Note that if $\mathcal{K}_{ij}=\emptyset$ then $Q_{ij}=0$. Moreover the matrix $Q$ is bi-infinite when $|\mathcal{E}|=\infty$, which is often the case in examples of interest.

Suppose $(X(t))_{t\geq 0}$ is the CTMC with transition rate matrix $Q$ and some initial state $X(0)\in\mathcal{E}$. Let $p(t)=(p_{0}(t),p_{1}(t),\dots)$ where

[TABLE]

is the probability that the dynamics is in state $x_{i}$ at time $t$. The time-evolution of the probability vector is given by the well-known Chemical Master Equation (CME) which can be expressed as

[TABLE]

Even though the CME is a first-order linear system of ODEs, solving it analytically is infeasible when the state-space $\mathcal{E}$ is large or infinite in size. For this reason the Finite State Projection (FSP) [4] was developed to approximately solve the CME by projecting it on a truncated state-space $\mathcal{E}_{n}\subset\mathcal{E}$.

2.2 Ergodicity of the CTMC

Often one is interested in the steady-state behavior of the CTMC $(X(t))_{t\geq 0}$, which is characterized by a stationary distribution $\pi=(\pi_{0},\pi_{1},\dots)$ which acts as a fixed-point for the CME (2.11),

[TABLE]

The CTMC is called ergodic if the fixed point is unique and globally attracting in the sense that

[TABLE]

In this paper we work under the assumption of exponential ergodicity [24] which holds if the above convergence is exponentially fast, i.e. there exists a constant $\rho>0$ and another constant $C>0$, that may depend on the initial distribution $p(0)$, such that for any $t>0$

[TABLE]

Verifying ergodicity of the CTMC corresponding to a reaction network is a challenging problem when the state-space is infinite. One way to address this challenge is to first establish that the underlying state-space $\mathcal{E}$ is irreducible, i.e. all states in $\mathcal{E}$ are accessible from each other via a sequence of positive-propensity reactions, and then construct a Foster-Lyapunov function on the state-space to show that the dynamics has a tendency to be attracted to a compact set within the state-space [24]. In particular, using Theorem 7.1 in [24] we can prove exponential ergodicity of the CTMC by showing that there exists a norm-like function222A positive function is called norm-like if all its sub-level sets are compact. $V:\mathcal{E}\to[1,\infty)$ such that for some $C_{1},C_{2}>0$, we have

[TABLE]

where $\mathbb{Q}$ (2.7) is the generator of the CTMC. In [5] and [25] computational procedures have been developed that systematically check state-space irreducibility and construct the required Foster-Lyapunov function $V$ satisfying (2.14) under some conditions. In fact, for many systems biology networks, a linear Foster-Lyapunov function of the form

[TABLE]

can be constructed [5]. Here $\langle\cdot,\cdot\rangle$ denotes the standard inner product on $\mathbb{R}^{d}$ and $v\in\mathbb{R}^{d}$ is a component-wise positive vector that can be found via linear programming. Note that for this linear Foster-Lyapunov function, condition (2.14) can be equivalently expressed as

[TABLE]

Often for typical reaction networks, in addition to this condition we also have that for some $C_{3},C_{4}>0$

[TABLE]

which ensures that all statistical moments of the stationary distribution $\pi$ are finite (see Theorem 5 in [5]).

2.3 The Stationary Finite State Projection Algorithm

The classical FSP algorithm [4] introduces an absorbing state into the dynamics, in order to restrict it to a finite truncated subset. Therefore even though FSP is very successful in solving the CME over finite time-intervals, it cannot be used to estimate the stationary distribution $\pi$, as all the probability-mass will flow into the absorbing state at steady-state (see [6] for more details).

Motivated by this problem, recently the stationary Finite State Projection (sFSP) algorithm [6] was proposed, that instead of using an absorbing state to capture all the outgoing transitions from the truncated state-space, forces all the outgoing transitions to flow into a designated state within the truncated state-space (see Figure 1). More formally, suppose $\mathcal{E}_{n}=\{x_{j_{1}},\dots,x_{j_{n}}\}\subset\mathcal{E}$ is the truncated state-space of size $n$, and let $Q_{n}$ be the matrix formed by considering only the rows and columns of $Q$ that are indexed by $j_{1},\dots,j_{n}$. Furthermore, let $c_{n}$ be the $n$-dimensional nonnegative column vector whose $i$-th component is the sum of the rates of all the reactions that originate at state $x_{i}$ and terminate at a state outside $\mathcal{E}_{n}$,

[TABLE]

Assuming that the $\ell$-th state $x_{j_{\ell}}$ is the designated state, in the sFSP algorithm the CTMC dynamics on the truncated state-space $\mathcal{E}_{n}$ has the transition rate matrix given by

[TABLE]

where $e_{\ell}$ is the $n\times 1$ vector whose $\ell$-th component is $1$ and the rest are zeros. The stationary distribution $\pi$ is then estimated by computing the finite-dimensional stationary distribution $\pi_{n}$ on $\mathcal{E}_{n}$ for the transition rate matrix $\overline{Q}_{n}$. This can be done by solving the linear system (2.12) with $Q$ and $\pi$ replaced by $\overline{Q}_{n}$ and $\pi_{n}$ respectively.

We now make some assumptions which are required to certify the accuracy of sFSP.

Assumption 2.1

The state-space $\mathcal{E}$ is irreducible for the original CTMC with transition rate matrix $Q$.

There exists a Foster-Lyapunov function $V:\mathcal{E}\to[1,\infty)$ satisfying (2.14).

$\{\mathcal{E}_{n}:n=1,2,\dots\}$* is a sequence of increasing finite sets (i.e. $n\leq m$ implies $\mathcal{E}_{n}\subset\mathcal{E}_{m}$) which converge to the full state-space $\mathcal{E}$ in the limit $n\to\infty$.*

Under these assumptions it can be shown (see Theorem 3.1 in [6]) that the stationary distribution $\pi_{n}$ for the projected CTMC on truncated state-space $\mathcal{E}_{n}$ exists uniquely and it converges to the stationary distribution $\pi$ for the original CTMC in the $\ell_{1}$ metric

[TABLE]

This shows that for a large enough truncated state-space $\mathcal{E}_{n}$, the sFSP-estimated stationary distribution $\pi_{n}$ is likely to be “close” to the true stationary distribution $\pi$. More precise estimates on the approximation error can be found in [6] and it involves the outflow rate, measured as $\langle c_{n},\pi_{n}\rangle$, at the estimated stationary distribution.

3 Steady-State Sensitivity Estimation

Consider the stochastic model of a reaction network, as described in Section 2.1, and suppose that in addition to the state $x$ the propensity functions $\lambda_{k}$-s depend on a scalar parameter $\theta$ and the mapping $\theta\mapsto\lambda_{k}(x,\theta)$ is differentiable for each fixed $x$. Replacing each $\lambda_{k}(x)$ by $\lambda_{k}(x,\theta)$ in (2.7), we obtain the generator $\mathbb{Q}_{\theta}$ of the $\theta$-dependent CTMC $(X_{\theta}(t))_{t\geq 0}$.

In this section, we present the main contributions of this paper, which is developing a method to estimate the steady-state sensitivity $S_{\theta}(f)$ defined by the limit (1.2). In the next section we prove that this limit exists under certain assumptions, and connect it with the solution of a Poisson Equation for the generator $\mathbb{Q}_{\theta}$.

3.1 An explicit expression for steady-state sensitivity

Suppose parts (A) and (B) of Assumption 2.1 hold for the original CTMC with generator $\mathbb{Q}_{\theta}$, and hence the state-space $\mathcal{E}$ is irreducible and suppose $V$ is the Foster-Lyapunov function satisfying (2.14). Under these assumptions this CTMC is exponentially ergodic with a unique stationary distribution $\pi_{\theta}$. Let $f:\mathcal{E}\to\mathbb{R}$ be the output function for which the steady-state sensitivity $S_{\theta}(f)$ needs to be computed. The Poisson equation for the pair $(\mathbb{Q}_{\theta},f)$ is given by

[TABLE]

where the inner product $\langle f,\pi_{\theta}\rangle$ is defined as in (1.4). It is known that the solution $g$ of this Poisson equation, if it exists, is unique up to addition by a constant function.

We next present our main result in this paper which connects the steady-state sensitivity $S_{\theta}(f)$ to the solution of the Poisson equation (3.20). Henceforth, we call a function $f:\mathcal{E}\to\mathbb{R}$ polynomially growing, if there exist constants $C,r>0$ such that

[TABLE]

where $\|\cdot\|$ denotes the standard norm on $\mathbb{R}^{d}$.

Theorem 3.1

Suppose parts (A) and (B) of Assumption 2.1 hold for the CTMC with generator $\mathbb{Q}_{\theta}$, and with a linear Foster-Lyapunov function $V$ (2.15) satisfying (2.16) and (2.17). Furthermore assume that the output function $f$ and $\theta$-derivatives of the propensity functions (i.e. $\partial\lambda_{k}/\partial\theta$ for $k=1,\dots,K$) are polynomially growing. Then there exists a solution $g$ of the Poisson equation (3.20) and the steady-state sensitivity is given by

[TABLE]

where $\pi_{\theta}$ is the stationary distribution of the CTMC.

Remark 3.2

Note that the sensitivity value $S_{\theta}(f)$ is independent of the choice of the solution $g$ of the Poisson equation (3.20), as all such solutions differ by a constant function and hence the differences $g(x+\zeta_{k})-g(x)$ are the same for all the solutions.

Proof. Let $V$ be the linear Foster-Lyapunov function (2.15) satisfying (2.16) and (2.17). Corresponding to $V$ we can define the $V$-norm of any function $h:\mathcal{E}\to\mathbb{R}$ as

[TABLE]

It is known that if function $f$ has finite $V$-norm (i.e. $\|f\|_{V}<\infty$), then there exists a solution $g$ of the Poisson equation (3.20) which also has finite $V$-norm (see Theorem 2.3333This result has been proved for discrete-time Markov chains but it can be easily transformed to the continuous-time setting using the resolvent operator as mentioned in the proof of Proposition 3.3. in [38]).

In our case the output function $f$ is polynomially growing and so $\|f\|_{V}$ may not be finite. However for any positive integer $m$, if we define a function $V_{m}:\mathcal{E}\to[1,\infty)$ by

[TABLE]

then it can be seen from the proof of Theorem 5 in [5] that, under conditions (2.16) and (2.17), $V_{m}$ also serves as a Foster-Lyapunov function for the dynamics (i.e. it satisfies (2.14) for generator $\mathbb{Q}_{\theta}$). Moreover from Theorem 4.2 in [24] it can be concluded that the stationary distribution $\pi_{\theta}$ is such that

[TABLE]

As our function $f$ is polynomially growing, we have $\|f\|_{V_{m}}<\infty$ for some positive integer $m$, and hence there is a solution $g$ of the Poisson equation (3.20) which also satisfies $\|g\|_{V_{m}}<\infty$. Similarly as $\theta$-derivatives of the propensity functions are polynomially growing as well, they also have finite $V_{m}$-norm for some $m$. As (3.24) holds for any $m$, it can be seen that the r.h.s. of (3.21) is finite even though the state-space $\mathcal{E}$ may be countably infinite.

We now prove that the steady-state sensitivity $S_{\theta}(f)$, defined by (1.2), is equal to the r.h.s. of (3.21). Let $(Y_{\theta}(t))_{t\geq 0}$ be a CTMC with generator $\mathbb{Q}_{\theta}$ and for any state $y\in\mathcal{E}$ and polynomially growing function $h:\mathcal{E}\to\mathbb{R}$ let

[TABLE]

Under conditions (2.16) and (2.17), this expectation is finite. Furthermore the limit

[TABLE]

holds (see Proposition S2.2 in [5]), where the inner product between function $h$ and stationary distribution $\pi_{\theta}$ is defined as in (1.4). Moreover since $h$ is polynomially growing, there exist constants $c,\beta>0$ and positive integer $m$ such that

[TABLE]

where $V_{m}$ is the Foster-Lyapunov function given by (3.23). This relation follows from Theorem 6.1 in [24].

Let $S_{\theta}(f,t)$ be the finite-time sensitivity value defined by (1.1) for the CTMC $(X_{\theta}(t))_{t\geq 0}$ with generator $\mathbb{Q}_{\theta}$ and initial state $x_{0}$. From Theorem 3.3 in [39], we can express $S_{\theta}(f,t)$ as

[TABLE]

Using (3.25) we know that for some constants $c,\beta>0$ and $m$

[TABLE]

This inequality implies that the integrand in (3.26) can be bounded by the product of a polynomially growing function of the state $X_{\theta}(s)$ and an exponential decay function. Since the moments of the CTMC $(X_{\theta}(t))_{t\geq 0}$ can be uniformly bounded over the infinite time-interval $[0,\infty)$ (see Theorem 2 in [5]), the integral on the r.h.s. of (3.1) exists as $t\to\infty$ and so we can write the steady-state sensitivity as

[TABLE]

Notice that due to (3.1) and uniform moment boundedness, for any $u\geq 0$ we have

[TABLE]

and therefore by exploiting the linearity of the integration and expectation operators we can represent the steady-state sensitivity as

[TABLE]

Letting $p_{\theta}(s,x)=\mathbb{P}(X_{\theta}(s)=x)$ we obtain

[TABLE]

where the last relation follows because (3.25) holds for any $m$, and $(t-s)\geq u$ for $s\in[0,t-u]$ and hence if $u$ is large $p_{\theta}(t-s,x)$ is close to the stationary distribution $\pi_{\theta}$ in the norm $\|\cdot\|_{V_{m}}$ defined by (3.24).

Now to prove this theorem it suffices to show that

[TABLE]

where $g$ is a solution of the Poisson equation (3.20). Let $(Y_{1}(t))_{t\geq 0}$ and $(Y_{2}(t))_{t\geq 0}$ be two CTMCs with generator $\mathbb{Q}_{\theta}$ and initial states $x$ and $(x+\zeta_{k})$ respectively. Then applying Dynkin’s formula (see [40]) on function $g$ we obtain

[TABLE]

and

[TABLE]

Subtracting (3.1) from (3.1) and rearranging we get

[TABLE]

Due to ergodicity $\lim_{t\to\infty}\mathbb{E}(g(Y_{2}(t)))=\lim_{t\to\infty}\mathbb{E}(g(Y_{1}(t)))=\langle g,\pi_{\theta}\rangle$, and hence letting $t\to\infty$ in (3.32) proves (3.29) and completes the proof of this theorem. $\Box$

3.2 Approximation of steady-state sensitivity

We mentioned before, that in most examples of interest the state-space $\mathcal{E}$ is infinite in size. Hence the stationary distribution $\pi_{\theta}$ cannot be exactly computed, and the problem of solving the Poisson equation (3.20) becomes infinite-dimensional. Due to these issues the exact formula (3.21) cannot be directly used to estimate the steady-state sensitivity value, and we need to replace the stationary distribution $\pi_{\theta}$ and the Poisson equation solution $g$ with its approximations. Let $\mathcal{E}_{n}=\{x_{j_{1}},\dots,x_{j_{n}}\}$ be the truncated state-space as in part (C) of Assumption 2.1, and let $\pi_{n,\theta}$ be the sFSP estimated approximation (see Section 2.3) of $\pi_{\theta}$ with $x_{j_{\ell}}$ as the designated state. We can view $\pi_{n,\theta}$ as a column-vector which sums to one and is in the left null-space of the matrix $\overline{Q}_{n,\theta}$ defined by (2.18).

We define the finite-dimensional Poisson equation, analogous to (3.20) as

[TABLE]

where ${\bf 1}$ is the $n\times 1$ vector of all ones, and $f_{n}$ denotes the $n\times 1$ vector given by $f_{n}=(f(x_{j_{1}}),\dots,f(x_{j_{n}}))$. A solution $g_{n}$ of this Poisson equation exists because $\pi_{n,\theta}$ is orthogonal to the r.h.s. of (3.33) and it spans the one-dimensional left null-space of matrix $\overline{Q}_{n,\theta}$. This null-space is one-dimensional because the CTMC specified by transition-rate matrix $\overline{Q}_{n,\theta}$ has a unique stationary distribution over the state-space $\mathcal{E}_{n}$. This uniqueness follows from the fact that the construction of $\overline{Q}_{n,\theta}$ ensures that the truncated state-space $\mathcal{E}_{n}$ is irreducible for the corresponding CTMC dynamics (for details see Section III.A in [6]). As for the original Poisson equation (3.20), a solution $g_{n}$ of the finite-dimensional Poisson equation (3.33), is unique up to addition by a constant vector. Using a solution $g_{n}$ we define the approximate steady-state sensitivity as

[TABLE]

where we set $(x+\zeta_{k})$ to the designated state $x_{j_{\ell}}$ whenever state $x$ is in the boundary of $\mathcal{E}_{n}$ denoted by

[TABLE]

Note that in this formula we regard $g_{n}$ and $\pi_{n,\theta}$ as real-valued functions over $\mathcal{E}_{n}$. The next proposition shows that for a large enough truncated state-space $\mathcal{E}_{n}$, $S_{n,\theta}(f)$ is an accurate approximation of the true steady-state sensitivity value $S_{\theta}(f)$.

Proposition 3.3

Suppose part (C) of Assumption 2.1 holds along with all the conditions of Theorem 3.1. Furthermore assume that the sequence of truncated state-spaces $\{\mathcal{E}_{n}\}$ grows uniformly w.r.t. to the Foster-Lyapunov function $V$ in the sense that there exists a constant $\theta\in(0,1)$ such that

[TABLE]

Then letting $S_{n,\theta}(f)$ be defined by (3.34), we have

[TABLE]

Proof. In this proof we drop the subscript $\theta$ for convenience. Let $Q$ be the bi-infinite transition rate matrix (see (2.10)) of the original CTMC with stationary distribution $\pi$ over state-space $\mathcal{E}=\{x_{0},x_{1},\dots\}$. Similarly let $\overline{Q}_{n}$ be the finite transition rate matrix (see (2.18)) of the CTMC projected over the truncated state-space $\mathcal{E}_{n}$ which has stationary distribution $\pi_{n}$. Without loss of generality we can assume that the truncated state-space comprises the first $n$ states in $\mathcal{E}$, i.e. $\mathcal{E}=\{x_{0},\dots,x_{n-1}\}$, and the designated state is $x_{\ell}=x_{0}$. Henceforth let $g$ and $g_{n}$ be the unique solutions of the Poisson equations (3.20) and (3.33) respectively, with $\langle g,\pi\rangle=\langle g_{n},\pi_{n}\rangle=0$ (see Remark 3.2).

As explained in the proof of Theorem 3.1, for any positive integer $m$, the function $V_{m}$ (see (3.23)) serves as a Foster-Lyapunov function for the original CTMC dynamics. The same holds true for the projected CTMC dynamics with transition rate matrix $\overline{Q}_{n}$, and the constants $C_{1}$ and $C_{2}$ that appear in the corresponding Foster-Lyapunov condition (2.14) are independent of $n$ (see [41]). Appealing to this fact and Theorem 2.3 in [38], we can conclude that the family of solutions $\{g_{n}\}$ to the Poisson equation (3.33) is uniformly polynomially growing, i.e. there exist positive constants $C$ and $m$ such that for each $n$ and $x\in\mathcal{E}$

[TABLE]

Moreover from Theorem 3.1 in [6] we know that for any $m$

[TABLE]

where the norm $\|\cdot\|_{V_{m}}$ is as defined in (3.24). This convergence along with (3.36) and the dominated convergence theorem, implies that in order to prove (3.35) it suffices to show point-wise convergence of $g_{n}$ to $g$ as $n\to\infty$, i.e.

[TABLE]

We shall establish this convergence with the help of the $\beta$-resolvent matrix (see [41]), defined as

[TABLE]

where $\beta>0$ and ${\bf I}$ is the identity matrix. It can be seen that $R_{\beta}$ is a positive matrix satisfying $R_{\beta}{\bf 1}={\bf 1}$, $\pi^{T}R_{\beta}=\pi^{T}$ and

[TABLE]

We can view $R_{\beta}$ as the probability transition matrix of a discrete time Markov chain over $\mathcal{E}$. The resolvent $\overline{R}_{\beta,n}$ corresponding to the projected CTMC with transition rate matrix $\overline{Q}_{n}$ and state-space $\mathcal{E}_{n}$, can be defined in a similar way.

Now select a state $x=x_{i}\in\mathcal{E}$ and let $\mu_{x}$ be the vector whose $i$-th coordinate is $1$ and the rest are all zeros. From the analysis in [6] we can conclude that for any positive integer $m$, there exist constants $C>0$ and $\rho\in(0,1)$ such that

[TABLE]

where $A^{j}$ denotes the $j$-th power of matrix $A$. Viewing functions $f$ and $g$ as vectors, we can rewrite Poisson equation (3.20) as

[TABLE]

Multiplying this equation by $R_{\beta}$ on the left and using (3.39) we obtain

[TABLE]

which allows us to express $g(x)$ as

[TABLE]

This expression is well-defined due to (3.40) and the fact that $f$ is polynomially growing. Moreover one can verify that $\langle g,\pi\rangle=0$. Similar to (3.41) we can represent the solution $g_{n}$ of the Poisson equation (3.33) as

[TABLE]

Since (3.37) holds, the point-wise convergence of $f_{n}$ to $f$ as $n\to\infty$ is immediate and the convergence of the resolvent operator $\overline{R}_{\beta,n}$ to $R_{\beta}$ is proved in [41]. Also notice that for some positive integer $m$, we have $\|f\|_{V_{m}}<\infty$ and $\sup_{n}\|f_{n}\|_{V_{m}}<\infty$, where the norm $\|\cdot\|_{V_{m}}$ is defined as in (3.22). Hence using (3.40) and the dominated convergence theorem we can conclude that (3.38) holds for each state $x\in\mathcal{E}$ and this concludes the proof of this proposition. $\Box$

3.3 Solving the Poisson Equation (3.33)

We now discuss how one can efficiently compute the steady-state sensitivity estimate $S_{n,\theta}(f)$ (3.34). We can obtain the approximation $\pi_{n,\theta}$ of the stationary distribution via sFSP [6]. Thereafter to compute $S_{n,\theta}(f)$ we need to find a solution to the finite-dimensional Poisson equation (3.33). Note that this equation has a unique solution $g_{n}$, if we impose that $\langle g_{n},{\bf 1}\rangle=0$. We now focus on how this unique solution can be computed. For convenience, we drop the subscript $\theta$ from $\overline{Q}_{n,\theta}$ and denote this matrix as $\overline{Q}_{n}$ in this section. Observe that matrix $\overline{Q}_{n}$ is singular (as $\textnormal{Rank}(\overline{Q}_{n})=n-1$) and usually very large in size. This poses difficulties in solving (3.33) with standard large-scale linear solvers that are currently in use (see Section 3.4). To overcome these issues, we solve the Poisson equation (3.33) using the Basis Function Method (BFM) which we describe next. This method is approximate in nature, but it is computationally very feasible and it often yields a very accurate estimate of the solution of the Poisson equation.

The idea of the BFM is very simple. We start by choosing a basis of $m$ vectors (with $m\ll n$) $\mathcal{B}=\{g_{1},\dots,g_{m}\}$ in $\mathbb{R}^{n}$. Then we construct the image of this basis under matrix $\overline{Q}_{n}$, by computing

[TABLE]

Letting $\overline{Q}_{n}(\mathcal{B}):=\{f_{1},\dots,f_{m}\}$, our goal is to project the r.h.s. of (3.33), given by

[TABLE]

on the linear space spanned by the vectors in $\overline{Q}_{n}(\mathcal{B})$. If this projection is given by

[TABLE]

for some coefficient vector $c=(c_{1},\dots,c_{m})\in\mathbb{R}^{m}$, then our solution approximation is given by

[TABLE]

To obtain the best solution approximation, we construct the optimal projection $\widehat{f}(c)$ by choosing the coefficient vector $c$ that minimizes the following error function

[TABLE]

where $\|\cdot\|_{\ell_{2}}$ denotes the standard $\ell_{2}$-norm and $\odot$ is the Hadamard (element-wise) vector product. In this error function, the discrepancy between $\overline{f}_{n}$ and its projection $\widehat{f}(c)$ is weighted by the approximate stationary distribution $\pi_{n,\theta}$. This ensures that our projection is biased towards minimizing discrepancies in states with high stationary probabilities, which is important because contribution from these states is more dominant in the expression (3.34) for the steady-state sensitivity $S_{n,\theta}(f)$.

Let ${\bf F}=\textnormal{Col}(f_{1},\dots,f_{m})$ be the $n\times m$ matrix whose columns are vectors in $\overline{Q}_{n}(\mathcal{B})$. Also let ${\bf\Pi}$ be the $n\times n$ diagonal matrix whose entries are given by the approximate stationary distribution $\pi_{n,\theta}$. The minimization of the error function $\epsilon(c)$ is essentially a weighted least-square problem whose optimal solution $c^{*}$ satisfies

[TABLE]

where ${\bf A}$ is the $m\times m$ matrix and ${\bf b}$ is the $m\times 1$ vector given by

[TABLE]

respectively. Since $m\ll n$, the $m$-dimensional linear-system (3.47) is much easier to solve than the Poisson equation (3.33). If the optimal error $\epsilon(c^{*})$ is below a small tolerance level, we can expect $\widehat{g}(c^{*})$ (see (3.45)) to be an accurate solution of the Poisson equation (3.33). On the other hand, if the error $\epsilon(c^{*})$ is above the tolerance level, it can be reduced by expanding the set of basis functions $\mathcal{B}$ by including more vectors.

Observe that the $n\times n$ matrix $\overline{Q}_{n}$ is extremely sparse, as each row can have up to $(K+1)$ non-zero entries, where $K$ is the total number of reactions. The BFM essentially substitutes the problem of solving the Poisson equation (3.33) for a large sparse matrix, with the problem of solving a much smaller (but dense) linear system (3.47). For this approach to work well, it is important that the set of basis vectors $\mathcal{B}=\{g_{1},\dots,g_{m}\}$ is chosen in such a way that the linear span of the image set $\overline{Q}_{n}(\mathcal{B})=\{f_{1},\dots,f_{m}\}$, contains the vector $\overline{f}_{n}$ (3.43). While it is difficult to determine the basis vectors that meet this criterion, we find that for most biological examples of interest, this containment condition is fulfilled with a very small approximation error when we choose these basis vectors as representatives of monomials of the form

[TABLE]

where exponents $e_{1},\dots,e_{m}$ are non-negative integers. Given such a monomial, its representative on the state-space $\mathcal{E}_{n}=\{x_{j_{1}},\dots,x_{j_{n}}\}$ is the $n$-dimensional vector $g=(G(x_{j_{1}}),\dots,G(x_{j_{n}}))$. Hence for any positive integer $d_{\textnormal{max}}$ specifying the maximal monomial degree, let $\mathcal{B}(d_{\textnormal{max}})$ be the set of all basis vectors representing monomials with exponents $e_{1},\dots,e_{M}$ satisfying $\sum_{i=1}^{M}e_{i}\leq d_{\textnormal{max}}$. Letting $\binom{a}{b}:=\frac{a!}{b!(a-b)!}$ denote the binomial coefficient, the number of basis vectors in $\mathcal{B}(d_{\textnormal{max}})$ is equal to

[TABLE]

which grows rapidly as the maximal degree $d_{\textnormal{max}}$ and the number of species $M$ increases.

We refer to our steady-state sensitivity estimation method as the Poisson Estimator (PE) as it is based on solving the Poisson Equation (3.33) associated with the stochastic reaction network. The main steps required by our method are summarized in Algorithm 1.

3.4 Software Implementation

We have computationally implemented our method (PE) as a C++ library, called cossmosLib, to which one can specify an arbitrary reaction network (as a script in Systems Biology Markup Language (SMBL) [26]) and output estimates of the desired steady-state parametric sensitivities. In our software implementation, the most computationally intensive tasks are parallelized via support of the Message Passing Interface (MPI). The cossmosLib can be downloaded from our GitLab repository [27]. This page contains full documentation of the library along with detailed installation instructions.

Our library has several dependencies on (mostly) open source software. The only non-open source software we employ, are the Intel MPI and Intel MKL (Math Kernel Library) libraries. Note that these libraries could in theory be interchanged with corresponding open source libraries (some implementations of the BLAS and LAPACK libraries as well as some implementation, like Open MPI, of the MPI library), but in our experience the Intel libraries performed the best. Furthermore we employ the ParMETIS library [28] from the Karypis Lab, the SuperLU_dist library [29] from the Computational Research Division of the Berkeley Lab Computing Sciences, the libSBML library [30] from the SBML community and the Trilinos library [31, 32, 33] from the Sandia National Laboratories.

4 Examples

In this section we test our method for steady-state sensitivity estimation on several examples, and compare it with the simulation method IntCLR [21] for estimating the steady-state sensitivity. In the preceding sections we discussed sensitivity computation w.r.t. a scalar parameter $\theta$, but the extension to a vector-valued parameter is straightforward. For vector-valued $\theta=(\theta_{1},\theta_{2},\dots)$, we view the steady-state sensitivity value $S_{\theta}(f)$ as the vector

[TABLE]

where each $S_{\theta_{i}}(f)$ is defined as (1.2). Suppose we obtain an estimate $\widehat{S}_{\theta_{i}}(f)$ of the true sensitivity value $S_{\theta_{i}}(f)$. Then the percentage Relative Error (RE) is defined as

[TABLE]

Moreover we define the percentage Total Relative Error (TRE) as

[TABLE]

From now on we denote the steady-state expectation of any output function $f$ of our $\theta$-dependent stochastic reaction network as $\mathbb{E}_{\theta}(f)$. In all our examples $f$ takes the form $f(x)=x_{i}$, i.e. the output is simply the copy-number of species $\mathcal{S}_{i}$. Therefore we shall express the steady-state expectation as $\mathbb{E}_{\theta}(x_{i})$ and the corresponding $\theta$-sensitivity vector as $S_{\theta}(x_{i})=(S_{\theta_{1}}(x_{i}),S_{\theta_{2}}(x_{i}),\dots)$.

In our first three examples, the irreducible state-space is $\mathbb{N}^{d}_{0}$, where $d$ is the number of species. For our computations, we shall employ trapezoidal truncations of the state-space (see [6]) which are defined using two cut-off values $C_{l}$ and $C_{r}$ (with $C_{l}\leq C_{r}$) as

[TABLE]

Here we follow the convention that $\binom{a}{b}=0$ if $a<b$, and the function $\phi_{d}$ is given by

[TABLE]

In the last example we consider, the irreducible state-space is a finite set and hence no truncations are necessary for the computations.

Finally we would like to point out that all the computations were performed on the Euler cluster of ETH Zurich[34]. This allowed us to test the parallel implementation of our library cossmosLib.

4.1 Gene Expression Network

Our first example is the standard gene-expression network introduced in [35] (see Figure 2(A)). It consists of two species $\mathcal{S}_{1}$ (mRNA) and $\mathcal{S}_{2}$ (Protein) which participate in four reactions whose propensities are parametrized by $\theta=(\theta_{1},\theta_{2},\theta_{3},\theta_{4})$. These reactions are described in Table A.1 and we set the parameters as $\theta_{1}=90$, $\theta_{2}=4$, $\theta_{3}=0.5$ and $\theta_{4}=0.2$.

We are interested in estimating the sensitivity-vectors $S_{\theta}(x_{1})$ and $S_{\theta}(x_{2})$, which correspond to parametric sensitivities of the steady-state expectation of the mRNA copy-numbers ($\mathbb{E}_{\theta}(x_{1})$) and the protein copy-numbers ($\mathbb{E}_{\theta}(x_{2})$) respectively. In this example all propensity functions are linear, and hence the steady-state expectations $\mathbb{E}_{\theta}(x_{1})$ and $\mathbb{E}_{\theta}(x_{2})$ can be explicitly computed as functions of $\theta$. Consequently the sensitivity vectors can also be exactly evaluated and used for estimating the RE.

We next estimate the sensitivity-vectors with our PE method given by Algorithm 1. We employ our method with five trapezoidal state-space truncations mentioned in Table A.2. This table also mentions the designated state chosen for sFSP. For estimating the solution of the finite-dimensional Poisson equation (3.33) via BFM we use $65$ basis vectors, which corresponds to $d_{\textnormal{max}}=10$ (see Section 3.3). The sensitivity-vectors estimated by PE, with the largest truncated state-space, are reported in Table A.3 and depicted as a bar chart in Figure 2(B) along with their percentage Relative Error (RE). Note that the RE$\%$ is always below $3\cdot 10^{-7}$ which demonstrates the accuracy of PE method. In Figure 2(C) we illustrate how the Total Relative Error (TRE) decreases, while the memory requirements get larger as the truncated state-space expands in size. Even though the memory requirements increase linearly with the state-space size, the decrease in the TRE$\%$ is almost exponential, which is consistent with the results reported in [6]. In Figure 2(D) we highlight potential drawbacks of the simulation-based approach IntCLR in estimating steady-state sensitivity $S_{\theta_{1}}(x_{2})$. Here the total central processing time is fixed to be $120$ hours, and the final time $T$ denotes the time at which steady-state is assumed to be reached and at which the sensitivity is computed. One can see that as $T$ increases, the estimation accuracy improves (i.e. RE decreases) but the number of simulation samples that are generated within 120 hours gets lower and hence the standard deviation of the estimator increases. Such a trade-off is natural to expect for any simulation-based method for estimating steady-state sensitivity. Note that even though our PE method has high memory requirements, even with the largest truncated state-space we consider, PE yields all the estimates within just $19$ minutes (with four processors on ETH Zurich’s Euler cluster) which is much smaller than the processing time for IntCLR (120 hours). Moreover even the most accurate estimate obtained by IntCLR has a RE$\%$ 10,000 times higher than the RE$\%$ for the corresponding PE-estimated value.

4.2 Toggle-Switch Network

Our next example is the synthetic toggle-switch proposed by Gardner et al.[36]. Here two species $\mathcal{S}_{1}$ and $\mathcal{S}_{2}$ are simply repressing the production of each other. This network consists of four reactions whose propensities are parametrized by $\theta=(\theta_{1},\dots,\theta_{6})$ (see Table A.4). We set the parameters as $\theta_{1}=500$, $\theta_{2}=3$, $\theta_{3}=200$, $\theta_{4}=0.4$, $\theta_{5}=1.5$ and $\theta_{6}=1$. As in the previous example, we are interested in estimating the sensitivity-vectors $S_{\theta}(x_{1})$ and $S_{\theta}(x_{2})$, which correspond to parametric sensitivities of the steady-state expectation of the copy-numbers of species $\mathcal{S}_{1}$ and $\mathcal{S}_{2}$ respectively. It is well-known that the stationary distribution for the CTMC model of the toggle-switch is bimodal, with each mode corresponding to one of the species being dominant [6]. This bimodality causes problems in estimating the stationary distribution and it arises due to nonlinearities in the propensity functions. Moreover these nonlinearities prevent us from determining the steady-state sensitivities exactly. Hence we cannot compute the RE$\%$ for this example.

We now estimate the sensitivity-vectors with our PE method (Algorithm 1). We employ six trapezoidal state-space truncations (see Table A.5), with the designated state set as $x_{\ell}=(235,115)$. We choose $d_{\textnormal{max}}=10$ to obtain $65$ basis vectors for the BFM. The sensitivity-vectors estimated by PE, with the largest truncated state-space, are reported in Table A.6 and depicted as a bar chart in Figure 3(B). To obtain these PE-estimated sensitivity values, overall 13 minutes were required with a single processor on ETH Zurich’s Euler cluster. Figure 3(C) shows that as expected, the increase in memory requirements is linear in the truncated state-space size. Next we estimated a single steady-state sensitivity $\mathcal{S}_{\theta_{5}}(x_{2})$ using IntCLR with total processing time of $120$ hours. As the true sensitivity value is unknown, we estimate the accuracy of the simulation-based estimate by computing eRE$\%$ obtained by substituting the true sensitivity value with the PE-estimated value in (4.51). Like in the previous example, as the simulation time-period increases, the estimate becomes more accurate (i.e. eRE$\%$ decreases) but the statistical accuracy of the estimate (measured as standard deviation of the estimator) deteriorates as the number of simulation samples obtained in the limited processing time gets smaller (see Figure 3(D)).

4.3 A Deficiency Zero Network

We now consider a larger network where the stationary distribution is explicitly known and hence the steady-state sensitivities can be computed explicitly. Consider the four-species network depicted in Figure 4(A) whose reactions are described in Table A.7. The parameter vector for this network is $\theta=(\theta_{1},\dots,\theta_{6})=(4.5,0.8,5,1,0.6,11,3,80)$. In this example we are interested in estimating the sensitivity vectors for all the species (i.e. $S_{\theta}(x_{i})$ for $i=1,2,3,4$). Note that each reaction $k$ of the form (2.5), can be represented as $\nu_{k}\longrightarrow\nu^{\prime}_{k}$, where $\nu_{k}=(\nu_{1k},\nu_{2k},\dots)$ and $\nu^{\prime}_{k}=(\nu^{\prime}_{1k},\nu^{\prime}_{2k},\dots)$ are non-negative integer vectors representing the source complex and the product complex, respectively. The network under consideration is weakly-reversible (i.e. for each reaction $\nu_{k}\longrightarrow\nu^{\prime}_{k}$ there is a sequence of reactions realizing $\nu^{\prime}_{k}\longrightarrow\nu_{k}$) and it has deficiency zero 444The deficiency of a reaction network is defined as $\delta=C-l-s$, where $C$ is the number of complexes, $l$ is the number of connected components in the complex reaction graph (see Figure 4(A) for example) and $s$ is the dimension of the subspace spanned by the stoichiometry vectors $\zeta_{k}$-s.. Moreover all the propensities have mass-action form (2.6). From Theorem 4.2 in [37] it can be seen that the stationary distribution for this network has the product-form Poisson distribution

[TABLE]

where $c(\theta)=(c_{1}(\theta),c_{2}(\theta),c_{3}(\theta),c_{4}(\theta))$ is the mean vector given by

[TABLE]

Using this explicit form of the stationary distribution we can exactly compute the sensitivity vectors at our chosen value of parameter-vector $\theta$. These exact values allow us to compute RE$\%$ and assess the accuracy of our sensitivity estimates.

We now employ our PE method (Algorithm 1) to estimate the sensitivity-vectors. We use nine trapezoidal state-space truncations (see Table A.8) and set the designated state as $x_{\ell}=(10,0,0,0)$. We choose $d_{\textnormal{max}}=10$ to obtain $1000$ basis vectors for the BFM. The sensitivity-vectors estimated by PE, with the largest truncated state-space, are reported in Table A.9 and depicted as a bar chart in Figure 4(B). Note that for all these sensitivity estimates the RE$\%$ is less than $10^{-6}$ certifying the accuracy of our PE method. All these estimations took 73 minutes with a single processor on ETH Zurich’s Euler cluster. Like in the first example, as the size of the truncated state-spaces increases, the overall error (TRE$\%$) decreases exponentially but the memory requirements increase linearly (see Figure 4(C)). Next a single steady-state sensitivity $S_{\theta_{5}}(x_{4})$ was estimated using IntCLR with total processing time of $120$ hours. As before, increasing the simulation time-period reduces the error but worsens the statistical accuracy of the estimator (see Figure 4(D)). Remarkably, even the most accurate estimate has a RE$\%$ which is 1’000 times higher than the RE$\%$ for the corresponding PE-estimated value.

4.4 A Simple Linear Network

As our last example, we consider the simple linear network given in [21] involving three species and four conversion reactions (see Figure 5(A) and Table A.10). The vector of reaction-rate constants is $\theta=(\theta_{1},\theta_{2},\theta_{3},\theta_{4})$ and we set the parameters as $\theta_{1}=10,\theta_{2}=20,\theta_{3}=0.03$ and $\theta_{4}=0.02$. We are interested in estimating all sensitivities of the form $S_{\theta_{i}}(x_{j})$ for $i=1,\dots,4$ and $j=1,\dots,3$. Exploiting the linearity of the propensity function, we can analytically compute the steady-state expectations and the exact sensitivity values. These values are used for estimating the RE for our sensitivity estimates.

Notice that in this reaction network there is no inflow or outflow of mass and the total number of molecules of all the species is conserved, i.e. the CTMC dynamics for this reaction network evolves on the finite state-space

[TABLE]

where the constant $c$ is determined by the initial condition. We choose $c=10$ and for this value the state-space $\mathcal{E}$ is reasonably small and sFSP can yield the exact stationary distribution $\pi_{\theta}$ without requiring any state-space truncation. Hence when we estimate the sensitivity values using our PE method (Algorithm 1) the main source of error is the error incurred by BFM in solving the Poisson equation. It can be seen from Figure 5(B) and Table A.11, that this error is very small (RE$\%\leq 2\cdot 10^{-5}$ for all estimates) and the PE-based estimation is highly accurate. Here only two basis vectors were used for BFM which corresponds to $d_{\textnormal{max}}=1$. We then estimate one steady-state sensitivity value $\mathcal{S}_{\theta_{3}}(x_{1})$ with the simulation-based approach IntCLR with the total central processing time fixed to $120$ hours. The results are shown in Figure 5(C) for different values of the final time $T$ where the steady-state is assumed to be reached. As in the previous examples, when $T$ increases the accuracy of the IntCLR estimate improves (i.e. RE$\%$ decreases), but the estimator standard deviation also increases due to lower number of simulation samples being generated in the limited processing time of 120 hours. However for this simple example, PE method yields all the sensitivity estimates within just 0.2 seconds (with a single processor on ETH Zurich’s Euler cluster) and the accuracy of the PE estimate (measured by RE$\%$) was 900 times better than the best estimate provided by IntCLR.

5 Conclusion

The aim of this paper is to develop a numerical method for the estimation of the steady-state parameter sensitivity for stochastic reaction network models where the dynamics is described by a continuous-time Markov chain (CTMC) [3]. Unlike other sensitivity estimation approaches [11, 12, 13, 14, 15, 16, 17, 21], our method does not require time-consuming simulations of the stochastic trajectories, but instead it relies on approximately solving a Poisson equation associated with the generator of the CTMC as well as on estimation of the stationary distribution obtained with the recently proposed stationary Finite State Projection (sFSP) method [6]. We call our method the Poisson Estimator (PE) and we mathematically prove the accuracy of this method under certain conditions. Moreover using many examples we demonstrate that it compares favorably to simulation-based methods, like IntCLR [21], that exist for steady-state sensitivity estimation. In particular we found that estimation with PE requires far less computational time than with IntCLR and the estimates provided by PE are much more accurate than the corresponding IntCLR estimates. However PE has much higher memory requirements than IntCLR, and to circumvent this problem parallel computations may be necessary for even moderately sized networks.

In order to help the usability of our method we have developed a C++ library, called cossmosLib [27], that provides parallel implementation of the PE method for efficient estimation of both the stationary distribution (via sFSP) and the steady-state parameter sensitivities. Currently the parallelization in the cossmosLib library is achieved by incorporating the Message Passing Interface (MPI) and existing sparse matrix factorization methods that are part of the Trilinos [31, 32, 33] package. Further improvements in the performance of cossmosLib can be obtained by customizing the parallel computations in such a way that the structure and the transition kernel of the CTMC are properly exploited in finding the stationary distribution and solving the Poisson equation.

Appendix: Supplementary Tables

This section contains the tables describing the reactions networks, sensitivity analysis and truncated state-spaces for all the examples.