Universally valid reduction of multiscale stochastic biochemical systems using simple non-elementary propensities

Abstract

Biochemical systems consist of numerous elementary reactions governed by the law of mass action. However, experimentally characterizing all the elementary reactions is nearly impossible. Thus, over a century, their deterministic models that typically contain rapid reversible bindings have been simplified with non-elementary reaction functions (e.g., Michaelis-Menten and Morrison equations). Although the non-elementary reaction functions are derived by applying the quasi-steady-state approximation (QSSA) to deterministic systems, they have also been widely used to derive propensities for stochastic simulations due to computational efficiency and simplicity. However, the validity condition for this heuristic approach has not been identified even for the reversible binding between molecules, such as protein-DNA, enzyme-substrate, and receptor-ligand, which is the basis for living cells. Here, we find that the non-elementary propensities based on the deterministic total QSSA can accurately capture the stochastic dynamics of the reversible binding in general. However, serious errors occur when reactant molecules with similar levels tightly bind, unlike deterministic systems. In that case, the non-elementary propensities distort the stochastic dynamics of a bistable switch in the cell cycle and an oscillator in the circadian clock. Accordingly, we derive alternative non-elementary propensities with the stochastic low-state QSSA, developed in this study. This provides a universally valid framework for simplifying multiscale stochastic biochemical systems with rapid reversible bindings, critical for efficient stochastic simulations of cell signaling and gene regulation. To facilitate the framework, we provide a user-friendly open-source computational package, ASSISTER, that automatically performs the present framework.

Introduction

To understand the complex dynamics of numerous molecular interactions in living cells, quantitative analysis using mathematical models is essential [1]. While elementary reactions in living cells can be modeled by the law of mass action, characterizing all their kinetics is challenging. Thus, over a century, the combined effect of a set of elementary reactions such as rapid reversible bindings has been described with non-elementary reaction functions (e.g., Michaelis-Menten and Morrison equations) to simplify deterministic models [2–7]. Since the early 2000s, these deterministically driven non-elementary reaction functions have also been widely used to derive propensity functions for stochastic simulations, which greatly reduces the computational cost [8–33]. This heuristic approach for efficient stochastic simulations was believed to be valid as long as the non-elementary reaction functions are accurate in the deterministic sense. However, unfortunately, this was not the case [33–40]. The reason for the discrepancy between the deterministic and stochastic simulations has been recently identified for some cases [37–40], but not for all [41]. Currently, guidelines for this popular but heuristic method for efficient stochastic simulations with non-elementary propensity functions are absent.

The non-elementary reaction functions are mainly the result of the reduction of deterministic models with the following reversible binding reactions:(1)The reversible binding between molecules, such as enzyme-substrate, receptor-ligand, and protein-DNA, is the first step for nearly all biological functions of living cells [42]. However, rather than the reversible binding itself, its outcome is usually our major interest. For instance, rather than the binding between a transcription factor and DNA, we are interested in its outcome, the transcription. Furthermore, the transcription factor binding to DNA takes at most one second while transcription takes about 30 minutes in a mammalian gene [43], which causes stiffness in numerical simulations [44].

Fortunately, such rapid reversible binding reactions can be eliminated from models using the property: the levels of the species (A, B, and C) regulated by the reversible binding more quickly equilibriate to their quasi-steady-states (QSSs) compared with the total levels of the bound and unbound species, which are not affected by the reversible binding. In deterministic models, their quasi-steady-state approximations (QSSAs), which are non-elementary reaction functions, can be obtained by finding the steady-state solution of the associated differential equation in terms of the total variables. Because the QSSAs are determined by the total variables, they are known as the “total” QSSA (tQSSA). After replacing the variables that represent the levels of A, B, and C with their tQSSAs, rapid reversible bindings have been successfully eliminated from various deterministic models describing enzyme catalysis, gene regulation, and cell cycle regulation [5–7, 23, 45–51]. Note that adopting the total variables leads to timescale separation among variables, while the sole rapidity of the reversible binding reactions does not guarantee timescale separation between the original variables, A, B, and C [7] (see Discussion for details).

In stochastic models, the QSSAs for the numbers of A, B, and C are their stationary average numbers (i.e., the first moment) conditioned on the total numbers of the bound and unbound species for uni- or bi-molecular reactions [27–30] (see S1 Appendix for details). These stochastic QSSAs can be obtained by finding the steady-state solution of the chemical master equation (CME). However, unlike the deterministic tQSSA, the stochastic QSSA has a complex form (Eq (4)), which does not provide any intuition, and importantly, increases computational cost. Thus, its approximation has been derived with the deterministic tQSSA. This approximation, often referred to as the stochastic tQSSA (stQSSA) [7, 31, 38, 39], leads to non-elementary propensity functions for stochastic simulations using the Gillespie algorithm [52]. In this way, the stochastic dynamics of various systems have been accurately captured with low computational cost [7, 30–33, 38, 39, 53, 54]. However, a recent study reported that the stQSSA can be inaccurate [41], which raises the question of validity conditions for the stQSSA.

Here, we identify the complete validity condition for using the stQSSA to simplify stochastic models containing rapid reversible bindings. Specifically, we find that the stQSSA is accurate for a wide range of conditions. However, when two species whose molar ratio is ∼1:1 tightly bind, the stQSSA highly overestimates the number of unbound species. In this case, using the stQSSA to simplify stochastic models distorts the stochastic dynamics of the transcriptional repression, the transcriptional negative feedback loop of the circadian clock, and the bistable switch for mitosis. Importantly, by using the fact that the number of the unbound species is low due to the tight binding when the stQSSA is inaccurate, we develop an alternative approach, stochastic “low-state” QSSA (slQSSA). In this way, when reversible bindings are tight and not tight, slQSSA and stQSSA can be used, respectively, which enables one to obtain accurately reduced stochastic models for any case. This proposes a complete and straightforward strategy for efficiently simulating multiscale stochastic biochemical systems containing the fundamental elementary reaction, i.e., rapid reversible binding. To facilitate this framework, we provide a user-friendly open-source computational package, ASSISTER (Adaptive Simplification of StochastIc SystEm with Reversible binding).

Results

stQSSA can overestimate the number of the unbound species

In the reversible binding reaction (Eq (1)), the concentration of A, denoted by , is governed by the following ODE:(2)where and are the total concentrations of the bound and unbound species. By solving in terms of and , the tQSSA for can be obtained as follows:(3)where the is the dissociation constant. Note that if the reversible binding (Eq (1)) is embedded in a larger system, there could be other reactions affecting the dynamics of and thus additional terms in Eq (2). However, as long as the reversible binding is fast (i.e., kf and kb are much larger than the other reaction rates), is still an accurate tQSSA for . Similarly, by solving and , the tQSSAs for and can be obtained. These tQSSAs, also known as the Morrison equations [6], are generally valid, unlike the Michaelis-Menten type equations which are valid only when the enzyme concentration is negligible [7, 47, 48, 50]. Thus, the tQSSAs have been used to simplify models containing not only interactions between metabolites but also proteins whose concentrations are typically comparable [7].

