STOCHASTIC CHEMICAL KINETICS

14 Pages·2007·0.52 MB·English

STOCHASTIC CHEMICAL KINETICS
Dan Gillespie
Dan T Gillespie Consulting
[email protected]

Current Support: University of California at Santa Barbara (DOE)
Caltech (Beckman Institute / BNMC)
Caltech (NIGMS)
Past Support: Molecular Sciences Institute (Sandia / DOE)
Caltech (DARPA / AFOSR)
ONR

Rethinking Chemical Kinetics
A Chemically Reacting System consists of …
• Molecules of N chemical species S ,…,S .
1 N
- Inside a volume Ω, at some temperature T.
• M “elemental” reaction channels R,…,R .
1 M
- R describes a single instantaneous physical event, which changes the
j
population of at least one species.
- “Elemental” means that R is one of two types:
j
S → something else (unimolecular),
i
or
S +S → something else (bimolecular).
i i′
- All other types (trimolecular, reversible, etc.) are made up of a series
of two or more elemental reactions.
1
How does a chemically reacting system evolve in time?
The traditional answer, for spatially homogeneous systems:
“According to the reaction rate equation (RRE).”
• A set of coupled, first-order ODEs.
• Derived using ad hoc, phenomenological reasoning.
o Is more than the “mass action equations” of thermodynamics,
which apply only to systems in equilibrium.
• Implies the system evolves continuously and deterministically, even
though molecules come in integer numbers and react stochastically.
• Is empirically accurate for large (test tube size) systems.
• But is often not adequate for very small (cell-size) systems.
* * *
Let’s take a fresh look at this question.
Doing it “right”: Molecular Dynamics
• The most exact way of describing the system’s evolution.
• The “motion picture” approach: Tracks the position and velocity of
every molecule in the system.
• Simulates every collision, non-reactive as well as reactive.
• Shows changes in species populations and their spatial distributions.
• But . . . it’s unfeasibly slow for nearly all realistic systems.
2
A great simplification occurs if successive reactive collisions tend
to be separated in time by very many non-reactive collisions.
• The overall effect of the non-reactive collisions is to randomize
- the velocities of the molecules (Maxwell-Boltzmann distribution).
- the positions of the molecules (spatially uniform or well-stirred),
• Then, instead of having to describe the system’s state as the
position, velocity and species of each molecule, we need only give
X(t)(cid:2)(X (t),…,X (t)),
1 N
X (t)(cid:2) the number of S molecules at time t.
i i
But this well-stirred simplification, which . . .
• ignores the non-reactive collisions,
• truncates the definition of the system’s state,
. . . comes at a price:
X(t) must now be viewed as a stochastic process.
(cid:1) But in fact, the system was never deterministic to begin with.
Even if molecules moved according to classical mechanics . . .
- Unimolecular reactions always involve randomness (QM).
- Bimolecular reactions usually do too.
- A system of many colliding molecules is so sensitive to initial
conditions that, for all practical purposes, it evolves “randomly”.
- The system is not isolated. It’s in a heat bath, which keeps it “at
temperature T” – via essentially random interactions.
3
For well-stirred systems, each R is completely characterized by …
j
• a propensity function a (x): Given the system is in state x,
j
a (x)dt (cid:2) probability that one R event will occur in the next dt.
j j
- The existence and form of a (x) follow from molecular physics.
j
• a state change vector ν ≡(ν ,…,v ):
j 1j Nj
ν (cid:2) the change in X caused by one R event.
ij i j
- R induces x→x+ν . {ν }≡ the “stoichiometric matrix.”
j j ij
Examples:
a (x)=c x x , ν =(+1,−1,0,…,0)
S +S (cid:5)(cid:3)(cid:3)(cid:3)c1(cid:3)(cid:3)(cid:4)(cid:3)2S :  1 1 x1 (2x −1) 1
1 2 c2 1 a2(x)=c2 1 21 , ν2 =(−1,+1,0,…,0)
Prob{v -collision in dt}=(πr122)(v12dt). Prob{R |v -collision}(cid:2) p (v ).
12 Ω j 12 j 12
(πr2)(v dt) πr2 
 12 12 ×p (v ) ×xx =  12 v p (v ) xx dt=c xx dt
 Ω  j 12 1 2  Ω 12 j 12 v12 1 2 (cid:6)j 1 2
