Transcript Introduction to ODE Modeling

Introduction to Stochastic
Models
Shlomo Ta’asan
Carnegie Mellon University
Download StochasticModels.zip
tsb.mssm.edu/summerschool/index.php/Introduction_to_stochastic_modeling
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
100*rand 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 < P
cell dies
end
The models
Reactions
Random Walks
Chemotaxis
Macrophage Looking for Bacteria
Stochastic Models for Reactions
Modeling Reactions
A -> 0
per event A count goes down by 1.
A -> B
per event A count goes down by 1, B count goes up by 1
A + B -> 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
Our first real ABM: A-> 0
A cell dies with probability p during time interval dt.
P(A-> 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
Efficient Probabilistic Model: A-> 0
A cell dies with probability p during time interval dt.
P(A-> 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 -> Na -1 (when one of the agents dies)
P( Na -> Na – 1) = p*Na*dt
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
Stochastic-SIR
Population has three groups:
Susceptible, Infected and Recovered
The dynamics is expressed in the reactions
S + I -> I + I (probability: r*dt) event 1
I -> R
(probability: a*dt)
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( S+I -> I + I) = r*Ns*Ni*dt
P(I -> R)
= a*dt*Ni
event 1 happens just once during dt
event 2 happens just once during dt
Since probabilities are < 1, dt should be small enough. More precisely,
We want the sum of all probabilities to be less than 1. !!!
Stochastic-SIR cont.
Each event changes the variables of the model Ns, Ni, Nr.
Event 1:
Event 2:
S + I -> I + I
I -> R
 Ns -> Ns -1 and Ni -> Ni + 1
 Ni -> Ni – 1 and Nr -> Nr + 1
P( Ns -> Ns -1 & Ni -> Ni + 1) = r*Ns*Ni*dt
P( Ni -> Ni -1 & Nr -> Nr + 1) = a*dt*Ni
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.
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!!
Modeling cell movement
Random Walks
A basic random walk in 2D:
P( x-> x+1, y->y) = ¼
P( x-> x-1, y->y) = ¼
P( x-> x, y->y-1) = ¼
P( x-> x, y->y+1) = ¼
¼
¼
¼
¼
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
Stochastic Modeling of Chemotaxis
1st Ingredient:
Bacterium performs a biased random walk based on
the gradient (gX,gY) of the chemoattractant molecules.
P( x-> x+1, y->y) = ¼ + gX
P( x-> x-1, y->y) = ¼ - gX
P( x-> x, y->y+1) = ¼ + gY
P( x-> x, y->y-1) = ¼ - gY
Simulation similar to previous random walk..
¼ +gY
¼ -gX
¼ +gX
¼ -gY
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.6*Cold(i,j) + 0.1*(Cold(i+1,j)+Cold(i-1,j)+Cold(I,j+1)+Cold(i,j-1))
Matlab code: Chemotaxis.m
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.6*Cold(i,j) + 0.1*(Cold(i+1,j)+Cold(i-1,j)+Cold(I,j+1)+Cold(i,j-1))
3rd Ingredient: Macrophage performs a random walk based on
chemoattractant molecules and kills bacteria as they are reached.
¼ +gY
¼ -gX
¼ +gX
Matlab code: Phagocytosis.m
¼ -gY
The End