Effective computational methods for hybrid stochastic gene networks

TL;DR

This paper introduces improved computational methods for analyzing hybrid stochastic gene networks modeled by PDMPs, enabling more accurate and efficient calculation of their probability distributions.

Contribution

The authors develop novel algorithms combining simulation and analytic techniques to better compute probability distributions in hybrid gene network models.

Findings

Monte-Carlo method effectively combines simulation with analytic ODE solutions.

Push-forward method accurately computes probability measures using analytic semigroup expressions.

Enhanced methods outperform earlier approaches in handling non-linear regulation functions.

Abstract

At the scale of the individual cell, protein production is a stochastic process with multiple time scales, combining quick and slow random steps with discontinuous and smooth variation. Hybrid stochastic processes, in particular piecewise-deterministic Markov processes (PDMP), are well adapted for describing such situations. PDMPs approximate the jump Markov processes traditionally used as models for stochastic chemical reaction networks. Although hybrid modelling is now well established in biology, these models remain computationally challenging. We propose several improved methods for computing time dependent multivariate probability distributions (MPD) of PDMP models of gene networks. In these models, the promoter dynamics is described by a finite state, continuous time Markov process, whereas the mRNA and protein levels follow ordinary differential equations (ODEs). The Monte-Carlo method combines direct simulation of the PDMP with analytic solutions of the ODEs. The push-forward method numerically computes the probability measure advected by the deterministic ODE flow, through the use of analytic expressions of the corresponding semigroup. Compared to earlier versions of this method, the probability of the promoter states sequence is computed beyond the naive mean field theory and adapted for non-linear regulation functions.

Figures

Tables

Table 1. Table 1: Execution times for different methods. All the methods were implemented in Matlab R2013b running on a single core (multi-threading inactivated) of a Intel i5-700u 2.5 GHz processor. The Monte-Carlo method computed the next jump waiting time using the analytical solution of Eq. 14 for M 1 subscript 𝑀 1 M_{1} and the numerical solution of Eq. 14 for M 2 subscript 𝑀 2 M_{2} . The push-forward method used analytic solutions for mRNA and protein trajectories from ( 15 ),( 16 ), and numerical computation of the integrals in Eq. 22 , for both models.

Model	Monte-Carlo [min]	Push-forward [s]
$M1$ slow-slow	45	20
$M1$ fast-fast	74	30
$M2$ slow-slow	447	20
$M2$ fast-fast	758	30

Equations

$$\frac{\mathrm{d}{\vec{\Phi }}_{\vec{s}}(t,\vec{x},\vec{y})}{\mathrm{d}t}={\vec{F}}_{\vec{s}}({\vec{\Phi }}_{\vec{s}}(t,\vec{x},\vec{y})),{\vec{\Phi }}_{\vec{s}}(0,\vec{x},\vec{y})=(\vec{x},\vec{y}).$$

$$\frac{\mathrm{d}{\Phi }_i^x}{\mathrm{d}t}$$

$$\frac{\mathrm{d}{\Phi }_i^y}{\mathrm{d}t}$$

$$\lambda (\vec{x},\vec{y},\vec{s})=\sum _{\vec{r}\ne \vec{s}}{H}_{\vec{r},\vec{s}}(\vec{x},\vec{y})=-H_{\vec{s},\vec{s}}(\vec{x},\vec{y}).$$

$${\vec{x}}_t(\omega ):={\vec{\Phi }}_{{\vec{s}}_0}(t,{\vec{x}}_0,{\vec{y}}_0)\text{ for }0\le t\le T_1(\omega ),$$

$$F(t)=\mathbb{P}\left[T_1>t\right]=\exp (-{\int }_0^t\lambda ({\vec{\Phi }}_{{\vec{s}}_0}(u,{\vec{x}}_0,{\vec{y}}_0))du),t>0,$$

$$\mathbb{P}\left[\vec{s}=\vec{r}\right]=\frac{H_{\vec{r},\overrightarrow{s_0}}({\vec{\Phi }}_{{\vec{s}}_0}(T_1(\omega ),{\vec{x}}_0,{\vec{y}}_0),\overrightarrow{s_0})}{\lambda ({\vec{\Phi }}_{{\vec{s}}_0}(T_1(\omega ),{\vec{x}}_0,{\vec{y}}_0),\overrightarrow{s_0})},\text{ for all }\vec{r}\ne {\vec{s}}_0.$$

$$\frac{\mathrm{d}{\Phi }_i^x}{\mathrm{d}t}$$

$$\frac{\mathrm{d}{\Phi }_i^y}{\mathrm{d}t}$$

$$\frac{\mathrm{d}\log F}{\mathrm{d}t}$$

$$\frac{\partial p_{\vec{s}}(t,\vec{x},\vec{y})}{\partial t}=-{\nabla }_{\vec{x},\vec{y}}.({\vec{F}}_{\vec{s}}(\vec{x},\vec{y})p_{\vec{s}}(t,\vec{x},\vec{y}))+\sum _{\vec{r}}{H}_{\vec{s},\vec{r}}(\vec{s},\vec{x},\vec{y})p_{\vec{r}}(t,\vec{x},\vec{y}).$$

$$\left[\begin{array}{cccc}-(f_1+f_2{x}_1)& h_1& h_2& 0\\ f_1& -(h_1+f_2{x}_1)& 0& h_2\\ f_2{x}_1& 0& -(f_1+h_2)& h_1\\ 0& f_2{x}_1& f_1& -(h_1+h_2)\end{array}\right].$$

$$\frac{\partial p_1}{\partial t}$$

$$\frac{\partial p_2}{\partial t}$$

$$\frac{\partial p_3}{\partial t}$$

$$\frac{\partial p_4}{\partial t}$$

