Introduction to Stochastic Simulation with the Gillespie ...

1 Introduction to Stochastic Simulation with the Gillespie MethodDavid KarigApril 18, 2005 Stochastic Systems Many systems driven by random, discrete interactions Traditional deterministic models may not accurately describe such systems$Example: The Lambda SwitchVirus Decision Dictated by NoiseArkin, Ross, McAdams, Genetics(1998)Example: The Lambda SwitchVirus Decision Dictated by NoiseArkin, Ross, McAdams, Genetics(1998)Example: The Lambda SwitchVirus Decision Dictated by NoiseArkin, Ross, McAdams, Genetics(1998)Outline Deterministic rate reaction model Gillespie method Examples Lambda phage Epidemiology OptimizationsDeterministic ModelReactions:A + B C C A + BA BC ODE s:dA/dt = -k1*A*B + k2*C k3*AdB/dt = -k1*A*B + k2*C + k3*AdC/dt = k1*A*B k2*C - k4*C Given initial conditions, integrate the coupled equations for some period of timek1k4k2k3 Problem StatementIf we start with N species which can interact through one of M reactions at a given time,what will be the population levels of species after a given period of time?

2 P** p**p**Deterministic Solution Continuous Average kinetic rates represent reaction probabilitiesdA/dt = -k1*A*B + k2*C k3*AdB/dt = -k1*A*B + k2*C + k3*AdC/dt = k1*A*B k2*C k4*CA + B CC A + B A B C 05100500100015002000tABCV alidity of Deterministic SolutionA + B C024681002004006008001000tABCV alidity of Deterministic Solution02004006008000246810tABC02004006 008000246810tABCR eaction Probabilities c dt= average probability that a particularcombination of reactants will react according to R in the next time interval dt h = number of reactant combinations h c dt= a dt= average probability that an R reaction will occur somewhere inside Vin the next time interval dtReaction Probabilities c dt= average probability that a particularcombination of reactants will react according to R in the next time interval dt h = number of reactant combinations h c dt= a dt= average probability that an R reaction will occur somewhere inside Vin the next time interval dtX molecules of AY molecules of BXYc1dt= probability that an R1reaction will occur somewhereinside Vin the next time interval dtR1:A + B Ch1= XY reactant combinationsExact Stochastic Simulation Avoid averaging assumptions Probabilistic formulation When does next reaction occur?