Unlike the deterministic QSSA (Eq (3)), the stochastic QSSA, which is the stationary average number conditioned on the total numbers of the bound and unbound species, has a complex form [41, 55, 56]. For instance, the stochastic QSSA for the number of A (⟨A⟩) can be expressed in terms of the dimensionless variables and parameters, , where Ω is the volume of a system (e.g., , ) as follows (see Methods for details):(4)where A0 = max{0, AT − BT}. This complex form of the stochastic QSSA does not provide any intuition and importantly increases computational cost. Thus, as an alternative to the stochastic QSSA, its approximation, the stQSSA was derived with the deterministic tQSSA [7, 22–26, 31]. Specifically, the stQSSA for A (Atq) can be derived from (Eq (3)) after replacing the concentration-based variables and parameters () with dimensionless variables and parameters (X) as follows:(5)Similarly, the stQSSA for B and C (Btq and Ctq) can be obtained from their deterministic tQSSAs as follows:(6)

To identify the validity conditions for these stQSSAs, we calculated the relative error (, X = A, B, C) of the stQSSA (Xtq) to the stochastic QSSA (⟨X⟩) (Fig 1a–1c). The errors are nearly zero in most of the parameter regions, which explains why various stochastic models reduced with the stQSSA have been accurate in most previous studies [7, 30–33, 38, 39, 53, 54]. However, the relative errors of the unbound species (RA and RB) are high when AT ≈ BT. Specifically, the relative error of the bound species (RC) is at most ∼0.2 but that of the unbound species (RA, RB) can be ∼100.

[Figure omitted. See PDF.]

Fig 1. stQSSA overestimates the number of the unbound species when their molar ratio is ∼1:1 and binding is tight.

(a-c) Heat maps of the relative errors () of the stQSSA (Xtq) to the stochastic QSSA (⟨X⟩) for X = A, B, C in the reversible binding reaction (Eq (1)). Color in the heat maps represents the maximum value of RX calculated by varying Kd from 10−4 to 102 for each total number of the bound and unbound species (AT = A + C and BT = B + C). RA and RB can be extremely large when AT ≈ BT while RC is always small. (d) RA calculated over BT/AT between 0 and 2 (gray arrow in a) for three fixed Kd values (10−4, 10−3 and 10−2). RA becomes larger as BT/AT is similar to 1 and the Kd becomes smaller (i.e., the binding becomes tighter). (e) RA mainly depends on the relative sensitivity of Atq (i.e., 2SA), which can be derived in a simple form, unlike RA (Eq (8)). The maximum value of 2SA is given by , which is achieved when BT/AT is similar to 1 as in the case of RA. (f) A trajectory (left) and the stationary probability distribution (right) of A for a parameter set where RA is large (green triangle in d, AT = BT = 100, kf/Ω = 104 s−1, kb = 1s−1), simulated using the Gillespie algorithm. Since AT = BT and A binds with B tightly, A = 0 (i.e., every A is bound) most of the time, and it rarely becomes 1 by the weak unbinding reaction (solid arrow) and immediately comes back to 0 by the strong binding reaction (dotted arrow). As a result, when Kd = 10−4, the probability that A = 1 is ∼0.01, but the stQSSA for A overestimates it as ∼0.1, which is 10 times larger (i.e., a 10-fold error).

https://doi.org/10.1371/journal.pcbi.1008952.g001

To investigate why RA is high when AT ≈ BT, we derived the exact upper and lower bounds for RA (see Methods for details):(7)where FA is the Fano factor of A (i.e., ), and SA is the relative sensitivity of Atq with respect to BT (i.e., ). Furthermore, we proved that the Fano factor (FA) is less than 1 (i.e., A has a sub-Poissonian stationary distribution; see S1 Appendix for details). Therefore, RA, especially its upper bound, mainly depends on SA (Fig 1d and 1e, and A in S1 Appendix) whose formula can be derived in the following simple form, unlike RA:(8)Because SA attains the maximum value at BT = AT − Kd, SA has a large maximum value when Kd ≪ 1 at AT = BT + Kd ≈ BT. This explains why RA, whose upper bound is mainly determined by 2SA, is large when the binding is tight (Kd ≪ 1) and the total numbers of the bound and unbound species are similar (AT ≈ BT) (Fig 1d). In this case, the majority of A is bound with B, and thus A = 0 most of the time (Fig 1f left). That is, A rarely becomes 1 by the weak unbinding reaction and then immediately A becomes 0 by the strong binding reaction. As a result, the probability that A = 1 is approximately 1% (i.e., ⟨A⟩ ≈ 0.01), but the stQSSA for A (Atq) overestimates it as 10%, which is 10 times larger (Fig 1f right). Since A and B are symmetric, the above analysis can be applied to B, analogously.

stQSSA can overestimate the transcriptional activity

We found that the stQSSA for the number of the unbound species is inaccurate if their molar ratio is ∼1:1 and their binding is tight (Fig 1d–1f). Thus, we expected that in such cases, using the stQSSA to eliminate a rapid reversible binding in a stochastic model can distort its dynamics. To illustrate this, we constructed a simple gene regulatory network where gene expressions are determined by a reversible binding between transcription factors and genes (Fig 2a left, Table A in S1 Appendix); DNA (D) and a transcription factor (P) reversibly bind to form a complex (D:P). As P acts as a repressor of MR transcription, the transcription rate of MR is proportional to the number of the unbound DNA (D). On the other hand, as P acts as an activator of MA transcription, the transcription rate of MA is proportional to the number of the bound DNA (D:P). Note that the number of unbound and bound DNA can be interpreted as the number of unbound and bound DNA binding sites. In this model, because the reversible binding reaction between D and P is much faster than the other reactions (i.e., the production and the decay of MR and MA), the variables (D and D:P) rapidly reach their QSSs. Thus, by replacing them with their stQSSAs (Dtq and D:Ptq), we can obtain a reduced model (Fig 2a right, Table B in S1 Appendix). The reduced model consists of only the slow variables, MR and MA, because Dtq and D:Ptq are fully determined by the conserved total number of the DNA (DT = D + D:P) and the conserved total number of the transcription factor (PT = P + D:P), as illustrated in Table B in S1 Appendix. This elimination of the fast variables, which are the major source of computational cost, greatly reduces the computation time of stochastic simulations [27–29].

[Figure omitted. See PDF.]

Fig 2. When DNA and a transcription factor bind tightly and their levels are similar, the stQSSA overestimates the number of the unbound DNA.