$$\left[\begin{array}{cccc}-(f_1+f_2{x}_1/(K_1+x_1))& h_1& h_2& 0\\ f_1& -(h_1+f_2{x}_1/(K_1+x_1))& 0& h_2\\ f_2{x}_1/(K_1+x_1)& 0& -(f_1+h_2)& h_1\\ 0& f_2{x}_1/(K_1+x_1)& f_1& -(h_1+h_2)\end{array}\right].$$

$$\frac{\partial p_1}{\partial t}$$

$$\frac{\partial p_2}{\partial t}$$

$$\frac{\partial p_3}{\partial t}$$

$$\frac{\partial p_4}{\partial t}$$

$${\Phi }_i^x(t,x_0,y_0)$$

$${\Phi }_i^y(t,x_0,y_0)$$

$$\lambda (\vec{x},\vec{y},\vec{s})=c_0(\vec{s})+\sum _i^N{c}_i(\vec{s})x_i+\sum _i^N{d}_i(\vec{s})f_i(x_i),$$

$$\lambda (\vec{x},\vec{y},\vec{s})=\left(f_1+f_2\frac{x_1}{K_1+x_1}\right){\delta }_{s,1}+\left(h_1+f_2\frac{x_1}{K_1+x_1}\right){\delta }_{s,2}+(h_2+f_1){\delta }_{s,3}+(h_2+h_1){\delta }_{s,4},$$

$$\begin{split}&-\log(U)=\bigg{[}(f_{1}+f_{2})T_{1}+f_{2}\int_{0}^{T_{1}}\frac{1}{K_{1}+\Phi_{1}^{x}(t^{\prime},x_{0},y_{0})}\,dt^{\prime}\bigg{]}\delta_{s_{0},1}+\bigg{[}(h_{1}+f_{2})T_{1}+\\ &f_{2}\int_{0}^{T_{1}}\frac{1}{K_{1}+\Phi_{1}^{x}(t^{\prime},x_{0},y_{0})}\,dt^{\prime}\bigg{]}\delta_{s_{0},2}+(h_{2}+f_{1})T_{1}\delta_{s_{0},3}+(h_{2}+h_{1})T_{1}\delta_{s_{0},4},\end{split}$$

$$y_i(\tau )$$

$$x_i(\tau )$$

$$p^{l_i,m_i}=\sum _{S_M}\sum _{\psi (l_i^{\prime },m_i^{\prime })=(l_i,m_i)}{p}_0^{l_i^{\prime },m_i^{\prime }}\mathbb{P}\left[S_M\right].$$

$$\mathbb{P}\left[S_M\right]={\Pi }_{s_{N-1},s_{N-2}}({\tau }_{N-2},{\tau }_{N-1})\dots {\Pi }_{s_1,s_0}({\tau }_0,{\tau }_1)P_0^S(s_0),$$

$$\vec{\Pi }({\tau }_j,{\tau }_{j+1})=\exp \left({\int }_{{\tau }_j}^{{\tau }_{j+1}}\vec{H}(\vec{x}(t))dt\right),$$

$$\mathbb{E}\left[\frac{f_2{x}_1(t)}{K_1+x_1(t)}\right]\approx \frac{f_2\mathbb{E}\left[x_1(t)\right]}{(K_1+\mathbb{E}\left[x_1(t)\right])}-\frac{f_2}{(K_1+\mathbb{E}\left[x_1(t)\right]{)}^3}Var(x_1(t)),$$

$${\vec{\Pi }}_1(\tau ,{\tau }^{\prime })=\left[\begin{array}{cc}(1-p_{1,on})+p_{1,on}{e}^{-{ϵ}_1({\tau }^{\prime }-\tau )},& (1-p_{1,on})-(1-p_{1,on})e^{-{ϵ}_1({\tau }^{\prime }-\tau )}\\ p_{1,on}-p_{1,on}{e}^{-{ϵ}_1({\tau }^{\prime }-\tau )},& p_{1,on}+(1-p_{1,on})e^{-{ϵ}_1({\tau }^{\prime }-\tau )}\end{array}\right],$$

$${\vec{\Pi }}_2(\tau ,{\tau }^{\prime })=$$

$$\left[\begin{array}{cc}{e}^{-{\int }_{\tau }^{{\tau }^{\prime }}(h_2+F_2(t))dt}+h_2{\int }_{\tau }^{{\tau }^{\prime }}{e}^{-{\int }_t^{{\tau }^{\prime }}(h_2+F_2(t^{\prime }))d{t}^{\prime }}dt,& h_2{\int }_{\tau }^{{\tau }^{\prime }}{e}^{-{\int }_t^{{\tau }^{\prime }}(h_2+F_2(t^{\prime }))d{t}^{\prime }}dt\\ 1-e^{-{\int }_{\tau }^{{\tau }^{\prime }}(h_2+F_2(t))dt}-h_2{\int }_{\tau }^{{\tau }^{\prime }}{e}^{-{\int }_t^{{\tau }^{\prime }}(h_2+F_2(t^{\prime }))d{t}^{\prime }}dt,& 1-h_2{\int }_{\tau }^{{\tau }^{\prime }}{e}^{-{\int }_t^{{\tau }^{\prime }}(h_2+F_2(t^{\prime }))d{t}^{\prime }}dt\end{array}\right],$$

$$d=\int |p(x)-\tilde{p}(x)|dx.$$

Full text

Availability. The algorithms described in this paper were implemented in MATLAB. The code is available on demand.

1 Introduction