(cid:7)(cid:8)(cid:8)(cid:8)(cid:8)(cid:8)(cid:9)(cid:8)(cid:8)(cid:8)(cid:8)(cid:8)(cid:10)v12 (cid:7)(cid:8)(cid:8)(cid:8)(cid:9)c (cid:8)(cid:8)(cid:8)(cid:10) aj(x)
Prob that a randomly chosen j
S1-S2 pair does an Rj in next dt
8k T  E 
Rj iff "collisional K.E." > Ej ⇒ v12pj(v12) v12 =(cid:7)(cid:8)π(cid:9)mB(cid:8)(cid:10)12 e(cid:7)xp(cid:8)(cid:8)(cid:9)−k(cid:8)BT(cid:8)j(cid:10)
v Arrhenius
12
4
Two exact, rigorously derivable consequences . . .

(cid:1) 1. The chemical master equation (CME):
∂P(x,t|x ,t ) M
0 0 =∑a (x−ν )P(x−ν ,t|x ,t )−a (x)P(x,t|x ,t ).
∂t  j j j 0 0 j 0 0 
j=1
• P(x,t|x ,t )(cid:2)Prob{X(t)=x, given X(t )=x } for t≥t .
0 0 0 0 0
• Follows from the probability statement
 M 
P(x,t+dt|x ,t )=P(x,t|x ,t )×1−∑(a (x)dt)
0 0 0 0 j
 j=1 

M
+∑P(x−ν ,t|x ,t )×(a (x−ν )dt).
j 0 0 j j
j=1
• But the CME is usually too hard to solve.
• Averages: f (X(t)) (cid:2)∑ f(x)P(x,t|x ,t ).
0 0
x
If we multiply the CME through by x and then sum over x, we find
d X(t) M
=∑ν a (X(t)) .
dt j j
j=1
• If there were no fluctuations,
a (X(t)) =a ( X(t) )=a (X(t)),
j j j
and the above would reduce to:
dX(t) M
=∑ν a (X(t)).
dt j j
j=1
- This is the reaction-rate equation (RRE).
- It’s usually written in terms of the concentration Z(t)(cid:2)X(t) Ω.

