Python Programs for Modelling Infectious Diseases book:Chapter 5:Program 5.4

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

Program 5.4 is a Rabbit Hemorrhagic Disease model (page 186 of the book). These are the equations and the code of the model:

Equations

 \frac{dX}{dt} = \alpha(t)*N-\beta(t)*X*Y- \Big( \mu+\frac{N}{K} \Big) *X

 \frac{dY}{dt} = \beta(t)*X*Y-\Big( \gamma+m+\mu+\frac{N}{K} \Big) *Y

 \frac{dN}{dt} = \Big( \alpha_{t}-\gamma-\mu-\frac{N}{K} \Big) *N-m*Y

 \beta(t) = \beta_{0}*(1+\beta_{1}* \sin(\omega*t))

 \alpha(t) = \alpha_{0}*(1+\alpha_{1}* \sin(\omega*t))

Code

Program 5.4: Rabbit Hemorrhagic Disease model
#!/usr/bin/env python

####################################################################
###    This is the PYTHON version of program 5.4 from page 186 of  #
### "Modeling Infectious Disease in humans and animals"            #
### by Keeling & Rohani.					   #
###								   #
### It is a model for Rabit Hemorrhagic Disease, in which both     #
### transmission rate and birth rates can be seasonally forced.    #
### Note that we are using numbers of individuals in this formulation.#
### Bifurcation plots are not possible with this code.             #
####################################################################

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

beta0=0.936;
beta1=0.1;
gamma=0.025;
alpha0=0.02;
alpha1=0.1;
mu=0.01;
m=0.475;
K=10000;

X0=0.5;
Y0=0.01;
N0=0.6;
Years = 60
ND=MaxTime=Years*365;
TS=1.0;

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

def diff_eqs(INP,t):  
	'''The main set of equations'''
	Y=np.zeros((3))
	V = INP   
	t=np.mod(t,365)
	alpha = alpha0*(1+alpha1*np.sin(2*np.pi*t/365))
	beta = beta0*(1+beta1*np.sin(2*np.pi*t/365))
	Y[0] = alpha * V[2]- beta*V[0]*V[1] - (mu+V[2]/K) * V[0] # dX/dt
	Y[1] = beta*V[0]*V[1] - (mu + m + gamma + V[2]/K) * V[1] # dY/dt
	Y[2] = (alpha - mu - V[2]/K) * V[2] - m * V[1]
	return Y   # For odeint


t_start = 0.0; t_end = ND; t_inc = TS
t_range = np.arange(t_start, t_end+t_inc, t_inc)
RES = spi.odeint(diff_eqs,INPUT,t_range)

print RES

t=(np.arange(ND)/365.)
#Ploting
pl.subplot(311)
pl.plot(t,RES[1:,0], 'g', label='S')
pl.xlabel ('Time (years)')
pl.ylabel ('Susceptibles')
pl.subplot(312)
pl.plot(t,RES[1:,1], 'r', label='I')
pl.xlabel ('Time (years)')
pl.ylabel ('Infectious')
pl.subplot(313)
pl.plot(t,1-(RES[1:,0]+RES[1:,1]), 'k', label='R')
pl.xlabel ('Time (years)')
pl.ylabel ('Recovereds')

pl.show()

--Ilias.soumpasis 07:49, 12 October 2008 (UTC)