Multi-level methods and approximating distribution functions Open Access

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.,² 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.²⁰ Deterministic approaches can also fail to accurately replicate the dynamics of systems with large molecular numbers, such as systems operating near bifurcation points.⁴ 

To include stochasticity into models of biochemical reaction networks one can consider the framework of the chemical master equation.¹⁰ 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.¹⁰ 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 Higham¹ 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.¹³ 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 method⁷ (Gillespie’s DM) and the tau-leap method.⁸ 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 system¹⁶ and a dimerisation model¹ 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 S₁, …, S_N represent N distinct populations or species that may be involved in M reaction channels, R₁, …, R_M, 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 S_i at time t as X_i(t). We define the state space vector at time t as

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

where ν_j is the stoichiometric vector, defined as the j^th column of the stoichiometric matrix,

Each reaction R_j has a propensity function a_j(X(t)) at time t, this quantity determines the rate at which the j^th reaction occurs. The propensity functions, a_j(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.⁵ 

The Kurtz representation¹² allows us to describe evolution of the state space vector as follows:

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

Gillespie’s DM⁷ 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, a_j = a_j(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/a₀log(1/r₁).

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 S_i we generate n sample paths to calculate

where $X i ( r )$ represents the copy number of S_i from the r^th 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.⁸ 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 approximation¹⁹ to Equation (4),

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, a_j = a_j(X(t)), for each reaction channel, j.

4. Generate Poisson random variates, p_j, 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.⁸ 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 Higham¹ and Giles⁶, as well as Lester et al.¹³ 

Suppose we are interested in calculating an estimate of X_i(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 τ = τ₀. If n₀ samples are generated, then we have the biased estimate,

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

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

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

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

to give

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 n_L+1 samples. We can then calculate

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

For the rest of this paper we shall choose τ_ℓ = τ₀/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.¹⁴ 

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, Z_c(0) and Z_f(0), and the stoichiometric matrix, ν. Choose a terminal time, T, base level leap size, τ₀, 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.¹³ 

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

subject to

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

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).¹ 

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 X_i at terminal time T using a finite set of M moments, {μ₁, …, μ_M}, and the method of maximising entropy¹¹ to first estimate the corresponding probability distribution function (PDF). For a distribution on the discrete state space Ω = {x₁, …, x_m} with PDF $p= p 1 , … , p m$⁠, the entropy function, ℍ(p), is defined as

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

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

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

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

for j ∈ {0, …, M}. We need to employ a numerical method to calculate λ = (λ₀, …, λ_M), the one used in this paper is adapted from work by Mohammad-Djafari.¹⁵ 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.

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

We fix the rates to be [k₁, k₂, k₃, k₄] = [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 10⁶ 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.

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.¹⁷ 

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

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 c_i, 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 P_i 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, 𝔼(X^k), 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 P₁ be the CDF generated using the maximum entropy method using m − 1 moments, and P₂ 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.

As our next example we consider gene expression for a single gene¹ (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:¹ 

We fix the rates to be [k₁, k₂, k₃, k₄, k₅] = [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.

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 c_i, where c_i = 𝔼[1{D(1) ≤ i}] for i ∈ Ω, calculated using 10⁶ 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.

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,

for i ∈ Ω. Although we refer to the collection of estimates c_i, 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. c_i ≤ c_i+1 for {i, i + 1} ∈ Ω. For two solutions to this problem see Section IV.

These statistics, c_i, can be estimated using the multi-level method for which we control the variance of each of the c_i, requiring, as before, V ar(1{X ≤ i}) < ϵ 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 c_i. 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, c_i, 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 J ⊂ I where each element j ∈ J corresponds to the estimate c_j that required the largest number of paths on some level, ℓ, of the multi-level method. Then, for each j ∈ J, 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 c_i. 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 i ∈ I 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 P₁ be the approximate CDF generated by estimating c_i for each i ∈ Ω.

6. Construct the maximal set J such that for every j ∈ J, $n ℓ j = N ℓ ˆ$ for some level ℓ. Then for every j ∈ J 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 i ∈ I 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 P₂ be the approximate CDF created by calculating c_i for each i ∈ Ω.

10. If 𝔻_P1,P2 < δ, terminate the algorithm, otherwise let P₁ = P₂ and return to step 6.

We note that c_i ∈ [0, 1] for every i ∈ Ω and so we do not scale the variance tolerance of each estimate by c_i, 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.

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.

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.

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 c_i in regions of Ω where the random variables 1{X ≤ i} 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 c_i 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 c_i 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.