Python Programs for Modelling Infectious Diseases book:Chapter 6:Program 6.4

From DeductiveThinking Wiki
(Redirected from Program 6.4)
Jump to: navigation, search

Program 6.4 is a SIR model with demographic stochasticity (page 203 of the book). These are the equations and the code of the model:

Equations

 \frac{dX}{dt} = \mu*N-\beta*X*Y/N -\mu*X

 \frac{dY}{dt} = \beta*X*Y/N-\gamma*Y-\mu*Y

 \frac{dZ}{dt} = \gamma*Y-\mu*Z

 \mbox{Infection }\ \ \ \mbox{Rate: }\beta*X*Y/N \ \ \ \ \mbox{Event: } X \mapsto X-1  \ Y \mapsto Y+1

 \mbox{Recovery }\ \ \ \mbox{Rate: }\gamma*Y \ \ \ \ \ \ \ \ \ \ \ \ \ \mbox{Event: }Y \mapsto Y-1 \ Z \mapsto Z+1

 \mbox{Birth }\ \ \ \ \ \ \ \ \mbox{Rate: } \mu*N \ \ \ \ \ \ \ \ \ \ \ \ \ \mbox{Event: } X \mapsto X+1

 \mbox{Death Susceptible } \mbox{Rate: } \mu*X\ \ \ \ \mbox{Event: } X \mapsto X-1

 \mbox{Death Infected }\ \ \ \ \mbox{Rate: } \mu*Y\ \ \ \ \ \mbox{Event: } Y \mapsto Y-1

 \mbox{Death Recovered } \ \mbox{Rate: } \mu*Z  \ \ \ \ \ \mbox{Event: } Z \mapsto Z-1

Code

Program 6.4: SIR model with demographic stochasticity (one run - Population size N=5000)
Program 6.4: SIR model with demographic stochasticity (one run - Population size N=50)
#!/usr/bin/env python

####################################################################
###    This is the PYTHON version of program 6.4 from page 203 of  #
### "Modeling Infectious Disease in humans and animals"            #
### by Keeling & Rohani.					   #
###								   #
### It is the SIR model (including births and deaths) with full    #
### (event-driven) demographic stochasticity.			   #
###								   #
### This is a more complex stochastic model as 6 events are	   #
### possible: infection, recovery, birth, death of susceptible,    #
### death of infected, death of recovered.			   #
### Note: by default we are using a very small population size 	   #
### to highlight the stochasticity.				   #
####################################################################

##########################################################################
### Copyright (C) <2008> Ilias Soumpasis                                 #
### ilias.soumpasis@deductivethinking.com                                #
### ilias.soumpasis@gmail.com	                                         #
###                                                                      #
### This program is free software: you can redistribute it and/or modify #
### it under the terms of the GNU General Public License as published by #
### the Free Software Foundation, version 3.                             #
###                                                                      #
### This program is distributed in the hope that it will be useful,      #
### but WITHOUT ANY WARRANTY; without even the implied warranty of       #
### MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the        #
### GNU General Public License for more details.                         #
###                                                                      #
### You should find a copy of the GNU General Public License at          #
###the Copyrights section or, see http://www.gnu.org/licenses.           #
##########################################################################


import numpy as np
import pylab as pl

beta=1.0
gamma=1/10.0
mu=5e-4
N0=5000.0
### You may want to try with popylation size of 50 (small) to see the events
### In this case uncomment the next line
#N0=50.0
ND=MaxTime=2*365.0
Y0=pl.ceil(mu*N0/gamma)
X0=pl.floor(gamma*N0/beta)
Z0=N0-X0-Y0

INPUT = np.array((X0,Y0,Z0))

timestep=0.0

def stoc_eqs(INP,ts): 
	V = INP
	Rate=np.zeros((6))
	Change=np.zeros((6,3))
	N=np.sum(V[range(3)])
	Rate[0] = beta*V[0]*V[1]/N; Change[0,:]=([-1, +1, 0]);
	Rate[1] = gamma*V[1];  Change[1,:]=([0, -1, +1]);
	Rate[2] = mu*N;  Change[2,:]=([+1, 0, 0]);
	Rate[3] = mu*V[0];  Change[3,:]=([-1, 0, 0]);
	Rate[4] = mu*V[1];  Change[4,:]=([0, -1, 0]);
	Rate[5] = mu*V[2];  Change[5,:]=([0, 0, -1]);
	R1=pl.rand();
	R2=pl.rand();
	ts = -np.log(R2)/(np.sum(Rate));
	m=min(pl.find(pl.cumsum(Rate)>=R1*pl.sum(Rate)));
	V[range(3)]=V[range(3)]+Change[m,:]
	return [V,ts]

def Stoch_Iteration(INPUT):
	lop=0
	ts=0
	T=[0]
	S=[0]
	I=[0]
	R=[0]
	while T[lop] < ND:
		lop=lop+1
		T.append(T[lop-1]+ts)
		S.append(INPUT[0])
		I.append(INPUT[1])
		R.append(INPUT[2])
		[res,ts] = stoc_eqs(INPUT,ts)
		lop=lop+1
		T.append(T[lop-1])
		S.append(INPUT[0])
		I.append(INPUT[1])
		R.append(INPUT[2])
	return [T,S,I,R]

[T,S,I,R]=Stoch_Iteration(INPUT)

tT=np.array(T)[1:,]/365.
tS=np.array(S)[1:,]
tI=np.array(I)[1:,]
tR=np.array(R)[1:,]

pl.subplot(311)
pl.plot(tT, tS, 'g')
#pl.xlabel ('Time (years)')
pl.ylabel ('Susceptible')
pl.subplot(312)
pl.plot(tT, tI, 'r')
#pl.xlabel ('Time (years)')
pl.ylabel ('Infectious')
pl.subplot(313)
pl.plot(tT, tR, 'k')
pl.xlabel ('Time (years)')
pl.ylabel ('Recovered')
pl.show()

--Ilias.soumpasis 09:46, 12 October 2008 (UTC)