A General Path–Integral Monte Carlo Method for Exact Simulations of Chemical Reaction Networks

Chemical Reaction Networks (CRNs) are stochastic many-body systems used to model real-world chemical systems through a differential Master Equation (ME); analytical solutions to these equations are only known for the simplest systems. In this paper, we construct a path-integral inspirited framework for studying CRNs. Under this scheme, the time-evolution of a reaction network can be encoded in a Hamiltonian-like operator. This operator yields a probability distribution which can be sampled, using Monte Carlo Methods, to generate exact numerical simulations of a reaction network. We recover the grand probability function used in the Gillespie Algorithm as an approximation to our probability distribution, which motivates the addition of a leapfrog correction step. To assess the utility of our method in forecasting real-world phenomena, and to contrast it with the Gillespie Algorithm, we simulated a COVID-19 epidemiological model using parameters from the United States for the Original Strain and the Alpha, Delta and Omicron Variants. By comparing the results of these simulations with official data, we found that our model closely agrees with the measured population dynamics, and given the generality of this framework it can also be applied to study the spread dynamics of other contagious diseases.