A Hybrid of the Chemical Master Equation and the Gillespie Algorithm for Efficient Stochastic Simulations of Sub-Networks

Abstract

Modeling stochastic behavior of chemical reaction networks is an important endeavor in many aspects of chemistry and systems biology. The chemical master equation (CME) and the Gillespie algorithm (GA) are the two most fundamental approaches to such modeling; however, each of them has its own limitations: the GA may require long computing times, while the CME may demand unrealistic memory storage capacity. We propose a method that combines the CME and the GA that allows one to simulate stochastically a part of a reaction network. First, a reaction network is divided into two parts. The first part is simulated via the GA, while the solution of the CME for the second part is fed into the GA in order to update its propensities. The advantage of this method is that it avoids the need to solve the CME or stochastically simulate the entire network, which makes it highly efficient. One of its drawbacks, however, is that most of the information about the second part of the network is lost in the process. Therefore, this method is most useful when only partial information about a reaction network is needed. We tested this method against the GA on two systems of interest in biology--the gene switch and the Griffith model of a genetic oscillator--and have shown it to be highly accurate. Comparing this method to four different stochastic algorithms revealed it to be at least an order of magnitude faster than the fastest among them.

Fig 1. The genetic switch: reactions and dynamics.

A) A schematic of the genetic switch. B) Dynamics of average mRNA, m, protein, n, and the three states of promoter, S₀, S₁ and S₂. The chosen reaction frequencies in inverse minutes were: α₁ = α₂ = 0.001, β₁ = β₂ = 1, r₀ = 0.1, r = 10, K = 1, k = 0.05, q = 0.01.

Fig 2. Comparison of the CME-GA with the GA—Genetic switch.

A) A superposition of 100 realizations generated by the GA. The black solid curve represents the average of 500 realizations, while the white curve is the solution of Eq (30) for $\overline{m}$. B) Probability distributions for m at t = 150min and t = 300min, showing the match between the CME-GA (asterisk) and the GA (bar) constructed from an ensemble size of 10000. C) Comparison between the CME-GA and the GA of the averages and standard deviations of m, S₀, S₁ and S₂. The ensemble size was 1000.

Fig 3. The Griffith model: reactions and dynamics.

A) A schematic of the Griffith model. B) Dynamics of average mRNA, m, protein, n, and the five states of promoter, S₀, S₁, S₂, S₃ and S₄. The chosen reaction frequencies in inverse minutes were: α₁ = α₂ = α₃ = α₄ = 0.01, β₁ = β₂ = β₃ = β₄ = 1, r = 10, K = 1, k = 0.05, q = 0.05, a = 0.1. The number of protein conformations, d was set to 10.

Fig 4. Comparison of the CME-GA with the GA—the Griffith model.

A) Graphs 1–3 show individual realizations generated by the GA. The forth graph shows a superposition of 50 realizations. B) Probability distributions for m at t = 350min and t = 500min, showing the match between the CME-GA (asterisk) and the GA (bar) constructed from an ensemble size of 10000. C) Comparison between the CME-GA and the GA of the averages and standard deviations of m, S₀, S₁, S₂, S₃ and S₄. The ensemble size was 1000.