Multi-level Monte Carlo for continuous time Markov chains, with applications in biochemical kinetics

1 Introduction

• Exact Monte Carlo costs O(Nϵ^-2) for O(ϵ) accuracy, where N is the order of the work per exact path.Exact paths are costly when many jump events must be simulated.Source: §1 · pp. 2–3
• For a first-order approximate method, O(ϵ^-3) work achieves O(ϵ) accuracy; a second-order method lowers the cost to O(ϵ^-2.5).The approximate estimator is biased, and its step size must balance discretization error and sampling error.Source: §1 · p. 3
• The paper extends multilevel Monte Carlo to continuous-time Markov chains using an easily simulated coupling that gives the estimator small variance.The authors identify the coupling as a central contribution.Source: §1 · p. 3
• A second multilevel algorithm uses the random time change representation to produce an unbiased estimator at significantly lower cost than exact simulation with crude Monte Carlo.The authors propose it for expected values in continuous-time Markov chains, particularly chemical reaction networks.Source: §1 · pp. 1, 4
• The complexity analysis accounts for system scaling without taking infinite-system-size or vanishing-time-step limits.The methods are presented as generic, with possible further tuning using system structure or analytical results.Source: §1 · p. 4

2 The basic stochastic model for chemical reaction networks

• A reaction network changes its state by a reaction vector ζk whenever reaction k occurs, and its intensity λk is the associated propensity.The count of each reaction channel is modeled through a Poisson-process random time change.Source: §2 · p. 4
• The random time change representation uses independent unit-rate Poisson processes and is equivalent to the chemical master equation representation.The state process is a continuous-time Markov chain.Source: §2 · pp. 4–5
• The paper adopts mass-action propensities, but its core ideas do not depend on that choice; the assumption is for analytical convenience and historical consistency.Under mass action, the propensity depends on the required molecule counts for a reaction.Source: §2 · p. 5
• For the example network, the generator combines reaction-state changes weighted by their propensities, including a quadratic propensity for the third reaction.The example specifies the reaction vectors and generator terms.Source: §2 · p. 5

3 Scaled models

• Scaling is used to account for natural model magnitudes when comparing computational costs, while simulations remain on the unscaled model.The scaled and unscaled models are equivalent representations, and scaling may inform approximate-method step size.Source: §3 · pp. 5–7
• The model is normalized so each abundance and propensity is order one, with the reaction-vector and propensity magnitudes described by scaling exponents.The exact-path work scale is defined from the model's parameters and reaction rates.Source: §3 · p. 6
• The time-scale parameter γ relates to the model's natural time scale; the analysis is most applicable for γ ≤ 0 because error bounds grow rapidly otherwise.The examples can still perform well for γ > 0, so the analysis does not fully capture method behavior there.Source: §3 · p. 6
• For mass-action networks, species abundances and rate constants are rescaled using N and exponents αi and βk to obtain normalized deterministic mass-action propensities.The rate-constant exponents βk are chosen so the rescaled constants are O(1).Source: §3 · pp. 6–7
• The scaled propensities are assumed globally Lipschitz, and classical scaling is the special case ck ≡ ρk ≡ 1 and γ = 0.The paper notes that the scaled propensities and their derivatives are better behaved in magnitude than those of the original model.Source: §3 · pp. 7–8

4 A review of Euler’s method in the current setting

• Tau-leaping holds each reaction intensity fixed over a time step and computes the number of reactions in that interval using Poisson samples.The interval is advanced by h, shortened at the terminal time T, and the process is updated before the next step.Source: §4 · p. 8
• Existing modifications include adaptive step-size selection and safeguards against negative population values.The paper uses h for the time step to avoid confusion with stopping-time notation.Source: §4 · pp. 8–9
• The path-wise representation treats the tau-leaped process as constant between step times, linking the discretized process to its random time change form.The scaled and unscaled tau-leaping models represent normalized abundances and counts, respectively.Source: §4 · pp. 8–9

