Biochemical reaction networks are often modelled using discrete-state, continuous-time Markov chains. System statistics of these Markov chains usually cannot be calculated analytically and therefore estimates must be generated via simulation techniques. There is a well documented class of simulation techniques known as exact stochastic simulation algorithms, an example of which is Gillespie’s direct method. These algorithms often come with high computational costs, therefore approximate stochastic simulation algorithms such as the tau-leap method are used. However, in order to minimise the bias in the estimates generated using them, a relatively small value of tau is needed, rendering the computational costs comparable to Gillespie’s direct method. The multi-level Monte Carlo method (Anderson and Higham, Multiscale Model. Simul. 10:146–179, 2012) provides a reduction in computational costs whilst minimising or even eliminating the bias in the estimates of system statistics. This is achieved by first crudely approximating required statistics with many sample paths of low accuracy. Then correction terms are added until a required level of accuracy is reached. Recent literature has primarily focussed on implementing the multi-level method efficiently to estimate a single system statistic. However, it is clearly also of interest to be able to approximate entire probability distributions of species counts. We present two novel methods that combine known techniques for distribution reconstruction with the multi-level method. We demonstrate the potential of our methods using a number of examples.

Biochemical reaction networks describe interactions between the molecular populations of biological and physiological systems. A classical deterministic approach to modelling these networks is to use systems of ordinary or partial differential equations to describe the evolution of the concentrations of each species. However, experimental researchers, such as Elowitz et al.,2 have shown that stochasticity can be observed in a variety of biological phenomena including gene expression within the cell. Thus, deterministic modelling can neglect important characteristics of a system. For example, in systems with low molecular populations using a deterministic model to describe the molecular concentrations can fail to predict a phenomenon known as stochastic focussing in signalling networks.20 Deterministic approaches can also fail to accurately replicate the dynamics of systems with large molecular numbers, such as systems operating near bifurcation points.4 

To include stochasticity into models of biochemical reaction networks one can consider the framework of the chemical master equation.10 This assumes spatial homogeneity in the system, which is the case we shall consider in this paper. The chemical master equation describes the temporal evolution of molecular species in terms of the number of molecules present rather than their concentrations. This description provides a system of ordinary differential equations (ODEs, one for each possible state of the system) that, for simple cases involving only zeroth and first order reactions, can be solved analytically.10 For more complicated systems that include higher order reactions analytical solutions are not generally feasible so one may try and use a numerical scheme to approximate the solution to the chemical master equation.3,9,10 However, this is also often infeasible due to the high (possibly infinite) number of distinct states in the system, so the only option to explore system behaviour is the use of simulation techniques.

There is a class of simulation techniques known as exact stochastic simulation algorithms (eSSAs) which model systems precisely, in the sense that the sample paths generated using them are mathematically equivalent to the chemical master equation. However, the more complex the system, the longer these simulations take to run, and for practical purposes approximate stochastic simulation algorithms (aSSAs) that favour speed over accuracy are often used. Unfortunately, obtaining accurate summary statistics for a system using aSSAs often renders the computational costs of those algorithms on a par with those of eSSAs.

Recently proposed by Anderson and Higham1 is the discrete-state multi-level method. This method vastly reduces the computational cost associated with generating summary statistics for a given system compared with traditional eSSAs and aSSAs. It does so by generating many sample paths of poor accuracy with very little computational effort and then reduces the error of summary statistics generated using these sample paths by adding correction terms from a smaller number of sample paths of higher accuracy. The paper by Lester et al.13 provides a clear introduction to the multi-level method with the focus being on implementation. In this work, we show how to exploit the multi-level method to accurately and efficiently estimate cumulative distribution functions (CDFs) of molecular species at a terminal time, T.

In Section II we briefly describe the basic simulation techniques that will be used throughout the paper, namely Gillespie’s direct method7 (Gillespie’s DM) and the tau-leap method.8 We also introduce the multi-level method, stating how to optimally choose the number of sample paths on each level and simulate these paths with low sample variance. Section III focusses on two methods of approximating the CDF of a chosen molecular species at time T, where we use the Schlögl system16 and a dimerisation model1 as example cases. Then we combine these two approximation methods with the multi-level method. We investigate how to optimise implementation of our proposed methods, resulting in two adaptive algorithms designed to use the multi-level method efficiently without any prior knowledge of the system. We conclude, in Section IV, with a short discussion.

Let S1, …, SN represent N distinct populations or species that may be involved in M reaction channels, R1, …, RM, in a system with volume ν. We assume that the system is well stirred in order that we may ignore any spatial effects in our model. Denote the size of population Si at time t as Xi(t). We define the state space vector at time t as

X ( t ) = ( X 1 ( t ) , , X N ( t ) ) T .
(1)

We also introduce the stoichiometric or state-change matrix, ν = {νij}, where νij is the change in the copy number of Si after the reaction channel Rj fires. Thus ν is an N × M matrix. If reaction Rj fires at time t + s given no other reactions in [t, t + s), we update the state space vector as follows

X ( t + s ) = X ( t ) + ν j ,
(2)