(a) Full model diagram of a gene regulatory network containing a rapid reversible binding between DNA (D) and a transcription factor (P) to form a complex (D:P) (left, Table A in S1 Appendix). The transcription rates of MR and MA are proportional to D and D:P, respectively. By replacing D and D:P with their stQSSAs (Dtq and D:Ptq), we can obtain a reduced model which consists of only slowly varing MR and MA (right, Table B in S1 Appendix). (b-d) Trajectories of MR (top) and MA (bottom) from the full model (red) and the reduced model (blue) simulated using the Gillespie algorithm (see Tables A and B in S1 Appendix for propensity functions). The lines with colored ranges and the histograms represent the mean ± standard deviation and the stationary distribution of 104 trajectories, respectively. When DT and PT are the same (DT = PT = 10) and the binding is tight (Kd = 10−2), the MR trajectories simulated with the reduced model largely exceed those simulated with the full model (b top) because Dtq overestimates the stochastic QSSA for D (⟨D⟩). On the other hand, D:Ptq accurately approximates the stochastic QSSA for D:P (⟨D:P⟩), and thus the reduced model accurately captures the dynamics of MA (b bottom). If DT is not similar to PT (DT = 15, PT = 10) (c) or the binding is weak (Kd = 10) (d), Dtq and D:Ptq accurately approximate ⟨D⟩ and ⟨D:P⟩, respectively, so that the reduced model accurately captures the dynamics of both MR and MA of the full model.

https://doi.org/10.1371/journal.pcbi.1008952.g002

To test whether the reduced model accurately captures the dynamics of the full model, we compared their stochastic simulations with the Gillespie algorithm (see Tables A and B in S1 Appendix for propensity functions) [52]. When DT and PT are the same and the binding between D and P is tight, MR simulated with the reduced model largely exceeds MR simulated with the full model (Fig 2b top) because the stQSSA (Dtq) overestimates the stochastic QSSA for the number of the unbound DNA (⟨D⟩) which determines the transcription rate of MR (Fig 2a), as seen in Fig 1f. On the other hand, when DT is not similar to PT (Fig 2c top) or the binding is weak (Fig 2d top), Dtq accurately approximates ⟨D⟩ as seen in Fig 1d, and thus the reduced model accurately captures the dynamics of MR in the full model.

Unlike MR (Fig 2b top), the stochastic dynamics of MA of the reduced model and the full model are identical (Fig 2b–2d bottom) because the stQSSA for D:P (D:Ptq) always accurately approximates the stochastic QSSA for the number of the bound DNA (⟨D:P⟩) which determines the transcription of MA (Fig 1c). Taken together, the stQSSA can be used to describe transcriptional activation depending on bound DNA under any conditions (Fig 2b–2d bottom). On the other hand, it needs to be restrictively used to describe transcriptional repression depending on unbound DNA (Fig 2b–2d top).

stQSSA can distort oscillatory dynamics

To illustrate how the stQSSA distorts the dynamics when the molar ratio between tightly binding species is ∼1:1, we investigated the simple model where the molar ratio is conserved (Fig 2). However, the total copy numbers of binding species and thus their molar ratio can be varied (e.g., oscillate) in a living cell due to other reactions in a larger system. This raises the question of whether the model reduction based on the stQSSA is accurate or not if the molar ratio is temporarily ∼1:1. To investigate this, we used a modified Kim-Forger model, which describes the transcriptional negative feedback loop of the mammalian circadian clock [24, 49, 51]. In this model (Fig 3a top, Table C in S1 Appendix), free activator (A) promotes the transcription of mRNA (M), and the protein translated from M produces repressor (R) passing through several steps (Pi, i = 1, 2, 3). Then R reversibly binds with A to form a complex (A:R) which no longer promotes the transcription, and thus represses its own transcription. In this model, the reversible binding between R and A is much faster than the other reactions (i.e., production and decay). Thus, by replacing the fast variable A, which determines the transcription rate of M, with its stQSSA (Atq), we can obtain a reduced model (Fig 3a bottom, Table D in S1 Appendix). The reduced model consists of only the slow variables, RT, M and Pi, because Atq is fully determined by the conserved total number of the activator (AT = A + A:R) and the slowly varying total number of the repressor (RT = R + A:R), as illustrated in Table D in S1 Appendix.

[Figure omitted. See PDF.]

Fig 3. stQSSA can distort the dynamics of a biological oscillator.

(a) Full model diagram of an oscillatory transcriptional negative feedback loop (top, Table C in S1 Appendix). Unbound activator (A) promotes the transcription of mRNA (M), and the protein translated from M produces repressor (R) passing through several steps (Pi, i = 1, 2, 3). Then R binds with A to form a complex (A:R) which is transcriptionally inactive, and thus represses its own transcription. As the reversible binding between R and A is rapid, by replacing A with its stQSSA (Atq), we can obtain a reduced model which consists of only slowly varying RT, M, and Pi (bottom, Table D in S1 Appendix). (b-c) Oscillatory trajectories of M (green) and RT/AT (orange) simulated with the full model (b top) and the reduced model (b bottom), using the Gillespie algorithm (see Tables C and D in S1 Appendix for propensity functions). When R binds with A tightly (Kd = 10−4) both the full model and the reduced model show the oscillatory behaviors. However, when the trajectory of RT/AT stays near 1 (dashed lines in b), Atq overestimates the stochastic QSSA for A (⟨A⟩), and thus the transcription more frequently occurs in the reduced model (b bottom) compared to the full model (b top). As a result, the reduced model predicts a shorter period than the full model (c). (d-e) On the other hand, when the degradation rate of R increases and thus the trajectory of RT/AT stays near 1 for a short time (d; dashed lines), the reduced model accurately captures the dynamics of the full model (e).

https://doi.org/10.1371/journal.pcbi.1008952.g003

In the model, because R tightly binds with A, when RT/AT ≈ 1, Atq overestimates the stochastic QSSA for A (⟨A⟩) and thus the transcription rate of M. As a result, when the trajectory of RT/AT reaches close to 1 (dashed lines in Fig 3b), the transcription more frequently occurs in the reduced model (Fig 3b bottom) compared to the full model (Fig 3b top). This overestimated transcriptional activity leads to the shorter peak-to-peak periods of the reduced model compared to the full model (Fig 3c). On the other hand, when the degradation rate of R increases and thus the trajectory of RT/AT stays near 1 for an extremely short time (Fig 3d dashed lines), the reduced model accurately captures the dynamics of the full model (Fig 3e). Taken together, if ∼1:1 molar ratio between the tightly binding activator and repressor of the transcriptional negative feedback loop persists for a considerable time, using the stQSSA overestimates the transitional activity and thus the frequency of oscillation.

stQSSA can distort bistable dynamics

To investigate how the misuse of the stQSSA distorts the dynamics of a bistable switch, we used a previously developed bistable switch model for the maturation promoting factor, cyclin B/Cdc2, whose activation promotes mitosis (Fig 4a top, Table E in S1 Appendix) [46, 57]. In the model, the inactive form of cyclin B/Cdc2 (P) is converted to an active form (M) by Cdc25 (D). Furthermore, as M activates D which converts P to M, M promotes its own activation (i.e., form a positive feedback loop; see [46, 57] for details). The positive feedback loop is suppressed by Suc1 protein (B) as it binds with M to form a complex (M:B) which no longer activates D. The total activated cyclin B/Cdc2 (M and M:B) become P with the same constant rate. In this model, the reversible binding between M and B is much faster than the other reactions. Thus, by replacing the fast variable M with its stQSSA (Mtq) a reduced model can be derived (Fig 4a bottom, Table F in S1 Appendix). The reduced model consists of only the slow variables, MT and P, because Mtq is fully determined by the conserved total number of Suc1 (BT = B + M:B) and the slowly varying total number of the activated cyclin B/Cdc2 (MT = M + M:B), as illustrated in Table F in S1 Appendix.

