[Paper Review] matLeap: A fast adaptive Matlab-ready tau-leaping implementation suitable for Bayesian inference

Justin Feigelman, Stefan Ganscha+1 more|arXiv (Cornell University)|Aug 25, 2016

Gene Regulatory Network Analysis8 references3 citations

TL;DR

matLeap is a fast, adaptive, Matlab-compatible C++ implementation of the explicit-implicit tau-leaping algorithm that dynamically switches between stochastic simulation algorithm (SSA), explicit, and implicit tau-leaping based on system stiffness, ensuring non-negativity via critical reactions. It uniquely provides essential summary statistics—reaction counts $ r_i $ and integrated propensity functions $ G_i $—enabling efficient Bayesian inference, outperforming StochKit and SimBiology in speed (up to 350× faster) while maintaining high accuracy.

ABSTRACT

Background: Species abundance distributions in chemical reaction network models cannot usually be computed analytically. Instead, stochas- tic simulation algorithms allow sample from the the system configuration. Although many algorithms have been described, no fast implementation has been provided for τ-leaping which i) is Matlab-compatible, ii) adap- tively alternates between SSA, implicit and explicit τ-leaping, and iii) provides summary statistics necessary for Bayesian inference. Results: We provide a Matlab-compatible implementation of the adap- tive explicit-implicit τ-leaping algorithm to address the above-mentioned deficits. matLeap provides equal or substantially faster results compared to two widely used simulation packages while maintaining accuracy. Lastly, matLeap yields summary statistics of the stochastic process unavailable with other methods, which are indispensable for Bayesian inference. Conclusions: matLeap addresses shortcomings in existing Matlab-compatible stochastic simulation software, providing significant speedups and sum- mary statistics that are especially useful for researchers utilizing particle- filter based methods for Bayesian inference. Code is available for download at https://github.com/claassengroup/matLeap. Contact: [e-mail]

Motivation & Objective

• To address the lack of a fast, adaptive, and Matlab-compatible tau-leaping implementation suitable for Bayesian inference in chemical reaction networks.
• To enable efficient parameter inference by providing reaction counts $ r_i $ and integrated propensity functions $ G_i $, critical for conjugate Bayesian updating under gamma priors.
• To dynamically switch between SSA, explicit, and implicit tau-leaping based on system stiffness, improving accuracy and performance in stiff and non-stiff regimes.
• To reduce computational overhead and memory usage compared to existing tools like SimBiology, which store full trajectories and are prohibitively slow.
• To offer a minimalistic, direct-to-Matlab interface using SBML models and mex-file integration, avoiding complex configuration workflows.

Proposed method

• The algorithm adaptively selects between explicit and implicit tau-leaping based on the relative change in reaction propensities, switching to SSA when critical reactions threaten non-negativity or when leap intervals become too small.
• Non-negativity is enforced using critical reactions that trigger a switch to SSA when their waiting time becomes comparable to the next leap interval.
• The implicit tau-leaping variant uses the inverse Jacobian of reaction propensities, which can be computed symbolically in Matlab for speed or estimated numerically.
• Summary statistics $ r_i $ (number of firings) and $ G_i = \int g_i(X_s) ds $ (integrated propensity) are computed using trapezoidal approximation during tau-leaping and exact integration during SSA.
• The implementation is compiled as a mex-file using C++ with Eigen and Boost libraries, enabling direct calling from Matlab with matrix inputs for initial conditions and parameters.
• Models are loaded from SBML files using libSBML, and the interface supports batch simulations with varying initial conditions and parameters.

Experimental results

Research questions

• RQ1 Can a fast, adaptive tau-leaping algorithm be implemented in a way that is natively compatible with Matlab and suitable for Bayesian inference?
• RQ2 How does the performance and accuracy of an adaptive explicit-implicit tau-leaping method compare to established tools like StochKit and SimBiology?
• RQ3 To what extent do the summary statistics $ r_i $ and $ G_i $ improve the efficiency of Bayesian parameter inference in stochastic reaction networks?
• RQ4 Does the use of symbolic Jacobians in implicit tau-leaping provide a measurable speedup in practice?
• RQ5 Can the algorithm maintain accuracy in stiff systems, particularly when variance underestimation is a known issue in implicit tau-leaping?

Key findings

• matLeap achieves at least 6× to 350× speedup over Matlab SimBiology and up to 100× speedup over StochKit in benchmark models, with runtimes as low as 0.0042 seconds per simulation on the stiff dimerizing model.
• The Kullback-Leibler divergence between matLeap and SSA is comparable to StochKit (mean $ \approx 0.0037 $), with only minor variance underestimation in the dimerizing model, a known limitation of implicit tau-leaping.
• Using symbolic Jacobians provides a ~14% speedup on the dimerizing model, demonstrating the benefit of symbolic pre-computation in high-performance settings.
• matLeap uniquely provides the summary statistics $ r_i $ and $ G_i $, which are essential for conjugate Bayesian updating under gamma priors, enabling efficient posterior computation.
• The adaptive switching mechanism prevents numerical instability and maintains accuracy by detecting stiffness and switching to SSA when necessary.
• The minimalistic mex-file interface allows direct integration into Matlab workflows without configuration overhead, unlike SimBiology or StochKit, which require complex setup or file I/O.