where νj is the stoichiometric vector, defined as the jth column of the stoichiometric matrix,

ν j = ( ν 1 j , , ν N j ) T .
(3)

Each reaction Rj has a propensity function aj(X(t)) at time t, this quantity determines the rate at which the jth reaction occurs. The propensity functions, aj(X(t)), can be calculated by considering the number of ways in which a reaction can occur given the copy numbers at time t. A detailed description can be found in the paper by Erban et al.5 

The Kurtz representation12 allows us to describe evolution of the state space vector as follows:

X ( t ) = X ( 0 ) + j = 1 M Y 1 j 0 t a j ( X ( s ) ) d s ν j ,
(4)

where the Y 1 j , j = 1, …, M, are Poisson processes of unit rate.

Gillespie’s DM7 is a standard example of an eSSA for generating sample paths from (4) We include a step by step algorithm below:

  1. Initialise the copy numbers, X(0), and the stoichiometric matrix, ν. Choose a terminal time, T, and set t = 0.

  2. Calculate the propensity functions, aj = aj(X(t)), for each reaction channel, j, and a 0 = j = 1 M a j .

  3. Generate r 1 U 0 , 1 and set Δ = 1/a0log(1/r1).

  4. If t + Δ > T then terminate the algorithm, otherwise continue to step 5.

  5. Choose the next reaction to occur by finding the minimal k such that j = 1 k a j > a 0 × r 2 , where r 2 U 0 , 1 .

  6. Update the copy numbers, X(t + Δ) = X(t) + νk.

  7. Let t = t + Δ and return to step 2.

The main advantage of this algorithm is that it can be used to produce unbiased estimates of system statistics at a terminal time, T. As an example, for the mean copy number of Si we generate n sample paths to calculate

E ( X i ( T ) ) = 1 n r = 1 n X i ( r ) ( T ) ,
(5)

where X i ( r ) represents the copy number of Si from the rth sample path. The disadvantage of this eSSA is that it simulates every reaction individually, which can render it very computationally expensive. In order to reduce computation time, we now introduce a more efficient algorithm, known as the tau-leap method, which provides sample paths that approximately follow (4).

The tau-leap method is an aSSA, first proposed by Gillespie.8 It is designed to reduce computation time by temporarily fixing the propensity functions and firing several reactions during each time interval of length τ. The method is equivalent to taking a Forward Euler approximation19 to Equation (4),

X ( t + τ ) = X ( t ) + j = 1 M Y 1 j a j X t τ ν j .
(6)

A step by step algorithm is again presented below:

  1. Initialise the copy numbers, X(0), and the stoichiometric matrix, ν. Choose a terminal time, T, and a time step, τ, such that T/τ is an integer. Set t = 0.

  2. If t + τ > T then terminate the algorithm, otherwise continue to step 3.

  3. Calculate the propensity functions, aj = aj(X(t)), for each reaction channel, j.

  4. Generate Poisson random variates, pj, from Y 1 j ( a j τ ) for j = 1, …, M.

  5. Update the copy numbers X ( t + τ ) = X ( t ) + j = 1 M p j ν j .

  6. Let t = t + τ and return to step 2.

An attractive feature of the algorithm is that the accuracy of the approximation is determined by the choice of τ. The smaller the time step, τ, the greater the accuracy of the method, and indeed as τ↓0 the method is equivalent to Gillespie’s DM.8 However, often the value of τ needed to gain the required level of accuracy is small enough that the CPU time becomes comparable to that of Gillespie’s DM, and so for further computational savings we utilise the multi-level Monte Carlo method.1,6

In this section we provide a very brief introduction to the multi-level Monte Carlo method.1,6 Heuristically, the method combines the best attributes of Gillespie’s DM and the tau-leap algorithm, to produce fast and unbiased estimates of system statistics. Further information can be found in the work by Anderson and Higham1 and Giles6, as well as Lester et al.13 

Suppose we are interested in calculating an estimate of Xi(T). We start at the base level, otherwise known as level 0. The estimate at this level is calculated using the tau-leap algorithm, as discussed in the previous section, with τ = τ0. If n0 samples are generated, then we have the biased estimate,

Q 0 = E Z τ 0 = 1 n 0 r = 1 n 0 Z τ 0 ( r ) ,
(7)

where Z τ 0 ( r ) is the copy number of population i at time T from the rth sample using the tau-leap algorithm with τ = τ0. Typically, a large τ0 is chosen so that Q0 can be calculated quickly. The purpose of the other levels is to reduce the bias of this crude approximation, Q0, by adding correction terms.

The first level correction is calculated by generating two sets of n1 sample paths using the tau-leap aSSA. One set has τ = τ0 and the other τ = τ1 < τ0. From the two sets of paths, we calculate the correction term as follows:

Q 1 = E Z τ 1 Z τ 0 1 n 1 r = 1 n 1 Z τ 1 ( r ) Z τ 0 ( r ) .
(8)

This correction is then added to the initial estimate from the base level:

Q 0 + Q 1 = E Z τ 0 + E Z τ 1 Z τ 0 = E Z τ 1 .
(9)

