Abstract
Probability reigns in biology, with random molecular events dictating the fate of individual organisms, and propelling populations of species through evolution. In principle, the master probability equation provides the most complete model of probabilistic behavior in biomolecular networks. In practice, master equations describing complex reaction networks have remained unsolved for over 70 years. This practical challenge is a reason why master equations, for all their potential, have not inspired biological discovery. Herein, we present a closure scheme that solves the master probability equation of networks of chemical or biochemical reactions. We cast the master equation in terms of ordinary differential equations that describe the time evolution of probability distribution moments. We postulate that a finite number of moments capture all of the necessary information, and compute the probability distribution and higher-order moments by maximizing the information entropy of the system. An accurate order closure is selected, and the dynamic evolution of molecular populations is simulated. Comparison with kinetic Monte Carlo simulations, which merely sample the probability distribution, demonstrates this closure scheme is accurate for several small reaction networks. The importance of this result notwithstanding, a most striking finding is that the steady state of stochastic reaction networks can now be readily computed in a single-step calculation, without the need to simulate the evolution of the probability distribution in time.