Weinan E. Principles of Multiscale Modeling

Подождите немного. Документ загружается.

2.3. AVERAGING METHODS 47

equilibrium distribution of Y :

= F (X, t), F (x, t) =

f(x, y, t)ρ

∗

(y; x)dy (2.3.28)

For the example given above, the limiting equation is easily computed explicitly:

= −X

−

X + sin(πt) + cos(

√

2πt), X(0) = x. (2.3.29)

Here the new term −

X, is due to the noise in (2.3.25). The solution of the limiting

equation is also shown in Figure 2.3.

The limiting system (2.3.28) can be derived using perturbation analysis on the back-

ward Kolmogorov equation. Let h = h(x, y) be a smooth function and let

u(x, y, t) = Eh(X(t), Y (t)) (2.3.30)

where x, y are the initial values of X and Y : X(0) = x, Y (0) = y. u satisﬁes the backward

Kolmogorov equation [?]

∂

u = L

u +

u, (2.3.31)

with initial condition u(x, y, 0) = h(x, y), where

= f(x, y, t)∂

, L

= −(y − x)∂

∂

(2.3.32)

We look for solutions of (2.3.31) in the form:

u ∼ u

+ εu

+ ε

+ ··· (2.3.33)

Substituting this ansatz into (2.3.31) and collecting terms of equal orders in ε, we get

= 0, (2.3.34)

∂u

∂t

− L

, (2.3.35)

= ··· (2.3.36)

The ﬁrst equation implies that u

is a function of x and t only. The second equation

requires a solvability condition:

ρ(y)(∂

− L

)(y)dy = 0 (2.3.37)

48 CHAPTER 2. ANALYTICAL METHODS

where ρ is any solution of the adjoint equation

∗

ρ = ∂

((y − x)ρ) +

∂

ρ = 0 (2.3.38)

It is easy to see that solutions to this equation have the form

ρ(y; x) = Ce

−(y−x)

(2.3.39)

where C is a normalization constant. Therefore, u

(x, t) must satisfy

∂

−

= 0 (2.3.40)

where

L = F (x, t)∂

(2.3.41)

and F is given by (2.3.28). This is the backward equation for the dynamics described by

(2.3.28).

This analysis can be generalized to the case when the equation for the slow variable

also contains noise:

= f(X, Y, t) + σ(X, Y, t)

, X(0) = x (2.3.42)

= −

(Y − X) +

√

, Y (0) = y. (2.3.43)

Here

and

are independent white noises. In this case, the limiting dynamics is

given by:

= F (X, t) + Σ(X, t)

W , (2.3.44)

where F is given in (2.3.28) and Σ satisﬁes

Σ(x, t)

σ(x, y, t)

∗

(y; x)dy (2.3.45)

Results of this type have been obtained in many papers, including [?, ?, ?, ?].

It may happen that the eﬀective dynamics of the slow variable is the trivial equation

dx/dt = 0. In this case, we should consider the dynamics at longer time scale, say

t = O(1/ε). To do this we introduce the new variable τ = tε, and proceed as before

using the new variable τ.

2.3. AVERAGING METHODS 49

Let us look at some simple examples (see [?]). Consider

= Y +

, (2.3.46)

= −

Y +

√

, (2.3.47)

where

and

are independent white noises. Notice that in this case, the fast dynam-

ics is independent of the slow variable, and is actually a standard Ornstein-Uhlenbeck

process. Therefore, the invariant measure for the fast variable is N(0, 1/2), a normal

distribution with mean 0 and variance 1/2. Obviously, F = 0 in this case. Therefore the

eﬀective dynamics for the slow variable X is given by

X(t) = W

(t) (2.3.48)

One can indeed show that E|X(t) −

X(t)|

≤ Cε.

Next, consider

= Y

(2.3.49)

= −

Y +

√

. (2.3.50)

The fast dynamics is still the Ornstein-Uhlenbeck process. According to (2.3.44), the

eﬀective dynamics for the slow variable is given by:

X(t) =

√

(t). (2.3.51)

However, it is easy to see that

E|X(t) −

X(t)|

= E|

(Y (s) −

)dW

(s)|

= E

|Y (s) −

ds = O(1).