This estimator has bias equivalent to just using the tau-leap aSSA with τ = τ1. To further decrease the bias of the estimate we introduce time steps τ < τℓ−1 < ⋯ < τ0, and then add additional correction levels,

Q = E Z τ Z τ 1 = 1 n r = 1 n Z τ ( r ) Z τ 1 ( r ) ,
(10)

to give

Q b = 0 L Q = E Z τ 0 + = 1 L E Z τ Z τ 1 = E Z τ L .
(11)

This means that the estimator ℚb has bias equivalent to that generated using the tau-leap aSSA with τ = τL. The number of levels, L, can be chosen such that a desired level of accuracy is reached. However, a final correction level may be used in order to gain an unbiased estimate, if required. Again we generate two sets of paths, one set from a tau-leap algorithm where τ = τL and one set from an eSSA such as Gillespie’s DM, using nL+1 samples. We can then calculate

Q L + 1 = E X i Z τ L = 1 n L + 1 r = 1 n L + 1 X i ( r ) T Z τ L ( r ) ,
(12)

where X i is an identical and independently distributed (i.i.d.) copy of Xi. This last correction can then be added to the previous estimate and, using the linearity of expectation, we have

E X i T = E Z τ 0 + = 1 L E Z τ Z τ 1 + E X i Z τ L = = 0 L Q + Q L + 1 .
(13)

For the rest of this paper we shall choose τ = τ0/K where K ∈ ℕ+. This means that the time steps used on consecutive levels are nested and renders the algorithm simple to implement. However, we note that it is not always the most appropriate choice: when dealing with systems that exhibit behaviours on a range of timescales a more sophisticated choice is necessary.14 

The reduction in computational cost for the multi-level method is achieved by introducing two coupling techniques. A tau-leap coupling technique for the sets of sample paths on the first L levels, and an exact coupling technique for the L + 1 s t level. The coupling is designed to create a positive correlation between coarse, Z τ 1 ( r ) , and fine, Z τ ( r ) , sample paths so that V ar Z τ Z τ 1 is reduced and fewer sample paths need to be generated to obtain an estimator with a desired accuracy. The following algorithm is the tau-leap coupling technique:

  1. Initialise the copy numbers for the coarse and fine paths, Zc(0) and Zf(0), and the stoichiometric matrix, ν. Choose a terminal time, T, base level leap size, τ0, and scaling factor, K. Set t = 0.

  2. For α = 1 to T/τℓ−1:

    • Calculate the propensity functions at the current time t for the coarse system, a j c = a j ( Z c ( t ) ) , for each reaction channel, j.

    • For β = 1 to K:

      • Calculate the propensity functions at the current time t for the fine system, a j f = a j ( Z f ( t ) ) , for each reaction channel, j.

      • Define the following three propensity functions:
        b j 1 = min ( a j c , a j f ) ;
        (14)
        b j 2 = a j c b j 1 ;
        (15)
        b j 3 = a j f b j 1 ;
        (16)
        for each reaction channel, j.
      • For j = 1, …, M generate Poisson random variates, Y j r , each with rate b j r τ for r ∈ {1, 2, 3}.

      • Update the copy numbers of the two systems as follows:
        Z c ( t + τ ) = Z c ( t ) + j = 1 M ( Y j 1 + Y j 2 ) ν j ;
        (17)
        Z f ( t + τ ) = Z f ( t ) + j = 1 M ( Y j 1 + Y j 3 ) ν j .
        (18)

        Let t = t + τ.

Further details, including those detailing the exact coupling technique, can be found in the paper by Lester et al.13 

Finally, we need an algorithm to evaluate how many paths, n, to generate on each level, ℓ, such that we minimise the expected CPU time of the method whilst controlling the variance, V = = 0 L V (the unbiased case can be considered analogously), where V is the estimator variance on the ℓth level, of the summary statistic of interest (here the first moment E X i ) to within a user prescribed tolerance, ϵ.

If a pair of paths on the ℓth level takes c units of CPU time to be generated, we wish to minimise

= 0 L n c ,
(19)

subject to

= 0 L V < ϵ .
(20)

This optimisation problem can be solved using the method of Lagrange multipliers13 and suggests that we should take

n = V / c ϵ = 0 L c V ,
(21)

where ⌈.⌉ represents the ceiling function. One means to estimate c and V involves generating an initial number of sample paths on each level (usually 100 or 1000).1 

In this section we introduce two methods to approximate CDFs of discrete state systems and then how to implement these methods efficiently alongside the multi-level method. We demonstrate the utility of our approach using two explanatory examples.

We first consider approximating the CDF of species Xi at terminal time T using a finite set of M moments, {μ1, …, μM}, and the method of maximising entropy11 to first estimate the corresponding probability distribution function (PDF). For a distribution on the discrete state space Ω = {x1, …, xm} with PDF p = p 1 , , p m , the entropy function, ℍ(p), is defined as

H ( p ) i = 1 m p i log 1 p i .
(22)

