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

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

Program 6.3 is a SIS model with demographic stochasticity (page 202 of the book). These are the equations and the code of the model:

Equations

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

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

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

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

Code

Program 6.3: SIS model with demographic stochasticity (one run)
#!/usr/bin/env python

####################################################################
###    This is the PYTHON version of program 6.3 from page 202 of  #
### "Modeling Infectious Disease in humans and animals"            #
### by Keeling & Rohani.					   #
###								   #
### This is a relatively simple stochastic model as only 2 events  #
### are possible: infection or recovery.			   #
### 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=0.03
gamma=1/100.0
Y0=70.0
N0=100.0
ND=MaxTime=10*365.0;


INPUT = Y0

timestep=0.0

def stoc_eqs(INP,ts):  
	Z=INP
	Rate1 = beta*(N0-Z)*Z/N0
	Rate2 = gamma*Z
	R1=pl.rand()
	R2=pl.rand()
	ts = -np.log(R2)/(Rate1+Rate2)
	if R1<(Rate1/(Rate1+Rate2)):
		Z += 1;  # do infection
	else:
		Z -= 1;  # do recovery
	return [Z,ts]

def Stoch_Iteration(INPUT):
	lop=0
	ts=0
	T=[0]
	RES=[0]
	while T[lop] < ND and INPUT > 0:
		[res,ts] = stoc_eqs(INPUT,ts)
		lop=lop+1
		T.append(T[lop-1]+ts)
		RES.append(INPUT)
		lop=lop+1
		T.append(T[lop-1])
		RES.append(res)
		INPUT=res
	return [RES, T]

[RES,t]=Stoch_Iteration(INPUT)

t=np.array(t)
RES=np.array(RES)
### plotting
pl.plot(t/365., RES, 'r')
pl.show()

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