rm(list=ls())
require(deSolve)
## Loading required package: deSolve
require(lhs)
## Loading required package: lhs
## Warning: package 'lhs' was built under R version 4.0.2
SARSCOV2Model <- function (t, y, params) {
  
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, y),
  {
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 and low risk groups from total population size 100000
}
peak.size <- function (out)
{
  max(out$I.h)
}
SIRR0<-function(b,k,g,p,beta,delta,c,S0=100000,omega,aH,gammaH,h,q){
    
with(
as.list(params),
{
R0<-(b*k*(g*(-1+ p)^2+p)*beta*q*delta+b*c*(-g*(-1+p)^3+p)*S0*omega+(-1+g)*(-1+p)^2*(aH*h+gammah-h*gammah)*(-k*beta*q*delta+c*(-1+p)*S0*omega))/(b*k*(b*g-(-1+g)*(aH*h+gammah-h*gammah))*delta)
} 
)
}
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)

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
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
params.set<-cbind(
g=lhs[,1]*(g.max-g.min)+g.min,
#q=lhs[,2]*(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,
#S0 = lhs[,9]*(S0.max-S0.min)+S0.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
)
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)
qseq<-seq(0,.1, by=.005)
levels<-length(qseq)
j <- 1
 data <- data.frame(matrix(rep(NA,levels*h2*19),nrow=levels*h2))
 for(i in 1:h2){
   for (q in qseq){
data[j,1:16] <- params <- as.list(c(params.set[i,], q=q, S0=100000)) #16 parameters
out <- (as.data.frame(ode(ystart,times, SARSCOV2Model, params)))
data[j,17] <- get.size(out)
data[j,18] <- SIRR0(params)
data[j,19] <- peak.size(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)
plot(data$q, data$R0, pch=19, cex=0.3, col='blue')

plot(data$q, data$outbreak.size/100000, pch=19, cex=0.3, col='blue') #divide epidemic final size by 100000 to get proportion

boxplot(data$outbreak.size/100000~data$q,  border='blue', log='y', pch='.')

plot(data$q, data$peak.size/100000, pch=19, cex=0.3, col='blue')

boxplot(data$peak.size/100000~data$q,  border='blue', log='y', pch='.')

boxplot(data$R0~data$q, border='blue', log='y', pch='.')

library(sensitivity)
## Warning: package 'sensitivity' was built under R version 4.0.2
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')
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.0008305276 -1.178127e-03 0.023973693 -0.07157297  0.081469332
## p      -0.2037078640 -8.087350e-05 0.020817064 -0.26405936 -0.143693112
## delta  -0.5850526364 -6.617860e-04 0.012882507 -0.62230724 -0.541584400
## b      -0.8171743342 -5.176645e-04 0.010332405 -0.84855912 -0.787269255
## k      -0.6189582594 -6.200367e-04 0.011911064 -0.65372323 -0.586026021
## beta    0.0936773940  1.324700e-03 0.023613190  0.03038058  0.157523916
## c       0.5959301927  8.337972e-04 0.013074036  0.56040254  0.632739247
## omega   0.5445632785  6.842973e-05 0.015858582  0.49644421  0.597275262
## aH     -0.0245303207  3.180153e-04 0.022378504 -0.10003868  0.039274301
## aQ      0.0309021448  9.718998e-04 0.018906465 -0.02454690  0.088253019
## h      -0.0367333953 -6.544705e-05 0.021326649 -0.10825535  0.027335705
## gammah  0.0729871818 -6.825706e-04 0.021037136  0.01498127  0.145458291
## gammaQ  0.2155707006 -8.160364e-04 0.021791102  0.14964040  0.286571883
## lambda -0.0806578252 -5.621312e-04 0.023850987 -0.15192654 -0.007319971
## q       0.9117831746  3.911217e-04 0.004648819  0.89806129  0.928553044
## S0     -0.0370946765  3.677662e-02 0.021412776 -0.18444484 -0.010559858
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)

```