Simulation of Chemical Reactions Gonzalo Mateos Dept. of Electrical - PowerPoint PPT Presentation

Simulation of Chemical Reactions Gonzalo Mateos Dept. of Electrical and Computer Engineering University of Rochester gmateosb@ece.rochester.edu http://www.ece.rochester.edu/~gmateosb/ November 16, 2014 Introduction to Random Processes Simulation of Chemical Reactions 1

Gillespie’s algorithm Gillespie’s algorithm Dimerization Kinetics Enzymatic Reactions Lactose digestion (lac operon) Introduction to Random Processes Simulation of Chemical Reactions 2

Simulation of chemical reactions ◮ Chemical system with m reactant types and n possible reactions ◮ Reactant quantities change over time as reactions occur ◮ Nr. of type j reactants at time t denoted as X j ( t ) ◮ System’s state ⇒ vector X ( t ) := [ X 1 ( t ) , X 2 ( t ) , . . . , X j ( t )] T ◮ To specify i -th reaction ⇒ reactants, products and rates h i ( X ) R i : s l i 1 X 1 + s l i 2 X 2 + . . . + s l → s r i 1 X 1 + s r i 2 X 2 + . . . + s r im X m im X m ◮ ( s l i 1 molecules of type 1) + . . . + ( s l im molecules of type m ) react ... ... to yield ( s r i 1 of type 1) + . . . + ( s r im of type m ) ◮ Rate of reaction h i ( X ) depends on number of molecules present ◮ Let T i ( X ) denote the time until the i -th reaction when state is X Introduction to Random Processes Simulation of Chemical Reactions 3

Stoichiometry matrices ◮ Can be more conveniently written using matrices ⇒ Define vector of rates h ( X ) = [ h 1 ( X ) , h 2 ( X ) , . . . , h n ( X )] T ⇒ Define stoichiometry left matrix S ( l ) with elements s l ij ⇒ Define stoichiometry right matrix S ( r ) with elements s r ij h ( X ) ◮ Write system of chemical reactions as ⇒ S ( l ) X → S ( r ) X     X 1 X 1     X 1 s l X 2 X 1 s r X 2     = X = X i 1 i 1 X 2 s r     X 2 s l     i 2 · i 2 ·     X m X m X m s l X m s r im im         s l 11 s l 12 · s l s l 11 X 1 + . . . + s l s r 11 s r 12 · s r s r 11 X 1 + . . . + s r 1 m X m 1 m X m 1 m 1 m             · · · · ·     · · · · ·         = s r s r i 2 · s r s r i 1 X 1 + . . . + s r  s l s l i 2 · s l   s l i 1 X 1 + . . . + s l  im X m im X m         i 1 im i 1 im                  · · · ·   ·  · · · · ·         s r n 2 s r n 2 · s r s r n 1 X 1 + . . . + s r     nm X m s l n 2 s l n 2 · s l s l n 1 X 1 + . . . + s l nm X m nm nm S ( r ) S ( r ) X S ( l ) S ( l ) X Introduction to Random Processes Simulation of Chemical Reactions 4

Example 1: Dimerization kinetics ◮ Molecule can exist in simple form P and as a dimer D ◮ Define vector X := [ P , D ] T ◮ Possible reactions are dimerization and dissociation h 1 ( X ) R 1 (Dimerization): 2 P → D h 2 ( X ) R 2 (Dissociation): D → 2 P ◮ Rates and stoichiometry matrices S ( l ) and S ( r ) given by � 2 � 0 � h 1 ( X ) � � � 0 1 S ( l ) = S ( r ) = , , h ( X ) = 0 1 2 0 h 2 ( X ) h ( X ) ◮ Rewrite equations more compactly as ⇒ S ( l ) X → S ( r ) X Introduction to Random Processes Simulation of Chemical Reactions 5

Example 2: Enzymatic reaction ◮ Substrate S converted to product P . Enzyme E catalyzes conversion ◮ Converting S into P directly requires significant energy ◮ Enzyme E reacts with S to form intermediate molecule SE (binding) ◮ Molecule SE then separates into product P liberating E (conversion) ◮ This cycle requires less energy than direct conversion ◮ SE may also separate back into S and E (dissociation) ◮ Possible reactions are binding, conversion and dissociation, then h 1 ( X ) R 1 (Binding): S + E → SE h 2 ( X ) R 2 (Dissociation): SE → S + E h 3 ( X ) → E + P R 3 (Conversion): SE Introduction to Random Processes Simulation of Chemical Reactions 6

Example 2: Enzymatic reaction (continued) ◮ System state represented by vector X := [ S , E , SE , P ] T ◮ Stoichiometry matrices S ( l ) and S ( r ) given by S E SE P S E SE P     1 1 0 0 R 1 0 0 1 0 R 1 S ( l ) = S ( r ) = 0 0 1 0 R 2 1 1 0 0 R 2     0 0 1 0 R 3 0 1 0 1 R 3 ◮ Reaction rate vector h ( X ) = [ h 1 ( X ) , h 2 ( X ) , h 3 ( X )] T h ( X ) ◮ Rewrite equations more compactly as ⇒ S ( l ) X → S ( r ) X Introduction to Random Processes Simulation of Chemical Reactions 7

