SARS-CoV-2 Model (R0)

Clearing environment

rm(list = ls())

require(deSolve)
require(lhs)

R0 function

SIRR0 <- function(b,K,g,p,beta,delta,c,S0,omega,aH,h,gammaH){
  R0 <- (b*K*(g*(p-1)^2 + p)*beta*delta + b*c*(-g*(p-1)^3 + p)*S0*omega + (g-1)*(p-1)^2*(aH*h + gammaH - h*gammaH)*(-K*beta*delta + c*(p-1)*S0*omega)) / (b*K*(b*g - (g-1)*aH*h + gammaH - h*gammaH)*delta)
}

Initial Values

R0 <- SIRR0(b=1/0.77,K=100000,g=0,p=0.2,beta=13,delta=1,c=0,S0=10000,omega=0,aH=1/2.7,h=0.00082,gammaH=1/2.7)

Plotting all parameters in R0

par(mfrow=c(2,3))

## b ##
vals.b <- seq(0,1,by=0.01)
R0vals.b <- SIRR0(b=vals.b,K=100000,g=0,p=0.2,beta=13,delta=1,c=0,S0=10000,omega=0,aH=1/2.7,h=0.00082,gammaH=1/2.7)
plot(R0vals.b~vals.b,pch=19,main="R0 of Varying b",xlab="b",ylab="R0")

## K ##
vals.K <- seq(0,100000,by=1000)
R0vals.K <- SIRR0(b=1/0.77,K=vals.K,g=0,p=0.2,beta=13,delta=1,c=0,S0=10000,omega=0,aH=1/2.7,h=0.00082,gammaH=1/2.7)
plot(R0vals.K~vals.K,pch=19,main="R0 of Varying K",xlab="K",ylab="R0")

## g ##
vals.g <- seq(0,1,by=0.01)
R0vals.g <- SIRR0(b=1/0.77,K=100000,g=vals.g,p=0.2,beta=13,delta=1,c=0,S0=10000,omega=0,aH=1/2.7,h=0.00082,gammaH=1/2.7)
plot(R0vals.g~vals.g,pch=19,main="R0 of Varying g",xlab="g",ylab="R0")

## p ##
vals.p <- seq(0,1,by=0.01)
R0vals.p <- SIRR0(b=1/0.77,K=100000,g=0,p=vals.p,beta=13,delta=1,c=0,S0=10000,omega=0,aH=1/2.7,h=0.00082,gammaH=1/2.7)
plot(R0vals.p~vals.p,pch=19,main="R0 of Varying p",xlab="p",ylab="R0")

## beta ##
vals.beta <- seq(0,20,by=0.2)
R0vals.beta <- SIRR0(b=1/0.77,K=100000,g=0,p=0.2,beta=vals.beta,delta=1,c=0,S0=10000,omega=0,aH=1/2.7,h=0.00082,gammaH=1/2.7)
plot(R0vals.beta~vals.beta,pch=19,main="R0 of Varying beta",xlab="beta",ylab="R0")

## delta ##
vals.delta <- seq(0,1,by=0.01)
R0vals.delta <- SIRR0(b=1/0.77,K=100000,g=0,p=0.2,beta=13,delta=vals.delta,c=0,S0=10000,omega=0,aH=1/2.7,h=0.00082,gammaH=1/2.7)
plot(R0vals.delta~vals.delta,pch=19,main="R0 of Varying delta",xlab="delta",ylab="R0")

## c ##
vals.c <- seq(0,1,by=0.01)
R0vals.c <- SIRR0(b=1/0.77,K=100000,g=0,p=0.2,beta=13,delta=1,c=vals.c,S0=10000,omega=0,aH=1/2.7,h=0.00082,gammaH=1/2.7)
plot(R0vals.c~vals.c,pch=19,main="R0 of Varying c",xlab="c",ylab="R0")

