-
Notifications
You must be signed in to change notification settings - Fork 5
/
sir_func.R
135 lines (125 loc) · 5.58 KB
/
sir_func.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
##################################################################################
# An R script to solve ODE's of an SIR model
# http://www.sherrytowers.com/sir_func.R
#
# Author: Sherry Towers
#
# Created: Dec 1st, 2012
# Copyright Sherry Towers, 2012
#
# This script is not guaranteed to be free of bugs and/or errors.
#
# This script can be freely used and shared as long as the author and
# copyright information in this header remain intact.
#
# In order to use this script in R, you must first have the deSolve library
# installed. To do this, type in the R command line:
# install.packages("deSolve")
# then pick a mirror site located close to you.
##################################################################################
require("deSolve")
##################################################################################
##################################################################################
# This is a function which, given a value of S, I, and R at time t
# calculates the time derivatives of S I and R
#
# The function gets called over and over again by functions in the R deSolve
# library to do the numerical integration over many time steps to solve
# the system of ODE's
#
# t is the current time at a particular step
# The vector x contains the current values in each compartment
# The list object vparameters contains the parameters of the model, like the
# recovery period, gamma, and the transmission rate, beta (in this case)
# In the main program, the vparameters list object is filled with the parameter
# values and their names.
#
# This function gets passed to the functions in the deSolve package
#
##################################################################################
derivative_calc_func=function(t, x, vparameters){
###############################################################################
# The vector x is the same length as the number of compartments. In your main
# program you identify which compartment corresponds to which element of x.
# You need to make sure that these are in the same order.
###############################################################################
S = x[1]
I = x[2]
R = x[3]
###############################################################################
# vparameters is a list object, filled in the main program, and passed
# Keep this next line the same when you are writing your own function to
# solve a system of ODE's
###############################################################################
with(as.list(vparameters),{
############################################################################
# calculate the population size, which for a simple SIR model is
# npop = S+I+R
# we will need this to calculate our SIR model derivatives below
############################################################################
npop = S+I+R
############################################################################
# Now give the expressions for the derivatives of S, I, and R wrt time
# these come from the model equations. The following equations are for
# an SIR model. When you write your own function, replace these with
# your model equations
############################################################################
dS = -beta*S*I/npop
dI = +beta*S*I/npop - gamma*I
dR = +gamma*I
############################################################################
# vout is an output vector that contains the values of the derivates
# the compartments in the vector on the RHS need to be in the same order as the
# compartments used to fill the x vector!
############################################################################
vout = c(dS,dI,dR)
list(vout)
})
}
##################################################################################
##################################################################################
##################################################################################
# this is the same as the above function, except now it includes births and
# deaths (both with rate mu) in the model equations
##################################################################################
derivative_calc_func_with_demographics=function(t, x, vparameters){
S = x[1]
I = x[2]
R = x[3]
with(as.list(vparameters),{
npop = S+I+R
dS = -beta*S*I/npop - mu*S + npop*mu
dI = +beta*S*I/npop - gamma*I - mu*I
dR = +gamma*I - mu*R
out = c(dS,dI,dR)
list(out)
})
}
##################################################################################
##################################################################################
# this is the same as the derivative_calc_function, except now this involves
# calculating the derivatives of an SIR model with vaccination
# Rvac is the vaccinated (and now immune) compartment
# The vaccination begins at time_vaccination_begins, and ends at
# time_vaccination_ends
##################################################################################
derivative_calc_func_with_vaccination=function(t, x, vparameters){
S = x[1]
I = x[2]
R = x[3]
Rvac = x[4]
with(as.list(vparameters),{
npop = S+I+R+Rvac
dS = -beta*S*I/npop
dI = +beta*S*I/npop - gamma*I
dR = +gamma*I
dRvac = 0
if (t>=time_vaccination_begins&t<=time_vaccination_ends){
dS = dS - rho*S
dRvac = +rho*S
}
out = c(dS,dI,dR,dRvac)
list(out)
})
}