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