## S0 ##
vals.S0 <- seq(0,10000,by=100)
R0vals.S0 <- SIRR0(b=1/0.77,K=100000,g=0,p=0.2,beta=13,delta=1,c=0,S0=vals.S0,omega=0,aH=1/2.7,h=0.00082,gammaH=1/2.7)
plot(R0vals.S0~vals.S0,pch=19,main="R0 of Varying S0",xlab="S0",ylab="R0")

## omega ##
vals.omega <- seq(0,1,by=0.01)
R0vals.omega <- SIRR0(b=1/0.77,K=100000,g=0,p=0.2,beta=13,delta=1,c=0,S0=10000,omega=vals.omega,aH=1/2.7,h=0.00082,gammaH=1/2.7)
plot(R0vals.omega~vals.omega,pch=19,main="R0 of Varying omega",xlab="omega",ylab="R0")

## aH ##
vals.aH <- seq(0,1,by=0.01)
R0vals.aH <- SIRR0(b=1/0.77,K=100000,g=0,p=0.2,beta=13,delta=1,c=0,S0=10000,omega=0,aH=vals.aH,h=0.00082,gammaH=1/2.7)
plot(R0vals.aH~vals.aH,pch=19,main="R0 of Varying aH",xlab="aH",ylab="R0")

## h ##
vals.h <- seq(0,0.001,by=0.00001)
R0vals.h <- SIRR0(b=1/0.77,K=100000,g=0,p=0.2,beta=13,delta=1,c=0,S0=10000,omega=0,aH=1/2.7,h=vals.h,gammaH=1/2.7)
plot(R0vals.h~vals.h,pch=19,main="R0 of Varying h",xlab="h",ylab="R0")

## gammaH ##
vals.gammaH <- seq(0,1,by=0.01)
R0vals.gammaH <- SIRR0(b=1/0.77,K=100000,g=0,p=0.2,beta=13,delta=1,c=0,S0=10000,omega=0,aH=1/2.7,h=0.00082,gammaH=vals.gammaH)
plot(R0vals.gammaH~vals.gammaH,pch=19,main="R0 of Varying gammaH",xlab="gammaH",ylab="R0")

Results of plots

# columns are parameters; rows are highest, lowest, and ending

results.R0 <- data.frame(
  b = c(max(R0vals.b), min(R0vals.b), tail(R0vals.b,n=1)),
  K = c(max(R0vals.K), min(R0vals.K), tail(R0vals.K,n=1)),
  g = c(max(R0vals.g), min(R0vals.g), tail(R0vals.g,n=1)),
  p = c(max(R0vals.p), min(R0vals.g), tail(R0vals.p,n=1)),
  beta = c(max(R0vals.beta), min(R0vals.beta), tail(R0vals.beta,n=1)),
  delta = c(max(R0vals.delta), min(R0vals.delta), tail(R0vals.delta,n=1)),
  c = c(max(R0vals.c), min(R0vals.c), tail(R0vals.c,n=1)),
  S0 = c(max(R0vals.S0), min(R0vals.S0), tail(R0vals.S0,n=1)),
  omega = c(max(R0vals.omega), min(R0vals.omega), tail(R0vals.omega,n=1)),
  aH = c(max(R0vals.aH), min(R0vals.aH), tail(R0vals.aH,n=1)),
  h = c(max(R0vals.h), min(R0vals.h), tail(R0vals.h,n=1)),
  gammaH = c(max(R0vals.gammaH), min(R0vals.gammaH), tail(R0vals.gammaH,n=1))
)

`.rowNamesDF<-`(results.R0,make.names=FALSE,c('Highest','Lowest','Ending'))

Sensitivity Analysis

x <- runif(50)
y <- runif(50)
h1 <- 50
lhs <- maximinLHS(h1,2)
par(mfrow=c(1,2))
plot(x,y,type='p',main='Random Uniform')
plot(lhs,type='p',main='LH Sampling')