[Figure omitted. See PDF.]

Fig 4. stQSSA can distort the dynamics of a bistable switch.

(a) Full model diagram of a bistable switch for mitosis (top, Table E in S1 Appendix). The inactive form of cyclin B/Cdc2 (P) becomes an active form (M) by Cdc25 (D). In this process, M enhances its own activation by activating D, and thus forms a positive feedback loop (see [46, 57] for details). The positive feedback loop is suppressed as Suc1 protein (B) binds with M to form a complex (M:B) which does not activate D. The total activated cyclin B/Cdc2, M and M:B, becomes P with the same constant rate. As the reversible binding between M and B is rapid, by replacing M with its stQSSA (Mtq), we can obtain a reduced model which consists of only slowly varying MT and P (bottom, Table F in S1 Appendix). (b-c) Simulated trajectories (b) and the stationary distributions (c) of MT from the full model and the reduced model using the Gillespie algorithm (see Tables E and F in S1 Appendix for propensity functions). When M binds with B tightly (Kd = 10−3), both the full model and the reduced model show the bistable behaviors between the upper and lower modes, which are separated by MT/BT = 1 (dashed line in b). However, because Mtq overestimates the stochastic QSSA for M (⟨M⟩) when MT/BT is close to 1, the trajectory from the reduced model is more attracted to the upper mode compared to the full model (b). As a result, the bimodal distribution of MT from the reduced model is biased to the upper mode (c). (d-e) On the other hand, when the binding between M and B becomes weak (Kd = 10), Mtq accurately estimates ⟨M⟩, and thus the reduced model accurately captures the dynamics of the full model, which no longer shows the bistable behavior.

https://doi.org/10.1371/journal.pcbi.1008952.g004

When M and B tightly bind, both the full model and the reduced model show the bistable behaviors (i.e., bimodal stationary distributions) of MT (Fig 4b). However, the trajectory of the reduced model is more attracted to the upper mode of MT compared to the full model (Fig 4b and 4c). This dynamics biased to the upper mode occurs because Mtq overestimates the stochastic QSSA for M (⟨M⟩) near the MT/BT = 1 region (Fig 4b dashed line) which separates the upper and lower modes. On the other hand, when the binding between M and B becomes weak, Mtq accurately approximates the stochastic QSSA for M even when MT is similar to BT. Thus, the reduced model accurately captures the dynamics of the full model, which no longer shows bistable behavior (Fig 4d and 4e). Taken together, when the binding between activated Cyclin B/Cdc2 and Suc1 protein is tight, which is essential to generate the bistable switch, using the stQSSA overestimates the activation of Cyclin B/Cdc2 and distorts the dynamics of the bistable switch.

An alternative approach when the stQSSA is not applicable

In the presence of a rapid and tight reversible binding between species whose molar ratio is ∼1:1, the reduction of stochastic models with the stQSSA for the number of the unbound species can cause errors (Figs 2b, 3c, and 4c). In such cases, due to the tight binding, the two species tend to bind until no molecules of one species are left (Fig 1f). Specifically, if AT ≤ BT (AT ≥ BT), the majority of the A (B) will be bound. Thus, in the presence of tight binding, we can assume that the stationary distributions of A or B are concentrated on 0 and 1. This low-state assumption allows us to derive the simple approximation for the stochastic QSSA (⟨A⟩ in Eq (4)) (see Methods for details):(9)We will refer to this approximation as the stochastic “low-state” QSSA (slQSSA). Since this approximation relies on the property that the state space is restricted at the low level, the stochastic QSSA with singular perturbation analysis introduced in [58] could be used to derive an alternative approximation.

The accuracy of the slQSSA for A (Eq (9)) is expected to increase when AT Kd decreases because AT Kd is an approximated number of the unbound A. On the other hand, the accuracy of the stQSSA for A decreases as AT Kd decreases (Fig 1d). To investigate this, we calculated the maximum relative error of Atq () and Alq () to the stochastic QSSA for A (⟨A⟩) for each AT Kd and Kd (Fig 5a and 5b). As expected, when AT Kd is low and high, the slQSSA and the stQSSA are accurate, respectively. In particular, when AT Kd < 10−1 and AT Kd > 101, and RA are less than 0.1 (i.e., the relative errors are less than 10%), respectively.

[Figure omitted. See PDF.]

Fig 5. slQSSA can be used to reduce multiscale stochastic biochemical systems containing rapid reversible bindings when the stQSSA is not applicable.

(a-b) Heat maps of the relative errors ( and ) when the stQSSA (Atq) and the two-state slQSSA (Alq) approximate the stochastic QSSA for A (⟨A⟩) in the reversible binding reaction (Eq (1)). Color represents the maximum value of RA and for each AT Kd and Kd when BT varies, and the dashed lines represent when those values are 0.1. When AT Kd are high and low, the stQSSA and the slQSSA are accurate, respectively. The parameters used in Fig 2b (triangle), Fig 3b (circle), and Fig 4b (square) are located in the region where the slQSSA (b), but not the stQSSA (a), is accurate (the circle is actually located outside of the heat maps; AT Kd = 5 × 10−4 and Kd = 10−4). (c-e) As a result, the full models are successfully reduced with the slQSSA (c-e) but not the stQSSA (Figs 2b, 3c, and 4c). See Tables B, D, and F in S1 Appendix for the propensity functions used for the simulations and Fig C in S1 Appendix for a benchmark comparison with GillesPy2, one of the major, standard software suites for stochastic simulation [59]. (f) The adaptive use of the stQSSA and the slQSSA to approximate the stochastic QSSA for A when AT Kd > 101 and otherwise, respectively, guarantees the successful reduction of stochastic models containing rapid reversible bindings. Note that when 10−1 < AT Kd < 101, the slQSSAs with more than two states need to be used (see Fig B in S1 Appendix for details).

https://doi.org/10.1371/journal.pcbi.1008952.g005