In PDMP models of gene networks, each gene promoter is described as a finite state Markov process [3, 10, 13, 12]. The promoter triggers synthesis of gene products (mRNAs and proteins) with intensities depending on its state. The promoter can exhibit two state (ON-OFF) dynamics, but also dynamics with more than two states and arbitrarily complex transitions[11, 19]. The transition rates between the states of the promoter depend on the expression levels of proteins expressed by the same or by other promoters. In PDMP models, the gene products are considered in sufficiently large copy numbers and are represented as continuous variables following ordinary differential equations (ODEs). The sources of noise in these models are thus the discrete transitions between the promoter states.

In single cell experimental settings the quantities of mRNA [18, 14, 1, 17] and proteins [5, 6] can be determined for each cell. By double or multiple- fluorophore fluorescence techniques products from several genes can be quantified simultaneously and one can have access to multivariate probability distributions (MPD) of mRNA or proteins. The stochastic dynamics of promoters and gene networks can have important consequences for fundamental biology [4] but also for HIV [15] and cancer research [7]. For this reason we aim to develop effective methods for computing time-dependent MPDs for PDMP models. Our main objective is the reduction of computation time which is prerequisite for parameter scans and machine learning applications [8].

PDMPs already represent a gain with respect to the chemical Markov equation from which they are derived by various limit theorems [2]. A gene network PDMP model can be simulated by numerical integration of ODEs coupled with a driven inhomogeneous Poisson process for the successive transitions of the promoters [20, 3, 16, 13]. The simulation becomes particularly effective when analytic solutions of the ODEs are available [10].

However, very little has been done to further improve the computational power by optimising simulation and analysis of PDMP models.

Numerical integration of the PDE satisfied by MPD is an interesting option combining precision and speed for small models. Finite difference methods, however, are of limited use in this context as they can not cope with many RNA and protein variables (extant examples are restricted to the dimension 2, corresponding to a single promoter, with or without self-regulation see [10, 12]).

Another interesting method for computing time dependent MPDs is the push-forward method. For gene networks, this method has been first introduced in [9] and further adapted for continuous mRNA variables in [10]. It is based on the idea to compute the MPD as the push-forward measure of the semigroup defined by the ODEs. This method is approximate, as one has to consider that the discrete PDMP variables are piecewise constant on a deterministic time partition. Furthermore, the transition rates between promoter states were computed in a mean field approximation. In this paper we replace the mean field approximation by the next order approximation taking into account the moments of the protein distribution.

2 Methods

2.1 PDMP models of gene networks

The state of a PDMP gene network model takes values in $E={\mathbb{R}}^{2N}\times\{0,\ldots,s_{max}\}^{N}-1$, where $N$ is the number of genes and $s_{max}$ is a positive integer representing the maximum number of states of a gene promoter. It is a process $\vec{\zeta}_{t}=(\vec{x}_{t},\vec{y}_{t},\vec{s}_{t})$, determined by three characteristics:

1)

For all $\vec{s}\in\{0,\ldots,s_{max}\}^{N}$ a vector field $\vec{F}_{\vec{s}}:{\mathbb{R}}^{2N}\to{\mathbb{R}}^{2N}$ determining a unique global flow $\vec{\Phi}_{\vec{s}}(t,\vec{x},\vec{y})$ in ${\mathbb{R}}^{2N}$, the space of all protein ($\vec{x}\in{\mathbb{R}}^{N}$) and mRNA ($\vec{y}\in{\mathbb{R}}^{N}$) values such that, for $t>0$,

[TABLE]

On coordinates, this reads

[TABLE]

where $b_{i}$, $k_{i}$, $a_{i}$, $\rho_{i}$ are translation efficiencies, transcription rates, protein degradation coefficients and mRNA degradation coefficients of the $i^{\text{th}}$ gene, respectively. Note that transcription rates depend on the relevant promoter states.

The flow $\vec{\Phi}_{\vec{s}}(t,\vec{x},\vec{y})$ represents a one parameter semigroup fulfilling the properties

(i)

$\vec{\Phi}_{\vec{s}}(0,\vec{x}_{0},\vec{y}_{0})=(\vec{x}_{0},\vec{y}_{0})$,

(ii)

$\vec{\Phi}_{\vec{s}}(t+t^{\prime},\vec{x}_{0},\vec{y}_{0})=\vec{\Phi}_{\vec{s}}(t^{\prime},\vec{\Phi}_{\vec{s}}^{x}(t,\vec{x}_{0},\vec{y}_{0}),\vec{\Phi}_{\vec{s}}^{y}(t,\vec{x}_{0},\vec{y}_{0}))$.

2)

A transition rate matrix for the promoter states $\vec{H}:{\mathbb{R}}^{2N}\to M_{N\times N}({\mathbb{R}})$, such that $H_{\vec{s},\vec{r}}(\vec{x},\vec{y})\geq 0$ and $H_{\vec{s},\vec{s}}(\vec{x},\vec{y})=-\sum_{\vec{r}\neq\vec{s}}H_{\vec{r},\vec{s}}(\vec{x},\vec{y})$ for all $\vec{s},\vec{r}\in\{0,\ldots,s_{max}\}^{N},\vec{s}\neq\vec{r}$ and for all $(\vec{x},\vec{y})\in R^{2N}$.

3)

A jump rate $\lambda:E\to{\mathbb{R}}^{+}$. The jump rate can be obtained from the transition rate matrix

[TABLE]

From these characteristics, right-continuous sample paths $\{(\vec{x}_{t},\vec{y}_{t}):t>0\}$ starting at $\vec{\zeta}_{0}=(\vec{x}_{0},\vec{y}_{0},\vec{s}_{0})\in E$ can be constructed as follows. Define

[TABLE]

where $T_{1}(\omega)$ is a realisation of the first jump time of $\vec{s}$, with the distribution

[TABLE]

