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){
  
  #Input arguments are given as a list. This was modified from Sheridan's code on Google Drive.
  
  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[,1]*(q.max-q.min)+q.min,
  p=lhs[,2]*(p.max-p.min)+p.min,
  delta=lhs[,3]*(delta.max-delta.min)+delta.min,
  b = lhs[,4]*(b.max-b.min)+b.min,
  K = lhs[,5]*(K.max-K.min)+K.min,
  beta = lhs[,6]*(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)
gseq<-seq(0,1, by=0.025)
levels<-length(gseq)

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 (g in gseq){
data[j,1:16] <- params <- as.list(c(params.set[i,], g=g, 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$g, data$R0, pch=19, cex=0.3, col='blue', xlab='Proportion of high risk individuals that self-isolate (g)', ylab='R0')

plot(data$g, data$outbreak.size/100000, pch=19, cex=0.3, col='blue', xlab='Proportion of high risk individuals that self-isolate (g)', ylab='Proportion of Epidemic Size') #divide epidemic final size by 100000 to get proportion

plot(data$g, data$peak.size/100000, pch=19, cex=0.3, col='blue', xlab='Proportion of high risk individuals that self-isolate (g)', 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$g,  border='blue', log='y', pch='.', xlab='Proportion of high risk individuals that self-isolate (g)', ylab='Proportion of Infected (High)')

boxplot(data$outbreak.size/100000~data$g,  border='blue', log='y', pch='.', xlab='Proportion of high risk individuals that self-isolate (g)', ylab='Proportion of Epidemic Size')

boxplot(data$R0~data$g, border='blue', log='y', pch='.', xlab='Proportion of high risk individuals that self-isolate (g)', 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.
## q       0.909796539  1.305029e-04 0.002339298  0.90315261  0.9171303736
## p      -0.253050618  2.996178e-04 0.013850946 -0.29934448 -0.2136402642
## delta  -0.574005437 -3.165269e-06 0.010289606 -0.60577609 -0.5443051692
## b      -0.852968045 -2.248241e-04 0.004545314 -0.86856359 -0.8392683832
## K      -0.613337115 -1.502886e-04 0.009977978 -0.64479543 -0.5795147872
## beta    0.306295660  7.536269e-04 0.016728191  0.25292166  0.3540530209
## c       0.546333747 -1.155479e-04 0.010458298  0.51329588  0.5791052637
## omega   0.481357184  3.599428e-05 0.011056508  0.44518836  0.5141437168
## aH      0.078093407  1.303721e-04 0.015828482  0.02851366  0.1289840925
## aQ      0.097174833  1.865933e-04 0.015093586  0.04583040  0.1419640861
## h       0.068433860 -8.177493e-04 0.014969057  0.01559213  0.1126726296
## gammah -0.237236083 -4.563170e-04 0.016333087 -0.28414224 -0.1895628417
## gammaQ  0.085239239 -1.945806e-04 0.015992440  0.03224657  0.1370402868
## lambda -0.057354544 -5.711785e-04 0.017174340 -0.10976134  0.0001338556
## g      -0.819873456 -2.852738e-04 0.005512659 -0.83490646 -0.8024108600
## S0     -0.002235468  1.535936e-03 0.015465722 -0.06377899  0.0369372470
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)