The parameters used in Fig 2b (triangle), Fig 3b (circle), and Fig 4b (square) are located in the region where the stQSSA is inaccurate (Fig 5a) but the slQSSA is accurate (Fig 5b). Therefore, with these parameters, the reduced models obtained by using the slQSSAs accurately capture the dynamics of the full models for the simple gene regulatory network (Fig 5c, Table B in S1 Appendix), the transcriptional negative feedback loop (Fig 5d, Table D in S1 Appendix), and the bistable switch for mitosis (Fig 5e, Table F in S1 Appendix), unlike the stQSSA (Figs 2b, 3c, and 4c). Furthermore, by allowing A or B to reach more than two states (e.g., 0, 1, and 2), more accurate slQSSAs can be derived (see Methods for details). In particular, the relative errors of the slQSSAs derived by allowing the 3/4/5 states are less than 0.1 when AT Kd is less than 2/5/10, respectively (Fig B in S1 Appendix). Consequently, for the error tolerance of 0.1, if AT Kd < 101 and thus the stQSSA is inaccurate, the slQSSA can be used to approximate the stochastic QSSA for A (Fig 5f). Taken together, by using either stQSSA or slQSSA depending on AT Kd, we can always accurately reduce multiscale stochastic biochemical systems with rapid reversible bindings. Of course, for a different error tolerance, we need a different threshold of AT Kd and the number of states for the slQSSA. To facilitate the calculation of such change depending on the error tolerance, we have developed a user-friendly open-source computational package, ASSISTER (Fig 6). In particular, the function Gillespie_Reduction in the package automatically constructs a reduced model adaptively using the more accurate one between the two approximation methods and performs accurate and efficient stochastic simulations (see S1 Appendix for the manual).

[Figure omitted. See PDF.]

Fig 6. Schematic diagram for the computational package, ASSISTER.

When a stochastic model containing a rapid reversible binding () and the error tolerance (ϵ) are given as inputs, the auxiliary function QSSA_Threshold determines the threshold of the total number of binding molecules, Dthres, and the number of states for the slQSSA, L. When the DT = D + D:P is less (larger) than Dthres, the L-state slQSSA (stQSSA) approximates the exact stochastic QSSA with a relative error less than ϵ. Based on the relationship between Dthres and DT, the more accurate one between the two models is adaptively chosen. In this way, Gillespie_Reduction performs efficient and accurate stochastic simulations, yielding the simulated trajectories as the final output. See the manual in S1 Appendix for a more detailed description of the input and output.

https://doi.org/10.1371/journal.pcbi.1008952.g006

Discussion

Reversible binding between molecules—for example, between DNA and a transcription factor, a ligand and a receptor, and an enzyme and a substrate—is a fundamental reaction for numerous biological functions [42]. As the reversible binding reactions occur typically on a timescale of 1∼1000ms, which is much faster than the other reactions (e.g., 30min for a mammalian mRNA transcription or a protein translation and 10h for their typical lifetimes) [43], a system containing the rapid reversible binding becomes a multi-timescale system. In such multi-timescale systems, the rapid reversible binding prohibitively increases the computational cost of stochastic simulations. Accordingly, to accelerate stochastic simulations, various methods have been developed [44, 60]. In particular, the model reduction using the stQSSA has successfully simplified various stochastic models in numerous studies [7, 30–33, 38, 39]. Thus, it has been commonly believed that the stQSSA is generally accurate for any conditions, until a recent counterexample was identified [41]. In this work, we rigorously derived the validity conditions for using the stQSSA to reduce stochastic models with a rapid reversible binding. Specifically, we showed that the relative error of the stQSSA for the number of unbound species (RA) mainly depends on the relative sensitivity of the stQSSA (SA, Eq (8)), which attains maximum value at AT = BT + Kd. This allowed us to find that the stQSSA for the number of the unbound species is inaccurate if their molar ratio is ∼1:1 and their binding is tight (Fig 1f). In that case, the stQSSA highly overestimates the number of the unbound species. Therefore, the reduced models obtained by using the stQSSA distort the dynamics of the gene regulatory model (Fig 2b), the transcriptional negative feedback loop model for circadian rhythms (Fig 3c), and the bistable switch model for mitosis (Fig 4c).

When the reversible binding reactions are sufficiently faster than the other reactions, the deterministic tQSSA is known to be accurate [7, 48, 50]. Indeed, for all examples considered in our work (Figs 2b, 3b and 4b), the deterministic simulations with the tQSSA are accurate, unlike the stochastic simulations. This indicates that it is risky to investigate the validity conditions of the stQSSA solely based on the validity conditions of the deterministic tQSSA. Instead, the direct derivation of the relative error of the stQSSA is needed, as demonstrated in this study (Eq (7)). It would be interesting in future work to perform such error analysis for more complex examples, such as coupled enzymatic networks with multiple rapid reversible bindings [26, 61, 62].

The rapid reversible reactions are typical conditions for model reductions using the “partial-equilibrium approximation,” which confines species concentrations to the equilibrium states. In general, this condition does not imply a timescale separation between the variables (A, B, and C), limiting the application of the QSSA. However, the rapid reversible binding guarantees that the total variables AT and BT always evolve more slowly than the variables A, B, and C. Therefore, the QSSAs in terms of the total variables can lead to accurate model reductions in both deterministic [7, 47, 48, 50] and stochastic [27–30] regimes. In this work, we investigated under which conditions the complex stochastic QSSA (Eq (4)) can be approximated by the corresponding simple deterministic QSSA (Eq (5)), referred to as the stQSSA. We found that such an approximation becomes inaccurate when two species with similar levels tightly bind. Thus, we derived an alternative approximation for the stochastic QSSA: slQSSA. Using the two approximations for the stochastic QSSA, we developed the universally valid reduction framework for stochastic models containing rapid reversible bindings. We also provided the user-friendly open-source computational package, ASSISTER, for this framework. On the other hand, in the absence of rapid reversible binding, one can reduce a model by assuming that ‘highly reactive species’ are in their QSSs [58]. Interestingly, the reduced stochastic models were often different from the heuristically reduced models obtained with the deterministic QSSA. It would be interesting in future work to investigate when such discrepancies occur.

While the deterministic tQSSA (Eq (3)) was used to approximate the stochastic QSSA for the number of reversibly binding species in this work, a simpler deterministic QSSA referred to as the “standard” QSSA (sQSSA) is more widely used to approximate the stochastic QSSA due to its simplicity [8–20, 27, 32]. For instance, the stochastic sQSSA for C in Eq (1), which has the Michaelis-Menten type form(10)has been widely used as a propensity function for the Gillespie algorithm. Though this sQSSA has been widely used, it is less accurate than the stQSSA (Eq (5)) [38, 39]. This is why many examples showing the inaccuracy of the stochastic sQSSA have been reported [33–40], whereas only one example showing the inaccuracy of the stQSSA has been reported [41]. Note that while Eq (10) is different from the typical “Michaelis-Menten” equation, which uses the Michaelis-Menten constant instead of the dissociation constant (Kd), they become nearly the same when the timescale of reversible binding is faster than the catalytic reaction. Importantly, our work also provides the validity condition for using the stochastic sQSSA (Eq (10)). That is, when BT + Kd ≫ AT, which is known as the low-enzyme concentration condition, Csq ≈ Ctq [7], indicating that the Michaelis-Menten type sQSSA for the bounded species (Eq (10)) can be used to reduce models containing the rapid reversible binding. Similarly, when BT + Kd ≫ AT, could also be used. This is consistent with the validity conditions for the stochastic sQSSA derived under the assumption of either a low fluctuation level [32] or a low copy number [40]. Furthermore, the “pre-factor” QSSA (pQSSA), which is more accurate than the sQSSA, has also been used for stochastic simulations [63, 64]. However, recent studies have shown that the stQSSA is more accurate than the stochastic pQSSA (see [38, 39] for details).

