SARS-CoV-2 Model (Initial Values)

Clearing environment

rm(list = ls())

Loading necessary package

require(deSolve)

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
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)
  }
)
}

Initial Values

times<-seq(0,180,by=1) 
covid.params<-c(q=0.5,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=100000)
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))

Plotting all classes with initial values

op1 <- par(mar=c(6,6,2,2))
plot(covid.out$S.h~covid.out$time,type="l", col="blue", xlab="Days", ylab="Individuals", ylim=c(0,80000))
lines(covid.out$E.h~covid.out$time, col="red")

lines(covid.out$I.h~covid.out$time, col="purple")
lines(covid.out$S.l~covid.out$time, col="green")

lines(covid.out$E.l~covid.out$time, col="orange")
lines(covid.out$I.l~covid.out$time, col="darkturquoise")

lines(covid.out$Q~covid.out$time, col="orchid2")
lines(covid.out$R~covid.out$time, col="deeppink1")

legend(100, 70000,legend = c("Susceptible (High)","Exposed (High)","Infected (High)","Susceptible (Low)","Exposed (Low)","Infected (Low)","Self-Isolating","Recovered"), col = c("blue","red","purple","green","orange","darkturquoise","orchid2","deeppink1"), lty=1, cex=0.8)
par(op1)

op2 <- par(mar=c(6,6,2,2))

plot(covid.out$V~covid.out$time, type="l", col="forestgreen", xlab="Days", ylab="Virus in Environment")
par(op2)

Results of Plot

# first dataframe lists highest values
# second dataframe lists lowest values
# third dataframe lists ending values

results <- data.frame(
  S.h = c(max(covid.out[ ,2]),min(covid.out[ ,2]),tail(covid.out[ ,2],n=1)),
  E.h = c(max(covid.out[ ,3]),min(covid.out[ ,3]),tail(covid.out[ ,3],n=1)),
  I.h = c(max(covid.out[ ,4]),min(covid.out[ ,4]),tail(covid.out[ ,4],n=1)),
  S.l = c(max(covid.out[ ,5]),min(covid.out[ ,5]),tail(covid.out[ ,5],n=1)),
  E.l = c(max(covid.out[ ,6]),min(covid.out[ ,6]),tail(covid.out[ ,6],n=1)),
  I.l = c(max(covid.out[ ,7]),min(covid.out[ ,7]),tail(covid.out[ ,7],n=1)),
  Q = c(max(covid.out[ ,8]),min(covid.out[ ,8]),tail(covid.out[ ,8],n=1)),
  R = c(max(covid.out[ ,9]),min(covid.out[ ,9]),tail(covid.out[ ,9],n=1)),
  V = c(max(covid.out[ ,10]),min(covid.out[ ,10]),tail(covid.out[ ,10],n=1))
)

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