Exact Stochastic Simulation of the Oregonator The Gillespie-Algorithm

(1)

Exact Stochastic Simulation of the Oregonator

The Gillespie-Algorithm

Marco Bussolera, Christoph Appel

Nonequilibrium Statistical Physics 21.03.2016

University of Augsburg

(2)

1

Introduction

2

Mathematical Basis

3

Stochastic Simulation Algorithm

4

Results for the Oregonator

(3)

Introduction

Two different approaches to calculate the time behaviour of a spatially homogeneous chemical system.

Deterministic approach: continuous time evolution for reaction-rate equations

→

set of coupled ordinary differential equations

Stochastic approach: discrete time evolution as random-walk governed by a single equation

→

Master-equation (mathematically complicated)

→

Exact Stochastic Simulation Algorithm (SSA) by Gillespie

(4)

Mathematical Basis

Description of the system

spatially homogeneous chemical system with fixed Volume V

N chemical species with X

_i

molecules (i = 1, 2, . . . , N)

M chemical reaction channels R

µ

with µ = 1, 2, . . . , M

thermal but no chemical equilibrium

(5)

Mathematical Basis

Stochastic formulation of chemical kinetics

Average probability that a combination of reactant molecules of R

µ

will react in the infinitesimal time interval (t, t +

dt)

= c

µdt

Probability that a reaction R

µ

occurs in (t, t +

dt)

= [number of distinct R

_µ

molecular combinations]

·

c

_µdt

≡

h

_µ

c

_µdt

≡

a

_µdt

(6)

Mathematical Basis

Reaction Probability Density Function

Probability that, given the state (X

1

, . . . , X

N

) at time t, the next reaction will occur in the time interval (t + τ, t + τ +

dτ

), and will be an R

_µ

reaction

= P (τ, µ)

dτ

= P

0

(τ )

·

a

µdτ

P

₀

(τ ) is the probability that no reaction occurs in (t, t + τ )

P

₀

(τ +

dτ

) = P

₀

(τ )

·

"

1

−

M

X

ν=1

a

_νdτ

#

resulting in

P

0

(τ ) =

exp −

M

X

ν=1

a

ν

τ

!

≡exp

(−a

₀

τ )

(7)

Stochastic Simulation Algorithm

Stochastic Simulation Algorithm (1)

Randomly generate a pair (τ, µ) according to the probability density function P (τ, µ)

Therefore generate a random pair (r

₁

, r

₂

) between 0 and 1 and assign:

τ = 1

a

0

·ln

1 r

1

µ−1

X

ν=1

a

ν

< r

2

a

0 ≤

µ

X

ν=1

a

ν

(8)

Stochastic Simulation Algorithm (2)

Step 0:

Input of population numbers X

i

and reaction constants c

µ

Step 1:

- Calculate h

₁

, . . . , h

_M

for current state - Calculate a

₁

= h

₁

c

₁

, . . . , a

_M

= h

_M

c

_M

- Calculate a

0

=

M

P

ν=1

a

ν

Step 2:

Generate a random pair (r

₁

, r

₂

) and calculate (τ, µ)

Step 3:

- Increase t by τ

- Update population numbers according to R

_µ Return to step 1

(9)

The Oregonator (1)

Chemical reactions:

X₁

+

Y₂ −→^c¹ Y₁ Y₁

+

Y₂ −→^c² Z₁ X₂

+

Y₁ −→^c³

2Y

1

+

Y₃

2Y

₁ −→^c⁴ Z₂ X₃

+

Y₃ −→^c⁵ Y₂

(10)

The Oregonator (2)

Reaction constants:

c

1

X

1

= ρ

1s

/Y

2s

c

2

= ρ

2s

/ (Y

1s

Y

2s

) c

₃

X

₂

= (ρ

₁_s

+ ρ

₂_s

) /Y

₁_s

c

4

= 2ρ

1s

/Y

1s

2

c

₅

X

₃

= (ρ

₁_s

+ ρ

₂_s

) /Y

₃_s

With the parameters:

ρ

1s ≡

c

1

X

1

Y

2s

ρ

2s ≡

c

2

Y

1s

Y

2s

(11)

Results for the Oregonator

Results: Oscillations, molecule numbers vs. time

0 2500 5000 7500 10000

0 1 2 3 4 5 6

numbersofmolecules

time (s)

Y1 Y2 Y3

0 2500 5000 7500 10000

1 1.5 2 2.5 3

numbersofmolecules

time (s)

Y1 Y2 Y3

Figure:These plots show Y₁(t),Y₂(t) andY₃(t). The initial values are equal to the steady-state values: Y₁_s =Y₁_i = 500, Y₂_s =Y₂_i = 1000, Y₃_s =Y₃_i = 2000, with the reaction rates: ρ1_s = 2000 andρ2_s = 50000. The simulation was done