The accuracy of the stQSSA for the number of the unbound species depends on both the molar ratio between reversibly binding species and the tightness of their binding (Fig 1d). However, as the molar ratio typically varies in larger models containing reversible binding, practically, the accuracy is mainly determined by the tightness of binding. Specifically, for the relative error of the stQSSA to be less than 0.1, AT Kd (≈ the number of the unbound A) should be larger than 10 (Fig 5a dashed line). This AT Kd value-based criteria explains the controversy about the accuracy of the stQSSA in previous studies. That is, AT Kd was less than 10 in a previous study where the reduced model obtained by using the stQSSA was inaccurate [41]. On the other hand, AT Kd were much greater than 10 in all of the examples investigated in previous studies reporting the accuracy of the stQSSA [7, 33, 38, 39, 53, 54]. Furthermore, the stQSSA always accurately approximates the stochastic QSSA for the number of the bound species (Fig 1c). This explains why the stQSSA was accurate in previous studies where the stQSSA was used to approximate the number of enzyme-substrate complex [30–33].

In real biological systems, the validity condition of the stQSSA (AT Kd > 10) is not always guaranteed. Specifically, the range of AT Kd can span approximately from 10−3 to 1010 in human cells (Ω = 10−15 ∼ 10−14 m3) since the protein-protein dissociation constant () is 10fM ∼ 1μM (i.e., 1012 ∼ 1020 m−3), and the numbers of molecules (AT) is 100 ∼ 104 [43, 65]. Moreover, in smaller cells like budding yeast (Ω = 10−17 ∼ 10−16 m3) or E. Coli (Ω = 10−19 ∼ 10−18 m3) cells, the range of AT Kd can span from 10−7 to 108, which contains the region in which the stQSSA can be extremely inaccurate. Accordingly, the slQSSA, which accurately approximates the stochastic QSSA when AT Kd is less than 10, is necessary. Specifically, the relative error of the slQSSA, unlike that of the stQSSA (Fig 5a and 5b), decreases as AT Kd decreases because the slQSSA relies on the assumption that the stationary distributions of the number of the unbound species (≈ AT Kd) are concentrated on the few lowest states. Taken together, by using the stQSSA and the slQSSA when the AT Kd value is greater and less than 10, respectively, one can always accurately simplify stochastic models containing rapid reversible binding reactions to accelerate simulation and also facilitate stochastic analysis (Fig 5f). This can be facilitated by the computational package, ASSISTER (Fig 6).

Methods

Exact bounds for the relative error of the stQSSA to the stochastic QSSA

In this section, we derive the exact upper and lower bounds for (Eq (7)) where Atq and ⟨A⟩ are the stQSSA and the stochastic QSSA for A, respectively. From the CME describing the reversible binding reaction (Eq (1)), the following steady-state moment equation can be derived:(11)where ⟨⋅⟩ is the stationary expectation. Eq (11) becomes ⟨A ⋅ (BT − AT + A)⟩ = Kd⟨AT − A⟩ by using the definitions AT = A + C, BT = B + C, and Kd = kbΩ/kf. Since AT and BT are invariant under the reversible binding reactions in Eq (1), we obtain ⟨A2⟩ − (AT − BT − Kd)⟨A⟩ − AT Kd = 0, and by using the relation ⟨A2⟩ = Var(A) + ⟨A⟩2, we get the following quadratic equation:(12)The non-negative root of this quadratic equation becomes ⟨A⟩:(13)By subtracting Eq (13) from Eq (5), we get(14)(15)Since 0 ≤ (AT − BT − Kd)2 + 4AT Kd − 4Var(A) ≤ (AT − BT − Kd)2 + 4AT Kd, we get the bounds for Atq − ⟨A⟩ from Eq (15):(16)By dividing Eq (16) by ⟨A⟩, we can get the bounds for the relative error, as follows:This can be re-expressed as FA SA ≤ RA ≤ 2FA SA (Eq (7)) because is the Fano factor of A (FA), and is the relative sensitivity of Atq, i.e., .

The relative sensitivity, SA, attains the maximum value when the term in the square root of the denominator has the minimum value, i.e., BT = AT − Kd (Eq (8)). In particular, SA has a large maximum value when Kd ≪ 1 at AT = BT + Kd ≈ BT. On the other hand, if AT ≪ BT, SA ≈ 0 because the majority of A presents in the bound state regardless of BT (i.e.,). When AT ≥ BT, because as BT decreases by one, approximately one A is released from the complex. In this case, if AT ≫ BT, the majority of A are free and thus , leading to SA ≈ 0. However, if AT ≈ BT, the majority of A is sequestered by B, Atq ≈ 0, leading to SA ≫ 1. When binding is weak (Kd ≫ 1), SA ≈ 0 because the number of A, which is approximated by Atq, changes little as BT changes (i.e., ). Taken together, SA is large only when the binding reaction is tight (Kd ≪ 1) and the binding species are present with 1:1 molar ratio (AT ≈ BT).

Since , SA = O(Ω−1) if the concentrations remain constant. This implies that when the volume Ω goes to infinity (i.e., thermodynamic limit), SA and thus RA become zero (i.e., the stochastic QSSA becomes nearly identical to the deterministic QSSA (tQSSA)). On the other hand, as Ω goes to zero (i.e., the volume of the system gets smaller), SA goes to infinity.

Derivation of the stochastic QSSA and the slQSSA

Here we derive the stochastic QSSA for A (⟨A⟩, Eq (4)). Let p(l) be the probability that A = l at its stationary distribution (i.e., the probability that A(∞) = l). Then the following recurrence relation of p(l) can be obtained from the steady-state CME:(17)Let A0 = max{AT − BT, 0}. Since A0 is the lowest state that A can reach, p(l) = 0 for l < A0. Then we can inductively prove that the following relation satisfies Eq (17):(18)where . Then, because ∑p(l) = 1, if A0 ≤ l ≤ AT, and p(l) = 0 otherwise by Eq (18). Therefore, we can obtain the stationary average number of A (Eq (4)) as

Next we derive the slQSSA, which is the approximation for Eq (4). In the presence of tight binding, we can assume that the stationary distributions of A and B are concentrated on the states {0, 1} when AT < BT and AT ≥ BT, respectively. Since when the distribution of B is concentrated on 0 and 1, the distribution of A is concentrated on AT − BT and AT − BT + 1, we can simply say that the distribution of A is concentrated on A0 and A0 + 1. Thus, by assuming that is approximately zero for l > A0 + 1 and , we can derive the two-state slQSSA for A (Eq (9)) as follows:In general, for any integer k ≥ 2, we can derive the k-state slQSSA as(19)