and $\omega$ is the element of the probability space for which the particular realisation of the process is given. The pre-jump state is $\vec{\zeta}_{T_{1}^{-}(\omega)}(\omega)=(\vec{\Phi}_{\vec{s}_{0}}(T_{1}(\omega),\vec{x}_{0},\vec{y}_{0}),\vec{s}_{0})$ and the post-jump state is $\vec{\zeta}_{T_{1}(\omega)}(\omega)=(\vec{\Phi}_{\vec{s}_{0}}(T_{1}(\omega),\vec{x}_{0},\vec{y}_{0}),\vec{s})$, where $\vec{s}$ has the distribution

[TABLE]

We then restart the process $\vec{\zeta}_{T_{1}(\omega)}$ and recursively apply the same procedure at jump times $T_{2}(\omega)$, etc..

Note that between each two consecutive jumps $(\vec{x}_{t},\vec{y}_{t})$ follow deterministic ODE dynamics defined by the vector field $\vec{F}_{\vec{s}}$. At the jumps, the protein and mRNA values $(\vec{x}_{t},\vec{y}_{t})$ are continuous.

The calculation of the flow between two jumps and of the jump time can be gathered in the same set of differential equations

[TABLE]

that has to be integrated with the stopping condition $F(T_{1})=U$, where $U$ is a random variable, uniformly distributed on $[0,1]$.

We define multivariate probability density functions $p_{\vec{s}}(t,\vec{x},\vec{y})$. These functions satisfy the Liouville-master equation which is a system of partial differential equations:

[TABLE]

2.2 ON/OFF gene networks

In this paper, for the purpose of illustration only, all the examples are constituted by ON/OFF gene networks.

For an ON/OFF gene each component $s_{i}$ has two possible values [math] for OFF and $1$ for ON.

As a first example that we denote as model $M1$, let us consider a two genes network; the expression of the first gene being constitutive and the expression of the second gene being activated by the first. We consider that the transcription activation rate of the second gene is proportional to the concentration of the first protein $f_{2}x_{1}$. All the other rates are constant $f_{1}$, $h_{1}$, $h_{2}$, representing the transcription activation rate of the first gene, and the transcription inactivation rates of gene one and gene two, respectively. For simplicity, we consider that the two genes have identical protein and mRNA parameters $b_{1}=b_{2}=b$, $a_{1}=a_{2}=a$, $\rho_{1}=\rho_{2}=\rho$. We further consider that $k_{i}(s_{i})=k_{0}$ if the gene $i$ is OFF and $k_{i}(s_{i})=k_{1}$ if the gene $i$ is ON.

The gene network has four discrete states, in order $(0,0)$, $(1,0)$, $(0,1)$, and $(1,1)$. Then, the transition rate matrix for the model $M_{1}$ is

[TABLE]

The Liouville-master equation for the model $M1$ reads

[TABLE]

The model $M_{2}$ differs from the model $M_{1}$ by the form of the activation function. Instead of a linear transcription rate $f_{2}x_{1}$ we use a Michaelis-Menten model $f_{2}x_{1}/(K_{1}+x_{1})$. This model is more realistic as it takes into account that the protein $x_{1}$ has to attach to specific promoter sites which become saturated when the concentration of this protein is high.

The transition rate matrix for the model $M_{2}$ is

[TABLE]

The Liouville-master equation for the model $M2$ reads

[TABLE]

2.3 Monte-Carlo method

The Monte-Carlo method utilizes the direct simulation of the PDMP based on Eq.7. A larger number $M$ of sample paths is generated and the values of $(\vec{x}_{t},\vec{y}_{t})$ are stored at selected times. Multivariate probability distributions are then estimated from this data.

The direct simulation of PDMPs needs the solutions of (7) which can be obtained by numerical integration. This is not always computationally easy. Problems may arise for fast switching promoters when the ODEs have to be integrated many times on small intervals between successive jumps. Alternatively, the numerical integration of the ODEs can be replaced by analytic solutions or quadratures. Analytic expressions are always available for the gene network flow (2) and read

[TABLE]

Let us consider the following general expression of the jump intensity function

[TABLE]

where $f_{i}$ are non-linear functions, for instance Michaelis-Menten $f_{i}(x_{i})=x_{i}/(K_{i}+x_{i})$ or Hill functions $f_{i}(x_{i})=x_{i}^{n_{i}}/(K_{i}^{n_{i}}+x_{i}^{n_{i}})$. If $d_{i}=0$ for all $1\leq i\leq N$, the cumulative distribution function of the waiting time $T_{1}$ can be solved analytically [10], otherwise it can be obtained by quadratures. For example, for the model $M_{2}f$ one has

[TABLE]

where $\delta_{i,j}$ is Kronecker’s delta. In this case the waiting time $T_{1}$ is obtained as the unique solution of the equation

[TABLE]

where $U$ is a random variable, uniformly distributed on $[0,1]$. In our implementation of the algorithm we solve (14) numerically, using the bisection method.

2.4 Push-forward method

This method allows one to compute the MPD of proteins and mRNAs at a time $\tau$ given the MPD of proteins and mRNAs at time [math].

In order to achieve this we use a deterministic partition $\tau_{0}=0<\tau_{1}<\ldots<\tau_{M}=\tau$ of the interval $[0,\tau]$ such that $\Delta_{M}=\max_{j\in[1,M]}(\tau_{j}-\tau_{j-1})$ is small. The main approximation of this method is to consider that $s(t)$ is piecewise constant on this partition, more precisely that $s(t)=s_{j}:=s(\tau_{j}),\,\text{for }t\in[\tau_{j},\tau_{j+1}),\,0\leq j\leq M-1$. This approximation is justified by Theorem 3.1 in Section 3.1.

