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