If m = ∞ we select an integer N at which to truncate the state space, so that we work with Ωt = {x1, …, xN}. The optimisation problem is now to find the PDF, p, which maximises the entropy function (22) subject to the moment constraints

μ j = i = 1 N ( x i ) j p i ,
(23)

for j ∈ {0, …, M}. Note we include μ0, and set μ0 = 1 as we require i = 1 N p i = 1 . The solution is obtained by introducing M + 1 Lagrange multipliers,18 λ = (λ0, …, λM), and the Lagrange functional

Γ = i = 1 N p i log ( 1 p i ) j = 1 M λ j ( i = 1 N ( x i ) j p i μ j ) c ( i = 1 N p i μ 0 ) ,
(24)

where c = λ0 − 1. By solving ∂Γ/∂pi = 0, we obtain the solution

p i = exp ( = 0 M λ ( x i ) ) ,
(25)

for i ∈ {1, …, N}. Now that we know the form of the solution, which has M + 1 unknowns, λ = (λ0, …, λM), we have a system of M + 1 nonlinear equations,

i = 1 N exp ( = 0 M λ ( x i ) ) ( x i ) j = μ j ,
(26)

for j ∈ {0, …, M}. We need to employ a numerical method to calculate λ = (λ0, …, λM), the one used in this paper is adapted from work by Mohammad-Djafari.15 The CDF can then be trivially calculated.

We now present possible ways of efficiently implementing this method in combination with the multi-level method in order to estimate CDFs of given species numbers. To aid us we use two example models, the Schlögl model and a dimerisation model.

1. The Schlögl model

The Schlögl model16,21 consists of a single species, A, and four possible reactions:

R 1 : k 1 A ; R 2 : A k 2 ; R 3 : 2 A k 3 3 A ; R 4 : 3 A k 4 2 A .
(27)

We fix the rates to be [k1, k2, k3, k4] = [2200, 37.5, 0.18, 0.00025] for the remainder of this paper. We chose the initial condition A(0) = 0 and terminal time T = 20. The first question we need to answer is how many moments are necessary to approximate the CDF of A to within a desired degree of accuracy. To gain insight into this problem, we generated 106 sample paths using Gillespie’s DM and calculated the first seven moments. With these moments we calculated approximate PDFs for A(20) using up to and including seven moments in the method of maximum entropy. The results are presented on the left of Figure 1.

FIG. 1.

Left: PDFs from the maximum entropy method using up to the first seven moments on the truncated state space Ωt = {0, 1, 2, …, 600}. Right: empirical PDF and the PDF generated using the maximum entropy method and the first seven moments.

FIG. 1.

Left: PDFs from the maximum entropy method using up to the first seven moments on the truncated state space Ωt = {0, 1, 2, …, 600}. Right: empirical PDF and the PDF generated using the maximum entropy method and the first seven moments.

Close modal

Figure 1 indicates that the peaks of the PDFs change significantly with the addition of the fourth, fifth and sixth moments into the entropy calculation. Thus it appears that for the Schlögl system we need at least six moments to accurately approximate the CDF of species A. The diagram on the right of Figure 1 compares the empirical PDF and the approximate PDF using seven moments. Qualitatively, the approximate PDF captures the behaviour of the empirical PDF very well. However, in order to have a quantitative measure of the accuracy of approximate PDFs generated with different numbers of moments, we use the Kolmogorov-Smirnov distance.17 

Given two cumulative distribution functions (CDFs), F(x) and G(x), defined on support Ω, the Kolmogorov-Smirnov distance is given by

D F , G sup x Ω | F ( x ) G ( x ) | .
(28)

The Kolmogorov-Smirnov distance can be interpreted as the largest vertical distance between the two CDFs. The Kolmogorov-Smirnov distance between the approximate distribution calculated using the first seven moments and the empirical distribution is 0.0037. As a comparison we define c i = E 1 { A ( 20 ) i } = P A ( 20 ) i for i ∈ Ω, where 1(⋅) is the indicator function, and calculate the 95% confidence intervals of each estimate ci, the largest of which was 0.0001. This means that the source of the error is in the maximum entropy method and is not a consequence of noise in the empirical CDF.

The multi-level method is used to generate the moments required by the maximum entropy method. In practice we do not want to rely on results from an eSSA to calculate the Kolmogorov-Smirnov distance so instead we assume that the PDFs on the left of Figure 1 converge to the empirical distribution. We can therefore choose the number of moments, i, to be the minimal i such that 𝔻Pi,Pi−1 < δ, where δ is chosen to provide the desired level of accuracy and Pi represents the CDF calculated using the first i moments in the method of maximum entropy.