Computational package for universally valid reduction of stochastic models containing rapid reversible binding reactions

We have developed a user-friendly computational package ASSISTER that contains three main codes implemented in MATLAB (Fig 6): LQSSA, QSSA_Threshold, and Gillespie_Reduction. LQSSA calculates the L-state slQSSA (Eq (19)) for given AT, BT, Kd, and L. QSSA_Threshold determines which of the stQSSA and the L-state slQSSA ensures a smaller error than a tolerance ϵ for a given Kd value. This allows the function Gillespie_Reduction to perform accurate stochastic simulations for any values of the parameters with the adaptive choice of the valid approximation method determined by using QSSA_Threshold (Fig 6). See S1 Appendix for details and the manual. ASSISTER can be found at https://github.com/Mathbiomed/ASSISTER.

Supporting information

S1 Appendix. Supplementary Methods, Tables A-F, and Figs A-C.

https://doi.org/10.1371/journal.pcbi.1008952.s001

(PDF)

Acknowledgments

We thank Krešimir Josić and John J. Tyson for valuable comments.

Citation: Song YM, Hong H, Kim JK (2021) Universally valid reduction of multiscale stochastic biochemical systems using simple non-elementary propensities. PLoS Comput Biol 17(10): e1008952. https://doi.org/10.1371/journal.pcbi.1008952

References

1. Tyson JJ, Novak B. A dynamical paradigm for molecular cell biology. Trends Cell Biol. 2020;. pmid:32362451

2. Hill AV. The possible effects of the aggregation of the molecules of haemoglobin on its dissociation curves. J Physiol. 1910;40:4–7.

3. Michaelis L, Menten ML. Die kinetik der invertinwirkung. Biochem Z. 1913;49(333-369):352.

4. Gunawardena J. Time-scale separation–Michaelis and Menten’s old idea, still bearing fruit. FEBS J. 2014;281(2):473–488. pmid:24103070

5. Laidler KJ. Theory of the transient phase in kinetics, with special reference to enzyme systems. Can J Chem. 1955;33(10):1614–1624.

6. Morrison JF. Kinetics of the reversible inhibition of enzyme-catalysed reactions by tight-binding inhibitors. Biochim Biophys Acta (BBA)-Enzymology. 1969;185(2):269–286. pmid:4980133

7. Kim JK, Tyson JJ. Misuse of the Michaelis–Menten rate law for protein interaction networks and its remedy. PLoS Comput Biol. 2020;16(10):e1008258. pmid:33090989

8. Gonze D, Halloy J, Gaspard P. Biochemical clocks and molecular noise: Theoretical study of robustness factors. J Chem Phys. 2002;116(24):10997–11010.

9. Gonze D, Halloy J, Goldbeter A. Deterministic versus stochastic models for circadian rhythms. J Biol Phys. 2002;28(4):637–653. pmid:23345804

10. Pedraza JM, van Oudenaarden A. Noise propagation in gene networks. Science. 2005;307(5717):1965–1969. pmid:15790857

11. Scott M, Ingalls B, Kærn M. Estimations of intrinsic and extrinsic noise in models of nonlinear genetic networks. Chaos. 2006;16(2):026107. pmid:16822039

12. Tian T, Burrage K. Stochastic models for regulatory networks of the genetic toggle switch. Proc Natl Acad Sci U S A. 2006;103(22):8372–8377. pmid:16714385

13. Çağatay T, Turcotte M, Elowitz MB, Garcia-Ojalvo J, Süel GM. Architecture-dependent noise discriminates functionally analogous differentiation circuits. Cell. 2009;139(3):512–522. pmid:19853288

14. Ouattara DA, Abou-Jaoudé W, Kaufman M. From structure to dynamics: Frequency tuning in the p53-Mdm2 network. II: Differential and stochastic approaches. J Theor Biol. 2010;264(4):1177–1189. pmid:20346959

15. Gonze D, Abou-Jaoudé W, Ouattara DA, Halloy J. How Molecular Should Your Molecular Model Be?: On the Level of Molecular Detail Required to Simulate Biological Networks in Systems and Synthetic Biology. Methods Enzymol. 2011;487:171–215. pmid:21187226

16. Smadbeck P, Kaznessis Y. Stochastic model reduction using a modified Hill-type kinetic rate law. J Chem Phys. 2012;137(23):234109. pmid:23267473

17. Riba A, Bosia C, El Baroudi M, Ollino L, Caselle M. A combination of transcriptional and microRNA regulation improves the stability of the relative concentrations of target genes. PLoS Comput Biol. 2014;10(2):e1003490. pmid:24586138

18. Dovzhenok AA, Baek M, Lim S, Hong CI. Mathematical modeling and validation of glucose compensation of the neurospora circadian clock. Biophys J. 2015;108(7):1830–1839. pmid:25863073

19. Zhang W, Tian T, Zou X. Negative feedback contributes to the stochastic expression of the interferon-β gene in virus-triggered type I interferon signaling pathways. Math Biosci. 2015;265:12–27. pmid:25892253

20. Schuh L, Saint-Antoine M, Sanford EM, Emert BL, Singh A, Marr C, et al. Gene networks with transcriptional bursting recapitulate rare transient coordinated high expression states in cancer. Cell Syst. 2020;10(4):363–378. pmid:32325034

21. Chen J, Grima R. Small protein number effects in stochastic models of autoregulated bursty gene expression. J Chem Phys. 2020;152(8):084115.

22. Bersani AM, Bersani E, Mastroeni L. Deterministic and stochastic models of enzymatic networks—applications to pharmaceutical research. Computers & Mathematics with Applications. 2008;55(5):879–888.

23. Choi B, Rempala GA, Kim JK. Beyond the Michaelis-Menten equation: Accurate and efficient estimation of enzyme kinetic parameters. Sci Rep. 2017;7(1):1–11. pmid:29208922

24. D’Alessandro M, Beesley S, Kim JK, Jones Z, Chen R, Wi J, et al. Stability of wake-sleep cycles requires robust degradation of the PERIOD protein. Curr Biol. 2017;27(22):3454–3467. pmid:29103939

25. Beesley S, Kim DW, D’Alessandro M, Jin Y, Lee K, Joo H, et al. Wake-sleep cycles are severely disrupted by diseases affecting cytoplasmic homeostasis. Proc Natl Acad Sci U S A. 2020;117(45):28402–28411. pmid:33106420

26. Bersani AM, Borri A, Carravetta F, Mavelli G, Palumbo P. On a stochastic approach to model the double phosphorylation/dephosphorylation cycle. Math Mech Complex Syst. 2020;8(4):261–285.

27. Rao CV, Arkin AP. Stochastic chemical kinetics and the quasi-steady-state assumption: Application to the Gillespie algorithm. J Chem Phys. 2003;118(11):4999–5010.

28. Cao Y, Gillespie DT, Petzold LR. The slow-scale stochastic simulation algorithm. J Chem Phys. 2005;122(1):014116. pmid:15638651

