### 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 shock### Background

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
*N*_{0} to 8×*N*_{0}molecules, and 0.0002 sec becomes
*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
*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 paper [9] included several simple examples of biochemical networks involving two to four species,
and derived limiting models in each time scale of interest. To apply this method,
more scaling exponents must be determined as the biochemical network involves more
species or reactions. Therefore, it is challenging to apply the method to complex
biochemical systems and to determine appropriate values for scaling exponents so that
the corresponding limiting models preserve important dynamical features of the full
system. One of the goals of this paper is to illustrate how to apply this method to
an example with significant complexity. In this paper, using a significantly complicated
biochemical network, we derive limiting models, show convergence of the scaled species
numbers to their limit, and estimate the error analytically between the scaled species
numbers and their limit. We analyze a heat shock response model of *Escherichia coli* (*E. coli*) developed by Srivastava, Peterson, and Bentley in [11]. The model involves 9 species and 18 reactions with significant complexity as shown
in Figure 1, and it has various time scales due to wide ranges of species numbers and reaction
rate constants. Because of various scales involved, this model has been used as an
example to show accuracy of the stochastic simulation algorithms which are developed
to increase computational efficiency using the multiscale nature of the chemical reaction
network [12,13]. Another version of a heat shock response model of *E. coli* is studied in [6] using an accelerated SSA that also exploits the multiscale nature of the system.

**Figure 1.** **A chemical reaction network in the heat shock response model of E. coli.** A dotted line represents the effect of the species acting as catalysts.

Applying the multiscale approximation method to the heat shock response model of *E. coli*, we derive limiting models in three time scales of our interests, which approximate
the full network given in Figure 1. Denote *∅* as species we are not interested in. Let *S*_{i} represent the *i*th species and *S*_{23}be addition of species *S*_{2}and *S*_{3}. *A*→*B* denotes a reaction where one molecule of species *A* is converted to one molecule of species *B*. In the early stage of time period of order 1 sec, we obtain the following reduced
network:

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}.

In the medium stage of time period of order 100 sec, the full network is reduced to

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.

In the late stage of time period of order 10,000 sec, we get a reduced network with more species involved than those in the previous time scales. However, the reduced network is still much simpler than the full network in Figure 1. At this time scale, we get

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

In the next several sections, we apply the multiscale approximation to the heat shock
response model of *E. coli* and derive the limiting models. The multiscale approximation method is described
in terms of the following steps so that the method can be applied to the general cases.

1. Write a chemical reaction network involving *s*_{0}species and *r*_{0} reactions in the form of

where *ν*_{ik} and
*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