In this situation, we see that the fast dynamics is ergodic but the strong version of the

“averaging theorem” fails. This is not a contradiction to the perturbation analysis since

the latter only suggests that weak convergence holds, i.e. only the averaged quantities

converge.

50 CHAPTER 2. ANALYTICAL METHODS

In general, if we have a system of the type

= a(X, Y ) + σ

(X)

, (2.3.52)

b(X, Y ) +

√

(X, Y )

i.e. if σ

does not depend on Y , then one can show that the averaging theorem holds in

a strong sense, namely strong convergence can be established. See [?, ?, ?]. Otherwise,

one can only expect weak converegence.

2.3.3 Stochastic simulation algorithms

We start with a simple example [?, ?] (see also a related example in [?]). Consider

a system that contains multiple copies of a molecule in four diﬀerent conformations,

denoted by S

, S

. The conformation changes are describ ed by:

⇋

, S

⇋

(2.3.53)

with transition rates a

, ··· , a

which depend on the number of molecules in each con-

formation. The state of the system is described by a vector x = (x

, x

) where x

denotes the number of molecules in the conformation S

. To be speciﬁc, let us assume:

= 10

, ν

= (−1, +1, 0, 0),

= 10

, ν

= (+1, −1, 0, 0),

= x

, ν

= ( 0, −1, +1, 0),

= x

, ν

= ( 0, +1, −1, 0),

= 10

, ν

= ( 0, 0, −1, +1),

= 10

, ν

= ( 0, 0, +1, −1).

(2.3.54)

The vectors {ν

} are the stoichiometric vectors (also called the state change vectors):

After the j-th reaction, the state of the system changes from x to x + ν

. For example,

if the ﬁrst reaction happens, then one molecule in the S

conformation changes to the S

conformation. Hence x

is decreased by 1, x

is increased by 1, and x

and x

remain the

same. Each reaction is speciﬁed by two objects: the transition (reaction) rate function

and the stoichiometric vector:

= (a

(x), ν

) (2.3.55)

2.3. AVERAGING METHODS 51

If we start the system with roughly the same number of molecules in each conformation

and observe how the system evolves, we will see fast conformation changes between S

and S

, and b etween S

and S

, with o ccasional, much more rare conformation changes

between S

and S

In many cases, we are not interested in the details of the fast processes, but only the

dynamics of the slow processes. We should be able to simplify the problem and obtain an

eﬀective description for the slow processes only. To see this, let us partition the reactions

into two sets: the slow ones and the fast ones. The fast reactions R



, ν

)



are:

= 10

, ν

= (−1, +1, 0, 0),

= 10

, ν

= (+1, −1, 0, 0),

= 10

, ν

= ( 0, 0, −1, +1),

= 10

, ν

= ( 0, 0, +1, −1).

(2.3.56)

The slow reactions R

= {(a

, ν

)} are

= x

, ν

= (0, −1, +1, 0),

= x

, ν

= (0, +1, −1, 0).

(2.3.57)

Observe that during the fast reactions, the variables

= x

+ x

, z

= x

+ x

(2.3.58)

do not change. They are the slow variables. Intuitively, one can think of the whole

process as follows. In the fast time scale, z

and z

are almost ﬁxed and there is fast

exchange between x

and x

, and between x

and x

. This process reaches equilibrium

during the fast time scale and the equilibrium is given by a Bernoulli distribution:

, x

) =

! z

! x









. (2.3.59)

At the slow time scale, the values of z

and z

do change. The eﬀective dynamics on

this time scale reduces to a dynamics on the space of z = (z

, z

) with eﬀective reaction

channels given by:

¯a

= hx

i =

+ x

, ¯ν

= (−1, +1),

¯a

= hx

i =

+ x

, ¯ν

= (+1, −1).

(2.3.60)

52 CHAPTER 2. ANALYTICAL METHODS

These expressions follow from the general theory that we describe now.

For the general case, let us assume that we have a total of N species of molecules,

denoted by S

, ··· , S

. The number of molecules of species S

is denoted by x

. The

state vector is then given by x = (x

, ··· , x

). We will denote by X the state space where

x lives in. Assume that there are M reaction channels, each described by its reaction

rate and stoichiometric vector:

= (a

, ν

), R = {R

, . . . , R

}. (2.3.61)