For each path realisation $S_{M}:=(s_{0},s_{1},\ldots,s_{M-1})\in\{0,1,\ldots,2^{N}-1\}^{M}$ of the promoter states, we can compute (see Appendix 2) the protein and mRNA levels $x(\tau),y(\tau)$ of all genes $i\in[1,N]$:

[TABLE]

where $w_{j}=\frac{e^{(a-\rho)\tau}-e^{(a-\rho)\tau_{j}}}{a-\rho}(e^{\rho\tau_{j}}-e^{\rho\tau_{j-1}})-\frac{e^{(a-\rho)\tau_{j}}-e^{(a-\rho)\tau_{j-1}}}{a-\rho}e^{\rho\tau_{j-1}}+\frac{e^{a\tau_{j}}-e^{a\tau_{j-1}}}{a}$ and $s_{j}^{i}:=0$ if promoter $i$ is OFF for $t\in[\tau_{j},\tau_{j+1})$ and $s_{j}^{i}:=1$ if promoter $i$ is ON for $t\in[\tau_{j},\tau_{j+1})$.

In order to compute the MDP at time $\tau$ one has to sum the contributions of all solutions (15),(16), obtained for the $2^{NM}$ realisations of promoter state paths with weights given by the probabilities of the paths.

Eqs.15,16 can straightforwardly be adapted to compute $x(t),y(t)$ for all $t\in[0,\tau]$. To this aim, $\tau$ should be replaced by $t$ and $M$ should be replaced by $M_{t}$ defined by the relation $t\in[\tau_{M_{t}},\tau_{{M_{t}}+1}]$.

Suppose that we want to estimate the MDP of all mRNAs and proteins of the gene network, using a multivariate histogram with bin centers $(x_{0}^{l_{i}},y_{0}^{m_{i}}),\,1\leq i\leq N,1\leq l_{i}\leq n_{x},1\leq m_{i}\leq n_{y}$ where $n_{x}$, $n_{y}$ are the numbers of bins in the protein and mRNA directions for each gene, respectively. Typically $x_{0}^{l_{i}}=b/(a\rho)(k0+(k1-k0)(l_{i}-1/2)),1\leq i\leq N,1\leq l_{i}\leq n_{x}$, $y_{0}^{m_{i}}=1/\rho(k0+(k1-k0)(m_{i}-1/2)),1\leq i\leq N,1\leq m_{i}\leq n_{y}$. The initial MDP at time $t=0$ is given by the bin probabilities $p_{0}^{l_{i},m_{i}},1\leq i\leq N,1\leq l_{i}\leq n_{x},1\leq m_{i}\leq n_{y}$. Let $(x^{l_{i},m_{i}},y^{l_{i},m_{i}})$ be the solutions , with $x_{i}(0)=x_{0}^{l_{i}}$ and $y_{i}(0)=y_{0}^{m_{i}}$. The many-to-one application $(l_{i}^{\prime},m_{i}^{\prime})=\psi(l_{i},m_{i})$ provides the histogram bin $(l_{i}^{\prime},m_{i}^{\prime})$ in which falls the vector $(x^{l_{i},m_{i}},y^{l_{i},m_{i}}))$. The push forward MDP at time $t=\tau$ is defined by the bin probabilities $p^{l_{i},m_{i}},1\leq i\leq N,1\leq l_{i}\leq n_{x},1\leq m_{i}\leq n_{y}$ that are computed as

[TABLE]

In order to compute ${\mathbb{P}}\!\left[S_{M}\right]$ we can use the fact that, given $x(t)$, $s(t)$ is a finite state Markov process, therefore

[TABLE]

where $P_{0}^{S}:\{0,1,\ldots,2^{N}-1\}\to[0,1]$ is the initial distribution of the promoter state,

[TABLE]

and $\vec{x}(t)$ are given by (16).

The push-forward method can be applied recursively to compute the MDP for times $\tau,\,2\tau,\ldots,n_{t}\tau$. The complexity of the calculation scales as $n_{t}(n_{x})^{N}(n_{y})^{N}2^{NM}$ which is exponential in the number of genes $N$. The exponential complexity comes from considering all the $2^{NM}$ possible paths $S_{M}$. However, many of these paths have almost the same probability and impose very similar trajectories to the variables $(\vec{x}(t),\vec{y}(t))$. In fact, a convenient approximation is to consider that different genes are switching between ON and OFF states according to Markov processes with rates given by the mean values of regulatory proteins (mean field approximation, [10]). This approximation consists in applying the push-forward procedure for each gene separately, using averaged transition probabilities. Thus, the $2^{N}$ states transition matrix $\vec{H}$ has to be replaced by $N$, $2\times 2$ state transition matrices for each gene. This approximation reduces the complexity of the calculations to $n_{t}n_{x}n_{y}N2^{M}$ which is linear in the number of genes.

In [10] we have replaced the regulation term $f_{2}x_{1}(t)$ occurring in the transition matrix by its mean $f_{2}{\mathbb{E}}\!\left[x_{1}(t)\right]$. In this case both $\vec{H}$ and $\vec{\Pi}$ can be computed analytically, which leads to a drastic reduction in the execution time. This approach is suitable for the model $M_{1}$, which contains only linear regulation terms. For non-linear regulation terms, $\vec{\Pi}$ can not generally be computed analytically. Furthermore, the mean field approximation introduces biases. For instance, in the case of the model $M_{2}$, the approximation $f_{2}x_{1}(t)/(K_{1}+x_{1}(t))\approx f_{2}{\mathbb{E}}\!\left[x_{1}(t)\right]/(K_{1}+{\mathbb{E}}\!\left[x_{1}(t)\right])$ is poor. A better approximation in this case is to replace $f_{2}x_{1}(t)/(K_{1}+x_{1}(t))$ by its mean and use

[TABLE]

in order to correct the bias. Here $Var$ indicates the variance.