The question we now need to answer is which moment should have its estimator variance controlled when choosing the number of paths on each level in the multi-level method. Recall in Section II we controlled the estimator variance of the first moment as this was the only system statistic of interest. However, now we have multiple system statistics. We choose to control the accuracy of every moment, 𝔼(Xk), by controlling its variance, V ar X k , to within ϵ E X k for some ϵ > 0. We now present an adaptive algorithm to generate an approximate CDF using the multi-level method together with the method of maximum entropy:

  1. Set the number of levels L and the number of initial moments m ≥ 2. Set the variance control parameter, ϵ, and the tolerance parameter for the Kolmogorov-Smirnov distance test, δ.

  2. Generate 1000 initial sample paths on each level using the multi-level algorithm.

  3. Calculate the number of additional paths, n j , necessary on each level, ℓ, to ensure V ar X j < ϵ E X j using equation (21) for each j ∈ {1, …, m}. Find n k = max i { 1 , , m } n i

  4. Generate n k sample paths on each level, ℓ, using the multi-level algorithm.

  5. Let P1 be the CDF generated using the maximum entropy method using m − 1 moments, and P2 the CDF using m moments. If 𝔻P1,P2 < δ terminate the algorithm, otherwise let m = m + 1 and return to step 3.

This adaptive algorithm systematically increases the number of moments used in the entropy calculation until the addition of a new moment has negligible effect on the resulting CDF, which is determined by the tolerance parameter, δ. The algorithm relies on a convergence assumption, i.e. increasing the number of moments indefinitely will yield CDFs that converge to the true CDF. We note that, however, in practice, the algorithm is limited by the robustness of the numerical scheme chosen to calculate the CDFs.

In Figure 2 we present for the Schlögl model: the CDF, P(x), generated using the adaptive algorithm with ϵ = 0.03 and δ = 0.02; and the empirical CDF, F(x), generated using Gillespie’s DM which took 4568 seconds and 38853 seconds of CPU time to compute, respectively. The adaptive algorithm selected six as the optimal number of moments.

FIG. 2.

Left: The blue curve is the CDF F(x) generated using Gillespie’s DM with 106 sample paths. The red curve is the CDF P(x) generated using the adaptive multi-level algorithm and the method of maximum entropy, with control parameters ϵ = 0.03 and δ = 0.02. The truncated state space is Ωt = {0, 1, 2, …, 600}. The Kolmogorov-Smirnov distance is 0.0125. Right: A plot of the error: |F(x) − P(x)|.

FIG. 2.

Left: The blue curve is the CDF F(x) generated using Gillespie’s DM with 106 sample paths. The red curve is the CDF P(x) generated using the adaptive multi-level algorithm and the method of maximum entropy, with control parameters ϵ = 0.03 and δ = 0.02. The truncated state space is Ωt = {0, 1, 2, …, 600}. The Kolmogorov-Smirnov distance is 0.0125. Right: A plot of the error: |F(x) − P(x)|.

Close modal

2. A dimerisation model

As our next example we consider gene expression for a single gene1 (G). It transcribes mRNA, known as messengers (M), which can, in turn, translate proteins (P). Pairs of these proteins may combine to form a dimer (D). The mRNA and proteins may also degrade. This biochemical network is described by the following set of reactions:1 

R 1 : G k 1 G + M ; R 2 : M k 2 M + P ; R 3 : P + P k 3 D ; R 4 : M k 4 ; R 5 : P k 5 .
(29)

We fix the rates to be [k1, k2, k3, k4, k5] = [25, 1000, 0.001, 0.1, 1] for the remainder of this paper. We wish to approximate the distribution of the dimer population at a fixed terminal time T = 1 and we take the initial conditions to be [G(0), M(0), P(0), D(0)]T = [1, 0, 0, 0]T.

FIG. 3.

Left: The blue curve is the CDF F(x) generated using Gillespie’s DM with 106 sample paths. The red curve is the CDF P(x) generated using the adaptive multi-level algorithm and the method of maximum entropy, with control parameters ϵ = 0.0025 and δ = 0.01. The truncated state space is Ωt = {0, 1, 2, …, 10000}. The Kolmogorov-Smirnov distance is 0.0029. Right: A plot of the error: |F(x) − P(x)|.

FIG. 3.

Left: The blue curve is the CDF F(x) generated using Gillespie’s DM with 106 sample paths. The red curve is the CDF P(x) generated using the adaptive multi-level algorithm and the method of maximum entropy, with control parameters ϵ = 0.0025 and δ = 0.01. The truncated state space is Ωt = {0, 1, 2, …, 10000}. The Kolmogorov-Smirnov distance is 0.0029. Right: A plot of the error: |F(x) − P(x)|.

Close modal

Figure 3 shows the CDF, P(x), produced using the adaptive multi-level algorithm, it took 321 seconds of CPU time to compute with parameters ϵ = 0.0025 and δ = 0.01. The figure also contains the empirical CDF, F(x), generated using Gillespie’s DM, which took 2598 seconds of CPU time to compute. The maximum 95% confidence interval for the estimates ci, where ci = 𝔼[1{D(1) ≤ i}] for i ∈ Ω, calculated using 106 realisations of Gillepsie’s DM was 0.0001. This confirms that the error lies in the maximum entropy method and is not a consequence of noise in the empirical distribution. The adaptive algorithm selected four as the optimal number of moments.