Given the state x, the occurrence of the reactions on an inﬁnitesimal time interval dt

is indep en dent for diﬀerent reactions and the probability for the reaction R

to happ en

during this time interval is given by a

(x)dt. After reaction R

, the state of the system

changes to x + ν

. In the chemistry and biology literature, this is often called the

stochastic simulation algorithm (SSA) or Gilliespie algorithm, named after an algorithm

that realizes this process exactly (see [?]).

Let X(t) be the state variable at time t, and denote by E

the expectation conditional

on X(0) = x. Consider the observable u(x, t) = E

f(X(t)). u(x, t) satisﬁes the following

backward Kolmogorov equation [?]:

∂u(x, t)

∂t

(x) (u(x + ν

, t) − u(x, t)) = (Lu)(x, t). (2.3.62)

The operator L is the inﬁnitesimal generator of the Markov process associated with the

chemical kinetic system we are considering.

Now we turn to chemical kinetic systems with two disparate time scales. Assume that

the rate functions have the following form

a(x) =



(x),

(x)



, (2.3.63)

where ε ≪ 1 represents the ratio of the fast and slow time scales of th e system. The

corresponding reactions and the associated stoichiometric vectors can be grouped accord-

ingly:

= {(a

, ν

)}, R

= {(

, ν

)}. (2.3.64)

We call R

the slow reactions and R

the fast reactions.

We now show how to derive the eﬀective dynamics on the O(1) time scale, following

the perturbation theory developed in [?, ?]. For this purpose it is helpful to introduce

2.3. AVERAGING METHODS 53

an auxiliary pr ocess, called the virtual fast process [?]. This auxiliary process retains

the fast reactions only, all slow reactions are turned oﬀ, assuming, for simplicity that

the set of fast and slow reactions do not change over time. Intuitively, each realization

of the SSA consists of a sequence of realizations of th e virtual fast process, punctuated

by occasional ﬁring of the slow reactions. Due to the time scale separation, with high

probability, the virtual fast process has enough time to relax to equilibrium before another

slow reaction takes place. Therefore it is crucial to characterize the ergodic components

and the quasi-equilibrium distributions associated with these ergodic components for the

virtual fast process. In particular, it is of interest to consider parametrizations of these

ergodic components. A set of variables that parametrizes these ergodic components can

be regarded as the slow variables.

Let v b e a function of the state variable x. We say that v(x) is a slow variable if it

does not change during the fast reactions, i.e. if for any x and any stoichiometric vector

associated with a fast reaction

v(x + ν

) = v(x). (2.3.65)

This is equivalent to saying that the slow variables are conserved quantities for the fast

reactions. It is easy to see that v(x) = b · x is a slow variable if and only if

b · ν

= 0 (2.3.66)

for all ν

. The set of such vectors form a linear subspace in R

. Let {b

, b

, . . . , b

} be

a set of basis vectors of this subspace, and let

= b

· x, for j = 1, . . . , J. (2.3.67)

Let

z = (z

, ··· , z

) (2.3.68)

and denote by Z the space over which z is deﬁned. The stochiometric vectors associated

with these slow variables are naturally deﬁned as:

¯ν

= (b

· ν

, ··· , b

· ν

), i = 1, . . . , N

. (2.3.69)

We now derive the eﬀective dynamics on the slow time scale using singular pertur-

bation theory. We assume, for the moment that z is a complete set of slow variables,

54 CHAPTER 2. ANALYTICAL METHODS

i.e. , for each ﬁxed value of the slow variable z, th e virtual fast process admits a unique

equilibrium distribution µ

(x) on the state space X.

Deﬁne the projection operator P

(Pv)(z) =

x∈X

(x)v(x). (2.3.70)

By this deﬁnition, for any v : X → R, Pv depends only on the slow variable z, i.e.

Pv : Z → R.

The backward Kolmogorov equation for the multiscale chemical kinetic system reads:

∂u

∂t

= L

u +

u (2.3.71)

where L

and L

are the inﬁnitesimal generators associated with the slow and fast reac-

tions, respectively: For any f : X → R,

f)(x) =

j=1

(x)(f(x + ν

) − f (x)),

f)(x) =

j=1