As in [10] we can use analytic expressions for ${\mathbb{E}}\!\left[x_{1}(t)\right]$, but also for $Var(x_{1}(t))$. These expressions can be found in Appendix 1. Although the elements of matrix $\vec{H}$ have analytic expressions, the elements of the matrix $\vec{\Pi}$ contain integrals that must be computed numerically. For the model $M_{2}$, we have

[TABLE]

for the transition rates of the first gene, where $p_{1,on}=f_{1}/(f_{1}+h_{1})$, $\epsilon_{1}=(f_{1}+h_{1})/\rho$, and

[TABLE]

for the transitions of the second gene, where $F_{2}(t)=f_{2}{\mathbb{E}}\!\left[\frac{x_{1}(t)}{K_{1}+x_{1}(t)}\right]$.

3 Results

3.1 Convergence of the push-forward method

The probability distribution obtained with the push-forward method converges to the exact PDMP distribution in the limit $M\to\infty$. This is a consequence of the following theorem

Theorem 3.1

Let $\Phi_{S_{M}}(t,\vec{x},\vec{y})$ be the flow defined by the formulas (16),(15), such that $(\vec{x}(t),\vec{y}(t))=\Phi_{S_{M}}(t,\vec{x}(0),\vec{y}(0))$ for $t\in[0,\tau]$, and let $\mu_{t}^{M}:\mathcal{B}({\mathbb{R}}^{2N})\to{\mathbb{R}}_{+}$ be the probability measure defined as $\mu_{t}^{M}(A)=\sum_{S_{M}}{\mathbb{P}}\!\left[S_{M}\right]\mu_{0}(\Phi_{S_{M}}^{-1}(t,A))$, where $\mu_{0}:\mathcal{B}({\mathbb{R}}^{2N})\to{\mathbb{R}}_{+}$ is the probability distribution of $(\vec{x},\vec{y})$ at $t=0$, ${\mathbb{P}}\!\left[S_{M}\right]$ are given by (18), and $\mathcal{B}({\mathbb{R}}^{2N})$ are the Borel sets on ${\mathbb{R}}^{2N}$. Let $\mu_{t}$, the exact distribution of $(\vec{x}(t),\vec{y}(t))$ for the PDMP defined by (1),(2),(3), with initial values $(\vec{x}_{0},\vec{y}_{0},\vec{s}_{0})$ distributed according to $\mu_{0}\times P_{0}^{S}$. Assume that $|\tau_{i}-\tau_{i-1}|<C/M$ for all $i\in[1,M]$, where $C$ is a positive constant. Then, for all $t\in[0,\tau]$, $\mu_{t}^{M}$ converges in distribution to $\mu_{t}$, when $M\to\infty$.

The proof of this theorem is given in the Appendix 3.

3.2 Testing the push-forward method

In order to test the push-forward method, we compared the resulting probability distributions with the ones obtained by the Monte Carlo method using the direct simulation of the PDMP. We considered the models $M_{1}$ and $M_{2}$ with the following parameters: $\rho=1$, $p1=1/2$, $a=1/5$, $b=4$, $k_{0}=4$, $k_{1}=40$ for the two genes. The parameter $\epsilon$ took two values $\epsilon=0.5$ for slow genes and $\epsilon=5.5$ for fast genes. We tested the slow-slow and the fast-fast combinations of parameters.

The initial distribution of the promoters states was $P_{0}^{S}((0,0))=1$ where the state $(0,0)$ means that both promoters are OFF. The initial probability measure $\mu_{0}$ was a delta Dirac distribution centered at $x_{1}=x_{2}=0$ and $y_{1}=y_{2}=0$. This is obtained by always starting the direct simulation of the PDMP from $x_{1}(0)=x_{2}(0)=0$, $y_{1}(0)=y_{2}(0)=0$, and $s_{1}(0)=s_{2}(0)=0$. The simulations were performed between $t_{0}=0$ and $t_{max}=20$ for fast genes and between $t_{0}=0$ and $t_{max}=90$ for slow genes. In order to estimate the distributions we have used $MC=50000$ samples.

The push-forward method was implemented with $M=10$ equal length sub-intervals of $[0,\tau]$. The time step $\tau$ was chosen $\tau=2$ for fast genes and $\tau=15$ for slow genes. The procedure was iterated $10$ times for fast genes (up to $t_{max}=20$) and $6$ times for slow genes (up to $t_{max}=90$).

The execution times are provided in the Table 1. The comparison of the probability distributions are illustrated in the Figures 1,2. In order to quantify the relative difference between methods we use the $L^{1}$ distance between distributions. More precisely, if $p(x)$ and $\tilde{p}(x)$ are probability density functions to be compared, the distance between distributions is

[TABLE]

4 Discussion and conclusion

Combining direct simulation of PDMP gene network models and analytic formulas for the ODE flow represents an effective, easy to implement method for computing time dependent MPD of these models. However, the precision of the Monte-Carlo estimates of the distributions increases like $\sqrt{MC}$, where $MC$ is the number of Monte-Carlo samples. For this reason, the execution time of this method, although smaller compared to PDMP simulation methods that implement numerical resolution of the ODEs such as reported in [13] (data not shown), is large compared to deterministic methods such as the push-forward method.

