rm(list = ls())
require(deSolve)
## Loading required package: deSolve
require(lhs)
## Loading required package: lhs

Required functions

SARSCOV2Model<- function (t, y, params) {
  #System of ODES
  
S.h <- y[1] #create local variable S, first element of y
E.h <- y[2] 
I.h <- y[3]
S.l <- y[4]
E.l <- y[5]
I.l <- y[6]
Q <- y[7]
R <- y[8]
V <- y[9]

with(
  as.list(params),
  {
dS.h<--q*beta*(I.h+I.l)*S.h/(S.h+E.h+I.h+S.l+E.l+I.l+Q+R)-c*(1-exp((-1/K)*V))*S.h
dE.h<-q*beta*(I.h+I.l)*S.h/(S.h+E.h+I.h+S.l+E.l+I.l+Q+R)+c*(1-exp((-1/K)*V))*S.h-lambda*E.h
dI.h<-lambda*E.h-b*g*I.h-aH*h*(1-g)*I.h-gammah*(1-h)*(1-g)*I.h
dS.l<--(1-p)*q*beta*(I.h+I.l)*S.l/(S.h+E.h+I.h+S.l+E.l+I.l+Q+R)-(1-p)*c*(1-exp((-1/K)*V))*S.l
dE.l<-(1-p)*q*beta*(I.h+I.l)/(S.h+E.h+I.h+S.l+E.l+I.l+Q+R)+(1-p)*c*(1-exp((-1/K)*V))*S.l-lambda*E.l
dI.l<-lambda*E.l-b*I.l
dQ<-b*I.l+g*b*I.h-aQ*h*Q-gammaQ*(1-h)*Q
dR<-gammah*(1-h)*(1-g)*I.h+gammaQ*(1-h)*Q
dV<-omega*I.h+(1-p)*omega*I.l-delta*V
  dy<-c(dS.h,dE.h,dI.h,dS.l,dE.l,dI.l,dQ,dR,dV) #combine results into one vector dy
list(dy)
  }
)
}
plot.SARSCOV2Model<-function(out)
  {
  par(mfrow=c(3,2))
  plot(out$time, out$S.h, xlab='Time', ylab='Susceptible(High Risk)',type='l',col='steelblue4', lwd=3)
  plot(out$time, out$E.h, xlab='Time', ylab='Exposed(High Risk)', type='l', col='steelblue4', lwd=3)
  plot(out$time, out$I.h, xlab='Time', ylab='Infectious(High Risk)', type='l', col='steelblue4', lwd=3)
  plot(out$time, out$S.l, xlab='Time', ylab='Susceptible(Low Risk)', type='l', col='steelblue4', lwd=3)
  plot(out$time, out$E.l, xlab='Time', ylab='Exposed(Low Risk)', type='l', col='steelblue4', lwd=3)
  plot(out$time, out$I.l, xlab='Time', ylab='Infectious(Low Risk)', type='l', col='steelblue4', lwd=3)
  plot(out$time, out$Q, xlab='Time', ylab='Quarantined', type='l', col='steelblue4', lwd=3)
  plot(out$time, out$R, xlab='Time', ylab='Recovered', type='l', col='steelblue4', lwd=3)
  plot(out$time, out$R, xlab='Time', ylab='V', type='l', col='steelblue4', lwd=3)
}
get.size <- function(out) {
  100000-(tail(out$S.h,1)+tail(out$S.l,1)) #subtract number of susceptibles in high risk and low risk groups from total population size 100000
}

Function to compute maximum of \(I_h\) (epidemic peak):

get.peak <- function (out)
{
  max(out$I.h)
}

R0 function for sensitivity analysis:

SIRR0list <- function(q, g, p, delta, b, K, beta, c, omega, aH, h, gammah, gammaQ, lambda, S0=100000){
  
  with(
  as.list(params),
  {
    R0 <- (b*K*(g*(p-1)^2 + p)*q*beta*delta + b*c*(-g*(p-1)^3 + p)*S0*omega + (g-1)*(p-1)^2*(aH*h + gammah - h*gammah)*(-K*q*beta*delta + c*(p-1)*S0*omega)) / (b*K*(b*g - (g-1)*aH*h + gammah - h*gammah)*delta)
  }
  )
  
}

