242 lines
11 KiB
Text
242 lines
11 KiB
Text
######################################################################
|
|
### SEIR Model Functions for OAW Measles Model ###
|
|
### Code by Paul Gastanaduy, Prabasaj Paul, Nina Masters ###
|
|
### Updated September 2 2022 ###
|
|
######################################################################
|
|
|
|
# Define parameters & states we will be using throughout these functions -------------------------------------------------------------
|
|
# Si : starting value for Susceptible type i hosts: numeric
|
|
# Svi : starting value for Susceptible type i (vaccine failures) - those who were vaccinated but didn't take
|
|
# Ei : starting value for Exposed type i hosts: numeric
|
|
# Ii : starting value for Infected type i hosts: numeric
|
|
# Ri : starting value for Recovered type i hosts: numeric
|
|
|
|
# bij : rate of transmission to susceptible type i host from infected type j host : numeric
|
|
# theta - same for all hosts: 1/duration of incubation period : numeric
|
|
# gamma - same for all hosts: 1/duration of infectiousness : numeric
|
|
|
|
# 1 Set up parameter definitions ------------------------------------------------------------------------------------------------------
|
|
#https://www.cdc.gov/vaccines/pubs/pinkbook/downloads/meas.pdf
|
|
|
|
R0 = 14 #Ro for Measles of 14, backtrack out the Ro based on WT, got estimates of Ro from 11.55 - 15.20
|
|
infectious.period = 5 # infectious period - mean = 5 days
|
|
preinfectious.period = 8 #preinfectious (incubation) period = 8 days
|
|
|
|
# betas = Ro/infectious period (set it up this way to give max flexibility to change to non homogeneous mixing)
|
|
b11=b12=b13=b14=b15 = R0/infectious.period
|
|
b21=b22=b23=b24=b25 = R0/infectious.period
|
|
b31=b32=b33=b34=b35 = R0/infectious.period
|
|
b41=b42=b43=b44=b45 = R0/infectious.period
|
|
b51=b52=b53=b54=b55 = R0/infectious.period
|
|
|
|
theta = 1/preinfectious.period
|
|
gamma = 1/infectious.period
|
|
|
|
# 2 Define SEIR model transitions with NO vaccination --------------------------------------------------------------------------------
|
|
SEIR_transitions <- list(
|
|
c(S1 = -1, E1 = 1), # infection of susceptible in group 1
|
|
c(S1 = -1, E1 = 1), # infection of susceptible in group 1
|
|
c(S1 = -1, E1 = 1), # infection of susceptible in group 1
|
|
c(S1 = -1, E1 = 1), # infection of susceptible in group 1
|
|
c(S1 = -1, E1 = 1), # infection of susceptible in group 1
|
|
|
|
c(S2 = -1, E2 = 1), # infection of susceptible in group 2
|
|
c(S2 = -1, E2 = 1), # infection of susceptible in group 2
|
|
c(S2 = -1, E2 = 1), # infection of susceptible in group 2
|
|
c(S2 = -1, E2 = 1), # infection of susceptible in group 2
|
|
c(S2 = -1, E2 = 1), # infection of susceptible in group 2
|
|
|
|
c(S3 = -1, E3 = 1), # infection of susceptible in group 3
|
|
c(S3 = -1, E3 = 1), # infection of susceptible in group 3
|
|
c(S3 = -1, E3 = 1), # infection of susceptible in group 3
|
|
c(S3 = -1, E3 = 1), # infection of susceptible in group 3
|
|
c(S3 = -1, E3 = 1), # infection of susceptible in group 3
|
|
|
|
c(S4 = -1, E4 = 1), # infection of susceptible in group 4
|
|
c(S4 = -1, E4 = 1), # infection of susceptible in group 4
|
|
c(S4 = -1, E4 = 1), # infection of susceptible in group 4
|
|
c(S4 = -1, E4 = 1), # infection of susceptible in group 4
|
|
c(S4 = -1, E4 = 1), # infection of susceptible in group 4
|
|
|
|
c(S5 = -1, E5 = 1), # infection of susceptible in group 5
|
|
c(S5 = -1, E5 = 1), # infection of susceptible in group 5
|
|
c(S5 = -1, E5 = 1), # infection of susceptible in group 5
|
|
c(S5 = -1, E5 = 1), # infection of susceptible in group 5
|
|
c(S5 = -1, E5 = 1), # infection of susceptible in group 5
|
|
|
|
c(E1 = -1, I1 = 1, Cases1 = 1), #from latency to active infection (also tallying cases) in group 1
|
|
c(E2 = -1, I2 = 1, Cases2 = 1), #from latency to active infection (also tallying cases) in group 2
|
|
c(E3 = -1, I3 = 1, Cases3 = 1), #from latency to active infection (also tallying cases) in group 3
|
|
c(E4 = -1, I4 = 1, Cases4 = 1), #from latency to active infection (also tallying cases) in group 4
|
|
c(E5 = -1, I5 = 1, Cases5 = 1), #from latency to active infection (also tallying cases) in group 5
|
|
|
|
c(I1 = -1, R1 = 1), # recovery from infection in group 1
|
|
c(I2 = -1, R2 = 1), # recovery from infection in group 2
|
|
c(I3 = -1, R3 = 1), # recovery from infection in group 3
|
|
c(I4 = -1, R4 = 1), # recovery from infection in group 4
|
|
c(I5 = -1, R5 = 1) # recovery from infection in group 5
|
|
)
|
|
|
|
# 3 Define SEIR Model Function --------------------------------------------------------------------------------
|
|
|
|
SEIR_Model <- function(x, parameters, t) {
|
|
|
|
S1 <- x["S1"]
|
|
S2 <- x["S2"]
|
|
S3 <- x["S3"]
|
|
S4 <- x["S4"]
|
|
S5 <- x["S5"]
|
|
|
|
E1 <- x["E1"]
|
|
E2 <- x["E2"]
|
|
E3 <- x["E3"]
|
|
E4 <- x["E4"]
|
|
E5 <- x["E5"]
|
|
|
|
I1 <- x["I1"]
|
|
I2 <- x["I2"]
|
|
I3 <- x["I3"]
|
|
I4 <- x["I4"]
|
|
I5 <- x["I5"]
|
|
|
|
R1 <- x["R1"]
|
|
R2 <- x["R2"]
|
|
R3 <- x["R3"]
|
|
R4 <- x["R4"]
|
|
R5 <- x["R5"]
|
|
|
|
Cases1 <- x["Cases1"]
|
|
Cases2 <- x["Cases2"]
|
|
Cases3 <- x["Cases3"]
|
|
Cases4 <- x["Cases4"]
|
|
Cases5 <- x["Cases5"]
|
|
|
|
N <- S1 + S2 + S3 + S4 + S5 + E1 + E2 + E3 + E4 + E5 + I1 + I2 + I3 + I4 + I5 + R1 + R2 + R3 + R4 +R5
|
|
|
|
#specify each rate/transition/reaction that can happen in the system
|
|
rates = c(b11*S1*I1/N, b12*S1*I2/N, b13*S1*I3/N, b14*S1*I4/N, b15*S1*I5/N, #infection in gp 1
|
|
b21*S2*I1/N, b22*S2*I2/N, b23*S2*I3/N, b24*S2*I4/N, b25*S2*I5/N, #infection in gp 2
|
|
b31*S3*I1/N, b32*S3*I2/N, b33*S3*I3/N, b34*S3*I4/N, b35*S3*I5/N, #infection in gp 3
|
|
b41*S4*I1/N, b42*S4*I2/N, b43*S4*I3/N, b44*S4*I4/N, b45*S4*I5/N, #infection in gp 4
|
|
b51*S5*I1/N, b52*S5*I2/N, b53*S5*I3/N, b54*S5*I4/N, b55*S5*I5/N, #infection in gp 5
|
|
theta*E1, theta*E2, theta*E3, theta*E4, theta*E5, #from pre-infectious to infectious (among exposed)
|
|
gamma*I1, gamma*I2, gamma*I3, gamma*I4, gamma*I5) #recovery in the infected in each group
|
|
return(rates)
|
|
}
|
|
|
|
|
|
# 4 Define SEIRv model transitions: from S --> E, Sv --> E, E--> I, I--> R from each model state ---------------------------------------------------
|
|
SEIRv_transitions <- list(
|
|
c(S1 = -1, E1 = 1), # infection of susceptible in group 1
|
|
c(S1 = -1, E1 = 1), # infection of susceptible in group 1
|
|
c(S1 = -1, E1 = 1), # infection of susceptible in group 1
|
|
c(S1 = -1, E1 = 1), # infection of susceptible in group 1
|
|
c(S1 = -1, E1 = 1), # infection of susceptible in group 1
|
|
|
|
c(S2 = -1, E2 = 1), # infection of susceptible in group 2
|
|
c(S2 = -1, E2 = 1), # infection of susceptible in group 2
|
|
c(S2 = -1, E2 = 1), # infection of susceptible in group 2
|
|
c(S2 = -1, E2 = 1), # infection of susceptible in group 2
|
|
c(S2 = -1, E2 = 1), # infection of susceptible in group 2
|
|
|
|
c(Sv2 = -1, E2 = 1), # infection of susceptible, vaccine failures in group 2
|
|
c(Sv2 = -1, E2 = 1), # infection of susceptible, vaccine failures in group 2
|
|
c(Sv2 = -1, E2 = 1), # infection of susceptible, vaccine failures in group 2
|
|
c(Sv2 = -1, E2 = 1), # infection of susceptible, vaccine failures in group 2
|
|
c(Sv2 = -1, E2 = 1), # infection of susceptible, vaccine failures in group 2
|
|
|
|
c(S3 = -1, E3 = 1), # infection of susceptible in group 3
|
|
c(S3 = -1, E3 = 1), # infection of susceptible in group 3
|
|
c(S3 = -1, E3 = 1), # infection of susceptible in group 3
|
|
c(S3 = -1, E3 = 1), # infection of susceptible in group 3
|
|
c(S3 = -1, E3 = 1), # infection of susceptible in group 3
|
|
|
|
c(Sv3 = -1, E3 = 1), # infection of susceptible, vaccine failures in group 3
|
|
c(Sv3 = -1, E3 = 1), # infection of susceptible, vaccine failures in group 3
|
|
c(Sv3 = -1, E3 = 1), # infection of susceptible, vaccine failures in group 3
|
|
c(Sv3 = -1, E3 = 1), # infection of susceptible, vaccine failures in group 3
|
|
c(Sv3 = -1, E3 = 1), # infection of susceptible, vaccine failures in group 3
|
|
|
|
c(S4 = -1, E4 = 1), # infection of susceptible in group 4
|
|
c(S4 = -1, E4 = 1), # infection of susceptible in group 4
|
|
c(S4 = -1, E4 = 1), # infection of susceptible in group 4
|
|
c(S4 = -1, E4 = 1), # infection of susceptible in group 4
|
|
c(S4 = -1, E4 = 1), # infection of susceptible in group 4
|
|
|
|
c(Sv4 = -1, E4 = 1), # infection of susceptible, vaccine failures in group 4
|
|
c(Sv4 = -1, E4 = 1), # infection of susceptible, vaccine failures in group 4
|
|
c(Sv4 = -1, E4 = 1), # infection of susceptible, vaccine failures in group 4
|
|
c(Sv4 = -1, E4 = 1), # infection of susceptible, vaccine failures in group 4
|
|
c(Sv4 = -1, E4 = 1), # infection of susceptible, vaccine failures in group 4
|
|
|
|
c(S5 = -1, E5 = 1), # infection of susceptible in group 5
|
|
c(S5 = -1, E5 = 1), # infection of susceptible in group 5
|
|
c(S5 = -1, E5 = 1), # infection of susceptible in group 5
|
|
c(S5 = -1, E5 = 1), # infection of susceptible in group 5
|
|
c(S5 = -1, E5 = 1), # infection of susceptible in group 5
|
|
|
|
c(E1 = -1, I1 = 1, Cases1 = 1), #from latency to active infection in group 1
|
|
c(E2 = -1, I2 = 1, Cases2 = 1), #from latency to active infection in group 2
|
|
c(E3 = -1, I3 = 1, Cases3 = 1), #from latency to active infection in group 3
|
|
c(E4 = -1, I4 = 1, Cases4 = 1), #from latency to active infection in group 4
|
|
c(E5 = -1, I5 = 1, Cases5 = 1), #from latency to active infection in group 5
|
|
|
|
c(I1 = -1, R1 = 1), # recovery from infection in group 1
|
|
c(I2 = -1, R2 = 1), # recovery from infection in group 2
|
|
c(I3 = -1, R3 = 1), # recovery from infection in group 3
|
|
c(I4 = -1, R4 = 1), # recovery from infection in group 4
|
|
c(I5 = -1, R5 = 1) # recovery from infection in group 5
|
|
)
|
|
|
|
# 5 Define SEIRv model rate function ------------------------------------------------------------------------------------------------------------
|
|
SEIRv_Model <- function(x, parameters, t) {
|
|
|
|
S1 <- x["S1"]
|
|
S2 <- x["S2"]
|
|
Sv2 <- x["Sv2"]
|
|
S3 <- x["S3"]
|
|
Sv3 <- x["Sv3"]
|
|
S4 <- x["S4"]
|
|
Sv4 <- x["Sv4"]
|
|
S5 <- x["S5"]
|
|
|
|
E1 <- x["E1"]
|
|
E2 <- x["E2"]
|
|
E3 <- x["E3"]
|
|
E4 <- x["E4"]
|
|
E5 <- x["E5"]
|
|
|
|
I1 <- x["I1"]
|
|
I2 <- x["I2"]
|
|
I3 <- x["I3"]
|
|
I4 <- x["I4"]
|
|
I5 <- x["I5"]
|
|
|
|
R1 <- x["R1"]
|
|
R2 <- x["R2"]
|
|
R3 <- x["R3"]
|
|
R4 <- x["R4"]
|
|
R5 <- x["R5"]
|
|
|
|
Cases1 <- x["Cases1"]
|
|
Cases2 <- x["Cases2"]
|
|
Cases3 <- x["Cases3"]
|
|
Cases4 <- x["Cases4"]
|
|
Cases5 <- x["Cases5"]
|
|
|
|
N <- S1 + S2 + Sv2 + S3 + Sv3 + S4 + Sv4 + S5 + E1 + E2 + E3 + E4 + E5 + I1 + I2 + I3 + I4 + I5 + R1 + R2 + R3 + R4 +R5
|
|
|
|
#specify each rate/transition/reaction that can happen in the system
|
|
rates = c(b11*S1*I1/N, b12*S1*I2/N, b13*S1*I3/N, b14*S1*I4/N, b15*S1*I5/N, #infection in gp 1
|
|
b21*S2*I1/N, b22*S2*I2/N, b23*S2*I3/N, b24*S2*I4/N, b25*S2*I5/N, #infection in gp 2
|
|
b21*Sv2*I1/N, b22*Sv2*I2/N, b23*Sv2*I3/N, b24*Sv2*I4/N, b25*Sv2*I5/N, #infection in gp 2 vaccine failures
|
|
b31*S3*I1/N, b32*S3*I2/N, b33*S3*I3/N, b34*S3*I4/N, b35*S3*I5/N, #infection in gp 3
|
|
b31*Sv3*I1/N, b32*Sv3*I2/N, b33*Sv3*I3/N, b34*Sv3*I4/N, b35*Sv3*I5/N, #infection in gp 3 vaccine failures
|
|
b41*S4*I1/N, b42*S4*I2/N, b43*S4*I3/N, b44*S4*I4/N, b45*S4*I5/N, #infection in gp 4
|
|
b41*Sv4*I1/N, b42*Sv4*I2/N, b43*Sv4*I3/N, b44*Sv4*I4/N, b45*Sv4*I5/N, #infection in gp 4 vaccine failures
|
|
b51*S5*I1/N, b52*S5*I2/N, b53*S5*I3/N, b54*S5*I4/N, b55*S5*I5/N, #infection in gp 5
|
|
theta*E1, theta*E2, theta*E3, theta*E4, theta*E5, #from pre-infectious to infectious (among exposed)
|
|
gamma*I1, gamma*I2, gamma*I3, gamma*I4, gamma*I5) #recovery in the infected in each group
|
|
return(rates)
|
|
}
|
|
|