A Stochastic Hybrid Approximation for Chemical Kinetics Based on the Linear Noise Approximation

Abstract

The Linear Noise Approximation (LNA) is a continuous approximation of the CME, which improves scalability and is accurate for those reactions satisfying the leap conditions. We formulate a novel stochastic hybrid approximation method for chemical reaction networks based on adaptive partitioning of the species and reactions according to leap conditions into two classes, one solved numerically via the CME and the other using the LNA. The leap criteria are more general than partitioning based on population thresholds, and the method can be combined with any numerical solution of the CME. We then use the hybrid model to derive a fast approximate model checking algorithm for Stochastic Evolution Logic (SEL). Experimental evaluation on several case studies demonstrates that the techniques are able to provide an accurate stochastic characterisation for a large class of systems, especially those presenting dynamical stiffness, resulting in significant improvement of computation time while still maintaining scalability.

Appendices

A Proofs

Proposition 1

Let $x^s \in S^s$ and $x^f \in S^f$. Then, for $t \in \mathbb {R}_{\ge 0}$

$$\begin{aligned} \frac{d P(x^s|t)}{dt}=\sum _{\tau \in R}\beta _{\tau }(x^s-\upsilon _{\tau },t)P(x^s-\upsilon _{\tau },t)- \beta _{\tau }(x^s,t)P(x^s,t) \end{aligned}$$

where $\beta _{\tau }(x^s,t)= \sum _{x^f \in S^f} \alpha _{\tau }(x^f,x^s)P(x^f|x^s,t)$.

Proof

By using the law of total probability we have

$$\begin{aligned} \frac{d P(x^s|t)}{dt}=\sum _{x^f \in S^f}\frac{d P(x^s,x^f,t)}{dt} \end{aligned}$$

Then, using Eq. (2), and rearranging terms we have

$$\begin{aligned} \sum _{x^f \in S^f}\frac{d P(x^s,x^f,t)}{dt}= \end{aligned}$$

$$\begin{aligned} \sum _{x^f \in S^f}\sum _{\tau \in R^f}\alpha _{\tau }(x^f-\upsilon _{\tau },x^s-\upsilon _{\tau })P(x^f-\upsilon _{\tau },x^s-\upsilon _{\tau },t)- \alpha _{\tau }(x^f,x^s)P(x^f,x^s,t)= \end{aligned}$$

$$\begin{aligned} \sum _{\tau \in R}\beta _{\tau }(x^s-\upsilon _{\tau },t)P(x^s-\upsilon _{\tau },t)- \beta _{\tau }(x^s,t)P(x^s,t) \end{aligned}$$

where $\beta _{\tau }(x^s,t)= \sum _{x^f \in S^f} \alpha _{\tau }(x^f,x^s)P(x^f|x^s,t)$, that is, the conditional expectation of the propensity rate of $\tau$ at time t given $X^s(t)=x^s$.

Theorem 2

Assume $\varLambda ^s_t$ is non-empty and $S^s$ is the state space of $X^s(t)$. Then, the stochastic process $Z^h:\varOmega \times \mathcal {R}_{\ge 0}\rightarrow \mathcal {S}$, with $\varOmega$ its sample space and $(\mathcal {S},\mathcal {B})$ a measurable space, is such that for $A\in \mathcal {B}$ and $t\in \mathbb {R}_{\ge 0}$

$$\begin{aligned} P(Z^h(t)\in A|X(0)=x_0)=\sum _{x^s \in S^s} P(Z_{x^s}(t)\in A)P(X^s(t)=x_s) \end{aligned}$$

where $Z_{x^s}(t)$ is a Gaussian random variable with expected value and variance

$$\begin{aligned} E[Z_{x^s}(t)]=B\cdot \begin{pmatrix} E[\bar{X}^f(t)]\\ x^s \end{pmatrix} \quad C[Z_{x^s}(t)]=B\cdot \begin{pmatrix} C[\bar{X}^f(t)]&{}0\\ 0&{}0 \end{pmatrix} \cdot B^T \end{aligned}$$

where $\bar{X}^f$ is the virtual fast process introduced in Sect. 3.

Proof

By the law of total probability we have

$$ P(Z(t)\in A|X(0)=x_0)=\sum _{x^s \in S^s}P(Z(t)\in A|X^s(t)\!=\!x^s,X(0)=x_0)P(X^s(t)\!=\!x^s|X(0)=x_0). $$

By application of the LNA it follows that $X^f(t)$ conditioned on the event $X^s(t)=x^s$ is a Gaussian random variable with expected value and variance

$$ E[X^f(t)|X^s(t)=x^s]=\begin{pmatrix} E[\bar{X}^f(t)]\\ x^s \end{pmatrix} $$

and covariance matrix

$$ C[X^f(t)|X^s(t)=x^s]=\begin{pmatrix} C[\bar{X}^f(t)]&{}0\\ 0&{}0 \end{pmatrix} $$

Given a multidimensional Gaussian distribution, each linear combination of its components is still Gaussian. As a consequence, $E[Z^h(t)|X^s(t)=x^s]=B\cdot E[X^f(t)|X^s(t)=x^s]$ and $C[Z^h(t)|X^s(t)=x^s]=B\cdot C[X^f(t)|X^s(t)=x^s]\cdot B^T$.

Theorem 3

Assume $\varLambda ^s_t$ is non-empty. Then, for $t\in \mathbb {R}_{\ge 0}$

$$\begin{aligned} E[Z^h(t)|X(0)=x_0]=\sum _{x^s \in S^s} E[Z_{x^s}(t)\in A]P(X^s(t)=x_s) \end{aligned}$$

$$\begin{aligned} C[Z^h(t)|X(0)=x_0]=\sum _{x^s \in S^s} C[Z_{x^s}(t)\in A]P(X^s(t)=x_s) \end{aligned}$$

Proof

The proof follows from the application of the law of total expectation for random variables with mutually exclusive and exhaustive events.

B Figures

Fig. 2.

Comparison of expected value and variance of mRNA in Example 2 in interval [0, 200] as calculated by direct solution of the CME (Fig. 2a) and by our algorithm (Fig. 2b).

Fig. 3.

Comparison of the probability distribution of RNA at time $t=200$ as calculated by numerical hybrid algorithm (Fig. 3a) and by the LNA (Fig. 3b).