A typical simulation:

times <- seq(0,1000,by=1) 
params <- c(q=0.2,beta=13,p=.2,c=0,lambda=1/4.43,b=1/0.77,g=0,gammaQ=0.1,gammah=1/2.7,aQ=1/1.93,aH=1/2.7,h=0.00082,omega=0,delta=1,K=10000)
ystart <- c(S.h=(.2)*100000,E.h=0,I.h=1,S.l=(1-0.2)*100000,E.l=0,I.l=0,Q=0,R=0,V=0)
out <- as.data.frame(lsoda(ystart,times,SARSCOV2Model,params))
plot.SARSCOV2Model(out)

Sensitivity analysis. Objective: To produce graphs of \(R_0\), final size of the epidemic and epidemic peak.

h1<-1000
set.seed(6242015) #this random seed can be changed
lhs<-maximinLHS(h1,14) #here 14 is number of parameters that we will generate

Min and max of all parameters except \(q\). Range of values for each parameter obtained from our parameters table on Google Drive; based on realistic ranges and confidence intervals where possible. Please check the values I have chosen below and feel free to modify.

g.min<-0
g.max<-1
q.min<-0
q.max<-0.2
p.min<-0.1
p.max<-0.2
delta.min<-1/7
delta.max<-24
b.max <- 1/1.29 # from Peak et al. 2020
b.min <- 0
K.max <- 1000000
K.min <- 10000
beta.max <- 18
beta.min <- 10
# c.max <- 5
# c.min <- 0.01
#S0.max <- 100000 #no need for S0. S0 = 100000
#S0.min <- 1
omega.max <- 20
omega.min <- 1
aH.max <- 1/2
aH.min <- 1/10
aQ.max <- 1/1
aQ.min <- 1/9
h.max <- 200/100000 #default: 82/100000 are hospitalized
h.min <- 50/100000
gammah.max <- 1/2.5
gammah.min <- 1/10
gammaQ.max <- 1/9
gammaQ.min <- 1/15
lambda.max <- 1/4
lambda.min <- 1/6

Create the 14 dimensional parameter space:

params.set<-cbind(
  g=lhs[,1]*(g.max-g.min)+g.min,
  q=lhs[,2]*(q.max-q.min)+q.min,
  p=lhs[,3]*(p.max-p.min)+p.min,
  delta=lhs[,4]*(delta.max-delta.min)+delta.min,
  b = lhs[,5]*(b.max-b.min)+b.min,
  K = lhs[,6]*(K.max-K.min)+K.min,
  beta = lhs[,7]*(beta.max-beta.min)+beta.min,
  # c = lhs[,7]*(c.max-c.min)+c.min,
  omega = lhs[,8]*(omega.max-omega.min)+omega.min,
  aH = lhs[,9]*(aH.max-aH.min)+aH.min,
  aQ = lhs[,10]*(aQ.max-aQ.min)+aQ.min,
  h = lhs[,11]*(h.max-h.min)+h.min,
  gammah = lhs[,12]*(gammah.max-gammah.min)+gammah.min,
  gammaQ = lhs[,13]*(gammaQ.max-gammaQ.min)+gammaQ.min,
  lambda = lhs[,14]*(lambda.max-lambda.min)+lambda.min
)
#q=lhs[,2]*(q.max-q.min)+q.min,
#S0 = lhs[,9]*(S0.max-S0.min)+S0.min,
h2<-100 #number of simulated parameter sets. We would like to increase this to at least 100 for our final plots
times <- seq(0,100,by=0.1)
cseq<-seq(0.01,5.01, by=0.25)
levels<-length(cseq)

For loop to compute \(R_0\), epidemic final size, epidemic peak as a function of \(q\):

j <- 1
 data <- data.frame(matrix(rep(NA,levels*h2*19),nrow=levels*h2))
 for(i in 1:h2){
   for (c in cseq){
data[j,1:16] <- params <- as.list(c(params.set[i,], c=c, S0=100000)) #16 parameters
out <- (as.data.frame(ode(ystart,times, SARSCOV2Model, params)))
data[j,17] <- get.size(out)
data[j,18] <- SIRR0list(params)
data[j,19] <- get.peak(out)
j <- j+1
}
}
names(data) <- c(names(params),'outbreak.size', 'R0','peak.size')
#save(data, file='data.Rdata') remove the comment symbol when you want to save a large matrix (e.g. 100 simulations)
#load('data.Rdata')remove the comment symbol when you want to load a large matrix (e.g. 100 simulations)
plot(data$c, data$R0, pch=19, cex=0.3, col='blue', xlab='Per-capita contact with environment(c)', ylab='R0')