× ⁵

(12)

Results: 3D-plot

0 2.5·10³ 5·10³ 7.5·10³ 1·10⁴

Y2

0 2.5·10³

5·10³ 7.5·10³

1·10⁴

0 Y1 2·10³ 4·10³ 6·10³ 8·10³ 1·10⁴

Y3

Figure:3D plot of the population numbersY₁,Y₂andY₃

(13)

Results: 2D-projections of the trajectory

0 2000 4000 6000 8000 10000

0 2000 4000

numberofY2molecules

0 2000 4000 6000 8000 10000

0 2000 4000

numberofY3molecules

0 2000 4000 6000 8000 10000

0 2000 4000 6000

numberofY3molecules

(14)

Results: 3D-plot for different initial values

0 2.5·10³ 5·10³ 7.5·10³ 1·10⁴

Y1

0 2·10³ 4·10³ 6·10³ 8·10³ 1·10⁴

Y3

0 2·10³

4·10³ 6·10³

8·10³ 1·10⁴ Y2

• Y_1i= 500, Y_2i= 1000,Y_3i= 2000

• Y_1i= 100, Y_2i= 500, Y_3i= 3000

• Y_1i= 1000,Y_2i= 1000,Y_3i= 1000

• Y1i= 1000,Y2i= 10, Y3i= 10

Figure:3D plot of the number of molecule of the species Y1,Y2andY3for four simulations with different initial values

(15)

Results: Closing the system to X

₁

-molecules

0 2500 5000 7500 10000

0 2 4 6 8

numbersofmolecules

time (s)

Y1 Y2 Y3 X1

Figure:This plot shows the behaviour of the system when the number of X₁

(16)

Results: Closing the system to X

₂

-molecules

0 2500 5000 7500 10000

0 2 4 6 8

0 25000 50000 75000 100000

numbersofmoleculesY1,Y2,Y3

time (s)

numbersofmoleculesX2

Y1 Y2 Y3 X2

Figure:This plot shows the behaviour of the system when X₂is not kept

constant, starting with an initial value ofX₂= 10⁵. In this case the amplitude of the oscillation decreases for decreasingX₂.

(17)

Results: Varying the reaction rates (1)

0 2500 5000 7500 10000

0 1 2 3 4 5 6

numbersofmolecules

time (s)

Y1 Y2 Y3

0 1000 2000 3000 4000 5000

Y2

0 1000 2000 3000 4000 5000

Y1 0

2·10³ 4·10³ 6·10³ 8·10³ 1·10⁴

Y3

Figure:Results after changing the reaction rates from ρ₁_s = 2000 and ρ₂_s = 50000 toρ₁_s = 2000 and ρ₂_s = 20000. The amplitude of the oscillation decreases by decreasingρ2.

(18)

Results: Varying the reaction rates (2)

0 2500 5000 7500 10000

0 1 2 3 4 5 6

numbersofmolecules

time (s)

Y1 Y2 Y3

0 1000 2000 3000 4000 5000

Y2

0 1000 2000 3000 4000 5000

Y1 0

2·10³ 4·10³ 6·10³ 8·10³ 1·10⁴

Y3

Figure:Results after changing the reaction rates to ρ₁_s = 4000 andρ₂_s = 20000.

The amplitude of the oscillation further decreases by increasing ρ₁_s.

(19)

Results: Varying the reaction rates (3)

0 2500 5000 7500 10000

0 1 2 3 4 5 6

numbersofmolecules

time (s)

Y1 Y2 Y3

0 1000 2000 3000 4000 5000

Y2

0 1000 2000 3000 4000 5000

Y1 0

2·10³ 4·10³ 6·10³ 8·10³ 1·10⁴

Y3

Figure:Results after changing the reaction rates to ρ1_s = 5000 andρ2_s = 5000.

By increasingρ1_s and decreasingρ2_s even further no oscillation occurs anymore and the molecule population numbers only fluctuate around a constant value.

(20)

Results: Varying the number of molecules

0 1·10⁴ 2·10⁴ 3·10⁴

4 6 8 10 12

numbersofmolecules

time (s)

Y1 Y2 Y3

0 5·10³ 1·10⁴ 1.5·10⁴ 2·10⁴ 2.5·10⁴

Y2

0 5·103 1·10⁴ 1.5·10⁴ 2·10⁴ Y1 0

1·10⁴ 2·10⁴ 3·10⁴ 4·10⁴

Y3

Figure:Results after increasing the steady state population levels by a factor of four: Y₁_s = 2000, Y₂_s = 4000, Y₃_s = 8000,with reaction ratesρ₁_s = 2000 and ρ₂_s = 50000. This leads to a decrease in the relative size of the stochastic fluctuations in the time-dependent oscillations. The “spreading” of the limit cycle path is not reduced.