Reduction of multiscale stochastic biochemical reaction networks using exact moment derivation

<div><p>Biochemical reaction networks (BRNs) in a cell frequently consist of reactions with disparate timescales. The stochastic simulations of such multiscale BRNs are prohibitively slow due to high computational cost for the simulations of fast reactions. One way to resolve this problem uses the fact that fast species regulated by fast reactions quickly equilibrate to their stationary distribution while slow species are unlikely to be changed. Thus, on a slow timescale, fast species can be replaced by their quasi-steady state (QSS): their stationary conditional expectation values for given slow species. As the QSS are determined solely by the state of slow species, such replacement leads to a reduced model, where fast species are eliminated. However, it is challenging to derive the QSS in the presence of nonlinear reactions. While various approximation schemes for the QSS have been developed, they often lead to considerable errors. Here, we propose two classes of multiscale BRNs which can be reduced by deriving an exact QSS rather than approximations. Specifically, if fast species constitute either a feedforward network or a complex balanced network, the reduced model based on the exact QSS can be derived. Such BRNs are frequently observed in a cell as the feedforward network is one of fundamental motifs of gene or protein regulatory networks. Furthermore, complex balanced networks also include various types of fast reversible bindings such as bindings between transcriptional factors and gene regulatory sites. The reduced models based on exact QSS, which can be calculated by the computational packages provided in this work, accurately approximate the slow scale dynamics of the original full model with much lower computational cost.</p></div

Model diagram of the feedforward network of the full model and the reduced model.

<p>In the diagram of the full model (left), red arrows indicate fast reactions. In the diagram of the reduced model (right), which consists of slow species <i>S</i>₂ and <i>S</i>₄ and slow reactions, 〈<i>S</i>₁〉 and 〈<i>S</i>₁<i>S</i>₃〉 represent conditional moments, 〈<i>X</i>₁|<i>X</i>₂, <i>X</i>₄〉 and 〈<i>X</i>₁<i>X</i>₃|<i>X</i>₂, <i>X</i>₄〉, respectively. In both of the diagrams, degradation reactions are not shown, for simplicity.</p

EMB model provides much more accurate approximation of the original feedforward network model than AMB model.

<p>(A-B) Trajectories of original full model with <i>ϵ</i> = 0.01 and the EMB model (A) and the AMB model (B). The lines and colored ranges indicate the mean of <i>X</i>₄ and standard deviations of <i>X</i>₄ from their mean, respectively. Histograms represent distributions of <i>X</i>₄ at the steady state. Here, <i>X</i>₁(0) = <i>α</i>₁/<i>β</i>₁, <i>X</i>₂(0) = <i>X</i>₃(0) = <i>X</i>₄(0) = 0. Here, 10⁴ stochastic simulations were performed. (C-D) Mean (C) and standard deviation (D) of stationary distribution (t = 8) simulated with the full model with varying <i>ϵ</i> = 0.1, 0.05, 0.01, the EMB model and the AMB model. (E-F) As <i>β</i>₁ = <i>α</i>₁ increases, the EMB model, but not the AMB model, predicts that the mean and the CV of <i>S</i>₄ decrease, which is consistent with the simulations of the original full model.</p

Model diagram of a positive feedback loop with decoys.

<p>(A) In the diagram of the full model, red arrows indicate the fast reversible binding and unbinding reactions between transcriptional factor (<i>P</i>) and either a regulatory promoter site (<i>A</i>) or one of <i>N</i> identical nonregulatory decoy binding sites (<i>D</i>). These two fast reversible bindings form a complex balanced network with species <i>P</i>, <i>G</i>₀, <i>G</i>_<i>A</i>, <i>D</i>, and <i>P</i>: <i>D</i> with conservations, and <i>X</i>_<i>D</i> + <i>X</i>_{<i>P</i>: <i>D</i>} = <i>N</i>. (B) The reduced model consists solely of a slow species, which is the total number of transcription factors <i>T</i> (). 〈<i>G</i>_<i>A</i>〉 represents stationary conditional moment, . In both of the diagrams (A and B), degradation reactions are not shown, for simplicity. (C) The exact moment , which is used for the EMB model and its approximation, which is used for the AMB model. When <i>X</i>_<i>T</i> < <i>N</i> = 10, they show a discrepancy.</p

Model diagram of a transcriptional negative feedback loop with a fast feedfoward subnetwork and a slow dimerization.

<p>In the diagram of the full model (top), red arrows indicate fast reactions. Note that the subnetwork consisting of the fast species <i>E</i>, <i>Q</i>, <i>F</i>, and <i>R</i> with fast reactions (red arrows) is feedforward. In the diagram of the reduced model (bottom), which solely consists of slow species and slow reactions, 〈<i>R</i>(<i>R</i> − 1)〉 represents a conditional moment, 〈<i>X</i>_<i>R</i>(<i>X</i>_<i>R</i> − 1)|<i>X</i>_<i>P</i>〉. In both of the diagrams, all degradation reactions are not shown, for simplicity.</p

Reactions and propensity functions in the transcriptional negative feedback loop with a dimerization (Fig 3).

<p>Reactions and propensity functions in the transcriptional negative feedback loop with a dimerization (<a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1005571#pcbi.1005571.g003" target="_blank">Fig 3</a>).</p

Reactions and propensity functions in the positive feedback loop with decoys (Fig 7A).

<p>Reactions and propensity functions in the positive feedback loop with decoys (<a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1005571#pcbi.1005571.g007" target="_blank">Fig 7A</a>).</p

Reactions and propensity functions in the genetic oscillator (Fig 5).

<p>Reactions and propensity functions in the genetic oscillator (<a href="http://www.ploscompbiol.org/article/info:doi/10.1371/journal.pcbi.1005571#pcbi.1005571.g005" target="_blank">Fig 5</a>).</p

Trajectories of the full model, the EMB model and the AMB model of the positive feedback loop with decoys.

<p>(A-B) The lines and colored ranges represent <i>E</i>(<i>X</i>_<i>T</i>) and <i>E</i>(<i>X</i>_<i>T</i>) ± <i>SD</i>(<i>X</i>_<i>T</i>)/2 of 10⁴ stochastic simulations when <i>α</i>_<i>A</i> = 10 and <i>α</i>₀ = 4. Here <i>X</i>_<i>i</i>(0) = 0. (C-D) When <i>α</i>_<i>A</i> = 20 and <i>α</i>₀ = 8, so that overall <i>X</i>_<i>T</i> is greater than the total number of decoy sites (i.e <i>N</i> = 10), the AMB model also becomes accurate.</p

The EMB model provides more accurate approximation of the transcriptional negative feedback loop model than AMB model.

<p>(A-B) Trajectories of the full model with <i>ϵ</i> = 0.01 and the EMB model (A) and the AMB model (B). The lines and colored ranges represent <i>E</i>(<i>X</i>_{<i>R</i>: <i>R</i>}) and <i>E</i>(<i>X</i>_{<i>R</i>: <i>R</i>}) ± <i>SD</i>(<i>X</i>_{<i>R</i>: <i>R</i>})/2 of 10⁴ stochastic simulations. Here <i>X</i>_<i>i</i>(0) = 0. (C-D) Mean (C) and standard deviation (D) of steady state distribution of the full model with varying <i>ϵ</i> = 0.1, 0.05, 0.01, the EMB model and the AMB model.</p