• The estimator uses a telescoping sum across levels, estimating the coarsest expectation and successive differences between adjacent step sizes.Source: §5 · p. 9
• Adjacent-level paths share the same randomness, making the difference f(Zℓ(T)) − f(Zℓ−1(T)) small in variance under the stated strong-error and Lipschitz conditions.The variance is assumed to be O(hℓ) when the method has strong order 1/2 and f is Lipschitz.Source: §5 · pp. 9–10
• MLMC achieves O(ϵ) accuracy with total computational complexity O(ϵ−2(log ϵ)2), compared with O(ϵ−3) for crude Monte Carlo.The finest step size is chosen as hL = O(ϵ).Source: §5 · pp. 9–10
• Coordinated nested simulations reduce cost because coarse paths are cheaper and small adjacent-level differences require fewer expensive fine paths.Source: §5 · p. 10
• The estimator targets E[f(ZL(T))], not E[f(X(T))], so the finest level is selected to keep the resulting bias within the desired tolerance.The analysis also uses both weak accuracy to choose hL and strong accuracy to control level-difference variance.Source: §5 · p. 10

6 Multi-level Monte Carlo for continuous time Markov chains

• The method couples fine and coarse tau-leap paths through shared Poisson processes, allowing simultaneous simulation with the usual marginal distributions and work comparable to a single Euler tau-leap path.The coupling gives paths that jump together when their reaction intensities overlap; the algorithm requires at most two Poisson random variables per reaction channel per step.Source: §6 · pp. 11–13
• The analytical bounds have constants uniform in N when gamma <= 0, including the classical chemical-kinetics scaling gamma = 0; performance is also reported for gamma > 0.For gamma > 0, the bounds’ constants could be large and of limited practical use, although the reported performance motivates further analysis.Source: §6 · pp. 13, 15
• For most nontrivial examples, biased MLMC complexity is dominated by epsilon^-2 h_l0^-1, with logarithmic terms potentially reducible through optimized sample counts.The usual exact-algorithm plus crude Monte Carlo approach has complexity O(N epsilon^-2).Source: §6 · pp. 10, 16

6.1 An unbiased MLMC

• The unbiased MLMC estimator couples the exact process X^N with the finest tau-leap approximation Z_L and is unbiased for the target expectation.Because weak error does not enter its analysis, h_L can instead be chosen small enough to control the exact–approximate variance, using preliminary simulations if needed.Source: §6.1 · pp. 16–17, §6.2 · p. 18

6.2 Some observations

• For the unbiased estimator, non-negativity enforcement is not essential if intensities are defined reasonably on all of Z^d and the required bounds hold, but poor definitions can increase sampling cost or cause instability.The section notes that implementation can help keep approximate paths stable with high probability.Source: §6.2 · p. 18
• The unbiased estimator’s complexity bound adds O(N A(h_L) epsilon^-2) relative to biased MLMC, which the authors judge unlikely to outweigh its benefit unless A(h_L) is substantial.The added term is the stated difference between the unbiased and biased bounds.Source: §6.2 · p. 18
• For γ ≤ 0, unbiased MLMC is recommended over exact simulation with crude Monte Carlo, and the analysis reports substantial gains for MLMC coupled via representations (18) and (22).Under classical scaling (ρ = 1, γ = 0, c_k ≡ 1), the passage assumes h_L > 1/N; it gives A(h_L) = O(h_L^2), but the supplied text truncates before stating the predicted speedup order.Source: §6.2 · pp. 18–19
• The unbiased method shifts computational burden from the exact process to an approximate process run with a crude time step, as an example in Section 9 demonstrates.That example uses a system not satisfying classical scaling.Source: §6.2 · p. 18
• Even when one approximate path costs only hundreds or thousands of steps, unbiased MLMC can substantially reduce complexity when high accuracy makes ϵ^-2 very large.The passage notes that these moderate-N cases typically have γ ≤ 0, so the analysis applies.Source: §6.2 · p. 19
• Further analysis is needed for γ > 0, which the authors identify as a focus for future work.Source: §6.2 · p. 19

7 Delayed proofs of Theorems 1 and 2

