A multiscale approximation in a heat shock response model of E. coli
- Hye-Won Kang^{1}Email author
DOI: 10.1186/1752-0509-6-143
© Kang; licensee BioMed Central Ltd. 2012
Received: 17 November 2011
Accepted: 7 November 2012
Published: 21 November 2012
Abstract
Background
A heat shock response model of Escherichia coli developed by Srivastava, Peterson, and Bentley (2001) has multiscale nature due to its species numbers and reaction rate constants varying over wide ranges. Applying the method of separation of time-scales and model reduction for stochastic reaction networks extended by Kang and Kurtz (2012), we approximate the chemical network in the heat shock response model.
Results
Scaling the species numbers and the rate constants by powers of the scaling parameter, we embed the model into a one-parameter family of models, each of which is a continuous-time Markov chain. Choosing an appropriate set of scaling exponents for the species numbers and for the rate constants satisfying balance conditions, the behavior of the full network in the time scales of interest is approximated by limiting models in three time scales. Due to the subset of species whose numbers are either approximated as constants or are averaged in terms of other species numbers, the limiting models are located on lower dimensional spaces than the full model and have a simpler structure than the full model does.
Conclusions
The goal of this paper is to illustrate how to apply the multiscale approximation method to the biological model with significant complexity. We applied the method to the heat shock response model involving 9 species and 18 reactions and derived simplified models in three time scales which capture the dynamics of the full model. Convergence of the scaled species numbers to their limit is obtained and errors between the scaled species numbers and their limit are estimated using the central limit theorem.
Keywords
Multiscale Markov chains Chemical reaction Reaction networks Heat shockBackground
Stochasticity may play an important role in biochemical systems. For example, stochasticity may be beneficial to give variability in gene expression, to produce population heterogeneity, and to adjust or respond to fluctuations in environment [1]. We are interested in local dynamics of biochemical networks involving some species with a small number of molecules so that the system is assumed to be well-mixed and relative fluctuations of small species numbers may play a role in the system dynamics.
The conventional stochastic model for the well-stirred biochemical network is based on the chemical master equation. The chemical master equation governs the evolution of the probability density of species numbers and is expressed as the balanced equation between influx and outflux of the probability density. When the biochemical network involves many species or bimolecular reactions, it is rarely possible to obtain an exact solution of the master equation in a closed form. Instead of searching for the solution of the master equation, stochastic simulation algorithms are used to obtain the temporal evolution of the species numbers. For example, Gillespie’s Stochastic Simulation Algorithm (SSA, or the direct method) is well known [2, 3] and provides a realization of the exact trajectory of the sample path for the species numbers. As the biochemical network has more species and reactions, SSA becomes computationally expensive and more efficient algorithms were suggested by many authors [4–6]. The detailed review of stochastic simulation methods, stochastic approximations, and hybrid simulation methods is given in [7]. For models with well-separated time scales, numerous authors suggested stochastic simulation algorithms for biochemical reaction networks by assuming that “fast” subnetworks have reached a “partial equilibrium” [6] or a “quasi-steady state” [4]. Using these assumptions, the approximate stochastic simulation algorithms involve a reduced number of species or reactions.
On the other hand, Ball et al. [8] described the state of the biochemical reaction network in the well-stirred system directly using stochastic equations for species numbers, and suggested an approximation of the reaction network via limiting models derived using different scalings for the species numbers and for the reaction rate constants. Kang and Kurtz [9] extended this multiscale approximation method and gave a systematic way to obtain limiting models in the time scales of interest. Conditions are given to help identify appropriate values for a set of scaling exponents which determine the time scale of each species and reaction. Using this method, nonstationary behavior of biochemical systems can be analyzed. Moreover, application of the method is flexible in the sense that the method does not require the exact parameter values but gives approximations valid for a range of parameter values. More recently, Crude et al. [10] also proposed a reduction method to derive simplified models with preserving stochastic properties and with key parameters using averaging and hybrid simplification.
The multiscale approximation method in [9] requires consideration of magnitude of both species numbers and rate constants of the reactions involving the corresponding species. When a moderately fast reaction involves two species, one with a small number of molecules and the other with a large number of molecules, the effects of this reaction on these species are different. Net molecule changes of species with large numbers due to the reaction is less noticeable than those of species with small numbers. Therefore, though the same reaction governs these species, their time scales may be different from each other. Letting N_{0} be a fixed constant and choosing a large value for N_{0}, for example N_{0}=100, we express magnitudes of species numbers and reaction rate constants in terms of powers of N_{0} with different scaling exponents. For instance, 1 to 10molecules are expressed as $1\times {N}_{0}^{0}$ to $10\times {N}_{0}^{0}\text{molecules}$, 500 to 800molecules are rewritten as 5×N_{0} to 8×N_{0}molecules, and 0.0002 sec becomes $2\times {N}_{0}^{-2}sec$. Assuming N_{0} is large, we replace N_{0}by a large parameter N and stochastic equations for species numbers are expressed in terms of N. Then, N is an analogue of 1/ε where ε is a small parameter in perturbation theory.
A specific time scale of interest is expressed in terms of a power of N, and its exponent contributes to reaction rates due to change of variables in time. For each species (or linear combination of species), we compare a power of N for the species number and those for reaction rates involving this species. Consider a case when the power for the species number is larger than those for the rates of all reactions where the species is involved. Then net molecule changes due to the reactions are not large enough to be noticeable in this time scale, and the species number is approximated as constant. Next, consider a case when the power for the species number is smaller than those for some reaction rates involving the species. In this case, the species number fluctuates very rapidly due to the fast reactions in this time scale, and the averaged behavior of the species number can be described in terms of other species numbers. The method of averaging is similar to approximation of one variable in terms of others using a quasi-steady state assumption. Last, when the power for the species number is equal to those for the rates of reactions where the species is involved, the scaled species number is approximated by a nondegenerate limit describing nonstationary behavior of the species number in the specific time scale of interest. The limit could be described in various kinds of variables: a continuous time Markov chain, a deterministic model given by a system of ordinary differential equations, or a hybrid model with both discrete and continuous variables. Since some of the scaled species numbers are approximated as constants or the averaged behavior of some species numbers is expressed in terms of other variables, dimension of species in the approximation of the biochemical network is reduced.
In the multiscale approximation method, scaling exponents for species numbers and for reaction rate constants are not uniquely determined, since the choice of values for the exponents is flexible. For example, 0.005 sec can be expressed as $0.5\times {N}_{0}^{-1}$ or $5\times {N}_{0}^{-1.5}$ when N_{0}=100. The goal in this method is to find an appropriate set of scaling exponents to obtain a nondegenerate limit of the scaled species numbers. Orders of magnitude of species numbers in the propensities affect reaction rates, and reaction rates contribute to determining rates of net molecule changes of the species involved in the reactions. Since species numbers and reaction rates interact, it is not easy to determine scaling exponents for all species numbers and reaction rate constants so that the limits of the scaled species numbers become balanced.
Kang and Kurtz [9] introduced balance conditions for the scaling exponents, which help to determine values for a set of exponents. The key idea in these conditions is that for each species (or linear combination of species) the maximum of scaling exponents in the rates of the reactions where this species is produced should be the same as that in the rates of the reactions where this species is consumed, i.e. maximal production and consumption rates of the species should be balanced in the order of magnitude. In case the maximums of scaling exponents for productions and consumptions are not balanced for some species, an increase or decrease of the scaled species number can be described by its limit during a certain time period. However after this time period, the scaled species number will either become zero or blow up to infinity. Therefore, if some of the scaled species numbers are not balanced due to a difference between orders of magnitude of production and consumption rates, the chosen scaling is valid up to a certain time scale. After this time scale, we need to choose different values for scaling exponents. In each time scale of interest we derive a limiting model including a subset of species and reactions, which is used to approximate the state of the full reaction network. The multiscale approximation method is applicable in case some of reaction rates are not known accurately, since the chosen scaling is applicable in some ranges of the parameters. Therefore, based on the behavior of the limiting models, we may be able to estimate behavior for a range of parameter values without performing a huge number of stochastic simulations.
The reduced network in the early stage has very simple structure without any bimolecular reactions, and all reactions involved are either production from a source or conversion. Moreover, the reduced network is well separated into two due to independence of S_{8}from S_{2}and S_{3}.
where a species over the arrow accelerates or inhibits the corresponding reaction. The reaction does not change this species number, but the propensity of the corresponding reaction is a function of this species number. In this time scale, conversion between S_{2} and S_{3} occurs very frequently and S_{2}and S_{3}play a role as a single “virtual” species rather than separate species. The species numbers of S_{23} and S_{8}are described as two independent birth processes and the species number of S_{7} is governed by conversion. In this time scale, the species number of S_{8}is normalized and treated as a continuous variable. The interesting thing is that the behavior of the species S_{8} which rapidly increases in time is well approximated in both first and second time scales.
As we see in Figure 1, the full network involves reactions with more than two reactants or products. However, all reactions in the reduced network at the times of order 10,000 sec consist of either production or degradation of each species, though most of the species (6 species out of 9) are involved in the reduced model. As in the medium stage of time period, S_{2}and S_{3}play a role as a single species. In the early and medium stages of time period propensities are in a form following the law of mass action, while in the late stage of time period the propensity for degradation of S_{23} is a nonlinear function of the species numbers similar to the reaction rate appearing in the Michaelis-Menten approximation for an enzyme reaction. The nonlinear function involves the species numbers of S_{23}, S_{8}, and S_{9}, which come from averaging of the species numbers of S_{2}and S_{6}which fluctuate rapidly in the third time scale. Similarly, the propensity of catalytic degradation of S_{8} is not proportional to the number of molecules of S_{8}.
In the late stage of time period of order 10,000 sec, we study the error between the scaled species numbers and their limit analytically using the central limit theorem derived in [14] and show that the error is of order 10^{−1}.
Methods
- 1.Write a chemical reaction network involving s _{0}species and r _{0} reactions in the form of$\sum _{i=1}^{{s}_{0}}{\nu}_{\mathit{ik}}{S}_{i}\to \sum _{i=1}^{{s}_{0}}{\nu}_{\mathit{ik}}^{\prime}{S}_{i},\phantom{\rule{2em}{0ex}}k=1,\cdots \phantom{\rule{0.3em}{0ex}},{r}_{0},$
where ν_{ ik } and ${\nu}_{\mathit{ik}}^{\prime}$ are nonnegative integers. Rearrange the reactions so that the reaction rate constants are decreasing monotonically as k gets large.
- 2.
Derive a system of stochastic equations for species numbers.
- (a)Letting X _{ i }(t) be the number of molecules of species S _{ i }at time t, the corresponding stochastic equation is${X}_{i}\left(t\right)={X}_{i}\left(0\right)+\sum _{k=1}^{{r}_{0}}{R}_{k}^{t}\left({\lambda}_{k}\right(X\left)\right)({\nu}_{\mathit{ik}}^{\prime}-{\nu}_{\mathit{ik}}),\phantom{\rule{1em}{0ex}}i=1,\cdots {s}_{0},$
where${R}_{k}^{t}(\xb7)$ counts the number of times that the k th reaction occurs up to time t.
- (b)
λ _{ k }(x) is determined by a stochastic version of mass action kinetics, and is expressed as a product of the rate constant and the numbers of molecules of reactants. If the k th reaction is second-order ($\sum _{i=1}^{{s}_{0}}{\nu}_{\mathit{ik}}=2$) with different types of reactants, ${\lambda}_{k}\left(x\right)={\kappa}_{k}^{\prime}{x}_{p}{x}_{q}$. When the reactants are two molecules of the same species, ${\lambda}_{k}\left(x\right)={\kappa}_{k}^{\prime}{x}_{p}({x}_{p}-1)$.
- (a)
- 3.
Derive a system of stochastic equations for the normalized species numbers after a time change, Z ^{N,γ}(t).
- (a)
In the equation for X _{ i }(t) obtained in Step 2 (a), replace X _{ i }by ${Z}_{i}^{N,\gamma}$ and divide reaction terms by N ^{ α } _{ i }. In the k th reaction term, put N ^{γ + ρ} _{ k } in the propensity and replace λ _{ k }(X) by ${\widehat{\lambda}}_{k}\left({Z}^{N,\gamma}\right)$. Then, we have
- (b)
In the equation in Step 3 (a), ${\rho}_{k}={\beta}_{k}+\sum _{j=1}^{{s}_{0}}{\alpha}_{j}{\nu}_{\mathit{jk}}$.
- (c)In the most reactions, ${\widehat{\lambda}}_{k}$ is obtained by replacing ${\kappa}_{k}^{\prime}$ by κ _{ k }in λ _{ k }. In case the k th reaction is second-order with reactants of the same species, ${\lambda}_{k}\left(x\right)={\kappa}_{k}^{\prime}{x}_{p}({x}_{p}-1)$ is replaced by ${\widehat{\lambda}}_{k}\left(z\right)={\kappa}_{k}{z}_{p}({z}_{p}-{N}^{-{\alpha}_{p}})$.$\begin{array}{l}{Z}_{i}^{N,\gamma}\left(t\right)={Z}_{i}^{N,\gamma}\left(0\right)+{N}^{-{\alpha}_{i}}\sum _{k=1}^{{r}_{0}}{R}_{k}^{t}\left({N}^{\gamma +{\rho}_{k}}{\widehat{\lambda}}_{k}\right({Z}^{N,\gamma}\left)\right)\\ \phantom{\rule{5em}{0ex}}\times ({\nu}_{\mathit{ik}}^{\prime}-{\nu}_{\mathit{ik}}),\phantom{\rule{1em}{0ex}}i=1,\cdots {s}_{0}.\end{array}$
- (a)
- 4.
Write a set of species balance equations and their time-scale constraints.
- (a)
Define ${\Gamma}_{i}^{+}$ and ${\Gamma}_{i}^{-}$ as subsets of reactions where the species number of S _{ i }increases or decreases every time the reaction occurs. Comparing ρ _{ k }’s for $k\in {\Gamma}_{i}^{+}$ and those for $k\in {\Gamma}_{i}^{-}$, set the balance equations as
- (b)
Time-scale constraints are given as
$\underset{k\in {\Gamma}_{i}^{+}}{max}{\rho}_{k}=\underset{k\in {\Gamma}_{i}^{-}}{max}{\rho}_{k},\phantom{\rule{2em}{0ex}}i=1,\cdots \phantom{\rule{0.3em}{0ex}},{s}_{0}.$$\gamma \le \underset{k\in {\Gamma}_{i}^{+}\cup {\Gamma}_{i}^{-}}{max}{\rho}_{k},\phantom{\rule{2em}{0ex}}i=1,\cdots \phantom{\rule{0.3em}{0ex}},{s}_{0}.$
- (a)
- 5.
Find a minimum set of linear combinations of species whose maximum of collective production (or consumption) rates may be different from that of one of any species. We construct a minimum set of linear combinations of species by selecting a linear combination of species if any reaction term involving the species consisting of the linear combination is canceled in the equation for the linear combination of species.
- 6.
For each selected linear combination of species, write a collective species balance equation and its time-scale constraint. They are obtained similarly to the ones in Step 4 using subsets of reactions where the number of molecules of linear combinations of species either increases or decreases instead of using ${\Gamma}_{i}^{+}$ and ${\Gamma}_{i}^{-}$.
- 7.Select a large value for N _{0}and choose an appropriate set of α _{ i }’s and β _{ k }’s so that
- (a)
the species number X _{ i }and the reaction rate constant ${\kappa}_{k}^{\prime}$ are approximately of orders ${N}_{0}^{{\alpha}_{i}}$ and ${N}_{0}^{{\beta}_{k}}$;
- (b)
the normalized species number ${Z}_{i}^{N,\gamma}$ and the scaled reaction rate constant κ _{ k }are of order 1;
- (c)
most of the balance equations obtained in Steps 4 and 6 are satisfied;
- (d)
β _{ k }’s are monotone decreasing among each class of reactions which have the same number of molecules of reactants.
- (a)
- 8.
Plugging the chosen values for α _{ i }’s and β _{ k }’s in the time-scale constraints obtained in Steps 4 and 6, compute an upper bound (denoted as γ _{0}) for a time-scale exponent. Then, the chosen set of exponents α _{ i }’s and β _{ k }’s can be used for γ satisfying γ≤γ _{0}. For γ>γ _{0}, select another set of exponents α _{ i }’s and β _{ k }’s using Steps 7 and 8.
- 9.Using each set of values for α _{ i }’s and β _{ k }’s, identify a natural time scale exponent of each species (denoted as γ _{ i } for species S _{ i }) so that γ _{ i } satisfies$\underset{k\in {\Gamma}_{i}^{+}\cup {\Gamma}_{i}^{-}}{max}({\gamma}_{i}+{\rho}_{k})\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}=\phantom{\rule{0.3em}{0ex}}\phantom{\rule{0.3em}{0ex}}{\alpha}_{i},\phantom{\rule{2em}{0ex}}i=1,\cdots \phantom{\rule{0.3em}{0ex}},{s}_{0}.$
We collect γ_{ i }’s with the same values, whose species are in the same time scales in the approximation.
- 10.
Modify α _{ i }’s and β _{ k }’s so that the conditions in Step 7 are satisfied and that γ _{ i }’s are divided into appropriate number of values, which gives the number of time scales, N ^{ γ }=N ^{ γ } _{ i }, we are interested in.
- 11.For each chosen γ, derive a limiting equation for each species S _{ i }with γ _{ i }=γ. Using the stochastic equation obtained in Step 3 (a), we let N go to infinity.
- (a)
For $k\in {\Gamma}_{i}^{+}\cup {\Gamma}_{i}^{-}$, the k th reaction term converges to zero if α _{ i }>γ + ρ _{ k }.
- (b)
If α _{ i }=γ + ρ _{ k }, the k th reaction term appears as a limit in the limiting equation. The limit of the k th reaction term is discrete if α _{ i }=0, while it is a continuous variable with the limit of its propensity if α _{ i }>0.
- (c)
There is no k satisfying α _{ i }<γ + ρ _{ k }in the equation for species S _{ i }with γ=γ _{ i }due to the definition of γ _{ i }given in Step 9.
- (a)
- 12.In the limiting equation for each species S _{ i }with γ _{ i }=γ, we approximate propensities in the reaction terms. Suppose that the normalized species number for S _{ j }appears in the propensities.
- (a)
If γ _{ j }>γ, the limit of the normalized species number for S _{ j }is its initial value.
- (b)
If γ _{ j }=γ, the limit of the normalized species number for S _{ j }appears as a variable in the propensities in the limiting equation.
- (c)
If γ _{ j }<γ, the limit of the normalized species number for S _{ j }is expressed as a function of the limits of the normalized species numbers for S _{ i }with γ _{ i }=γ. The function for S _{ j }is obtained by dividing the equation for S _{ j }by ${N}^{\underset{k\in {\Gamma}_{j}^{+}\cup {\Gamma}_{j}^{-}}{max}(\gamma +{\rho}_{k})-{\alpha}_{j}}$ and letting N go to infinity.
- (a)
- 13.
If a limiting model is not closed, consider limiting equations for some linear combinations of species selected in Step 5 whose natural time scale exponents are equal to the chosen γ.
The method for multiscale approximation described above can be applied to general chemical reaction networks containing different scales in species numbers and reaction rate constants. We can apply the method in case the rates of chemical reactions are determined by law of mass action and when there is no species whose number is either zero or infinity at all times. As given in [9], in the reaction network involving ∅→S_{1}, ∅→S_{2}, ∅→S_{3}, S_{1} + S_{2}→∅, and S_{1} + S_{3}→∅, convergence of the limit for the scaled species numbers may not be guaranteed at some time scales. Suppose that production rate of S_{1} is larger than that of S_{2}but with the same order of magnitude, and that production rate of S_{3} is much smaller than those of S_{1}and S_{2}. Then, X_{1}(t) may blow up to infinity and X_{2}(t) may go to zero at some time scales. In this case, the method is not applicable.
Results and discussion
Model description
We analyze a heat shock response model of E. coli developed by Srivastava, Peterson, and Bentley [11]. The heat shock response model gives a simplified mechanism occurring in the E. coli to respond to high temperature. Heat causes unfolding, misfolding, or aggregation of proteins, and cells overcome the heat stress by producing heat shock proteins, which refold or degrade denatured proteins. In E. coli, σ^{32}factors play an important role in recovery from the stress under the high temperature. σ^{32}factors catalyze production of the heat shock proteins such as chaperon proteins and other proteases. In this model, J denotes a chaperon complex, FtsH represents a σ^{32}-regulated stress protein, and GroEL is a σ^{32}-mediated stress response protein.
σ^{32} factors are in three different forms, free σ^{32}protein, σ^{32} combined with RNA polymerase (E σ^{32}), and σ^{32} combined with a chaperon complex (σ^{32}-J). Under the normal situation without stress, most of the σ^{32} factors combine with chaperon complexes and form σ^{32}-J. A chaperon complex J keeps σ^{32}factors in an inactive form, and σ^{32}factors can directly respond to the stress by changing into different forms. When there exist σ^{32}factors combined with chaperon complexes, FtsH catalyzes degradation of σ^{32} factors. Thus, if enough σ^{32}-regulated stress proteins are produced, σ^{32}factors are degraded.
Species in the heat shock response model of E. coli and their initial values
X _{1} | = | # of S_{1} | σ^{32} mRNA | X_{1}(0) | = | 10 |
X _{2} | = | # of S_{2} | σ^{32} protein | X_{2}(0) | = | 1 |
X _{3} | = | # of S_{3} | E σ ^{32} | X_{3}(0) | = | 1 |
X _{4} | = | # of S_{4} | FtsH | X_{4}(0) | = | 93 |
X _{5} | = | # of S_{5} | GroEL | X_{5}(0) | = | 172 |
X _{6} | = | # of S_{6} | J | X_{6}(0) | = | 54 |
X _{7} | = | # of S_{7} | σ^{32}-J | X_{7}(0) | = | 7 |
X _{8} | = | # of S_{8} | Recombinant protein | X_{8}(0) | = | 50 |
X _{9} | = | # of S_{9} | J-Recombinant protein | X_{9}(0) | = | 0 |
Reactions in the heat shock response model of E. coli
Reaction | Transition | |
---|---|---|
R1 | $\varnothing \stackrel{\text{gene}}{\to}{S}_{8}$ | Recombinant protein synthesis |
R2 | S_{2}→S_{3} | Holoenzyme association |
R3 | S_{3}→S_{2} | Holoenzyme disassociation |
R4 | $\varnothing \stackrel{{S}_{1}}{\to}{S}_{2}$ | σ^{32} translation |
R5 | ${S}_{3}\stackrel{\text{gene}}{\to}{S}_{2}+{S}_{5}$ | GroEL synthesis |
R6 | ${S}_{3}\stackrel{\text{gene}}{\to}{S}_{2}+{S}_{4}$ | FtsH synthesis |
R7 | ${S}_{3}\stackrel{\text{gene}}{\to}{S}_{2}+{S}_{6}$ | J-production |
R8 | S_{7}→S_{2} + S_{6} | σ^{32}-J-disassociation |
R9 | S_{2} + S_{6}→S_{7} | σ^{32}-J-association |
R10 | S_{6} + S_{8}→S_{9} | Recombinant protein-J association |
R11 | S_{8}→∅ | Recombinant protein degradation |
R12 | S_{9}→S_{6} + S_{8} | Recombinant protein-J disassociation |
R13 | $\varnothing \stackrel{\text{gene}}{\to}{S}_{1}$ | σ^{32} transcription |
R14 | S_{1}→∅ | σ^{32} mRNA decay |
R15 | ${S}_{7}\stackrel{{S}_{4}}{\to}{S}_{6}$ | σ^{32} degradation |
R16 | ${S}_{5}\to \varnothing $ | GroEL degradation |
R17 | ${S}_{6}\to \varnothing $ | J-disassociation |
R18 | ${S}_{4}\to \varnothing $ | FtsH degradation |
Stochastic reaction rate constants in the heat shock response model of E. coli
Rates | Rates | ||
---|---|---|---|
${\kappa}_{1}^{\prime}$ | 4.00×10^{0} | ${\kappa}_{10}^{\prime}$ | 3.62×10^{−4} |
${\kappa}_{2}^{\prime}$ | 7.00×10^{−1} | ${\kappa}_{11}^{\prime}$ | 9.99×10^{−5} |
${\kappa}_{3}^{\prime}$ | 1.30×10^{−1} | ${\kappa}_{12}^{\prime}$ | 4.40×10^{−5} |
${\kappa}_{4}^{\prime}$ | 7.00×10^{−3} | ${\kappa}_{13}^{\prime}$ | 1.40×10^{−5} |
${\kappa}_{5}^{\prime}$ | 6.30×10^{−3} | ${\kappa}_{14}^{\prime}$ | 1.40×10^{−6} |
${\kappa}_{6}^{\prime}$ | 4.88×10^{−3} | ${\kappa}_{15}^{\prime}$ | 1.42×10^{−6} |
${\kappa}_{7}^{\prime}$ | 4.88×10^{−3} | ${\kappa}_{16}^{\prime}$ | 1.80×10^{−8} |
${\kappa}_{8}^{\prime}$ | 4.40×10^{−4} | ${\kappa}_{17}^{\prime}$ | 6.40×10^{−10} |
${\kappa}_{9}^{\prime}$ | 3.62×10^{−4} | ${\kappa}_{18}^{\prime}$ | 7.40×10^{−11} |
approximates the network at the times of order 10,000 sec. Detailed derivation is given in the later sections. Note that it is possible to identify different numbers of time scales depending on the scaling of the species numbers and reaction rate constants. In the heat shock response model of E. coli, it is possible to obtain approximate models with two or four time scales. However, if the number of time scales are too many, the limiting model in each time scale may involve one species and a few number of reactions and the model in this case may not be interesting to consider.
Derivation of the scaled models
so that ${X}_{i}^{{N}_{0}}\left(0\right)={X}_{i}\left(0\right)$ and $\underset{N\to \infty}{lim}{Z}_{i}^{N}\left(0\right)=\underset{N\to \infty}{lim}$${N}^{-{\alpha}_{i}}{X}_{i}^{N}\left(0\right)={N}_{0}^{-{\alpha}_{i}}{X}_{i}\left(0\right)$.
For each reaction, ρ_{ k }is given in terms of α_{ i }and β_{ k } in the Additional file 1: Table S1.
We are interested in dynamics of species numbers ${Z}_{2}^{N}\left(t\right)$ and ${Z}_{3}^{N}\left(t\right)$ in various stages of time period. In the early stage of time period, normalized species numbers of S_{2} and S_{3} are very close to their scaled initial values, since these species numbers have not changed yet. In the medium stage of time period, the normalized species numbers of S_{2}and S_{3} are asymptotically equal to non-constant limits. In the late stage of time period, the normalized species numbers of S_{2} and S_{3}fluctuate very rapidly and their averaged behavior is captured in terms of some function of other species numbers.
Then, ${Z}_{i}^{N,\gamma}\left(t\right)$ gives a normalized species number at the times of order N^{ γ }. A natural time scale of S_{ i }is the time when ${Z}_{i}^{N,\gamma}\left(t\right)$ has a nonzero finite limit which is not constant and of order 1.
where N^{ γ }in each propensity comes from the change of the time variable. Here, the initial values may depend on γ, since we can choose different values for α_{ i }for each γ due to changes in order of magnitude of species numbers in time. The stochastic equations after scaling and a time change for all species are given in the Additional file 1: Section 1.
Balance conditions
Inequalities in (14) mean that if maximal production and consumption rates are not balanced either for S_{2} or S_{3}, the chosen set of values for scaling exponents can be used to approximate the dynamics of the full network up to times of order N^{ u }_{2} or N^{ u }_{3}. For times later than those of order N^{ u }_{2}or N^{ u }_{3}, we need to choose another set of values for scaling exponents based on the balance equations. We call the balance equation and the time-scale constraint for each species as the species balance condition. If either (12a) or (??) is satisfied, we say that the species balance condition for S_{2} is satisfied.
Similarly to the time-scale constraint in the species balance condition, (18) implies that if maximal collective production and consumption rates for S_{23}are not balanced, our choice of values for scaling exponents are valid up to times of order N^{ u }_{23}.
Balance equations and time-scale constraints for each species and for each collective species chosen
Balance equations | Time-scale constraints | |
---|---|---|
S _{1} | ρ_{13}=ρ_{14} | $\gamma \le {\alpha}_{1}-max({\rho}_{13},{\rho}_{14})$ |
S _{2} | $max({\rho}_{3},{\rho}_{4},{\rho}_{5},{\rho}_{6},{\rho}_{7},{\rho}_{8})=max({\rho}_{2},{\rho}_{9})$ | $\gamma \le {\alpha}_{2}-max({\rho}_{2},{\rho}_{3},{\rho}_{4},{\rho}_{5},{\rho}_{6},{\rho}_{7},{\rho}_{8},{\rho}_{9})$ |
S _{3} | ${\rho}_{2}=max({\rho}_{3},{\rho}_{5},{\rho}_{6},{\rho}_{7})$ | $\gamma \le {\alpha}_{3}-max({\rho}_{2},{\rho}_{3},{\rho}_{5},{\rho}_{6},{\rho}_{7})$ |
S _{4} | ρ_{6}=ρ_{18} | $\gamma \le {\alpha}_{4}-max({\rho}_{6},{\rho}_{18})$ |
S _{5} | ρ_{5}=ρ_{16} | $\gamma \le {\alpha}_{5}-max({\rho}_{5},{\rho}_{16})$ |
S _{6} | $max({\rho}_{7},{\rho}_{8},{\rho}_{12},{\rho}_{15})=max({\rho}_{9},{\rho}_{10},{\rho}_{17})$ | $\gamma \le {\alpha}_{6}-max({\rho}_{7},{\rho}_{8},{\rho}_{9},{\rho}_{10},{\rho}_{12},{\rho}_{15},{\rho}_{17})$ |
S _{7} | ${\rho}_{9}=max({\rho}_{8},{\rho}_{15})$ | $\gamma \le {\alpha}_{7}-max({\rho}_{8},{\rho}_{9},{\rho}_{15})$ |
S _{8} | $max({\rho}_{1},{\rho}_{12})=max({\rho}_{10},{\rho}_{11})$ | $\gamma \le {\alpha}_{8}-max({\rho}_{1},{\rho}_{10},{\rho}_{11},{\rho}_{12})$ |
S _{9} | ρ_{10}=ρ_{12} | $\gamma \le {\alpha}_{9}-max({\rho}_{10},{\rho}_{12})$ |
S_{2} + S_{3} + S_{7} | ρ_{4}=ρ_{15} | $\gamma \le max({\alpha}_{2},{\alpha}_{3},{\alpha}_{7})-max({\rho}_{4},{\rho}_{15})$ |
S_{2} + S_{3} | $max({\rho}_{4},{\rho}_{8})={\rho}_{9}$ | $\gamma \le max({\alpha}_{2},{\alpha}_{3})-max({\rho}_{4},{\rho}_{8},{\rho}_{9})$ |
S_{2} + S_{7} | $max({\rho}_{3},{\rho}_{4},{\rho}_{5},{\rho}_{6},{\rho}_{7})=max({\rho}_{2},{\rho}_{15})$ | $\gamma \le max({\alpha}_{2},{\alpha}_{7})-max({\rho}_{2},{\rho}_{3},{\rho}_{4},{\rho}_{5},{\rho}_{6},{\rho}_{7},{\rho}_{15})$ |
S_{6} + S_{7} + S_{9} | ρ_{7}=ρ_{17} | $\gamma \le max({\alpha}_{6},{\alpha}_{7},{\alpha}_{9})-max({\rho}_{7},{\rho}_{17})$ |
S_{6} + S_{7} | $max({\rho}_{7},{\rho}_{12})=max({\rho}_{10},{\rho}_{17})$ | $\gamma \le max({\alpha}_{6},{\alpha}_{7})-max({\rho}_{7},{\rho}_{10},{\rho}_{12},{\rho}_{17})$ |
S_{6} + S_{9} | $max({\rho}_{7},{\rho}_{8},{\rho}_{15})=max({\rho}_{9},{\rho}_{17})$ | $\gamma \le max({\alpha}_{6},{\alpha}_{9})-max({\rho}_{7},{\rho}_{8},{\rho}_{9},{\rho}_{15},{\rho}_{17})$ |
S_{8} + S_{9} | ρ_{1}=ρ_{11} | $\gamma \le max({\alpha}_{8},{\alpha}_{9})-max({\rho}_{1},{\rho}_{11})$ |
Based on species and collective species balance equations in Table 4, we choose appropriate values for α_{ i }’s and β_{ k }’s so that most of the balance equations are satisfied. If some of the balance equations are not satisfied, corresponding time-scale constraints give a range of γ where the chosen α_{ i }’s and β_{ k }’s are valid. The time-scale constraint, γ≤γ_{0}, implies that the set of scaling exponents α_{ i }’s and β_{ k }’s chosen is appropriate only up to time whose order of magnitude is equal to N^{ γ }_{0}. For the times larger than O(N^{ γ }_{0}), we need to choose a different set of values for the scaling exponents, α_{ i }’s. Assuming that reaction rate constants do not change in time and that the species numbers vary in time, we in general use one set of β_{ k }’s for all time scales and may use several sets of α_{ i }’s. A large change of the species numbers in time requires different α_{ i }’s in different time scales. For the heat shock model we identify three different time scales as we will see in the section of limiting models in three time scales, and α_{1}, α_{2}, α_{3}, α_{8}, and α_{9} may depend on the time scale. α_{4}, α_{5}, α_{6}, and α_{7} are the same for all time scales.
Then, the first set of scaling exponents with α_{1}=1 and α_{2}=α_{3}=0 is valid only when γ≤0. Next, based on the fact that X_{2}(t)≈O(10) and X_{3}(t)≈O(10) in the medium stage of time period, we choose α_{2}=α_{3}=0 for γ>0. At this stage of time period, we set ${X}_{1}\left(t\right)=O\left(10\right)\approx O\left({N}_{0}^{{\alpha}_{1}}\right)$ with α_{1}=0. Then, (12a) and (12b) are satisfied but not (16). The condition (18) gives γ≤1, and the second set of scaling exponents with α_{1}=α_{2}=α_{3}=0 is valid when γ≤1. Finally, we set α_{1}=0 and α_{2}=α_{3}=1 for γ>1 based on the fact that the numbers of molecules of S_{2}and S_{3} grow in time and are of order 100. Then, (12a), (12b), and (16) are all satisfied, and the third set of scaling exponents with α_{1}=0 and α_{2}=α_{3}=1 can be used for γ>1.
The three sets of values for the scaling exponents chosen are given in the Additional file 1: Table S4. With chosen values for the scaling exponents, we check whether each balance equation is satisfied and give a time-scale constraint in the Additional file 1: Table S6 in case the balance equation is not satisfied. Different choices of α_{ i }’s and β_{ k }’s from the ones in the Additional file 1: Table S4 give different limiting models. As long as the chosen values for α_{ i }’s and β_{ k }’s satisfy balance conditions, the limiting model will describe nontrivial behavior of the species numbers which are nonzero and finite in the specific time of interest.
Limiting models in three time scales
In the heat shock response model of E. coli, we identify a time scale of interest using the chosen set of scaling exponents and derive a limiting model which approximates dynamics of the full chemical reaction network. Each limiting model involves a subset of species and reactions, and gives features of the full network during the time interval of interest.
where Γ i + denotes the collection of reactions where the species number of S_{ i } increases every time the reaction occurs. Similarly, Γ i− is the subset of reactions where the species number of S_{ i }decreases every time the reaction occurs. In (19), the left-side term is the maximal order of magnitude of rates of reactions involving S_{ i }and the right-side term is the order of magnitude of the species number for S_{ i }. If times are earlier than those of order N^{ γ }_{ i }(γ<γ_{ i }), fluctuations of species number of S_{ i } due to the reactions involving S_{ i }are not noticeable compared to magnitude of the species number of S_{ i }. Then, the species number of S_{ i } is approximated as its initial value. In the times of order N^{ γ }_{ i }(γ=γ_{ i }), changes of species number of S_{ i } due to the reactions and the species number of S_{ i } are similar in magnitude and behavior of the species number of S_{ i }is described by its nondegenerate limit. If times are later than those of order N^{ γ }_{ i }(γ>γ_{ i }), the species number of S_{ i } fluctuates very rapidly due to the reactions involving S_{ i } compared to the magnitude of the species number of S_{ i }. Then, the averaged behavior of the species number of S_{ i }is approximated by some function of other species numbers. Note that γ_{ i } depends on α_{ i }’s and β_{ k }’s, and the time scale of the i th species may change if we use several sets of α_{ i }’s.
All values of α_{ i }’s and ρ_{ k }’s for three scalings which are used to derive limiting models are given in the Additional file 1: Table S4. The equations for normalized species numbers and the equation for ${Z}_{23}^{N,\gamma}$ which are used later in this section are given in the Additional file 1: Section 1 and Section 2, respectively. When we derive limiting models in three time scales, boundedness of the normalized species numbers is required. For first two time scales, we define stopping times so that the normalized species numbers are bounded up to those times. For the last time scale, we proved stochastic boundedness of some normalized species numbers in a finite time interval. For more details, see Additional file 1: Section 5.
and we get γ_{2}=0. Similarly, we get γ_{3}=γ_{8}=0.
and we get γ_{1}=2. Similarly, we get γ_{ i }>0 for i=4,5,6,7,9. Among all natural time scale exponents of species, we choose the smallest one, γ=0, and set t∼O(N^{0})=O(1) as the first time scale we are interested in. Since γ_{1}>0, ${Z}_{1}^{N,0}\left(t\right)\to {Z}_{1}^{0}\left(0\right)$ as N→∞. Similarly, ${Z}_{i}^{N,0}\left(t\right)\to {Z}_{i}^{0}\left(0\right)={N}_{0}^{-{\alpha}_{i}}{X}_{i}\left(0\right)$ for i=4,5,6,7,9 as N→∞. To sum up, in this time scale with γ=0, the species numbers of S_{ i }’s for i=1,4,5,6,7,9 change more slowly than other species numbers, and the species numbers with slow time scales are approximated as constant.
Similarly, we get a limiting model with ${Z}_{2}^{0}$, ${Z}_{3}^{0}$, and ${Z}_{8}^{0}$ for γ=0 as given in (3).
and we get γ_{6}=1. Similarly, we get γ_{7}=γ_{8}=1, γ_{ i }<1 for i=2,3, and γ_{ i }>1 for i=1,4,5,9. We already get the temporal behavior of species numbers of S_{2}, S_{3}, and S_{8} through the limiting model when γ=0. Thus, we set t∼O(N^{1}) as the second time scale we are interested in, and derive a limiting model for S_{6}, S_{7}, and S_{8} when γ=1. Note that species S_{8} is involved in the limiting models for both γ=0 and γ=1, since we use different sets of scaling exponents in these models. For i=1,4,5,9${Z}_{i}^{N,1}\left(t\right)\to {Z}_{i}^{1}\left(0\right)$ as N→∞, since γ_{ i }>1. Thus, in the 12th and 15th reaction terms in (24), ${Z}_{9}^{N,1}\left(s\right)\to {Z}_{9}^{1}\left(0\right)$ and ${Z}_{4}^{N,1}\left(s\right)\to {Z}_{4}^{1}\left(0\right)$ as N→∞. Since the propensities of the 8th, 9th, and 17th reaction terms in (24) are of order N^{γ−2}=N^{−1} for γ=1 and the species number of S_{6} is of order 1, these reaction terms go to zero as N→∞. In the 10th and 15th reaction terms in (24), ${Z}_{6}^{N,1}\left(s\right)$, ${Z}_{7}^{N,1}\left(s\right)$, and ${Z}_{8}^{N,1}\left(s\right)$ are asymptotically O(1) and converge to ${Z}_{6}^{1}\left(s\right)$, ${Z}_{7}^{1}\left(s\right)$, and ${Z}_{8}^{1}\left(s\right)$ as N→∞ since γ_{6}=γ_{7}=γ_{8}=1.
as N→∞.
In (30), note that ${R}_{12}^{t}\left({\kappa}_{12}{Z}_{9}^{1}\right(0\left)\right)=0$ since X_{9}(0)=0 as given in Table 1. Limiting equations for ${Z}_{7}^{N,1}$ and ${Z}_{8}^{N,1}$ can be derived similarly, and a limiting model with ${Z}_{23}^{1}$, ${Z}_{6}^{1}$, ${Z}_{7}^{1}$, and ${Z}_{8}^{1}$ for γ=1 is given in (4).
uniformly as N→∞.
as N→∞.
In (41), note that ${Z}_{9}^{2}\left(0\right)=0$ since X_{9}(0)=0 as given in Table 1.