plot(data$c, data$outbreak.size/100000, pch=19, cex=0.3, col='blue', xlab='Per-capita contact with environment (c)', ylab='Proportion of Epidemic Size') #divide epidemic final size by 100000 to get proportion

plot(data$c, data$peak.size/100000, pch=19, cex=0.3, col='blue', xlab='Per-capita contact with environment (c)', ylab='Proportion of Infected (High)')

Boxplots of the above data are what we’d like in our report and poster:

boxplot(data$peak.size/100000~data$c,  border='blue', log='y', pch='.', xlab='Per-capita contact with environment (c)', ylab='Proportion of Infected (High)')

boxplot(data$outbreak.size/100000~data$c,  border='blue', log='y', pch='.', xlab='Per-capita contact with environment (c)', ylab='Proportion of Epidemic Size')

boxplot(data$R0~data$c, border='blue', log='y', pch='.', xlab='Per-capita contact with environment (c)', ylab='R0')

Partial rank correlation coefficients - how sensitive is the output of interest to each parameter?

library(sensitivity)
bonferroni.alpha <- 0.05/16
prcc <- pcc(data[,1:16], data[,18], nboot = 1000, rank=TRUE, conf=1-bonferroni.alpha) #how sensitive is R0 to each parameter?
save(prcc, file='prcc.Rdata')

The greater the magnitude of PRCC, the more influence that parameter has on \(R_0\). Above the line: positively correlated Below the line: negatively correlated

library(sensitivity)
load('prcc.Rdata')
summary <- print(prcc)
## 
## Call:
## pcc(X = data[, 1:16], y = data[, 18], rank = TRUE, nboot = 1000,     conf = 1 - bonferroni.alpha)
## 
## Partial Rank Correlation Coefficients (PRCC):
##            original          bias  std. error    min. c.i.   max. c.i.
## g      -0.587146313  2.568046e-04 0.015518140 -0.639478941 -0.54012166
## q       0.825142000 -2.884810e-04 0.007704553  0.804349933  0.85078024
## p       0.056450017 -2.095763e-04 0.020619000 -0.006783609  0.12474125
## delta  -0.478862233 -2.624720e-04 0.017910763 -0.532967935 -0.42114291
## b      -0.858152372  4.937260e-05 0.006517995 -0.877117343 -0.83925026
## K      -0.536048943 -2.876352e-04 0.015458858 -0.589280467 -0.49045803
## beta    0.387491358 -4.416468e-04 0.021719511  0.323510829  0.45405092
## omega   0.440127983 -4.969173e-05 0.019271687  0.383719655  0.50195224
## aH     -0.042717620 -3.444674e-05 0.023461041 -0.120944792  0.02179767
## aQ      0.130196697 -1.908086e-04 0.021222706  0.063408769  0.19921509
## h      -0.145297127 -6.329461e-04 0.022367799 -0.215223658 -0.07872360
## gammah -0.164366879 -2.909383e-04 0.021834058 -0.232107315 -0.08665109
## gammaQ  0.005940182 -3.368744e-04 0.021621717 -0.060206007  0.06936703
## lambda -0.004703382 -5.946155e-04 0.025302851 -0.093904864  0.06976988
## c       0.748928132  5.447422e-04 0.009254748  0.720093418  0.77449562
## S0     -0.027103806  2.611217e-02 0.021337699 -0.156897758  0.02458901
par(mar=c(7,4,4,2)+0.1)
plot(summary$original, main='Partial rank correlation coefficients', ylim=c(-1,1),
     xlab='', ylab='Coefficient',
     axes=FALSE)
axis(2)
axis(1, at=seq(1:16), labels=row.names(summary), las=2)
mtext(text='Parameter', side=1, line=4.5)
box()
for(i in 1:16) lines(c(i,i),c(summary[i,4], summary[i,5]))
abline(h=0)