• Lemma 3’s Nγ exponent gives a worst-case bound that is sharp for the network S1 → 2S1, which has exponential growth.The paper identifies classifying networks where this upper bound can be substantially decreased as future work.Source: §7 · p. 20
• The proof of Theorem 2 bounds the martingale terms using quadratic covariation and the Burkholder-Davis-Gundy inequality, then applies Gronwall’s inequality.Source: §7 · pp. 20–22
• Theorem 1’s proof can follow Theorem 2’s argument or use a cruder triangle-inequality bound; for γ ≤ 0, its C1(NγT) and C2(NγT) terms are uniformly bounded in N.Source: §7 · p. 22

8 Implementation issues

• The analysis gives only an order of magnitude for the path count at each level, not a precise or optimal choice, which depends on f and the model.Source: §8 · p. 22
• For a given tolerance, preliminary simulations estimate each level’s cost–variance constant and set target variances to minimize estimated computation time.Each level is then simulated until its estimator variance falls below its target; this resembles the strategy in.Source: §8 · pp. 22–23
• Optimizing pre-computation selects the finest and crudest levels before the full run, reducing the probability that approximate processes become negative or otherwise non-physical.The authors say this helps keep the approximate processes stable and increases efficiency.Source: §8 · p. 23
• If pre-computation finds no level and path-count choice significantly faster than exact simulation with crude Monte Carlo, the method reverts to the exact algorithm.The authors conclude it will therefore not be appreciably slower than that exact-algorithm baseline.Source: §8 · p. 23
• The examples generate approximate paths in batches to reduce separate calls to Matlab’s built-in Poisson random number generator.Batch sizes range from hundreds to tens of thousands.Source: §8 · p. 23

Gene transcription and translation

• The exact crude Monte Carlo baseline for the dimer mean required approximately five million paths, estimated at 155,000 CPU seconds using the modified next reaction method.The target was accuracy ±1 with 95% confidence for the expected dimer count at time 1.Source: §9 · pp. 23–24
• Biased MLMC reduced update counts 56-fold versus tau-leaping at hL = 3^-7, while its approximation bias remained apparent.Reported CPU times include the embedded optimization used to choose paths per level.Source: §9 · p. 24, Table 3 · p. 25
• Unbiased MLMC had no bias for any tested M or L and used approximately 80 times fewer updates and 140 times less CPU time than exact crude Monte Carlo.The speedups depend on the choice of M and L.Source: §9 · p. 25, Table 4 · p. 25
• For M = 3, L = 5, and ℓ0 = 2, unbiased MLMC needed 3,900 finest-level exact paths versus nearly five million for exact crude Monte Carlo.It used more than eight million paths at the coarsest level, where paths were cheap; most CPU time was spent there.Source: §9 · p. 25, Table 5 · p. 26

Parameter-dependent example

• As θ increased, unbiased MLMC progressively lost its advantage, yet remained more efficient through θ = 100; at θ = 1,000, it was slightly slower than exact crude Monte Carlo.For this example, the reported speedup followed approximately 54.6 θ^-0.62.Source: §9 · p. 27, Table 7 · p. 27, Figure 1 · p. 28

Viral kinetics

• In viral kinetics, a coupling-based control variate without multilevel discretization used 60 times fewer updates and 22 times less CPU time than crude Monte Carlo.The approximate process uses a Poisson equilibrium mean for the protein count when viral templates are present, and remains a valid continuous-time Markov chain.Source: §9 · pp. 29–30

10 Conclusions

• The new multilevel algorithm dramatically outperforms existing methods across a wide range of scaling regimes, and numerical results confirm the predicted practical benefits.Its two conceptual advances are using coupling previously introduced as a theoretical tool and exploiting exact sampling to construct an unbiased estimator.Source: §10 · p. 31
• For fixed systems, tau-leaping becomes infinitely more expensive than exact sampling as h → 0, motivating analysis that incorporates system scaling without asymptotic limits.The paper frames Monte Carlo expectation estimation as a tradeoff between accuracy and the cost of each sample.Source: §10 · pp. 30–31
• Future work includes quasi-Monte Carlo sampling and adapting the method to multiscale or hybrid models with fast/slow reactions or coupled diffusion and ODE models.These models may offer special structure to exploit.Source: §10 · p. 31
• The authors propose extending the complexity analysis to γ > 0 and using coupling without discretization for a structure-aware control variate.The latter approach is illustrated in the paper's third example.Source: §10 · p. 31