(cid:2) But as yet, we have no justification for ignoring fluctuations.
5
(cid:1) 2. The stochastic simulation algorithm (SSA):
A procedure for constructing sample paths or realizations of X(t).
Idea: Generate properly distributed random numbers for
- the time τ to the next reaction,
- the index j of that reaction.
• p(τ, j|x,t)dτ(cid:2) probability, given X(t)=x, that the next reaction
will occur in [t+τ,t+τ+dτ), and will be R .
j
=P(τ)×a (x)dτ, P(τ)(cid:2)Pr(no reactions in time τ).
0 j 0
P(τ+dτ)=P(τ)×(1−a (x)dτ), where a (x)(cid:2)∑Ma (x).
0 0 0 0 1 j′
dP(τ)
Implies 0 =−a (x)P(τ). Solution: P(τ)=e−a0(x)τ.
dτ 0 0 0
a (x)
∴ p(τ, j|x,t)=e−a0(x)τa (x)=a (x)e−a0(x)τ× j .
j 0 a (x)
0
Thus,
- τ is an exponential random variable with mean 1 a (x),
0
- j is an integer random variable with probabilities a (x) a (x).
j 0
The “Direct” Version of the SSA
M
1. In state x at time t, evaluate a (x),…,a (x), and a (x)≡∑a (x).
1 M 0 j′
j′=1
2. Draw two unit-interval uniform random numbers r and r , and
1 2
compute τ and j according to
1 1
• τ= ln ,
a0(x) r1
j
• j = the smallest integer satisfying ∑a (x)>r a (x).
k 2 0
k=1
3. Replace t←t+τ and x←x+ν .
j
4. Record (x,t). Return to Step 1, or else end the simulation.
6
A Simple Example: S c1→0.
1
a (x )=c x , ν =−1. Take X (0)=x0.
1 1 1 1 1 1 1
dX (t)
RRE: 1 =−c X (t). Solution is X (t)=x0e−c1t.
dt 1 1 1 1
∂P(x ,t|x0,0)
CME: 1 1 =c (x +1)P(x +1,t|x0,0)−xP(x ,t|x0,0).
∂t 1 1 1 1 1 1 1 
Solution: P(x ,t|x0,0)= x10! e−c1x1t(1−e−c1t)x10−x1 (x =0,1,…,x0)
1 1 x !(x0−x )! 1 1
1 1 1
which implies X (t) =x0e−c1t, sdev{X (t)}= x0e−c1t(1−e−c1t).
1 1 1 1
1 1
SSA: Given X1(t)=x1, generate τ= c x lnr, then update:
1 1
t←t+τ, x ←x −1.
1 1
100
S → 0
1
80 c = 1, X(0) = 100
1 1
60
40
20
0
0 1 2 3 4 5 6
7
100
S → 0
1
80 c = 1, X(0) = 100
1 1
60
40
20
0
0 1 2 3 4 5 6
The SSA . . .
• Is exact.
• Does not entail approximating “dt” by “∆t”.
• Is logically on par with the CME (but is not a method for numerically
solving the CME).
• Is procedurally simple, even when the CME is intractable.
• Comes in a variety of implementations …
- Direct Method (Gillespie, 1976)
- First Reaction Method (Gillespie, 1976)
- Next Reaction Method (Gibson & Bruck, 2000)
- First Family Method (Lok, 2003)
- Modified Direct Method (Cao, Li & Petzold, 2004)
- Sorting Direct Method (McCollum, et al. 2006)
• Remains too slow for most practical problems: Simulating every
reaction event one at a time just takes too much time if any reactants
are present in very large numbers.
8
We would be willing to sacrifice a little exactness . . .
. . . if that would buy us a faster simulation.
Tau-Leaping
• Approximately advances the process by a pre-selected time τ,
which may encompass more than one reaction event.
• Key: The definition of “the Poisson random variable with mean aτ”:
P(aτ)(cid:2) the number of events that will occur in a time τ,
given that the probability of an event in any dt is adt
where a can be any positive constant.
• With X(t)=x, let us choose τ small enough to satisfy the
Leap Condition: Each a (x)≈ constant in [t,t+τ].
j
• Then: The number of R firings in [t,t+τ] ≈ P(a (x)τ).
j j
M
X(t+τ)(cid:11)x+∑P (a (x)τ)ν
j j j
j=1
- Practical Implementation of Tau-Leaping -
• We have two control parameters, ε and n :
c
- To satisfy the Leap Condition, restrict τ so that ∆a a ≤ε,∀j.
τ j j
- To avoid populations <0, allow only one firing of all critical reactions
((cid:2) reactions that are within n firings of exhausting any reactant).
c
• We take τ=min(τ′,τ′′), where:
- τ′ maximally satisfies the Leap Condition for firings of the non-
critical reactions. (We have a fairly efficient way to estimate τ′.)
- τ′′ is the time to the next critical reaction . ( Generate τ′′ by
applying the SSA to the critical reactions.)
• For each non-critical R , generate k as a sample of P(a (x)τ).
j j j
• If τ′<τ′′: Set the k ’s for all the critical R ’s to 0.
j j
If τ′′≤τ′: Use the SSA to determine which critical reaction fires, set
its k to 1, and set all other critical k ’s to 0.
j j
M
• Leap: t←t+τ and x←x+∑k ν .
j j
j=1
(cid:1) Becomes the SSA if all reactions are critical (n →∞).
c
9
50000
45000
S1--> 0 c1 = 1
40000 S1 + S1 --> S2 c2 = 0.002
mer X3 SS22 ---->> SS13 + S 1 cc43 == 00..054
X1, unstable dimer X2, stable di22330505000000000000 X2 ---- EI55n01xt07iaa ,clr0leyt6 :aS7 c SXtrieA1oa=n Rc1stu 0ipo0nen,.r0s 0pto0lot;ta tXel.2d= dXo3t=0.
mer 15000
mono X3
10000
5000
X1
0
0 2 4 6 8 10 12 14 16 18 20 22 24
Time
50000
45000
S1--> 0 c1 = 1
40000 S1 + S1 --> S2 c2 = 0.002
mer X3 SS22 ---->> SS13 + S 1 cc43 == 00..054
X1, unstable dimer X2, stable di22330505000000000000 X2 ----- εεεε19RE 0 xul=e5pn a 0l litpe.ci0m aip4tpe e;sT rs a ntppoulecto aeL t=ltde.e ua1dpp0 d .ionovgt.e Rr SuSnA > 10X.
mer 15000
mono X3
10000
5000
X1
0
0 2 4 6 8 10 12 14 16 18 20 22 24
Time
10