Simulation of bimolecular reactions: Numerical challenges with the graph Laplacian

Abstract

An important framework for modelling and simulation of chemical reactions is a Markov process sometimes known as a master equation. Explicit solutions of master equations are rare; in general the explicit solution of the governing master equation for a bimolecular reaction remains an open question. We show that a solution is possible in special cases. One method of solution is diagonalization. The crucial class of matrices that describe this family of models are non-symmetric graph Laplacians. We illustrate how standard numerical algorithms for finding eigenvalues fail for the non-symmetric graph Laplacians that arise in master equations for models of chemical kinetics. We propose a novel way to explore the pseudospectra of the non-symmetric graph Laplacians that arise in this class of applications, and illustrate our proposal by Monte Carlo. Finally, we apply the Magnus expansion, which provides a method of simulation when rates change in time. Again the graph Laplacian structure presents some unique issues: standard numerical methods of more than second-order fail to preserve positivity. We therefore propose a method that achieves fourth-order accuracy, and maintain positivity.

References

• A. Basak, E. Paquette, and O. Zeitouni. Regularization of non-normal matrices by gaussian noise–-the banded toeplitz and twisted toeplitz cases. In Forum of Mathematics, Sigma, volume 7. Cambridge University Press, 2019. doi:10.1017/fms.2018.29.
• S. Blanes, F. Casas, J. A. Oteo, and J. Ros. The magnus expansion and some of its applications. Phys. Rep., 470(5-6):151–238, 2009. doi:10.1016/j.physrep.2008.11.001.
• B. A. Earnshaw and J. P. Keener. Invariant manifolds of binomial-like nonautonomous master equations. SIAM J. Appl. Dyn. Sys., 9(2):568–588, 2010. doi10.1137/090759689.
• J. Gunawardena. A linear framework for time-scale separation in nonlinear biochemical systems. PloS One, 7(5):e36321, 2012. doi:10.1371/journal.pone.0036321.
• A. Iserles and S. MacNamara. Applications of magnus expansions and pseudospectra to markov processes. Euro. J. Appl. Math., 30(2):400–425, 2019. doi:10.1017/S0956792518000177.
• S. MacNamara. Cauchy integrals for computational solutions of master equations. ANZIAM Journal, 56:32–51, 2015. doi:10.21914/anziamj.v56i0.9345.
• S. MacNamara, A. M. Bersani, K. Burrage, and R. B. Sidje. Stochastic chemical kinetics and the total quasi-steady-state assumption: Application to the stochastic simulation algorithm and chemical master equation. J. Chem. Phys., 129:095105, 2008. doi:10.1063/1.2971036.
• S. MacNamara and K. Burrage. Stochastic modeling of naive T cell homeostasis for competing clonotypes via the master equation. SIAM Multiscale Model. Sim., 8(4):1325–1347, 2010.
• S. MacNamara, K. Burrage, and R. B. Sidje. Multiscale modeling of chemical kinetics via the master equation. SIAM Multiscale Model. and Sim., 6(4):1146–1168, 2008. doi:10.1137/060678154.
• S. MacNamara, Wi. McLean, and K. Burrage. Wider contours and adaptive contours, pages 79–98. Springer International Publishing, 2019. doi:10.1007/978-3-030-04161-8_7.
• M. J. Shon. Trapping and manipulating single molecules of DNA. PhD thesis, Harvard University, 2014. http://nrs.harvard.edu/urn-3:HUL.InstRepos:11744428.
• M. J. Shon and A. E. Cohen. Mass action at the single-molecule level. J. Am. Chem. Soc., 134(35):14618–14623, 2012. doi:10.1021/ja3062425.
• C. Timm. Random transition-rate matrices for the master equation. Phys. Rev. E, 80(2):021140, 2009. doi:10.1103/PhysRevE.80.021140.
• L. N. Trefethen and M. Embree. Spectra and pseudospectra: The behavior of nonnormal matrices and operators. Princeton University Press, 2005. https://press.princeton.edu/books/hardcover/9780691119465/spectra-and-pseudospectra.