Stochastic Representations of Ion Channel Kinetics and Exact Stochastic Simulation of Neuronal Dynamics

In this paper we provide two representations for stochastic ion channel kinetics, and compare the performance of exact simulation with a commonly used numerical approximation strategy. The first representation we present is a random time change representation, popularized by Thomas Kurtz, with the second being analogous to a "Gillespie" representation. Exact stochastic algorithms are provided for the different representations, which are preferable to either (a) fixed time step or (b) piecewise constant propensity algorithms, which still appear in the literature. As examples, we provide versions of the exact algorithms for the Morris-Lecar conductance based model, and detail the error induced, both in a weak and a strong sense, by the use of approximate algorithms on this model. We include ready-to-use implementations of the random time change algorithm in both XPP and Matlab. Finally, through the consideration of parametric sensitivity analysis, we show how the representations presented here are useful in the development of further computational methods. The general representations and simulation strategies provided here are known in other parts of the sciences, but less so in the present setting.Comment: 39 pages, 6 figures, appendix with XPP and Matlab cod

Positivity-preserving methods for ordinary differential equations

[EN] Many important applications are modelled by differential equations with positive solutions. However, it remains an outstanding open problem to develop numerical methods that are both (i) of a high order of accuracy and (ii) capable of preserving positivity. It is known that the two main families of numerical methods, Runge-Kutta methods and multistep methods, face an order barrier. If they preserve positivity, then they are constrained to low accuracy: they cannot be better than first order. We propose novel methods that overcome this barrier: second order methods that preserve positivity unconditionally and a third order method that preserves positivity under very mild conditions. Our methods apply to a large class of differential equations that have a special graph Laplacian structure, which we elucidate. The equations need be neither linear nor autonomous and the graph Laplacian need not be symmetric. This algebraic structure arises naturally in many important applications where positivity is required. We showcase our new methods on applications where standard high order methods fail to preserve positivity, including infectious diseases, Markov processes, master equations and chemical reactions.The authors thank the Isaac Newton Institute for Mathematical Sciences for support and hospitality during the programme "Geometry, compatibility and structure preservation in computational differential equations" when work on this paper was undertaken. This work was supported by EPSRC grant EP/R014604/1. S.B. has been supported by project PID2019-104927GB-C21 (AEI/FEDER, UE).Blanes Zamora, S.; Iserles, A.; Macnamara, S. (2022). Positivity-preserving methods for ordinary differential equations. ESAIM Mathematical Modelling and Numerical Analysis. 56(6):1843-1870. https://doi.org/10.1051/m2an/20220421843187056

Simulation of bimolecular reactions: Numerical challenges with the graph Laplacian

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