The push-forward method represents an effective alternative to Monte-Carlo methods, ensuring reduced execution time. With respect to an earlier implementation of this method in [9] we used promoter states instead of mRNA copy numbers as discrete variables of the PDMP. As a consequence, the number of discrete states is lower and we can afford increasing the number $M$ of time subdivisions. Compared to the similar work in [10] we used second moments of the protein distribution which took into account the correlation of the promoter states and lead to increased accuracy in the case of nonlinear regulation. We proved rigorously the convergence of the distributions calculated with the push-forward method to the exact distributions of the PDMP. However, the push-forward method is an approximate method, and its accuracy relies on the careful choice of the time and space steps, namely of the integers $M$, $n_{t}$, $n_{x}$, $n_{y}$. We will present elsewhere error estimates allowing an optimal choice of these parameters. Although the protein moments and the exponential transition rate matrix $\vec{\Pi}$ can be computed numerically, the effectiveness of the push-forward method is increased when analytic expressions are available for these quantities. In this paper, these expressions were computed for particular cases. In the future, we will provide expressions, as well as symbolic computation tools to compute these quantities in more general cases. We situate our findings in the broader effort of the community to produce new effective tools for computational biology by combining numerical and symbolic methods.

Appendix1: mean and variance of the protein

We compute here the mean and the variance of the protein synthesized by a constitutive promoter (gene 1 of models $M_{1}$ and $M_{2}$).

We start with

[TABLE]

where $s(t)=0$ if the promoter is OFF and $s(t)=1$ if the promoter is ON at the time $t$.

Eq.24 leads to

[TABLE]

From (25) it follows

[TABLE]

The promoter state variable $s(t)$ follows the master equation

[TABLE]

that has the solution

[TABLE]

where $p10={\mathbb{P}}\!\left[s(0)=1\right]$, $\epsilon=(f+h)/\rho$, and $p_{1}=f/(f+h)$. Using straightforward algebra, we find

[TABLE]

where

[TABLE]

From (25) we find also

[TABLE]

We have considered here that $x(0)$, $y(0)$ are uncorrelated, but more general expressions can be obtained.