where
*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 (

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
*N*^{α}_{i}. In the *k*th reaction term, put *N*^{γ + ρ}_{k} in the propensity and replace *λ*_{k}(*X*) by

(b) In the equation in Step 3 (a),

(c) In the most reactions,
*κ*_{k}in *λ*_{k}. In case the *k*th reaction is second-order with reactants of the same species,

4. Write a set of species balance equations and their time-scale constraints.

(a) Define
*S*_{i}increases or decreases every time the reaction occurs. Comparing *ρ*_{k}’s for

(b) Time-scale constraints are given as

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

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

(b) the normalized species number
*κ*_{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.

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

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*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.

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* go to infinity.

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.

Not only *σ*^{32}factors, but recombinant proteins also require chaperon complexes to form a complex
so that denatured protein can be fixed. Therefore, *σ*^{32}factors and recombinant proteins compete to bind chaperon complexes, and different
levels of binding affinity of recombinant proteins to chaperon complexes change the
evolution of the system state. In the model, we assume that *σ*^{32} factors and recombinant proteins have the same affinity to bind to chaperon complexes.
The system is sensitive to the amount and forms of *σ*^{32} factors: a small decrease of *σ*^{32}factors causes a large reduction of production of chaperon complexes and *σ*^{32}-regulated stress proteins, and the ratio of three different forms of *σ*^{32}factors determines system dynamics in the stress response [11]. The total initial number of molecules of *σ*^{32} factors in each cell is small [11] (also see initial values for *S*_{2}, *S*_{3}, and *S*_{7} which are 1, 1, and 7 in Table 1), and the stochastic model is appropriate to be considered.

**Table 1.** **Species in the heat shock response model of E. coli and their initial values**

The model involves 9 species and 18 reactions. Denote *s*_{0} as the number of species and *r*_{0} as the number of reactions. Let *X*(*t*) be a state vector whose *i*th component represents the number of molecules of species *S*_{i} at time *t* for *i*=1,⋯,*s*_{0}. Define a random process which counts the number of times that the *k*th reaction occurs by time *t* as

where *λ*_{k}(*X*) is the propensity of the *k*th reaction and the *Y*_{k}’s are independent unit Poisson processes. Therefore,
*λ*_{k}(·) gets large, the moment when
*ν*_{ik}(
*S*_{i} that are consumed (produced) in the *k*th reaction. Define *ν*_{k}(
*s*_{0}-dimensional vector whose *i*th component is *ν*_{ik}(
*X*(*t*) is given as

That is, species numbers at time *t* are expressed in terms of their initial values and sum of the number of times that
each reaction occurs multiplied by net molecule changes in the corresponding reaction.
In our model, the system of equations are derived using a set of reactions in Table
2 as:

**Table 2.** **Reactions in the heat shock response model of E. coli**

*k*th reaction, and their values from [11] are given in Table 3.

**Table 3.** **Stochastic reaction rate constants in the heat shock response model of E. coli**

We derive the limiting models in three time scales, which approximate a full network
in a certain time period involving a subset of species and reactions. In what follows,
*S*_{i} at some time scales depending on *γ*, and as *γ* gets larger the times are in the later stage. Note that the exponent *γ*in
*Z*_{i})^{γ} but it shows dependence of
*γ*. Let *κ*_{k} be a scaled reaction rate constant for the *k*th reaction. In the first time scale (when the times are in the early stage), the
subnetwork governed by

approximates the network when the times are of order 1 sec. Denote
*S*_{2}and *S*_{3}. In the second time scale (when the times are in the medium stage), the subnetwork
governed by

approximates the network at the times of order 100 sec. In the third time scale, set
the limit of the averaged scaled species numbers of fast-fluctuating species *S*_{2}, *S*_{3}, and *S*_{6} as

When the times are in a late stage, the subnetwork governed by

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

The stochastic equations given in Equations (2) describe temporal evolution of the
species numbers. For example, the equations for species *S*_{2}and *S*_{3} are

In Equation (6), species numbers of *S*_{2}and *S*_{3} are determined by the times when reactions occur and by the number of times that
reactions happen. On the other hand, reaction time and frequency are determined by
propensities which are some functions of species numbers. Therefore, reaction rates
and species numbers interact one another. Reaction rates vary from *O*(10^{−11}) to *O*(1) as we see in Table 3, and species numbers in this model are from *O*(1) to *O*(10^{4}) as we see later in the simulation of the full network. We express each species number
and rate constant in terms of powers of a common number with different weights on
exponents. Define *N*_{0}=100 as a fixed unitless constant used to express the magnitude of the species numbers
and the reaction rate constants. Define *α*_{i} for *i*=1,⋯,*s*_{0} and *β*_{k} for *k*=1,⋯,*r*_{0} as the scaling exponent for species *S*_{i} and for the reaction rate constant
*κ*_{k} is of order 1 and is determined so that
*β*_{6}=−1 so that the reaction rate is expressed as
*N*_{0} is large, we replace *N*_{0} by *N* and express the process as
*N*. Note that {*X*^{N}(*t*)} is a family of processes depending on *N* and
*N*=*N*_{0}. Then, the equation for

where
*N*=*N*_{0}. Since the numbers of molecules of species are in different orders of magnitude,
we scale the number of molecules of the *i*th species by *N*^{α}_{i} and set a normalized species number as

The *i*th species number may have different orders of magnitude at different times so *α*_{i} may have different values for different time scales. Now, we set the initial values
as

so that

Next, we scale the propensities of reactions by replacing
*N*^{β}_{k}*κ*_{k}. For example, consider the 9th reaction term in (6a). Replacing
*N*^{β}_{9}*κ*_{9},

For simplicity, set *ρ*_{9}=*β*_{9} + *α*_{2} + *α*_{6} and define a scaling exponent in the propensity of the *k*th reaction term as

where
*ν*_{ik} gives the number of molecules of species *S*_{i}consumed in the *k*th reaction. Then, (8) is rewritten as

Dividing (6a) by *N*^{α}_{2} and (6b) by *N*^{α}_{3} and scaling the propensities, we get

For each reaction, *ρ*_{k}is given in terms of *α*_{i}and *β*_{k} in the Additional file 1: Table S1.

**Additional file 1.** **Supplementary material for “A multiscale approximation in a heat shock response model
of E. coli.”** This is a supplementary material of the paper including calculations and tables.

Format: PDF Size: 201KB Download file

This file can be viewed with: Adobe Acrobat Reader

We are interested in dynamics of species numbers
*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.

We want to express the time scale of each species in terms of power of *N*. First, we express order of magnitude of a specific time period of interest as a
power of *N* with a time scale exponent *γ*. Applying a time change by replacing *t* by *N*^{γ}*t*in

Then,
*N*^{γ}. A *natural time scale* of *S*_{i}is the time when

Changing a time variable by replacing *t* by *N*^{γ}*t* in (9a) and (9b), the normalized species numbers of *S*_{2}and *S*_{3}after a time change satisfy

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

Our goal is to approximate dynamics of the full network in the heat shock response
model of *E. coli* in specific times of interest in terms of simplified subnetworks preserving significant
biological features. In each time period of interest, we obtain a nondegenerate limiting
model which is not equal to zero and does not blow up to infinity. In this section,
we introduce balance conditions which help us to choose appropriate values for the
scaling exponents *α*_{i}’s and *β*_{k}’s so that the limit is nonzero finite. For each time period of interest of order
*N*_{0}=100, we choose values for scaling exponents so that orders of magnitude of the species
number for *S*_{i} and the *k*th reaction rate constant are about

It is natural to choose *β*_{k}’s in monotone decreasing manner in *k*, since
*β*_{k}’s holds in each class of reactions. In other words, we choose *β*_{k}’s so that

Next, in order to make the normalized specie number
*S*_{i}should be in the same order of magnitude. If the order of magnitude of production
rate is larger than that of consumption, the normalized species number asymptotically
goes to infinity. In the opposite case, the normalized species number asymptotically
becomes zero. Therefore, for each species *S*_{i}, we set the balance equation for *α*_{i}’s and *β*_{k}’s so that the maximal exponent in the propensities of the reactions producing *S*_{i} is equal to that in the propensities of the reactions consuming *S*_{i}. For example, to obtain a balance equation for species *S*_{2}, we compare the scaling exponents in propensities of reactions involving *S*_{2}using (11a), and set the maximal exponents of production and consumption of *S*_{2} equal. Similarly, using (11b), we set the maximal exponents in the production rates
and the consumption rates of *S*_{3} equal. Then, the balance equations for species *S*_{2}and *S*_{3} are

If the maximal orders of magnitudes of production and consumption rates for *S*_{2} are different from each other, the species number of *S*_{2}should be large enough so that a difference between production and consumption of
*S*_{i} is not noticeable. In other words if *α*_{i}’s and *β*_{k}’s do not satisfy (12a), *α*_{2}should be at least as large as the scaling exponents located in all reaction terms
in (11a) to prevent the limit becoming zero or blowing up to infinity. Similarly,
in case (12b) is not satisfied, *α*_{3} should be at least as large as the scaling exponents located in the reaction terms
in (11b) to prevent the limit becoming zero or blowing up to infinity.

Solving (13) for *γ*, we get the following time-scale constraints:

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.

Even though species balance conditions for *S*_{2}and *S*_{3} are satisfied, the limit of the normalized species numbers for *S*_{2}or *S*_{3} may become degenerate. Consider addition of species *S*_{2}and *S*_{3} as a single virtual species, and compare the collective rates of production and consumption
of this species. Recall that *S*_{23} denotes addition of species *S*_{2}and *S*_{3}. Since production of one species is canceled by consumption of the other species,
maximal production rate of *S*_{23} may be different from that of *S*_{2}or *S*_{3}. Suppose that the maximal collective rates of production or consumption of *S*_{23} are slower than the maximal production or consumption rates of *S*_{2}and *S*_{3}. Also, suppose that the maximal collective rates of production and consumption of
the complex have different orders of magnitude. Then, a limit of the normalized species
number of *S*_{23}can be zero or infinity, even though the species balance conditions for *S*_{2} and *S*_{3} are satisfied. Therefore, we need an additional condition to obtain balance between
collective production and consumption rates for *S*_{23}. To obtain a balance equation for *S*_{23}, we unnormalize (11a) and (11b) by multiplying *N*^{α}_{2} and *N*^{α}_{3}, respectively. Adding the unnormalized equations for species *S*_{2}and *S*_{3} and dividing it by

Comparing the maximal scaling exponents of production and consumption of *S*_{23} in (15), a balance equation for *S*_{23}is given as

In case (16) is not satisfied, the order of magnitude of the species number for *S*_{23} should be larger than those of collective production and consumption rates so that
a difference between production and consumption is not noticeable. This gives

Solving (??) for *γ*, we get

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}.

We call (16) and (18) the *collective species balance condition* for *S*_{23}, that is, either (16) or (18) must hold. The species balance conditions for all species
and the collective species balance conditions for all positive linear combinations
of species should be satisfied to obtain a nondegenerate limit of
*S*_{2} and *S*_{3} should be satisfied, since reactions producing *S*_{2}or *S*_{3} may not increase the species number of *S*_{23}. In Table 4, we choose linear combinations of species whose collective species balance conditions
may not be satisfied by the species balance conditions. For other linear combinations
of species, their collective species balance conditions are derived from the ones
in Table 4. Satisfying all balance conditions in Table 4 guarantees satisfying balance conditions for all positive linear combination of species,
and these conditions help to identify scaling exponents which give a nondegenerate
limit of the normalized species numbers in the heat shock response model of *E. coli*. In most cases satisfying balance conditions gives nondegenerate limiting models
in the times of interest, but we can still find counter examples as given in the last
paragraph in the section for methods.

**Table 4.** **Balance equations and time-scale constraints for each species and for each collective
species chosen**

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.

Before we determine scaling exponents for *S*_{1}, *S*_{2}, and *S*_{3}, we run one realization of stochastic simulation to find ranges of the species numbers
in time. Using initial values for species *S*_{1}, *S*_{2}, and *S*_{3}, *X*_{1}(0)=10 and *X*_{2}(0)=*X*_{3}(0)=1 as given in Table 1 and using *N*_{0}=100, we set
*α*_{1}=1 and *α*_{2}=*α*_{3}=0 in the early stage of time period. Plugging in *α*_{i}’s and *β*_{k}’s in the balance equations for *S*_{2}, *S*_{3}, and *S*_{23}, equality holds in (12a) and (12b) but not in (16). Therefore, (18) gives

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
*α*_{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.

To identify a time scale involving a limiting model with interesting dynamics (nondegenerate),
we first need to determine a natural time scale of each species. Recall that a natural
time scale of species *S*_{i} is the time period of order *N*^{γ}_{i} when
*γ*_{i}for species *S*_{i} is rigorously determined by

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

Consider a model with the first set of scaling exponents including *α*_{1}=1 and *α*_{2}=*α*_{3}=0. Note that the first set of scaling exponents is valid when *γ*≤0 based on the time scale constraints given in Table 4. Substituting *α*_{2}=0 and *ρ*_{k}’s for the first scaling to the equation for

When *γ*=*γ*_{2}, the maximal scaling exponent in the propensities of all reaction terms in (20) should
be equal to the scaling exponent for the species number of *S*_{2}. Therefore, *γ*_{2}satisfies

and we get *γ*_{2}=0. Similarly, we get *γ*_{3}=*γ*_{8}=0.

Next, we plug *α*_{1}=1 and *ρ*_{k}’s for the first scaling in the equation for

By comparing the maximal scaling exponent in the propensities of all reaction terms
in (22) and the scaling exponent for the species number of *S*_{1}, *γ*_{1} satisfies

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,
*N*→*∞*. Similarly,
*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.

To derive the limiting equation for *S*_{2}, we set *γ*=0 in (20). Since the 2nd, 3rd, and 4th reaction terms have propensities with *N*^{0}=1 and the species number of *S*_{2} is of order 1, these reaction terms converge to nonzero limits in the limiting equation.
On the other hand, the propensities of the 5th, 6th, 7th, 8th and 9th reaction terms
are of order *N*^{−1} or *N*^{−2} which are smaller than the species number for *S*_{2}of order 1. Therefore, these reaction terms converge to zero as *N*→*∞*at least in the finite time interval. In the 2nd and 3rd reaction terms in (20),
*N*→*∞* since *γ*_{2}=*γ*_{3}=0. Then, using
*N*→*∞*, the limit of

Similarly, we get a limiting model with
*γ*=0 as given in (3).

Next, consider a model with the second set of scaling exponents including *α*_{1}=*α*_{2}=*α*_{3}=0. Note that the second set of scaling exponents is valid when *γ*≤1 based on the time scale constraints given in Table 4. To determine the natural time scale of *S*_{6}, substitute *α*_{6}=0 and *ρ*_{k}’s for the second scaling in the equation for

Comparing the exponents inside and outside of the reaction terms in (24), *γ*_{6} satisfies

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
*N*→*∞*, since *γ*_{i}>1. Thus, in the 12th and 15th reaction terms in (24),
*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),
*O*(1) and converge to
*N*→*∞* since *γ*_{6}=*γ*_{7}=*γ*_{8}=1.

Now, consider the asymptotic behavior of the 7th reaction term in (24) when *γ*=1. Since *γ*_{3}<1,
*N*→*∞*. However,
*S*_{3}. To get the limit of

The law of large numbers of Poisson processes gives an asymptotic limit of the scaled reaction terms as

where the *Y*_{k}’s are unit Poisson processes and *α*_{i}>0. For example, the 2nd reaction term in (26) divided by *N* is approximated as

Dividing (26) by *N* and using the law of large numbers for Poisson processes, we get

as *N*→*∞*.

We introduce an auxiliary variable to make the limiting model closed and define

Plugging *α*_{2}=*α*_{3}=0 and *ρ*_{k}’s in the second scaling in the equation for

Since
*γ*_{23} denotes a natural time scale exponent of *S*_{23}, we compare the scaling exponents of *N* in the reaction terms in (29) and the scaling exponent of *N* outside the reaction terms. Then *γ*_{23}satisfies

and we get *γ*_{23}=1. Since these reaction terms have *N*^{γ−2}=*N*^{−1} in their propensities when *γ*=1, which is smaller than the species number for *S*_{23} of order 1, these reaction terms converge to zero as *N*→*∞*. Using

Adding and subtracting terms in (28) and dividing the equation by −(*κ*_{2} + *κ*_{3}), we get

as *N*→*∞*, and this is used to obtain the limit of the 7th reaction term in (24). Letting *N*→*∞*, the limiting equation for

In (30), note that
*X*_{9}(0)=0 as given in Table 1. Limiting equations for
*γ*=1 is given in (4).

Last, consider a model with the third scaling exponents with *α*_{1}=0 and *α*_{2}=*α*_{3}=1. To derive a limiting equation for
*ρ*_{k}’s and *α*_{2}=*α*_{3}=1 for the third scaling in the equation for

In (31), the 8th reaction term is asymptotically zero, since the term is of order
*N*^{−1}. Using the law of large numbers for Poisson processes in (27), the 4th and the 9th
terms in (31) are asymptotically equal to

Since *γ*_{1}=2,
*N*→*∞*. On the other hand, since *γ*_{2},*γ*_{6}<2, both
*γ*_{3}<2. We actually show convergence of the fast-fluctuating species numbers of *S*_{2} and *S*_{3} to some limits in the Additional file 1: Section 5.1. For any *ε*>0 and for any *t* such that

uniformly as *N*→*∞*.

On the other hand, since *γ*_{6}<2,
*S*_{6}. Using the equation for

Dividing (35) by *N*^{2}, using the law of large numbers for Poisson processes in (27), and using the stochastic
boundedness of the propensities of the 8th, 9th, 15th, and 17th reaction terms in
the finite time interval shown in the Additional file 1: Section 5.1, we get

as *N*→*∞*. Therefore, a difference between the 10th and 12th reaction terms is approximated
in terms of

which converges to

as *N*→*∞*.

Now, from the equations for

Using the law of large numbers of Poisson processes in (27), the reaction terms in (39a) and (39b) are asymptotically equal to

Using (40a), (40b), and (38), the limiting equations of (39a) and (39b) are given as

In (41), note that
*X*_{9}(0)=0 as given in Table 1.

Since
*γ*_{8}=2,

Using (38) and (42),

From (33) and (43), we get

For more details used in (43) and (44), see Lemma 1.5 and Theorem 2.1 in [15]. Finally, we get the limiting equation of

#### Theorem 1

For *γ*=0,
*γ*=1,
*γ*=2,