Second order reaction ◮ Consider second order reaction R i : X 1 + X 2 → . . . (two reactants) ◮ Let T i ( X 1 , X 2 ) be time until R occurs when there are X 1 type 1 and X 2 type 2 molecules ◮ Have seen that T i ( X 1 , X 2 ) is exponentially distributed with rate h i ( X ) = h i ( X 1 , X 2 ) = c i X 1 X 2 ◮ Constant c i measures reactivity of X 1 and X 2 ◮ Argument ⇒ T i (1 , 1) memoryless (depends on chance encounter) ⇒ Thus T i (1 , 1) is exponential with, say, parameter c i ⇒ T i ( X 1 , X 2 ) is the minimum of X 1 X 2 exponentials ⇒ T i ( X 1 , X 2 ) exponential with parameter c i X 1 X 2 Introduction to Random Processes Simulation of Chemical Reactions 8

Second order involving molecules of same type ◮ Second order reaction with two molecules of same type R i : X 1 + X 1 → . . . ◮ Hazard depends on the number of molecules X 1 , i.e. h i ( X ) = h i ( X 1 ) ◮ Reaction does not occur if there is a single molecule ◮ If there are 2 molecules T i (2) is exponential with parameter, say, c i ◮ For arbitrary X 1 there are X 1 ( X 1 − 1) / 2 possible encounters ◮ Then, T i ( X 1 ) is exponential with parameter h i ( X ) = h i ( X 1 ) = c i X 1 ( X 1 − 1) / 2 ◮ c i X 1 ( X 1 − 1) / 2 substantially different from c i X 2 1 / 2 for small X 1 Introduction to Random Processes Simulation of Chemical Reactions 9

Zero-th, first and higher order reactions ◮ Zero-th order reaction R i : ∅ → X 1 (spontaneous generation) ◮ Assume an exponential model with constant rate h i = c i ◮ Used to model exogenous factors (and biblical phenomena) ◮ First order reaction R i : X 1 → . . . (decay) ◮ Exponential with rate h i ( X ) = h i ( X 1 ) = c i X 1 ◮ Higher order reactions involving more than two reactants ◮ E.g., third order reaction R i : X 1 + X 2 + X 3 → X 4 ◮ Time until next R i reaction exponential. Hazard: h i ( X ) = c i X 1 X 2 X 3 ◮ Reactions of order more than 2 are rare ◮ Most likely, R i is encapsulating two second order reactions X 1 + X 2 → X 5 , X 5 + X 3 → X 4 Introduction to Random Processes Simulation of Chemical Reactions 10

The hazard function ◮ All reaction times are exponential RVs ⇒ CTMC with state X ◮ Hazards h i ( X ) determine transition rates of CTMC ◮ Hazards for zero-th, first and second order reactions (for reference) Order Reaction Rate c zero-th ∅ → · · · c c first X 1 → · · · cX 1 c second X 1 + X 2 → · · · cX 1 X 2 c second 2 X 1 → · · · cX 1 ( X 1 − 1) / 2 ◮ Probability of reaction R i happening in infinitesimal time ǫ is P ( T i ( X ) < ǫ ) = h i ( X ) ǫ + o ( ǫ ) ◮ That’s why the name hazard Introduction to Random Processes Simulation of Chemical Reactions 11

State transition for given reaction ◮ State is X ( t ) = X . Reaction R i occurs. Next state X ( t + dt ) = Y ? ◮ Number of reactants per type = = i -th row of left stoichiometry matrix s ( l ) = [ s l i 1 , s l i 2 , . . . , s l im ] T i h i ( X ) s l i 1 X 1 + s l i 2 X 2 + . . . + s l im X m → . . . ◮ Number of products per type = = i -th row of right stoichiometry matrix s ( r ) = [ s r i 1 , s r i 2 , . . . , s r im ] T i h i ( X ) → s r i 1 X 1 + s r i 2 X 2 + . . . + s r . . . im X m ◮ X decreases by nr. of reactants and increases by nr. of products ◮ Next sate is ⇒ Y = X − s ( l ) + s ( r ) (upon reaction R i ) i i Introduction to Random Processes Simulation of Chemical Reactions 12

Transition rates and probabilities ◮ q ( X , Y ) = transition rate from state X to state Y . Given by � � X , X − s ( l ) + s ( r ) q = h i ( X ) , i = 1 , . . . , n i i ◮ Transition from state X to X − s ( l ) + s ( r ) when reaction R i occurs i i ◮ ν ( X ) = Transition rate out of X into any state (any reaction occurs) n n � � � X , X − s ( l ) + s ( r ) � ν ( X ) = = h i ( X ) q i i i =1 i =1 ◮ P ( X , Y ) = Prob. of going into Y given transition out of X occurs � � X , X − s ( l ) + s ( r ) q = h i ( X ) i i � � X , X − s ( l ) + s ( r ) = P i i ν ( X ) ν ( X ) ◮ Probability that i -th reaction occurs given that a reaction occurred Introduction to Random Processes Simulation of Chemical Reactions 13

Gillespie’s algorithm Gillespie’s algorithm = Simulation of CTMC Input: Stoichiometry matrices S ( l ) and S ( r ) . Initial state X (0) Output: Molecule numbers as a function of time X ( t ) (1) Initialize time and CTMC’s state t = 0 , X = X (0) (2) Calculate all hazards ⇒ h i ( X ) (3) Calculate transition rate ⇒ ν ( X ) = � n i =1 h i ( X ) � � (4) Draw random time of next reaction ∆ t ∼ Exp ν ( X ) (5) Advance time to t = t + ∆ t (6) Draw reaction at time t + ∆ t ⇒ R i drawn with prob. h i ( X ) /ν ( X ) (7) Update state vector to account for this reaction ⇒ X − s ( l ) + s ( r ) i i (8) Repeat from (2) Introduction to Random Processes Simulation of Chemical Reactions 14