In order to compute the two times covariance ${\mathbb{E}}\!\left[s(t_{1})s(t_{3})\right]-{\mathbb{E}}\!\left[s(t_{1})\right]{\mathbb{E}}\!\left[s(t_{3})\right]$ we combine the tower property of the conditional expectation with the Markov property satisfied by $s(t)$. More precisely, for $t_{1}\geq t_{3}$ we find ${\mathbb{E}}\!\left[s(t_{1})s(t_{3})\right]={\mathbb{E}}\!\left[{\mathbb{E}}\!\left[s(t_{1})s(t_{3})|s(t_{3}\right]\right]={\mathbb{E}}\!\left[((s(t_{3})-p_{1})e^{-\rho\epsilon(t_{1}-t_{3})}+p_{1})s(t_{3})\right]$ and ${\mathbb{E}}\!\left[s(t_{1})\right]{\mathbb{E}}\!\left[s(t_{3})\right]={\mathbb{E}}\!\left[{\mathbb{E}}\!\left[s(t_{1})|s(t_{3}\right]\right]{\mathbb{E}}\!\left[s(t_{3})\right]=(({\mathbb{E}}\!\left[s(t_{3})\right]-p_{1})e^{-\rho\epsilon(t_{1}-t_{3})}+p_{1}){\mathbb{E}}\!\left[s(t_{3})\right]$. Then, it follows

[TABLE]

$s(t_{3})$ is a Bernoulli variable, therefore $Var[s(t_{3})]={\mathbb{E}}\!\left[s(t_{3})\right](1-{\mathbb{E}}\!\left[s(t_{3})\right])$. From (35) and (28) it follows

[TABLE]

Similarly, one gets

[TABLE]

The domain of the multiple integral in (34) should be split in two sub-domains corresponding to $t_{2}<t_{4}$ and to $t_{2}>t_{4}$. Each of these sub-domains should be subdivided into two smaller sub-domains corresponding to $t_{3}>t_{1}$ and $t_{1}<t_{3}$. Symmetry arguments imply that the integrals on $t_{2}<t_{4}$ and on $t_{4}<t_{2}$ are equal, which allows us to perform the calculation of the integral on only two sub-domains, instead of four. After some algebra we find

[TABLE]

where

[TABLE]

Appendix2: details of the derivation of (15),(16)

$x(t)$ and $y(t)$ satisfy the following system of equations

[TABLE]

For simplification, we rescale variables and parameters $t\to t\rho$, $k_{i}\to k_{i}/\rho$, $a\to a/\rho$, $b\to b/\rho$ and obtain

[TABLE]

From (41) it follows $y(\tau)=y(0)e^{-\tau}+\int_{0}^{\tau}d\tau^{\prime}e^{-(\tau-\tau^{\prime})}[k_{0}+(k_{1}-k_{0})s(\tau^{\prime})]=y(0)e^{-\tau}+\sum_{j=1}{M-1}\int_{\tau_{j}}^{\tau_{j+1}}d\tau^{\prime}e^{-(\tau-\tau^{\prime})}[k_{0}+(k_{1}-k_{0})s(\tau_{j})]=y(0)e^{-\tau}+k_{0}(1-e^{-\tau})+(k_{1}-k_{0})\sum_{j=0}{M-1}e^{-\tau}(e^{\tau_{j+1}}-e^{\tau_{j}})s(\tau_{j})$ and hence (15).

From (41) we also obtain $x(\tau)=x(0)e^{-a\tau}+b\int_{0}^{\tau}e^{a(\tau^{\prime}-\tau)}y(\tau^{\prime})d\tau^{\prime}=x(0)e^{-a\tau}+\frac{by(0)}{a-1}(e^{-\tau}-e^{-a\tau})+b\int_{0}^{\tau}d\tau^{\prime}e^{a(\tau^{\prime}-\tau)}\int_{0}^{\tau^{\prime}}d\tau^{\prime\prime}e^{-(\tau^{\prime}-\tau^{\prime\prime})}(k_{0}+(k_{1}-k_{0})s(\tau^{\prime\prime}))=x(0)e^{-a\tau}+\frac{by(0)}{a-1}(e^{-\tau}-e^{-a\tau})+bk_{0}\left(\frac{1-e^{-a\tau}}{a}+\frac{e^{-a(\tau-\tau_{0})}-e^{-(\tau-\tau_{0})}}{a-1}\right)+b(k_{1}-k_{0})I,$ where

[TABLE]

In order to compute the integral $I$ we decompose the triangular integration domain into $M-1$ rectangles and $M$ triangles on each of which $s(\tau^{\prime\prime})$ is constant, as in Figure3.

The contribution of each rectangle to the integral $I$ is $\int_{\tau_{i}}^{\tau}d\tau^{\prime}e^{a(\tau^{\prime}-\tau)}\int_{\tau_{i-1}}^{\tau_{i}}d\tau^{\prime\prime}e^{-(\tau^{\prime}-\tau^{\prime\prime})}s(\tau_{i-1})=e^{-a\tau}\frac{e^{(a-1)\tau}-e^{(a-1)\tau_{i}}}{a-1}(e^{\tau_{i}}-e^{\tau_{i-1}})s(\tau_{i-1})$.

The contribution of each triangle to the integral $I$ is $\int_{\tau_{i-1}}^{\tau_{i}}d\tau^{\prime}e^{a(\tau^{\prime}-\tau)}\int_{\tau_{i-1}}^{\tau^{\prime}}d\tau^{\prime\prime}e^{-(\tau^{\prime}-\tau^{\prime\prime})}s(\tau_{i-1})=e^{-a\tau}s(\tau_{i-1})\left(\frac{e^{a\tau_{i}}-e^{a\tau_{i-1}}}{a}-e^{\tau_{i-1}}\frac{e^{(a-1)\tau_{i}}-e^{(a-1)\tau_{i-1}}}{a}\right)$.

It follows that $I=\sum_{i=1}^{M-1}e^{-a\tau}\frac{e^{(a-1)\tau}-e^{(a-1)\tau_{i}}}{a-1}(e^{\tau_{i}}-e^{\tau_{i-1}})s(\tau_{i-1})+\sum_{i=1}^{M}e^{-a\tau}s(\tau_{i-1})\left(\frac{e^{a\tau_{i}}-e^{a\tau_{i-1}}}{a}-e^{\tau_{i-1}}\frac{e^{(a-1)\tau_{i}}-e^{(a-1)\tau_{i-1}}}{a}\right).$

Noting that the first sum in the expression of $I$ can go to $i=M$ (the M-th term is zero) we obtain (16).

Appendix3: proof of the Theorem 3.1

The proof the Theorem 1 relies on the following Lemma:

Lemma 4.1

Let $(x(t),y(t),s(t))$ be a realization of the PDMP such that $s(\tau_{j})=s_{j}$ for all $j\in[0,M-1]$ and let $(x_{M}(t),y_{M}(t))$ be the push-forward solutions computed with (15),(16) considering that $s(t)$ is piecewise constant on the intervals $[\tau_{i-1},\tau_{i}]$. Then, under the conditions of Theorem 1, ${\mathbb{P}}\!\left[|y(t)-y_{M}(t)|>\epsilon\right]\to 0$ and ${\mathbb{P}}\!\left[|x(t)-x_{M}(t)|>\epsilon\right]\to 0$ when $M\to\infty$, for any $\epsilon>0$ and for all $t\in[0,\tau]$.

We prove this Lemma for $y(t)$. The proof for $x(t)$ follows the same principles.

According to the constructive definitions of PDMP (see Section 2.1), and considering that $\lambda(\vec{x}(t),\vec{y}(t),s(t))<A$, for all $t\in[0,\tau]$ (this follows from the continuity of $\lambda$ and the boundedness of $\vec{x}(t)$, and $\vec{y}(t)$) then, with probability one, $s(t)$ has a finite number of jumps inside the interval $[0,\tau]$. Consider that for the $i^{\text{th}}$ gene there are $k$ jumps such that $s_{i}$ changes from ON to OFF or from OFF to ON. The positions of these jumps are $\tau^{*}_{j_{l}}\in[\tau_{j_{l}},\tau_{j_{l}+1})$, for $l\in[1,k]$. Using (15) it follows $|y_{i}(t)-y_{i,M}(t)|\leq e^{-\rho\tau}\frac{k_{1}-k_{0}}{\rho}\sum_{l=1}{k}|e^{\rho\tau^{*}_{j_{l}}}-e^{\rho\tau_{j_{l}}}|\leq\frac{k_{1}-k_{0}}{\rho}\frac{k_{1}-k_{0}}{\rho}\sum_{l=1}{k}(\tau_{j_{l+1}}-\tau_{j_{l}})<\frac{k_{1}-k_{0}}{\rho}\frac{kC}{M}$.

For constitutive promoters the number of jumps $k$ of the promoter of the gene $i$ inside $[0,\tau]$ has a mean ${\mathbb{E}}\!\left[k\right]=\frac{\tau}{h_{i}+f_{i}}$. Slightly more complex, but finite bounds, can be obtained for a regulated gene as well.

Using Markov’s inequality we find that ${\mathbb{P}}\!\left[\frac{k_{1}-k_{0}}{\rho}\frac{kC}{M}>\epsilon\right]\leq\frac{k_{1}-k_{0}}{\rho}\frac{C\tau}{\epsilon(h_{i}+f_{i})M}$. It follows that ${\mathbb{P}}\!\left[|y(t)-y_{M}(t)|>\epsilon\right]\to 0$ when $M\to\infty$, for any $\epsilon>0$.

The proof of the Theorem 1 follows from the Lemma 4.1 because by construction, the promoter states have the same distribution in the push-forward and PDMP schemes, and the convergence in probability of the mRNAs and of the proteins implies the convergence in distribution of these variables.