Figure 4 shows how the CPU times increase for both the Schlögl model and the dimerisation model as the variance control parameter, ϵ, is decreased. The set of ϵ parameters chosen for the Schlögl system were {0.01, 0.015, 0.02, 0.025, 0.03, 0.035, 0.04, 0.045, 0.05, 0.055} and the corresponding CPU times (in seconds) were {78978, 31794, 21975, 6973, 4568, 3565, 2359, 1583, 668}. For the dimerisation model the ϵ parameters chosen were {0.0005, 0.001, 0.0015, 0.002, 0.0025, 0.003, 0.0035, 0.004, 0.0045, 0.005} and the corresponding CPU times (in seconds) were {10542, 2414, 997, 526, 321, 216, 152, 117, 91, 76}. The CPU times for the adaptive algorithm using the method of maximum entropy for the Schlögl model are seen to jump at around ϵ = 0.02. This is because as ϵ is decreased the number of moments necessary to pass the Kolmogorov-Smirnov distance test increases, which results in a significant increase in the CPU time.

FIG. 4.

Left: CPU times for the adaptive algorithm using the method of maximum entropy for the Schlögl model where δ = 0.02. Right: CPU times for the adaptive algorithm using the method of maximum entropy for the dimerisation model where δ = 0.01. In both plots the red line shows the CPU for 106 sample paths from Gillespie’s DM.

FIG. 4.

Left: CPU times for the adaptive algorithm using the method of maximum entropy for the Schlögl model where δ = 0.02. Right: CPU times for the adaptive algorithm using the method of maximum entropy for the dimerisation model where δ = 0.01. In both plots the red line shows the CPU for 106 sample paths from Gillespie’s DM.

Close modal

We now consider another method to approximate CDFs. Let X be a distribution on a discrete state space Ω. We can approximate the CDF of X directly by estimating the following statistics,

c i = P X i = E 1 X i ,
(30)

for i ∈ Ω. Although we refer to the collection of estimates ci, for i ∈ Ω, as an approximate CDF it is important to note that this method does not always generate a true CDF because estimates are not necessarily monotonic increasing i.e. cici+1 for {i, i + 1} ∈ Ω. For two solutions to this problem see Section IV.

These statistics, ci, can be estimated using the multi-level method for which we control the variance of each of the ci, requiring, as before, V ar(1{Xi}) < ϵ for some ϵ > 0. One option is then to generate N = max i Ω { n i } paths on each level, where n i is the number of paths required on the ℓth level to control the variance of ci. However, this requires large CPU times. Instead we introduce I where I ⊂ Ω. We generate N ˆ = max i I n i on each level and then calculate the estimates, ci, for i ∈ Ω. Then we iteratively increase the size of the subset I until further increases have only negligible effects on the approximate CDFs produced. We update I by finding the maximal set JI where each element jJ corresponds to the estimate cj that required the largest number of paths on some level, ℓ, of the multi-level method. Then, for each jJ, we add new entries to I as the midpoints of j and the nearest neighbours of j that are already in I. We now present an algorithm for the adaptive method and the indicator function approach:

  1. Set the number of levels L and the variance tolerance parameter, ϵ, for all the estimates ci. Set the tolerance parameter, δ, for the Kolmogorov-Smirnov distance test. Set the initial subset I ⊂ Ω from which the number of paths to be generated will be decided.

  2. Generate 1000 initial sample paths on each level using the multi-level method.

  3. Calculate the number of additional paths, n i , necessary on each level, ℓ, to ensure V ar 1 X i < ϵ for each iI using equation (21).

  4. Generate N ˆ paths on each level, ℓ, where N ˆ is given by,
    N ˆ = max i I n i ,
    (31)
    using the multi-level method.
  5. Let P1 be the approximate CDF generated by estimating ci for each i ∈ Ω.

  6. Construct the maximal set J such that for every jJ, n j = N ˆ for some level ℓ. Then for every jJ add additional entries to I as the midpoints between j and its direct neighbours. If I cannot be updated further, i.e. there is no midpoint not already in I, then terminate the algorithm.

  7. Calculate the number of additional paths, n i , necessary on each level ℓ to ensure V ar 1 X i < ϵ , for each iI using equation (21).

  8. Generate N ˆ paths on each level, ℓ, where N ˆ is given by,
    N ˆ = max i I n i ,
    (32)
    using the multi-level method.
  9. Let P2 be the approximate CDF created by calculating ci for each i ∈ Ω.

  10. If 𝔻P1,P2 < δ, terminate the algorithm, otherwise let P1 = P2 and return to step 6.

We note that ci ∈ [0, 1] for every i ∈ Ω and so we do not scale the variance tolerance of each estimate by ci, instead we select a single variance tolerance, ϵ, for all estimates. We now use this adaptive algorithm on the two example systems we have seen previously.

1. The Schlögl model

We consider the molecular species A and the reactions (27). We truncate the state space to Ωt = {0, …, 600} and consider the initial condition A(0) = 0 and terminal time T = 20. The initial subspace used is I = {100, 200, 300, 400, 500}. In Figure 5 we present the approximate CDF P(x) generated using the adaptive multi-level algorithm with ϵ = 0.0025 and δ = 0.001 and the CDF F(x) generated using Gillespie’s DM which took 6762 seconds and 38853 seconds of CPU time, respectively. The CDFs are very similar and the exact error is given on the right of Figure 5.

