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

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

Program 6.2 is a SIR model with Scaled additive noise (page 197 of the book). These are the equations and the code of the model:

Equations

 \frac{dX}{dt} = [\nu*N-\sqrt{\nu*N}*\xi_{1}]-[\beta*X*Y/N+\sqrt{\beta*X*Y/N}*\xi_{2}]-[\mu*X+\sqrt{\mu*X}*\xi_{3}]

 \frac{dY}{dt} = [[\beta*X*Y/N+\sqrt{\beta*X*Y/N}*\xi_{2}]-[\gamma*Y+\sqrt{\gamma*Y}*\xi_{4}]-[\mu*Y+\sqrt{\mu*Y}*\xi_{5}]

 \frac{dZ}{dt} = [\gamma*Y+\sqrt{\gamma*Y}*\xi_{4}]-[\mu*Z+\sqrt{\mu*Z}*\xi_{6}]

Code

Program 6.2: SIR model with Scaled additive noise (one run)
#!/usr/bin/env python

####################################################################
###    This is the PYTHON version of program 6.2 from page 197 of  #
### "Modeling Infectious Disease in humans and animals"            #
### by Keeling & Rohani.					   #
###								   #
### % It is the SIR epidemic model with constant additive noise    #
### added all the various rates.				   #
### Given the difficulties in integrating the dynamics, the user   #
### is prompted for a integration time-step.			   #
####################################################################

##########################################################################
### 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 scipy.integrate as spi
import numpy as np
import pylab as pl


beta=1.0;
gamma=1/10.0;
mu=1/(50*365.0);
X0=1e5;
Y0=500;
N0=1e6;
Step=1;
ND=MaxTime=5*365.0;
TS=1.0

INPUT0=np.hstack((X0,Y0))

def diff_eqs(INP,t):  
	'''The main set of equations'''
	Y=np.zeros((2))
	V = INP     
	Y[0] = (mu*N0 + np.sqrt(mu*N0)*P[0]) - (beta*V[0]*V[1]/N0 \
	+ np.sqrt(beta*V[0]*V[1]/N0)*P[1])	- (mu*V[1] + np.sqrt(mu*V[1])*P[2]);
	Y[1] = (beta*V[0]*V[1]/N0 + np.sqrt(beta*V[0]*V[1]/N0)*P[1]) - \
	(gamma*V[1] + np.sqrt(gamma*V[1])*P[3]) - (mu*V[1] + np.sqrt(mu*V[1])*P[4]);
	return Y   # For odeint

T=np.zeros((np.ceil(ND/Step),1))
RES=np.zeros((np.ceil(ND/Step),2))
INPUT=INPUT0
t=0
loop=0
sqrtStep=np.sqrt(Step)

while t<ND and INPUT[0]>0 and INPUT[1]>0:
	t_start = 0.0; t_end = t_start+Step; t_inc = TS
	t_range = np.arange(t_start, t_end+t_inc, t_inc)
	P=np.random.normal(size=5)/sqrtStep
	PRES = spi.odeint(diff_eqs,INPUT,t_range)
	T[loop]=t=t+Step
	INPUT=PRES[-1]
	RES[loop]=PRES[-1]
	loop += 1

print(RES)

### plotting
pl.subplot(211)
pl.plot(T/365., RES[:,0], '.-g')
pl.xlabel('Time (Years)')
pl.ylabel('Susceptibles')
pl.subplot(212)
pl.plot(T/365., RES[:,1], '.-r')
pl.ylabel('Infected')
pl.xlabel('Time (Years)')

pl.show()

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