Simulation of bimolecular reactions: Numerical challenges with the graph Laplacian

Abstract

[EN] 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.