Introduction to Stochastic Models

1
Introduction to Stochastic Models

• Shlomo Taasan
• Carnegie Mellon University

3nd Biodefense Summer School, August 2008
2
How to work with Probabilities?
Simplest model Tossing a fair coin P (Head)
1/2 P (Tail) 1/2 How to toss a coin using
matlab? rand - generates a random number
between 0 and 1
rand(100,1) generate 100 random
numbers 100rand generate a random number
between 0 and 100 (decimal numbers) We also
need to know how to do something (kill the cell)
with probability P r rand if r lt P cell
dies end

3
The models

• Reactions
• Random Walks
• Chemotaxis
• Macrophage Looking for Bacteria

4

• Stochastic Models for Reactions

5
Modeling Reactions
A -gt 0 per event A count goes down by 1. A -gt
B per event A count goes down by 1, B count goes
up by 1 A B -gt C per event A,B count go down
by 1, C count goes up by 1 Model in two ways
1. ABM modeling of each molecule, 1 it
exists, 0 is it deleted. for each agent we
perform an action with some probability 2.
probabilistic model for total number of molecules
in each group. Na, Nb, etc the total number of
molecules of type A,B, we change Na, Nb etc
by one according to some probability

6
Our first real ABM A-gt 0
A cell dies with probability p during time
interval dt. P(A-gt 0 dt) p dt For
simulation we need a counter for time, and a
variable to monitor if the cell is dead or alive.
Use a variable x for this. x 1 cell is
alive x 0 cell is dead What if we have many
cells? What are the questions we can ask about
this model? population size as a function of
time, average time to extinction, fluctuation
in time to extinction ? Matlab code
realABM_A.m

7
Efficient Probabilistic Model A-gt 0
A cell dies with probability p during time
interval dt. P(A-gt 0 dt) p dt We model
the total number of agents, Na, instead of each
agent. We pick dt small enough such that only
one cell can die (i.e the probability that two
die is very small, so we ignore it) The only
event in this model Na -gt Na -1 (when one of
the agents dies) P( Na -gt Na 1)
pNadt This is much faster than the real ABM
model. Sometimes it is enough to do this type of
models. Sometimes we need a full ABM model Lets
do the simulation! Matlab Code ABM_A.m

8
Stochastic-SIR
Population has three groups Susceptible,
Infected and Recovered The dynamics is expressed
in the reactions S I -gt I I
(probability rdt) event 1 I -gt R
(probability adt) event 2 When to expect
epidemic? A relation between parameters Ns, Ni,
Nr population size for Susceptible, Infected and
Recovered. We pick dt small enough such that only
one event will happen, The probability of events
P( SI -gt I I) rNsNidt event 1 happens
just once during dt P(I -gt R)
adtNi event 2 happens just once during
dt Since probabilities are lt 1, dt should be
small enough. More precisely, We want the sum
of all probabilities to be less than 1. !!!

9
Stochastic-SIR cont.
Each event changes the variables of the model Ns,
Ni, Nr. Event 1 S I -gt I I ? Ns -gt
Ns -1 and Ni -gt Ni 1 Event 2 I -gt R
? Ni -gt Ni 1 and Nr -gt Nr 1 P( Ns -gt
Ns -1 Ni -gt Ni 1) rNsNidt P( Ni -gt Ni
-1 Nr -gt Nr 1) adtNi Now lets simulate
it. Matlab Code ABM_SIR.m How does it
compare with the ODE model? Notice the finite
time to eradicate infection in this model.
Compare to the ODE.

10
Stochastic vs. ODE
Reaction A 2 X ?? 3X A ?? X Multiple steady
states! ODE Stochastic model

Stochastic models do more than just adding noise
to results of an ODE!!
11
Modeling cell movement
12
Random Walks
A basic random walk in 2D P( x-gt x1, y-gty)
¼ P( x-gt x-1, y-gty) ¼ P( x-gt x, y-gty-1) ¼
P( x-gt x, y-gty1) ¼ Pick a random number, r.
if r between 0 and ¼ go right If r between ¼ and
½ go left, if r is between ½ and ¾ go down, if r
between ¾ and 1 go up. Questions If we start
at (0,0), where will we be after n steps? in
different places each time How to describe
this? Probability distribution Lets simulate
it Matlab code randomWalk.m - ABM for
bacterial movement

13
Stochasic Modeling of Chemotaxis
1st Ingredient Bacterium performs a biased
random walk based on the gradient (gX,gY) of the
chemoattractant molecules. P( x-gt x1, y-gty)
¼ gX P( x-gt x-1, y-gty) ¼ - gX P( x-gt x,
y-gty1) ¼ gY P( x-gt x, y-gty-1) ¼ -
gY Simulation similar to previous random walk..
2nd Ingredient Chemoattractant molecules
diffuse (spread) we model them using
concentration at every lattice point
C(I,j) Evolution an average between value at
(I,j) and values at nbrs Cnew(i,j)
0.6Cold(i,j) 0.1(Cold(i1,j)Cold(i-1,j)Cold(
I,j1)Cold(i,j-1)) Matlab code Chemotaxis.m

14
A Macrophage looking for Bacteria
1st Ingredient Bacterium performs a simple
random walk and releases chemoattractant
molecules. 2nd Ingredient Chemoattractant
molecules diffuse Cnew(i,j) 0.6Cold(i,j)
0.1(Cold(i1,j)Cold(i-1,j)Cold(I,j1)Cold(i,j-
1)) 3rd Ingredient Macrophage performs a random
walk based on chemoattractant molecules and kills
bacteria as they are reached. Matlab code
Phagocytosis.m

15
(No Transcript)
16
(No Transcript)
17
Plan
We will learn in this tutorial translating
biological knowledge to differential equations
models use matlab to simulate models generate
graphs, make predictions, ... In particular
we will model Reactions, Trafficking, Simple
infections Participants will also have an
opportunity to play with more complex
models Voltera-Lotka periodic solutions

Main material - Matlab programs http//tsb.mssm.ed
u/summerschool/images/c/cb/StochasticLab.zip
18
An HIV model
T target cells (CD4 T cels) I infected cells V
virus Model assumptions -gt T (lambda)
target cells production T -gt 0 (d)
target cells natural death T V -gt I V (k)
target cell becomes infected by virus I -gt 0
(delta) infected cells death I -gt I V (p)
virus replication in infected cells V -gt 0
(c) virus clearance Ex Construct a
probabilistic model. Follow the ABM_SIR
model Step 1 Define variables NT, NI, NV. Step
2 determine probabilities for each event. Step
3 Determine changes of these variables per each
event. Step 4 simulate. More involved models
NatureReviewsPereleson.pdf