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

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

Program 6.5 is a SIR model with tau leap method (page 204 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}  \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {P1=Poisson}(\tau*\beta*X*Y/N)\ \ \mbox{Event: } X \mapsto X-P1, \ \ Y \mapsto Y+P1

\mbox{Recovery}  \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \  {P2=Poisson}(\tau*\gamma*Y)\ \ \ \ \ \ \ \ \ \ \ \ \ \ \mbox{Event: } Y \mapsto Y-P2, \ \ Z \mapsto Z+P2

\mbox{Birth}\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \mbox{P3=Poisson}(\tau*\mu*N) \ \ \ \ \ \ \ \ \ \mbox{Event: } X \mapsto X+P3

\mbox{Death Susceptible} \ \mbox{P4=Poisson}(\tau*\mu*X)\ \ \ \ \ \ \ \ \ \mbox{Event: } X \mapsto X-P4

\mbox{Death Infected} \ \ \ \ \ \mbox{P5=Poisson}(\tau*\mu*Y)\ \ \ \ \ \ \ \ \ \ \mbox{Event: } Y \mapsto Y-P5

\mbox{Death Recovered} \ \ \mbox{P6=Poisson}(\tau*\mu*Z)\ \ \ \ \ \ \ \ \ \ \mbox{Event: } Z \mapsto Z-P6

Code

Program 6.5: SIR model with tau leap method (one run)
#!/usr/bin/env python

####################################################################
###    This is the PYTHON version of program 6.5 from page 204 of  #
### "Modeling Infectious Disease in humans and animals"            #
### by Keeling & Rohani.					   #
###								   #
### It is the SIR model (including births and deaths) with         #
### (event-driven) demographic stochasticity approximated using    #
### the tau-leap method and assuming Poisson distributions.	   #
###								   #
### 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
tau=1.0
INPUT = np.array((X0,Y0,Z0))

def stoc_eqs(INP): 
	V = INP
	Rate=np.zeros((6))
	Change=np.zeros((6,3))
	N=V[0]+V[1]+V[2]
	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]);
	for i in range(6):
		Num=np.random.poisson(Rate[i]*tau);
		## Make sure things don't go negative
		Use=min([Num, V[pl.find(Change[i,:]<0)]]);
		V=V+Change[i,:]*Use;
	return V

def Stoch_Iteration(INPUT):
	lop=0
	S=[0]
	I=[0]
	R=[0]
	for lop in T:
		res = stoc_eqs(INPUT)
		S.append(INPUT[0])
		I.append(INPUT[1])
		R.append(INPUT[2])
		INPUT=res
	return [S,I,R]

T=np.arange(0.0, ND, tau)
[S,I,R]=Stoch_Iteration(INPUT)

tT=np.array(T)/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:57, 12 October 2008 (UTC)