(x)(f(x + ν

) − f (x)),

(2.3.72)

where M

is the number of slow reactions in R

and M

is the number of fast reactions

in R

. We look for a solution of (2.3.71) in the form of

u ∼ u

+ εu

+ ε

+ ··· (2.3.73)

Inserting this into (2.3.71) and equating the coeﬃcients of equal powers of ε, we arrive

at the hierarchy of equations:

= 0, (2.3.74)

∂u

∂t

− L

, (2.3.75)

= ··· (2.3.76)

The ﬁrst equation implies that u

belongs to the null-space of L

, which by the ergodicity

assumption of the fast process, is equivalent to saying that

(x, t) =

U(b · x, t) ≡

U(z, t), (2.3.77)

2.3. AVERAGING METHODS 55

for some

U yet to be determined. Inserting (2.3.77) into the second equation in (2.3.74)

gives (using the explicit expression for L

)

(x, t) =

∂

U(b · x, t)

∂t

−

i=1

(x)



U(b · (x + ν

), t) −

U(b · x, t)



∂

U(z, t)

∂t

−

i=1

(x)



U(z + ¯ν

, t) −

U(z, t)



(2.3.78)

This equation requires a solvability condition, namely that the right-hand side b e orthog-

onal to the null-space of the adjoint operator of L

. By the ergodicity assumption of the

fast process, this amounts to requiring that

∂

∂t

(z, t) =

i=1

¯a

(z)



U(z + ¯ν

, t) −

U(z, t)



(2.3.79)

where

¯a

(z) = (Pa

)(z) =

x∈X

(x)µ

(x). (2.3.80)

(2.3.79) is the eﬀective rates on the slow time scale. The eﬀective reaction kinetics on

this time scale are given in terms of the slow variable by

R = (¯a

(z), ¯ν

). (2.3.81)

It is important to note that the limiting dynamics described in (2.3.79) can be re-

formulated on the original state space X. Indeed, it is easy to check that (2.3.79) is

equivalent to

∂u

∂t

(x, t) =

i=1

˜a

(x) (u

(x + ν

, t) − u

(x, t)) , (2.3.82)

where ˜a

(x) = ¯a

(b · x), in the sense that if u

(x, t = 0) = f(b · x), then

(x, t) =

U(b · x, t) (2.3.83)

where

U(z, t) solves (2.3.79) with the initial condition

U(z, 0) = f(z). This fact is useful

for developing numerical algorithms for this kind of problems [?, ?, ?]. In fact, (2.3.82)

provides an eﬀective model that is more general than the setting considered above. Its

validity can be understood intuitively as follows [?]. To ﬁnd the eﬀective dynamics in

56 CHAPTER 2. ANALYTICAL METHODS

which the fast reactions are averaged out, we ask what happens for the original SSA

during a time scale that is very short compared with the time scale for the slow reactions

but long compared with the time scale of the fast reactions. During this time scale,

the original SSA essentially functions like the virtual fast process. Therefore to obtain

the eﬀective slow reaction rates, we need to average over the states of the virtual fast

process. However, during this time scale the virtual fast process has had suﬃcient time

to relax to its quasi-equilibrium. Hence the eﬀective slow reaction rates can be obtained

by averaging the original slow reaction rates over the quasi-equilibrium distribution of

the virtual fast p rocess. Note that all these can be formulated over the original state

space using the original variables. This is the main content of (2.3.82).

In the large volume limit when there are a large number of molecules in each species,

this stochastic, discrete model can be approximated by a deterministic, continuous model,

in the spirit of the law of large numbers [?]:

(x), j = 1, ··· , N (2.3.84)

The disparity in the rates results in the stiﬀness of this system of ODEs. Formally, we

can write (2.3.84) as

j=1

(x)ν

j=1

(x)ν

. (2.3.85)

The averagining theorem suggests that, in the limit as ε → 0, the solution of (2.3.85)

converges to the solution of

j=1

¯a

(b ·

x)ν

, (2.3.86)

Here

¯a

(z) =

(x)dµ

(x), (2.3.87)

and dµ

(x) is the equilibrium of the virtual fast process

j=1

)ν

. (2.3.88)

For instance, in the large volume limit, the ODEs corresponding to the simple example