Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
IV121: Computer science applications in
biology
Quantitative Models in Biology
David ˇSafr´anek
March 5, 2012
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Obsah
Continuous mass action
Stochastic mass action
Beyond elementary reaction kinetics
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
What is continous? What does continuous mean?
ˆ think of physical motion
ˆ by means of classical mechanics
ˆ by means of classical electrodynamics
ˆ all are models. . .
ˆ compare with quantum mechanics – the scale of 10−8
m makes
the barrier between views. . .
ˆ think of a crowd of thousand people
ˆ what you observe when someone disappears?
ˆ what you observe when someone new appears?
ˆ think of molecules in a solution or in a cell . . .
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Continuous model of reaction kinetics
ˆ assume well-stirred solution
ˆ high amounts of all substances
ˆ fixed thermodynamics conditions (temperature, pressure, . . . )
ˆ fixed volume of the solution
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Continuous model of reaction kinetics
A
k
−→
ˆ consider a barrel with a substance A of molar volume [A] [M]
ˆ how much of substance A “flows out” per a single time unit?
ˆ value proportional to [A](t) in a given time t
−
d[A](t)
dt
= k · [A](t)
ˆ coefficient of proporcionality is denoted k [s−1
]
so-called kinetic constant (coefficient)
- determines the speed of mass decay (“outflow”)
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Continuous model of reaction kinetics
[A](t)
dt
= k · [A](t)
ˆ which functions has the same form as its derivation?
ˆ f (t) = 1 + t + t2
/2! + t3
/3! + t4
/4! + ...
f (t) = et
ˆ plat´ı
det
dt
= et
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Continuous model of reaction kinetics
A
k
−→
−d[A](t)
dt = k · [A](t)
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Continuous model of reaction kinetics
A
k
−→
−d[A](t)
dt = k · [A](t)⇔ [A](t) = [A](0) · e−kt
ˆ so-called first-order kinetics (a special case of mass action)
ˆ linear autonomous differential equation
ˆ unique solution
ˆ can be either analytically solved or numerically approximated
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Continuous model of reaction kinetics
ˆ state is a vector of actual amounts of all susbtances in the
system
ˆ continuous-time dynamics: the state change X(t) → X(t + dt)
updates all components of X (continuous concurrent flow of
all reactions)
ˆ we consider reaction rate as a function of time: for a
reaction R in time t we denote the actual rate vR(t)
reaction type rate function vR state update
→ A vR(t) = k dA
dt = −vR
A → B vR(t) = k · [A](t) dA
dt = −vR, dB
dt = vR
A + B → AB vR(t) = k · [A](t) · [B](t) dA
dt = dB
dt = −vR, dAB
dt = vR
2A → AA vR(t) = k · [A]2 dA
dt = −2vR, dAA
dt = vR
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
General Mass Action Kinetics
k
i=1
si · Ai →
l
j=1
pj · Bj
v =
k
i=1
Ai
si
∀1 ≤ i ≤ k. dAi
dt = −si · v
∀1 ≤ j ≤ l.
dBj
dt = pj · v
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Example: Michaelis-Menten
S + E
k1
k2
ES
k3
−→ P + E
d[S]
dt
= −k1[E][S] + k2[ES]
d[E]
dt
= −k1[E][S] + k2[ES] + k3[ES]
d[ES]
dt
= k1[E][S] − k2[ES] − k3[ES]
d[P]
dt
= k3[ES]
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Euler method
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Euler method
ˆ approximate solution y(t) (Euler):
y (t) = f (t, y(t))
y(0) = y0
ˆ exact solution ϕ(t):
ϕ (t) = f (t, ϕ(t))
ϕ(0) = y0
ˆ for each n ≥ 0, tn = n∆t:
yn ≈ ϕ(tn)
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Euler method
Exact solution ϕ(t) satisfies:
ϕ(tn+1) = ϕ(tn) +
tn+1
tn
ϕ (t)dt
= ϕ(tn) +
tn+1
tn
f (t, ϕ(t))dt
Numerical approximation:
yn+1 = yn + σ
where
σ ≈
tn+1
tn
f (t, ϕ(t))dt
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Euler method I
dy
dt
= f (t, y)
yn+1 = yn + ∆t · f (tn, yn)
1. init t0, y0, ∆t, n;
2. for j from 1 to n do
2.1 m := f (t0, y0);
2.2 y1 := y0 + ∆tm;
2.3 t1 := t0 + ∆t;
2.4 t0 := t1;
2.5 y0 := y1;
3. end
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Petri Net Analysis Framework
quantitative parameters ignored
quantitative parameters required
continuous model
qualitative model
stochastic model
variables
continuous
abstracted
modeled
time
discrete
approximation
abstraction abstraction
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Petri Net Analysis Framework
continuous model
qualitative model
stochastic model
variables
continuous
abstracted
modeled
time
discrete
Monte Carlo simulation
Static analysis
Behavioral analysis
Simulation analysis
Steady state analysis
Numerical simulation
Simulation analysis
approximation
abstraction abstraction
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Petri Net Representation of Models
Continuous Petri Net
ODE
Reaction Network
ˆ for mass action kinetics both transformations are direct
ˆ unique Petri net representation of ODEs always achievable
S. Soliman, M. Heiner (2010) “A Unique Transformation from Ordinary Differential Equations to Reaction
Networks.” PLoS ONE 5(12): e14284. doi:10.1371/journal.pone.0014284
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Continuous Petri Nets
Structure
Continuous Petri net is a quadruple N = S, R, f , v, m(0) where
ˆ S finite set of places (substances),
ˆ T finite set of transitions (reactions),
ˆ f : ((P × T) ∪ (T × P)) → N0 set of weighted edges,
ˆ x• = {y ∈ S ∪ R | f (x, y) = 0} denotes target of x
ˆ •x = {y ∈ S ∪ R | f (y, x) = 0} denotes source of x
ˆ weight represents stoichiometric coefficients
ˆ v is mapping which assigns each transition r ∈ R a function
hr : R|•r|
→ R
ˆ v represents transition (reaction) rate
ˆ m(0) : S → R+
0 is initial marking (initial condition).
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Continuous Petri Nets
Dynamics
Number of places denotes the dimension of the system, n = |S|.
Each place s ∈ S is marked by a value in R+
0 (representing
concentration of the respective species):
ˆ in Petri net terminology evaluation of places is called marking
and represented as an n-dimensional vector m ∈ Rn
ˆ marking evolves in time: m(t)
Dynamics of each place s ∈ S is defined by an ODE:
dms(t)
dt
=
r∈•s
f (r, s)v(r) −
r∈s•
f (s, r)v(r)
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Continuous Petri Nets
Michaelis-Menten Mass Action Kinetics Example
r1 : E + S → ES, r2 : ES → E + S, r3 : ES → P + E
dmS
dt = v(r2) − v(r1)
dmE
dt = v(r2) + v(r3) − v(r1)
dmES
dt = v(r1) − v(r2) − v(r3)
dmP
dt = v(r3)
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Continuous Petri Nets
Michaelis-Menten Mass Action Kinetics Example
E
ES
P
1*ES
1*ES
0.1*E*S
10
50
S
dmS
dt = k2mES − k1mE mS
dmE
dt = k2mES + k3mES − k1mE mS
dmES
dt = k1mE mS − k2mES − k3mES
dmP
dt = k3mES
k1 = 0.1 [M−1s−1]
k2 = 1 [s−1]
k3 = 1 [s−1]
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Parameter Estimation Problem
ˆ inverse problems – determine the model from measured data
ˆ quite easy for linear systems, but what for non-linear?
ˆ general steps in inverse problem solution:
1. identify relations among variables
2. identify functions describing relations (e.g., mass action)
3. estimate constants appearing in the functions – parameter
estimation
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Parameter Estimation Problem
ˆ parameter estimation is solved as optimization problem w.r.t.
measured data
ˆ the goal is to minimize average deviation of model from data
ˆ so-called least squares method
ˆ we seek for global minima
ˆ many heuristics for optimization procedure, many algorithms
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Parameter Landscape
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Parameter Landscape
Walking the landscape to find the global minimum
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Tool Support for Continuous Models
ˆ format transformation and models editing
ˆ Snoopy – Petri nets representation visual editing, SBML
export/import, Petri net variants transformation
ˆ CellDesigner – SBGN visual editing, SBML export/import
ˆ CellIllustartor – visual editing, hybrid Petri nets simulation
ˆ analysis
ˆ Octave, Matlab (simulation and SBML: SBMLToolBox,
SimBiology ToolBox)
ˆ COPASI (simulation, SBML export/import, other analysis
tasks)
ˆ BioCHAM (robustness analysis, model checking)
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Literature
M. Feinberg. Lectures on Chemical Reaction Networks.
http://www.che.eng.ohio-state.edu/~FEINBERG/
LecturesOnReactionNetworks/
M. Heiner, D. Gilbert & R. Donaldson. Petri Nets for Systems
and Synthetic Biology. SFM 2008: 215-264
Hoops S. et al. COPASI – a COmplex PAthway SImulator.,
Bioinformatics 22, 3067-74
T. Vejpustek. Parameter estimation v COPASI – tutotial.
http://anna.fi.muni.cz/~xsafran1/PV225/parameter_
estimation/copasi.html.
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Obsah
Continuous mass action
Stochastic mass action
Beyond elementary reaction kinetics
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Stochastic model of reaction kinetics
ˆ assume well-stirred solution
ˆ uniform distribution of molecules in solution
ˆ low amounts of substances
ˆ fixed thermodynamics conditions (temperature, pressure, . . . )
ˆ fixed volume of the solution
ˆ reactions (molecule collisions) viewed as discrete events
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Stochastic model of reaction kinetics
ˆ discrete events happen in continuous time
ˆ time between two events is a stochastic quantity
ˆ average probability (over the whole solution) of reaction
realization in given time
ˆ some reactions faster (more probable), some slower (less
probable)
ˆ probability depends on amounts of reactant molecules
ˆ stochasticity is a measure of uncertainty caused by other
(non-reactive) events happening in solution
⇒ approximation of the following aspects:
ˆ molecule position and rotation
ˆ molecule motion (speed)
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Stochastic model of reaction kinetics
Gillespie’s Hypothesis
D. T. Gillespie. Exact Stochastic Simulation of Coupled Chemical Reactions. In Journal of Physical Chemistry,
volume 81, No. 25, pages 2340-2381. 1977.
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Stochastic model of reaction kinetics
Gillespie’s Hypothesis
ˆ basic Newtonian physics and thermodynamics is assumed
ˆ realization probability for reaction Rj globally characterized by
the rate constant cj
ˆ depends on radii of colliding molecules and their average
relative velocities (considerred relatively to the solution
volume)
ˆ direct function of temperature and molecular structure
ˆ if a pair of two molecules has kinetic energy higher than
reaction energy then the reaction is realized
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Stochastic model of reaction kinetics
Comparison of models
ˆ continuous kinetics provides a macro-scale view
ˆ systems view abstracting from location (space)
ˆ continuous dynamics of large quantities – quantitity as a
population
ˆ single average evolution of averaged (well-stirred) events
ˆ stochastic kinetics provides a meso-scale view
ˆ systems view still abstracting from location (space)
ˆ discrete dynamics of low quantities
ˆ all possible evolutions of averaged (well-stirred) events
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Stochastic model of reaction kinetics
Grand probability function
ˆ Gillespie’s hypothesis enables stochastic formulation of
molecular (low population) dynamics
ˆ for time t the grand probability function Pr(X; t)
characterizes the probability that there will be present Xi
molecules of species Si , X = X1, ..., Xn is a vector quantity
denoting configuration of the population
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Stochastic model of reaction kinetics
Grand probability function
ˆ Gillespie’s hypothesis enables stochastic formulation of
molecular (low population) dynamics
ˆ for time t the grand probability function Pr(X; t)
characterizes the probability that there will be present Xi
molecules of species Si , X = X1, ..., Xn is a vector quantity
denoting configuration of the population
ˆ how to compute?
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Stochastic model of reaction kinetics
Grand probability function
ˆ considering reactions as discrete events leads to:
Pr(X; t + dt) = Pr(X; t) · Pr(no state change)
+ m
i=1 Pr(X − ui ; t) · Pr(state changed to X)
where
ˆ dt is a small time step in which at most 1 reaction occurs
ˆ ui is update caused by the effect of reaction Ri (X → X + ui )
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Stochastic model of reaction kinetics
Grand probability function
ˆ considering reactions as discrete events leads to:
Pr(X; t + dt) = Pr(X; t) · (1 − m
i=1 χi (X)dt)
+ m
i=1 Pr(X − ui ; t)χi (X − ui )dt
where
ˆ dt is a small time step in which at most 1 reaction occurs
ˆ ui is update caused by the effect of reaction Ri (X → X + ui )
ˆ χi is hazard function characterizing the probability of exactly
one occurrence of Ri in time interval dt
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Stochastic model of reaction kinetics
Hazard function
In a particular configuration, probability of reaction realization in
given time is characterized by hazard function.
Hazard function for reaction R is denoted χR(X) where X is a
current state (configuration, marking). Assume R is assigned a
stochastic rate constant cR. The table below shows the hazard
function for all forms of elementary reactions:
→ ∗ χR(X) = cR
A → ∗ χR(X) = cR · XA
A + B → ∗ χR(X) = cR · XA · XB
2A → ∗ χR(X) = cR · XA·(XA−1)
2
stochastic mass action
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Stochastic model of reaction kinetics
Stochastic Simulation Algorithm – SSA
ˆ single transition X(t) → X(t + dt) updates just one
component of X
ˆ realization of just one reaction R
ˆ reaction realization does not take time
ˆ in a state X(t), the time to next realization of reaction Ri is
characterized by distribution Exp(χRi
(X))
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Exponential distribution
E(X) =
1
λ
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Stochastic model of reaction kinetics
ˆ for transition X(t) → X(t + dt), dt is the time to earliest
reaction event
ˆ dt is sampled as minimal time over all n reactions:
dt ∼ Exp(χ(X)) χ(m) =
n
i=1
χRi
(X)
ˆ reaction Ri is chosen with probability: P(Ri ) =
χRi
(X)
χ(X)
ˆ formally this comes from the property of exponential
distribution
ˆ the model behind is continuous-time Markov process
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Gillespie direct method (SSA)
Output: a single trajectory realizing the grand probability distribution
1. initialize t = 0, X(0)
2. compute χRi
(X) ∀i ∈ {1, ..., n}
3. compute χ(X) ≡
n
i=1 χRi
(X)
4. sample τ ∈ Exp(χ(X))
5. t := t + τ
6. choose reaction Ri with probability
χRi
(X)
χ(X)
7. update: X(t) = X(t − τ) + uRi
8. while t < Tmax go to (2)
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Example
Consider reaction: A → B
A B
5
0
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Example
Consider reaction: A → B
A B
2s
T~Exp(5)
5
0
→
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Example
Consider reaction: A → B
A B
5
0
→ (2s) →
4
1
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Example
Consider reaction: A → B
A B
0.5s
T~Exp(4)
5
0
→ (2s) →
4
1
→
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Example
Consider reaction: A → B
A B
5
0
→ (2s) →
4
1
→ (0.5s) →
3
2
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Example
Consider reaction: A → B
A B
1s
T~Exp(3)
5
0
→ (2s) →
4
1
→ (0.5s) →
3
2
→
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Example
Consider reaction: A → B
A B
5
0
→ (2s) →
4
1
→ (0.5s) →
3
2
→ (1s) →
2
3
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Example
Consider reaction: A → B
A B
T~Exp(2)
0.8s
5
0
→ (2s) →
4
1
→ (0.5s) →
3
2
→ (1s) →
2
3
→
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Example
Consider reaction: A → B
A B
5
0
→ (2s) →
4
1
→ (0.5s) →
3
2
→ (1s) →
2
3
→
(0.8s) →
1
4
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Example
Consider reaction: A → B
2s
T~Exp(1)
A B
5
0
→ (2s) →
4
1
→ (0.5s) →
3
2
→ (1s) →
2
3
→
(0.8s) →
1
4
→
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Example
Consider reaction: A → B
A B
5
0
→ (2s) →
4
1
→ (0.5s) →
3
2
→ (1s) →
2
3
→
(0.8s) →
1
4
→ (2s) →
0
5
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Example
Consider reaction: A → B
5
0
→ (2s) →
4
1
→ (0.5s) →
3
2
→ (1s) →
2
3
→
(0.8s) →
1
4
→ (2s) →
0
5
ˆ hazard function considered: χ(X) = 1 · XA
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Example – positive autoregulation of gene expression
R1 : P + g → gP
R2 : gP → g + P
R3 : gP → gP + P
R4 : P →
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Example – positive autoregulation of gene expression
ˆ initial settings: g(0) = 5, P(0) = 2, gP(0) = 0;
c1 = c2 = 1, c3 = 0.1, c4 = 0.01
ˆ distribution in t = 1000 for 2000 simulations
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Petri Net Analysis Framework
Stochastic Model
continuous model
qualitative model
stochastic model
variables
continuous
abstracted
modeled
time
discrete
Monte Carlo simulation
Static analysis
Behavioral analysis
Simulation analysis
Steady state analysis
Numerical simulation
Simulation analysis
approximation
abstraction abstraction
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Discrete Approximation
ˆ notion of Petri Net token
ˆ token represent molecule or a certain concentration level
ˆ suppose bounded concentration for all substrates:
0, max) ⊂ R
ˆ uniform partitioning into N intervals:
0, (0, 1 ·
max
N
, (1 ·
max
N
, 2 ·
max
N
, . . . , (N − 1 ·
max
N
, N ·
max
N
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Discrete Approximation
Stochastic vs. continuous model
ˆ substance concentration [M]:
c =
n
V
where n substance quantity [mol], V solution volume [l]
ˆ expressed in terms of Avogadro constant (number of particles
in 1 mol):
c =
N
NA · V
where NA Avogadro constant [mol−1], V solution volume [l]
and N number of molecules.
ˆ transformation factor:
γ = NA · V [mol−1
l] ⇒ N = c · γ, c =
N
γ
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Discrete Approximation
Stochastic vs. continuous model
A continuous Petri net N = S, R, f , v, m(0) can be transformed
to a stochastic Petri net N = S, R, f , v , m(0) :
ˆ m(0) : S → N0 is initial marking
ˆ v assigns each transition a hazard:
reaction type r ∈ R v(r) → v (r) transformation
→ A v (r) = v(r)
A → B v (r) = v(r)
A + B → AB v (r) = v(r)
γ
A + A → AA v (r) = 2v(r)
γ
L. Cardelli (2008), “From Processes to ODEs by Chemistry”. In 5th International Conference on
Theoretical Computer Science, pages 261-281. DOI:10.1007/978-0-387-09680-3 18
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Stochastic Petri Nets
Michaelis-Menten Stochastic Mass Action Kinetics Example
E
ES
P
1*ES
1*ES
0.1*E*S
10
50
S
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Tool Support
ˆ Monte Carlo simulation
ˆ Dizzy, COPASI, SPiM
ˆ simulation analysis
ˆ BioNessie (statistical model checking)
ˆ BioCHAM
ˆ symbolic analysis
ˆ PRISM (strong transient and steady state analysis)
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Petri Net Analysis Framework
Tool Support Overview
continuous model
qualitative model
stochastic model
variables
continuous
abstracted
modeled
time
discrete
Monte Carlo simulators
COPASI, STOC2, SPIM
Monte Carlo analyzers
BioNessie, BioCHAM
Symbolic analyzers
PRISM
ODE solvers
BioNessie, BioCHAM
ODE dynamics analyzers
Static analyzers
Charlie, PIPE, COPASI
Dynamic analyzers
Charlie, BioCHAM
approximation
abstraction abstraction
Continuous mass action Stochastic mass action Beyond elementary reaction kinetics
Literature
M. Heiner, D. Gilbert & R. Donaldson. Petri Nets for Systems
and Synthetic Biology. SFM 2008: 215-264
D. Wilkinson. Stochastic Modelling for Systems Biology,
second edition., Chapman & Hall/CRC, 2011.
D. T. Gillespie. Exact Stochastic Simulation of Coupled
Chemical Reactions. In Journal of Physical Chemistry, volume
81, No. 25, pages 2340-2381. 1977