SARS-CoV-2 Model (Sensitivity Analysis)

Clearing environment

rm(list = ls())

Loading necessary package

require(deSolve)

Loading required package: deSolve
package 㤼㸱deSolve㤼㸲 was built under R version 3.6.3

require(lhs)

Loading required package: lhs
package 㤼㸱lhs㤼㸲 was built under R version 3.6.3

System of DEs

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

Initial Values

times <- seq(0,180,by=1) 
covid.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)
covid.out <- as.data.frame(lsoda(ystart,times,SARSCOV2Model,covid.params))

Creating get.size Function

get.size <- function(out) tail(out$R,1)

Latin Hypercube Sampling

h1 <- 100
set.seed(6242015)
lhs <- maximinLHS(h1,16)

Max’s & Min’s for all parameters

b.max <- 1/0.1
b.min <- 1/5
K.max <- 1000000
K.min <- 1
g.max <- 1
g.min <- 0
p.max <- 1
p.min <- 0
q.max <- 1
q.min <- 0
beta.max <- 20
beta.min <- 1
delta.max <- 1/0.1
delta.min <- 1/10
c.max <- 20
c.min <- 1
S0.max <- 100000
S0.min <- 1
omega.max <- 1
omega.min <- 0
aH.max <- 1/0.1
aH.min <- 1/10
aQ.max <- 1/0.1
aQ.min <- 1/10
h.max <- 1
h.min <- 0
gammaH.max <- 1/0.1
gammaH.min <- 1/10
gammaQ.max <- 1/0.1
gammaQ.min <- 1/15
lambda.max <- 1/0.1
lambda.min <- 1/10

Creating Parameter Sets

params.set <- cbind(
  b = lhs[,1]*(b.max-b.min)+b.min,
  K = lhs[,2]*(K.max-K.min)+K.min,
  g = lhs[,3]*(g.max-g.min)+g.min,
  p = lhs[,4]*(p.max-p.min)+p.min,
  q = lhs[,5]*(q.max-q.min)+q.min,
  beta = lhs[,6]*(beta.max-beta.min)+beta.min,
  delta = lhs[,7]*(delta.max-delta.min)+delta.min,
  c = lhs[,8]*(c.max-c.min)+c.min,
  S0 = lhs[,9]*(S0.max-S0.min)+S0.min,
  omega = lhs[,10]*(omega.max-omega.min)+omega.min,
  aH = lhs[,11]*(aH.max-aH.min)+aH.min,
  aQ = lhs[,12]*(aQ.max-aQ.min)+aQ.min,
  h = lhs[,13]*(h.max-h.min)+h.min,
  gammaH = lhs[,14]*(gammaH.max-gammaH.min)+gammaH.min,
  gammaQ = lhs[,15]*(gammaQ.max-gammaQ.min)+gammaQ.min,
  lambda = lhs[,16]*(lambda.max-lambda.min)+lambda.min
)

levels <- 15
h2 <- 25

Cycling through Different Values of T

j <- 1
data <- data.frame(matrix(rep(NA,levels*h2*18),nrow=levels*h2))
for(i in 1:h2){
  for(T in times){
    data[j,1:17] <- covid.params1 <- as.list(c(params.set[i,],T=T))
    covid.out1 <- as.data.frame(lsoda(ystart,times,SARSCOV2Model,covid.params1))
    data[j,18] <- get.size(covid.out1)
    j <- j+1
  }
}
names(data) <- c(names(covid.params1),'outbreak.size')
save(data,file='data.Rdata')

Plotting Data

load('data.Rdata')

## S.h ##
plot(times, covid.out1$S.h,type='l', ylim=c(0,2e6), lwd=3, xlab='Days', main='Susceptibles (High)', ylab='Susceptibles (High)')
points(data$T,data$outbreak.size,pch=19,cex=0.3,col='blue')

## E.h ##
plot(times, covid.out1$E.h,type='l', log = 'y', ylim=c(10,2e6), lwd=3,xlab='Days', main='Exposed (High)', ylab='Exposed (High)')

1 y value <= 0 omitted from logarithmic plot

points(data$T,data$outbreak.size,pch=19,cex=0.3,col='blue')

## I.h ##
plot(times, covid.out1$I.h,type='l', log = 'y', ylim=c(10,2e6), lwd=3,xlab='Days', main='Infected (High)', ylab='Infected (High)')
points(data$T,data$outbreak.size,pch=19,cex=0.3,col='blue')

## S.l ##
plot(times, covid.out1$S.l,type='l', log = 'y', ylim=c(10,2e6), lwd=3,xlab='Days', main='Susceptibles (Low)', ylab='Susceptibles (Low)')
points(data$T,data$outbreak.size,pch=19,cex=0.3,col='blue')

## E.l ##
plot(times, covid.out1$E.l,type='l', log = 'y', ylim=c(10,2e6), lwd=3,xlab='Days', main='Exposed (Low)', ylab='Exposed (Low)')

1 y value <= 0 omitted from logarithmic plot

points(data$T,data$outbreak.size,pch=19,cex=0.3,col='blue')

## I.l ##
plot(times, covid.out1$I.l,type='l', log = 'y', ylim=c(10,2e6), lwd=3,xlab='Days', main='Infected (Low)', ylab='Infected (Low)')

2 y values <= 0 omitted from logarithmic plot

points(data$T,data$outbreak.size,pch=19,cex=0.3,col='blue')

## Q ##
plot(times, covid.out1$Q,type='l', log = 'y', ylim=c(10,2e6), lwd=3,xlab='Days', main='Self-Isolating', ylab='Self-Isolating')

1 y value <= 0 omitted from logarithmic plot

points(data$T,data$outbreak.size,pch=19,cex=0.3,col='blue')

## R ##
plot(times, covid.out1$R,type='l', log = 'y', ylim=c(10,2e6), lwd=3,xlab='Days', main='Recovered', ylab='Recovered')

1 y value <= 0 omitted from logarithmic plot

points(data$T,data$outbreak.size,pch=19,cex=0.3,col='blue')

## V ##
plot(times, covid.out1$V,type='l', log = 'y', ylim=c(10,2e6), lwd=3,xlab='Days', main='Virus in Environment', ylab='Virus in Environment')

1 y value <= 0 omitted from logarithmic plot

points(data$T,data$outbreak.size,pch=19,cex=0.3,col='blue')

Experimenting with Random Uniform and LH Sampling