Introduction to Stochastic Models

1. Introduction to Stochastic Models Shlomo Ta’asan Carnegie Mellon University

2. Download StochasticModels.ziptsb.mssm.edu/summerschool/index.php/Introduction_to_stochastic_modeling

3. 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

4. The models Reactions Random Walks Chemotaxis Macrophage Looking for Bacteria

5. Stochastic Models for Reactions

6. 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

7. 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

8. 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

9. 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 event 1 happens just once during dt P(I -> R) = a*dt*Ni 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

10. Each event changes the variables of the model Ns, Ni, Nr. Event 1: S + I -> I + I Ns -> Ns -1 and Ni -> Ni + 1 Event 2: I -> R  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-SIR cont.

11. 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!!

12. Modeling cell movement

13. ¼ ¼ ¼ ¼ 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

14. ¼ +gY ¼ -gX ¼ +gX ¼ -gY 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.. 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

15. ¼ +gY ¼ -gX ¼ +gX ¼ ¼ -gY ¼ ¼ ¼ 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. Matlab code: Phagocytosis.m

16. The End