2. A dimerisation model

Consider the four species G, M, P and D and the reactions (29). We wish to approximate the distribution of the dimer population at a fixed time T = 1 where we take the initial conditions to be [G(0), M(0), P(0), D(0)]T = [1, 0, 0, 0]T. We truncate the state space to Ωt = {0, …, 10000} and the initial subspace is I = {1000, 2000, …, 8000, 9000}. In Figure 6 we present the approximate CDF P(x) generated using the adaptive multi-level algorithm with ϵ = 0.002 and δ = 0.001 and the CDF F(x) generated using Gillespie’s DM which took 249 seconds and 2598 seconds of CPU time, respectively. Again the CDFs are very similar and the exact error is given on the right of Figure 6.

Figure 7 shows how the CPU times increase for both the Schlögl model and the dimerisation model as the variance control parameter, ϵ, is decreased. The set of ϵ parameters chosen for the Schlögl system were {0.001, 0.0015, 0.002, 0.0025, 0.003, 0.0035, 0.004, 0.0045, 0.005, 0.0055} and the corresponding CPU times (in seconds) were {50782, 20614, 12772, 6093, 4582, 3342, 2383, 1696, 1394, 1138}. For the dimerisation model the ϵ parameters chosen were {0.0005, 0.001, 0.0015, 0.002, 0.0025, 0.003, 0.0035, 0.004, 0.0045, 0.005} and the corresponding CPU times (in seconds) were {5521, 1036, 518, 249, 176, 106, 76, 63, 52, 38}. As ϵ is decreased we see that the CPU time of the adaptive method and the indicator function approach becomes larger than the CPU time for Gillespie’s DM.

FIG. 5.

Left: The blue curve is the CDF F(x) generated using Gillespie’s DM with 106 sample paths. The red curve is the approximate CDF P(x) generated using the adaptive multi-level algorithm for the indicator function approach, with control parameters ϵ = 0.0025 and δ = 0.001. The truncated state space is Ωt = {0, 1, 2, …, 600}. The Kolmogorov-Smirnov distance is 0.0042. Right: A plot of the error: |F(x) − P(x)|.

FIG. 5.

Left: The blue curve is the CDF F(x) generated using Gillespie’s DM with 106 sample paths. The red curve is the approximate CDF P(x) generated using the adaptive multi-level algorithm for the indicator function approach, with control parameters ϵ = 0.0025 and δ = 0.001. The truncated state space is Ωt = {0, 1, 2, …, 600}. The Kolmogorov-Smirnov distance is 0.0042. Right: A plot of the error: |F(x) − P(x)|.

Close modal
FIG. 6.

Left: The blue curve is the CDF F(x) generated using Gillespie’s DM with 106 sample paths. The red curve is the approximate CDF P(x) generated using the adaptive multi-level algorithm for the indicator function approach, with control parameters ϵ = 0.002 and δ = 0.001. The truncated state space is Ωt = {0, 1, 2, …, 10000}. The Kolmogorov-Smirnov distance is 0.0028. Right: A plot of the error: |F(x) − P(x)|.

FIG. 6.

Left: The blue curve is the CDF F(x) generated using Gillespie’s DM with 106 sample paths. The red curve is the approximate CDF P(x) generated using the adaptive multi-level algorithm for the indicator function approach, with control parameters ϵ = 0.002 and δ = 0.001. The truncated state space is Ωt = {0, 1, 2, …, 10000}. The Kolmogorov-Smirnov distance is 0.0028. Right: A plot of the error: |F(x) − P(x)|.

Close modal
FIG. 7.

Left: CPU times for the adaptive algorithm using the indicator function approach for the Schlögl model where δ = 0.001. Right: CPU times for the adaptive algorithm using the indicator function approach for the dimerisation model where δ = 0.001. In both plots the red line shows the CPU for 106 sample paths from Gillespie’s DM.

FIG. 7.

Left: CPU times for the adaptive algorithm using the indicator function approach for the Schlögl model where δ = 0.001. Right: CPU times for the adaptive algorithm using the indicator function approach for the dimerisation model where δ = 0.001. In both plots the red line shows the CPU for 106 sample paths from Gillespie’s DM.

Close modal

We have presented two novel algorithms designed to employ the multi-level Monte Carlo method to generate approximate CDFs for molecular populations at a terminal time T. The first method was the method of maximum entropy which used the multi-level method to generate a finite set of moments and from these moments calculated an approximate distribution. We note that this method is favourable if an analytical approximation of the PDF is required. The second method was an indicator function approach. This adaptive method focussed on generating more paths to control the estimates of ci in regions of Ω where the random variables 1{Xi} have a higher variance. These regions are often found to be near the peaks of the corresponding PDFs which is where the method of maximum entropy returns the largest errors. Thus we see that the indicator function approach achieves smaller Kolmogorov-Smirnov distances than the method of maximum entropy when using the Schlögl model.