29. Gómez-Uribe CA, Verghese GC, Tzafriri AR. Enhanced identification and exploitation of time scales for model reduction in stochastic chemical kinetics. J Chem Phys. 2008;129(24):244112. pmid:19123500

30. Cao Y, Gillespie DT, Petzold LR. Accelerated stochastic simulation of the stiff enzyme-substrate reaction. J Chem Phys. 2005;123(14):144917. pmid:16238434

31. Barik D, Paul MR, Baumann WT, Cao Y, Tyson JJ. Stochastic simulation of enzyme-catalyzed reactions with disparate timescales. Biophys J. 2008;95(8):3563–3574. pmid:18621809

32. Sanft KR, Gillespie DT, Petzold LR. Legitimacy of the stochastic Michaelis–Menten approximation. IET Syst Biol. 2011;5(1):58–69. pmid:21261403

33. Kang HW, KhudaBukhsh WR, Koeppl H, Rempała GA. Quasi-steady-state approximations derived from the stochastic model of enzyme kinetics. Bull Math Biol. 2019;81(5):1303–1336. pmid:30756234

34. Thomas P, Straube AV, Grima R. Communication: limitations of the stochastic quasi-steady-state approximation in open biochemical reaction networks. J Chem Phys. 2011;135(18):181103. pmid:22088045

35. Thomas P, Straube AV, Grima R. The slow-scale linear noise approximation: an accurate, reduced stochastic description of biochemical networks under timescale separation conditions. BMC Syst Biol. 2012;6(1):1–23. pmid:22583770

36. Agarwal A, Adams R, Castellani GC, Shouval HZ. On the precision of quasi steady state assumptions in stochastic dynamics. J Chem Phys. 2012;137(4):044105. pmid:22852595

37. Lawson MJ, Petzold L, Hellander A. Accuracy of the Michaelis–Menten approximation when analysing effects of molecular noise. J R Soc Interface. 2015;12(106):20150054. pmid:25833240

38. Kim JK, Josić K, Bennett MR. The validity of quasi-steady-state approximations in discrete stochastic simulations. Biophys J. 2014;107(3):783–793. pmid:25099817

39. Kim JK, Josić K, Bennett MR. The relationship between stochastic and deterministic quasi-steady state approximations. BMC Syst Biol. 2015;9(1):1–13. pmid:26597159

40. Holehouse J, Grima R. Revisiting the reduction of stochastic models of genetic feedback loops with fast promoter switching. Biophys J. 2019;117(7):1311–1330. pmid:31540707

41. Kim JK, Sontag ED. Reduction of multiscale stochastic biochemical reaction networks using exact moment derivation. PLoS Comput Biol. 2017;13(6):e1005571. pmid:28582397

42. Williams M, Daviter T. Protein-Ligand Interactions. Springer; 2016.

43. Milo R, Phillips R. Cell biology by the numbers. Garland Science; 2015.

44. Marchetti L, Priami C, Thanh VH. Simulation algorithms for computational systems biology. Springer; 2017.

45. Borghans JA, De Boer RJ, Segel LA. Extending the quasi-steady state approximation by changing variables. Bull Math Biol. 1996;58(1):43–63. pmid:8819753

46. Thron C. A model for a bistable biochemical trigger of mitosis. Biophys Chem. 1996;57(2-3):239–251. pmid:8573678

47. Schnell S, Maini PK. Enzyme kinetics far from the standard quasi-steady-state and equilibrium approximations. Math Comput Model. 2002;35(1-2):137–144.

48. Tzafriri AR. Michaelis-Menten kinetics at high enzyme concentrations. Bull Math Biol. 2003;65(6):1111–1129. pmid:14607291

49. Kim JK, Forger DB. A mechanism for robust circadian timekeeping via stoichiometric balance. Mol Syst Biol. 2012;8(1):630. pmid:23212247

50. Bersani AM, Bersani E, Dell’Acqua G, Pedersen MG. New trends and perspectives in nonlinear intracellular dynamics: one century from Michaelis–Menten paper. Contin Mech Thermodyn. 2015;27(4):659–684.

51. Kim JK. Protein sequestration versus Hill-type repression in circadian clock models. IET Syst Biol. 2016;10(4):125–135. pmid:27444022

52. Gillespie DT. Exact stochastic simulation of coupled chemical reactions. J Chem Phys. 1977;81(25):2340–2361.

53. Kim JK, Rempala GA, Kang HW. Reduction for stochastic biochemical reaction networks with multiscale conservations. Multiscale Model Simul. 2017;15(4):1376–1403.

54. Herath N, Del Vecchio D. Reduced linear noise approximation for biochemical reaction networks with time-scale separation: The stochastic tQSSA+. J Chem Phys. 2018;148(9):094108.

55. Sontag ED, Zeilberger D. A symbolic computation approach to a problem involving multivariate Poisson distributions. Adv Appl Math. 2010;44(4):359–377.

56. Hong H, Kim J, Al-Radhawi MA, Sontag ED, Kim JK. Derivation of stationary distributions of biochemical reaction networks via structure transformation. Commun Biol. 2021;4(1):1–10. pmid:34031517

57. Thron CD. Bistable biochemical switching and the control of the events of the cell cycle. Oncogene. 1997;15(3):317–325. pmid:9233766

58. Mastny EA, Haseltine EL, Rawlings JB. Two classes of quasi-steady-state model reductions for stochastic kinetics. J Chem Phys. 2007;127(9):094106. pmid:17824731

59. Drawert B, Hellander A, Bales B, Banerjee D, Bellesia G, Daigle BJ Jr, et al. Stochastic simulation service: bridging the gap between the computational expert and the biologist. PLoS Comput Biol. 2016;12(12):e1005220. pmid:27930676

60. Schnoerr D, Sanguinetti G, Grima R. Approximation and inference methods for stochastic biochemical kinetics—a tutorial review. J Phys A Math Theor. 2017;50(9):093001.

61. Ciliberto A, Capuani F, Tyson JJ. Modeling networks of coupled enzymatic reactions using the total quasi-steady state approximation. PLoS Comput Biol. 2007;3(3):e45. pmid:17367203

62. Kumar A, Josić K. Reduced models of networks of coupled enzymatic reactions. J Theor Biol. 2011;278(1):87–106. pmid:21377474

63. Isaacs FJ, Hasty J, Cantor CR, Collins JJ. Prediction and measurement of an autoregulatory genetic module. Proc Natl Acad Sci U S A. 2003;100(13):7714–7719. pmid:12808135

64. Kepler TB, Elston TC. Stochasticity in transcriptional regulation: origins, consequences, and mathematical representations. Biophys J. 2001;81(6):3116–3136. pmid:11720979

65. Kastritis PL, Moal IH, Hwang H, Weng Z, Bates PA, Bonvin AM, et al. A structure-based benchmark for protein–protein binding affinity. Protein Science. 2011;20(3):482–491. pmid:21213247