Challenges in simulation of Chemical Reaction Networks

The study of deterministic models of Chemical Reaction Networks has been a central topic in Chemical Physics since the 1970s thanks to the work of M. Feinberg, F. Horn, and R. Jackson. The stochastic counterpart was formulated in the same period and investigated mostly through simulations. Stochastic simulation algorithms have been developed, the most known one being the exact simulation method usually associated to the name of Gillespie who made it popular in the field of Chemical Physics. A mathematical theory of stochastic models of Chemical Reaction Networks has been developed more recently, but efficient simulation methods are still needed to address most models that remains intractable with analytical tools. The interest in studying complicated models and the availability of an increased computer power made it possible to simulate larger examples and it stimulated the development of approximate simulation methods to speed up the computations further. The most known approximated simulation algorithm is the so-called tau-leaping method, again due to Gillespie. Recently, several authors proposed improvements to tau-leaping algorithm, including the MidPoint corrections and Post-Leap Checks both introduced by D. Anderson and collaborators on which we focus. In this thesis, we develop a few non-trivial numerical examples where the application of such improved tau-leaping methods are compared among each other and against exact simulations to evaluate their performance. The first setting is that of the stochastic Lotka-Volterra predator-prey model, in a parameter regime where the extinction of one of the populations is a rare event, but still possible when the initial population consistency is not extremely large. This is an interesting classical case where we expect that the trivial tau-leaping method can easily fail, while the improved ones should show their capabilities. We illustrate our findings showing that both improvements (MidPoint corrections and Post-Leap Checks) are important and that they can be successfully combined. In the second setting we simulate a stochastic model of nanoparticle growth that has been subject of recent studies. In this case the simulation algorithm needs a careful implementation due to several source of complexity, both in terms of speed and memory usage, and in terms of the peculiarity of the model. In this case MidPoint corrections are not applicable, while Post-Leap Checks are strictly necessary for the tau-leaping method to be applicable. However, in some interesting parameter regimes, the computational advantage provided by the approximated tau-leaping method with respect to the exact Gillespie algorithm partly fades away due to the high rate of failure of the Post-Leap checks.