An Alternative to Gillespie’s Algorithm for Simulating Chemical Reactions Barbuti An Alternative to Gillespie’s Algorithm for Maggiolo-Schettini Milazzo Troina Simulating Chemical Reactions Outline Introduction Roberto Barbuti Andrea Maggiolo–Schettini Simulating chemical reactions Paolo Milazzo Angelo Troina Gillespie’s algorithm Our algorithm A comparison Examples and Dipartimento di Informatica, Universit` a di Pisa, Italy applications Probabilistic MultiSet Rewriting Edinburgh – April 4, 2005 Simulation results Conclusions 1/19
An Alternative to Gillespie’s Algorithm for Simulating Introduction Chemical Reactions Barbuti Maggiolo-Schettini Simulating chemical reactions Milazzo Troina Gillespie’s algorithm Outline Our algorithm Introduction A comparison Simulating chemical reactions Gillespie’s algorithm Examples and applications Our algorithm A comparison Probabilistic MultiSet Rewriting Examples and applications Simulation results Probabilistic MultiSet Rewriting Simulation results Conclusions Conclusions 2/19
Introduction An Alternative to Gillespie’s Algorithm for Simulating Chemical reactions are described by the law of mass action Chemical Reactions ◮ the speed of a reaction is proportional to the Barbuti Maggiolo-Schettini concentrations of the individual reactants involved Milazzo Troina ◮ differential equations Outline Introduction Gillespie’s simulation algorithm Simulating ◮ stochastic method based on the theory of collisions chemical reactions Gillespie’s algorithm ◮ each reaction takes a (continuous) random time which Our algorithm A comparison is exponentially distributed Examples and applications Probabilistic MultiSet Rewriting The simulation algorithm we propose Simulation results ◮ performs discrete time steps of (fixed) lenght ∆ t ; Conclusions ◮ assume that at each step at most one reaction may occur (randomly chosen). 3/19
Background An Alternative to Gillespie’s Algorithm for Usual notation for chemical reactions: Simulating Chemical Reactions k k − 1 ℓ ′ 1 P 1 + . . . + ℓ ′ ℓ 1 S 1 + . . . + ℓ ρ S ρ γ P γ Barbuti ⇋ Maggiolo-Schettini Milazzo Troina where: Outline ◮ S i , P i are molecules Introduction ◮ ℓ i , ℓ ′ i are stoichiometric coefficients Simulating chemical reactions ◮ k , k − 1 are the kinetic constants Gillespie’s algorithm Our algorithm A comparison For the law of mass action , the forward rate of a reaction is: Examples and applications dP Probabilistic MultiSet dt = k [ S 1 ] ℓ 1 · · · [ S ρ ] ℓ ρ Rewriting Simulation results Conclusions and the backward rate is: dS dt = k − 1 [ P 1 ] ℓ ′ 1 · · · [ P γ ] ℓ ′ γ 4/19
Gillespie’s stochastic approach (1) An Alternative to Gillespie’s Algorithm for Simulating Chemical Reactions ◮ assumes a stochastic reaction constant c µ for each Barbuti Maggiolo-Schettini Milazzo chemical reaction R µ Troina ◮ c µ dt is the probability that a particular combination of Outline reactant molecules of R µ react in an infinitesimal time Introduction interval dt Simulating chemical reactions Gillespie’s algorithm Our algorithm The total probability (denoted a µ dt ) of R µ to occur in the A comparison infinitesimal time interval dt is Examples and applications Probabilistic MultiSet Rewriting a µ dt = h µ c µ dt Simulation results Conclusions where h µ is the number of distinct molecular reactant combinations. 5/19
Gillespie’s stochastic approach (2) An Alternative to Gillespie’s Algorithm for Simulating Chemical Reactions Example: Barbuti solution with X 1 molecules S 1 and X 2 molecules S 2 Maggiolo-Schettini Milazzo Troina reaction R 1 : S 1 + S 2 → 2 S 1 Outline ◮ h 1 = X 1 X 2 Introduction ◮ a 1 = X 1 X 2 c 1 Simulating chemical reactions ◮ a 1 dt = X 1 X 2 c 1 dt Gillespie’s algorithm Our algorithm A comparison Examples and applications reaction R 2 : 2 S 1 → S 1 + S 2 Probabilistic MultiSet ◮ h 2 = X 1 ( X 1 − 1) Rewriting Simulation results 2 Conclusions ◮ a 2 = X 1 ( X 1 − 1) c 2 2 ◮ a 2 dt = X 1 ( X 1 − 1) c 2 dt 2 6/19
Gillespie’s algorithm An Alternative to Gillespie’s Algorithm for Simulating Given a set of reactions { R 1 , . . . R M } and a current time t Chemical Reactions 1. The time t + τ at which the next reaction will occur is Barbuti Maggiolo-Schettini randomly chosen with τ exponentially distributed with Milazzo parameter � M ν =1 a ν ; Troina 2. The reaction R µ that has to occur at time t + τ is Outline randomly chosen with probability a µ dt . Introduction Simulating chemical reactions At each step t is incremented by τ and the chemical solution Gillespie’s algorithm is updated. Our algorithm A comparison Examples and At each step the probability density function applications Probabilistic MultiSet Rewriting � M � Simulation results � P g ( τ, µ ) = exp − a ν τ · a µ dt Conclusions ν =1 gives the probability that the next reaction will occur in the time interval ( t + τ, t + τ + dt ) and will be R µ . 7/19
Our algorithm (1) An Alternative to Gillespie’s Algorithm for Simulating Chemical Reactions Barbuti Maggiolo-Schettini ◮ assumes that in a very small (fixed) time interval ∆ t at Milazzo Troina most one reaction may occur ◮ ∆ t depends on the number and on the rates of the Outline Introduction chemical reactions Simulating chemical reactions Gillespie’s algorithm Our algorithm Basic idea: A comparison ◮ divide the rate of each reaction (given by the law of Examples and applications mass action) by an arbitrarily great integer value N Probabilistic MultiSet Rewriting Simulation results ◮ use the result as the probability of each reaction to Conclusions occur in ∆ t 8/19
Our algorithm (2) An Alternative to Gillespie’s Algorithm for Simulating Given a set of reactions { R 1 , . . . , R M } , and assuming a Chemical Reactions volume of 1 litre Barbuti Maggiolo-Schettini Milazzo 1 ◮ ∆ t has to be fixed to Troina MN ◮ N is such that Outline Introduction 0 < k µ [ S µ 1 ] ℓ µ 1 · · · [ S µρ ] ℓ µρ ≤ 1 Simulating N chemical reactions Gillespie’s algorithm Our algorithm for 1 ≤ µ ≤ M and for all the possible concentrations A comparison (assumed to be finite) of S µ 1 , . . . , S µρ Examples and applications Probabilistic MultiSet Rewriting Simulation results The probability of R µ is Conclusions k µ [ S µ 1 ] ℓµ 1 ··· [ S µρ ] ℓµρ � if R µ can occur N P ( R µ ) = 0 otherwise 9/19
Our algorithm (3) An Alternative to Gillespie’s Algorithm for Simulating Chemical The algorithm iterates the following steps: Reactions Barbuti 1. A reaction R µ is randomly chosen (all the reactions are Maggiolo-Schettini Milazzo equiprobable); Troina 2. The chosen R µ is performed with probability P ( R µ ). Outline Introduction The probability of choosing and performing R µ in ∆ t is Simulating chemical reactions Gillespie’s algorithm P ( µ ) = 1 Our algorithm M P ( R µ ) A comparison Examples and applications and the probability of performing no reactions in ∆ t is Probabilistic MultiSet Rewriting Simulation results Conclusions M � P 0 = 1 − P ( ν ) ν =1 10/19
Comparing the two algorithms (1) An Alternative to Gillespie’s Algorithm for Simulating Chemical Reactions Barbuti Assume ∆ t infinitesimal. Results: Maggiolo-Schettini Milazzo ◮ The probability of performing R µ in Gillespie’s algorithm Troina is equivalent (with an approximation) to the probability Outline of choosing and performing R µ in our algorithm, that is Introduction Simulating a µ dt ≈ P ( µ ) chemical reactions Gillespie’s algorithm Our algorithm ◮ A step in Gillespie’s algorithm can be simulated by a A comparison Examples and sequence of steps in our algorithm having applications Probabilistic MultiSet (approximatively) the same probability Rewriting Simulation results Conclusions Approximations are introduced by deriving c µ from k µ 11/19
Comparing the two algorithms (2) An Alternative to Gillespie’s Algorithm for Simulating Chemical Gillespie’s algorithm Our Algorithm Reactions Barbuti Maggiolo-Schettini Milazzo – Assumes c µ + Uses k µ Troina in general unknown usually well–known Outline derived by k µ (approx) Introduction Simulating – Considers reactions + Is based on the law chemical reactions Gillespie’s algorithm individually of mass action Our algorithm A comparison does not depend on the Examples and applications measure units (scalability) Probabilistic MultiSet Rewriting Simulation results + Precise reaction times – Assumes at most one Conclusions reaction in ∆ t 12/19
Recommend
More recommend