3 Whichreaction occurs next? Reaction probability density function:tt+ P( , )d = probability at time tthat the next reaction is R and occurs in interval (t+ ,t+ +d ) t+ +d R A + B CC A + B A B C R1:R2:R3:R4:System State:t+ t78C21B21AP( , )d = P0( ) h c d probability that no reaction occurs during (t,t+ )probability that reaction occurs during (t+ , t+ +d )Deriving Reaction PDFtt+ t+ +d R Deriving P0( ) for Reaction PDFK =tt+ Divide (t,t+ ) into K subintervals of width Probability that none of the reactions occur in any of the K subintervals:()[] = = == = = MiiichKMiiiKKMiiieKchchP11101lim1 Direct Method for Generating P( , )() ()())()(,000001 PPaaeaeaechPchPaachMiii = == == =when next reaction occurswhich ==Ma1 RMR3R2R1 Direct Method for Generating P( , )() ()())()(,000001 PPaaeaeaechPchPaachMiii = == == =when next reaction occurswhich reactionoccurs = =<< 1211iioiiaaranda is the integer for ==Miia1rand2a0 = (1/a0)ln(1/rand1)RMR3R2R1 Stochastic Simulation Algorithm1.

4 Initialization Set values of c for the M reactions. Set initial population sizes2. Calculate the M values a and a0 = a .3. Generate ( , ) based on P( , )4. Adjust population levels according to the reaction R , and increase t by 5. Return to Step 2 Lambda Phage Developmental PathwayArkin, Ross, McAdams, Genetics(1998) Regulatory circuit exploits Stochastic noise to produce different outcomes Stochastic model can predict statistics of regulatory outcomesLambda Phage: Regulatory Circuit CI and Crocompetitively bind OR1, OR2, OR3 Cro represses PRand PRM CI represses PR, can activate PRML ambda Phage: Regulatory Circuit CI and Crocompetitively bind OR1, OR2, OR3 Cro represses PRand PRM CI represses PR, can activate PRMTime Evolution of Two RunsEpidemiology ExampleSusceptibleInfectedRecoveredp** p**p**050100150200700800900# susceptibles050100150200100200# infected050100150200420044004600# recoveredt (weeks)050100150200700800900# susceptibles050100150200100200# infected050100150200420044004600# recoveredt (weeks)

5 Evaluation of Stochastic Simulation Advantages continuous time, discrete population changes captures effects of noise simple implementation small memory requirements Disadvantages CPU intensive typically must simulate many runs must use good random number generator periodicity affects size of Simulation resolution limits range of probabilitiesComputational Requirements Memory (N + 2M + 1) N species populations cand avalues for each of M reactions; a0 Total time scales with number of reactions that occur Operations per reaction:Calculate avaluesSearch Calculate Calculate a0 Generate 2 random numbersO(M)O(M)O(M)Optimized Direct Method (ODM)Y. Cao, H. Li, L. Petzold. 2004. J. Chem.

6 Phys. 121:4059-4067 = =<< 1211iioiiaaranda Reduce cost of searching for index Observation: Reactions are typically multiscale in a large system. Subset will frequently occur. Sort index of reactions based on how often they ==Miia1rand2a0 Optimized Direct Method (ODM) Reduce cost of calculating all a Dependency graph Reduce cost of summing all a to calculate a0 Modify a0by subtracting old values, adding newGibson and Bruck, 2000 Conclusion Provides means of studying role of noise in complex systems Can predict statistics (Lambda phage) Can depict behavior that deterministic simulations do not capture (epidemiology example) Enhancing performance an active area of investigationUseful References / Further Study Implementations STOCKS: BioNetS: Algorithms/optimizations Gillespie .

7 A General Method for Numerically Simulating the Stochastic Time Evolution of Coupled Chemical Reactions. 1976. J Comput Phys 22:403-434. Gillespie . Exact Stochastic Simulation of Coupled Chemical Reactions. 1977. J Phys Chem81:2340-2361 Gibson and J. Bruck. 2000. Efficient Exact Stochastic Simulation of Chemical Systems with Many Species and Many Channels. J Phys Chem104:1876-1889 Gillespie . 2001. Approximate accelerated Stochastic Simulation of chemically reacting systems. J. Chem. Phys. 115:1716-1733. Y. Cao, H. Li and L. Petzold, Efficient formulation of the Stochastic Simulation algorithm for chemically reacting system. 2004. J Chem Phys121:4059-4067 Y. Cao, D. Gillespie , L. Petzold, The slow-scale Stochastic Simulation algorithm.

8 2005. J Chem Phys122(1). A. Chatterjee and Vlachos. Binomial distribution based -leap accelerated Stochastic Simulation . 2005. J Chem Phys122:024112 Example Uses A. Arkin, J. Ross, H. McAdams. Stochastic Kinetic Analysis of Developmental Pathway Bifurcation in Phage-Infected Escherichia coliCells. 1998. Genetics 149:1633-1648 J. Dushoff, Plotkin, Levin, Earn. Dynamical resonance can account for seasonality of influenza epidemics. 2004 PNAS. S. Hooshangi, S. Thiberge, R. Weiss. Ultrasensitivity and noise propagation in a synthetic transcriptional cascade. 2005. PNAS 102:3581-3586 Acknowledgements Ron Weiss and the Weiss group PICASso Steven Kleinstein SinghComparing to Experimental Results Compare predicted percent lysogenization to experimental results at different infection levels Average phage input (API) = ratio of phage particles to cells at time of infection Multiplicity of infection (MOI) = phage particles per cell Poisson probability that a given cell will be infected with MOI=M when API=A Expected fraction of lysogens:()AMeAAAMP =!

9 ,()()() =MlysogensMFAMPAF,Comparing to Experimental ResultsExample (deterministic assumptions)Average rate of R1 in dtper unit volume:VcXXVcXX121121=Using concentrations xi= Xi/V and dividing by density of reactants:211211xxVcxxk=Deterministic assumption:Vck11=12112S SS :R +21122111xxktxxxktx = = Master equation B dt= probability that single R reaction brings us to state X1,..,XN a dt= h c dt, where h is the number of reactant combinations Often difficult to solve analytically and even numerically Instead simulate individual trajectories using Monte Carlo algorithm()()[] = = MNNtXXPaBtXXPt111.