Simulation of stochastic network dynamics via entropic matching

The simulation of complex stochastic network dynamics arising, for instance, from models of coupled biomolecular processes remains computationally challenging. Often, the necessity to scan a model's dynamics over a large parameter space renders full-fledged stochastic simulations impractical, motivating approximation schemes. Here we propose an approximation scheme which improves upon the standard linear noise approximation while retaining similar computational complexity. The underlying idea is to minimize, at each time step, the Kullback-Leibler divergence between the true time evolved probability distribution and a Gaussian approximation (entropic matching). This condition leads to ordinary differential equations for the mean and the covariance matrix of the Gaussian. For cases of weak nonlinearity, the method is more accurate than the linear method when both are compared to stochastic simulations.

Figure 1

Test application to a system of van der Pol oscillators in the regime of weak nonlinearity. Top panel: Stochastic simulation with 1000 runs (some individual runs are plotted) of system (19) and corresponding average value for the first species of a system with $\mu =0.05$ and $N=3$. Middle panel: Comparison of stochastic average value (solid) with prediction using Eq. (16) (dashed) and Eq. (7) (dotted). Bottom panel: Comparison of standard deviation for the stochastic simulation (solid) with prediction using Eq. (16) (dashed) and Eq. (7) (dotted).

Figure 2

Histogram of stochastic simulation with 1000 runs of system (19) with $\mu =0.05$ and $N=3$ and the corresponding Gaussian (solid line) predicted by our method (16) at the final time $t=T$ (in this case $T=10$). The first three panels show the distributions of $x_0$, $x_1$, and $x_2$, respectively. Bottom panel: histogram of chi-squared distribution for the stochastic simulation data compared to the ${\chi }^2$ distribution with three degrees of freedom.

Figure 3

Same as Fig. 1 but for moderately strong nonlinearity ($\mu =1.5$).

Figure 4

Same as Fig. 2 but for moderately strong nonlinearity ($\mu =1.5$).

Figure 5

Analogous to Fig. 2 but for the genetic circuit model (21) with $n=2$ and $N=3$.

Figure 6

Final histogram of stochastic simulation with 1000 runs of genetic circuit model (21) with $n=2$ and $N=3$.