However, we note that the indicator function approach does not necessarily provide a CDF as the estimates ci are not always monotonic increasing in i. Although the estimates can still be used to give an accurate and meaningful description of the distribution, if a physical CDF is required there are at least two solutions. Firstly one can vastly decrease the value of the variance control parameter ϵ to reduce the probability of non-montonicity, this however will significantly increase the CPU time of the method such that it is comparable with the tau-leap method and Gillespie’s DM. Secondly one can take estimates for ci from a coarse subset S ⊂ Ω, and generate the CDF using polynomial interpolation whilst enforcing a positive gradient over the range of Ω. This second approach uses a coarser subset S to reduce the probability of non-monotonicity without the further increase in the CPU time.

Future work will involve investigating the potential for improvements on the efficiency of the maximum entropy method, as well as improvements in the efficiency of using multi-level methods to produce approximate CDFs in general.

The authors would like to thank Chris Lester for many insightful discussions.

1.
D. F.
Anderson
and
D. J.
Higham
, “
Multi-level Monte Carlo for continuous time Markov chains with applications in biochemical kinetics
,”
SIAM Multiscale Modeling and Simulation
10
(
1
),
146
179
(
2012
).
2.
M. B.
Elowitz
,
A. J.
Levine
,
E. D.
Siggia
, and
P. S.
Swain
, “
Stochastic gene expression in a single cell
,”
Science
297
(
5584
),
1183
1186
(
2002
).
3.
S.
Engblom
, “
Spectral approximation of solutions to the chemical master equation
,”
Journal of Computational and Applied Mathematics
229
(
1
),
208:221
(
2009
).
4.
R.
Erban
,
J. S.
Chapman
,
I. G.
Kevrekidis
, and
T.
Vejchodský
, “
Analysis of a stochastic chemical system close to a sniper bifurcation of its mean-field model
,”
SIAM Journal of Applied Mathematics
70
(
3
),
984:1016
(
2009
).
5.
R.
Erban
,
J. S.
Chapman
, and
P. K.
Maini
, A practical guide to stochastic simulations of reaction-diffusion processes.arXiv, 0704(1908), 2007.
6.
M. B.
Giles
, “
Multilevel Monte Carlo path simulation
,”
Operations Research
56
(
3
),
607
617
(
2008
).
7.
D. T.
Gillespie
, “
Exact stochastic simulation of coupled chemical reactions
,”
The Journal of Physical Chemistry
81
(
25
),
2340
2361
(
1977
).
8.
D. T.
Gillespie
, “
Approximate accelerated stochastic simulation of chemically reacting systems
,”
Journal of Chemical Physics
115
(
4
),
1716
1733
(
2001
).
9.
T.
Jahnke
, “
On reduced models for the chemical master equation
,”
Multiscale Modelling and Simulation
9
(
4
),
1646
1676
(
2011
).
10.
T.
Jahnke
and
W.
Huisinga
, “
Solving the chemical master equation for monomolecular reaction systems analytically
,”
Journal of Mathematical Biology
54
(
1
),
1
26
(
2007
).
11.
E. T.
Jaynes
,
Probability Theory: The Logic of Science
(
Cambridge University Press
,
2003
).
12.
T. G.
Kurtz
and
D. F.
Anderson
, “
Continuous time Markov chain models for chemical reaction networks
,” in
Design and Analysis of Biomolecular Circuits
(
Springer
,
2011
), pp.
3
42
.
13.
C.
Lester
,
R. E.
Baker
,
M. B.
Giles
, and
C. A.
Yates
, A guide to efficient discrete-state multi-level simulation of stochastic biological systems.arXiv, 1412(4069), 2014.
14.
C.
Lester
,
C. A.
Yates
,
M. B.
Giles
, and
R. E.
Baker
, “
An adaptive multi-level simulation algorithm for stochastic biological systems
,”
The Journal of Chemical Physics
142
(
2
),
024113
(
2015
).
15.
A.
Mohammad-Djafari
, “
A Matlab program to calculate maximum entropy distributions
,”
Fundamental Theories of Physics
50
,
221
233
(
1992
).
16.
F.
Schlögl
, “
Chemical reaction models for non-equilibrium phase transitions
,”
Zeitschrift für Physik
253
(
2
),
147:161
(
1972
).
17.
D. J.
Sheskin
,
Handbook of Parametric and Non-parametric Statistical Procedures
, third ed. (
Chapman and Hall
,
2003
).
18.
J.
Stewart
,
Calculus
, sixth ed. (
Brooks/Cole
,
2009
).
19.
E.
Süli
and
D. F.
Mayers
,
An Introduction to Numerical Analysis
(
Cambridge University Press
,
2003
).
20.
T.
Székely
,
K.
Burrage
, and
R.
Erban
, “
A higher order numerical framework for stochastic simulation of chemical reaction systems
,”
BMC Systems Biology
6
(
1
),
85
(
2012
).
21.
T.
Wilhelm
, “
The smallest chemical reaction system with bistability
,”
BMC Systems Biology
3
(
1
),
90
(
2009
).