#install.packages("deSolve")
library(deSolve)
## 
## Attaching package: 'deSolve'
## 
## The following object is masked from 'package:graphics':
## 
##     matplot
seir_model = function (current_timepoint, state_values, parameters)
{
    # create state variables (local variables)
    S = state_values [1]           #Initial number of susceptible individuals
    E = state_values [2]           #Initial number of exposed individuals
    I = state_values [3]           #Initial number of infectious individuals
    R = state_values [4]           #Initial number of recovered individuals
    
    with (
    as.list (parameters),          # variable names within parameters can be used
  {
    #compute derivatives
    dS = (-beta * S * I)
    dE = (beta * S * I) - (delta * E)
    dI = (delta * E) - (gamma * I)
    dR = (gamma * I)
    
    # combine results
    results = c (dS, dE, dI, dR)
    list (results)
    }
  )
}

Defining parameters

infectious_period = 4.5                       # infectious period
latent_period = 1.9                           # latent period
transmission_rate = beta_value = 0.7          # Transmission rate 
vaccine_compliance = 0                        # Vaccine compliance at the base case
vaccine_efficacy = 0.5                        # Vaccine efficacy

vac.total = vaccine_compliance * vaccine_efficacy

Compute values for delta and gamma

beta_value = transmission_rate 
delta_value = 1 / latent_period              
gamma_value = 1 / infectious_period

Compute R0

R0 = beta_value / gamma_value

Disease dynamics parameters

parameter_list = c (beta = beta_value, gamma = gamma_value, delta = delta_value)

Defining the initial condition

X = 9999             #Initial number of susceptible individuals
V = 0                #Initial number of exposed inidviduals
Y = 1                #Initial number of infectious individuals
Z = 0                #Initial number of recovered individuals

Compute total population

N = X + V + Y + Z

Initial state values for the differential equations

initial_values = c (S=(X*(1-vac.total))/N,E=V/N,I=Y/N,R=(Z+X*vac.total)/N)

Output timepoints

timepoints = seq (1, 100, by=1)

Simulate the SEIR epidemic

output = lsoda (initial_values, timepoints, seir_model, parameter_list)

Base case plot

plot (S~time, data=output, type="l", lwd=2, col=1, ylab="", ylim=c(0,1), main="SEIR Epidemic: Base Case")
lines (E~time, data=output, type="l", lwd=2, col=2, ylab="")
lines (I~time, data=output, type="l", lwd=2, col=3, ylab="")
lines (R~time, data=output, type="l", lwd=2, col=4, ylab="")
legend ("topright", legend=c("S","E","I","R"), col=c(1,2,3,4), lwd=2, cex=0.7)

Attack rate and total health care cost of treatment (base case)

treatment_cost = 500       # per infected person
Attack_rate = sum (output[,4]) / (1/gamma_value)
basecase_cost = treatment_cost * Attack_rate * N

#Summary
print (paste("Base case attack rate=", round(Attack_rate,3)))
## [1] "Base case attack rate= 0.95"
print (paste("Base case cost = $", round(basecase_cost,3)))
## [1] "Base case cost = $ 4748966.747"

COMMENT

The attack rate of the base case in the outbreak is 0.95 (95%) and the total treatment cost is $4,748,966.747

Vaccine Intervention

vaccine_compliance = 0.4
vaccine_efficacy = seq (0.4, 0.6, 0.2)     # Vaccine efficacies
cost_per_vaccine = 20                      # average vaccination cost per person

par(mfrow=c(2,2))
for (i in vaccine_efficacy) {
  vaccine_efficacy = i
  vac.total = vaccine_efficacy * vaccine_compliance
  initial_values = c(S=(X*(1-vac.total))/N,E=V/N,I=Y/N,R=(Z+X*vac.total)/N)
  output = lsoda (initial_values, timepoints, seir_model, parameter_list)
  
  # plot intervention
  plot (S~time, data=output,type="l", lwd=2, col=1, ylab="", ylim=c(0,1), main=paste("vaccine efficacy=",i))
  lines (E~time, data=output, type="l", lwd=2, col=2, ylab="")
  lines (I~time, data=output, type="l", lwd=2, col=3, ylab="")
  lines (R~time, data=output, type="l", lwd=2, col=4, ylab="")
  legend ("topright", legend=c("S","E","I","R"), col=c(1,2,3,4), lwd=2,cex=0.7)

   # Cost
   vaccine_cost = cost_per_vaccine * vac.total * X
   new_Attack_rate = sum (output[,4])/(1/gamma_value)
 
   new_treatment_cost = treatment_cost * new_Attack_rate * N 
   total_cost = vaccine_cost + new_treatment_cost
  
   # Incremental cost effectiveness ratio
   ICER = (total_cost - basecase_cost) / {(new_Attack_rate * N) - (Attack_rate * N)}
    
   #Summary
   print (paste("Base case attack rate =", round(Attack_rate,3),"Base case cost =",round(basecase_cost,3)))
   print (paste("Vaccine compliance =",vaccine_compliance, "Vaccine efficacy =",i,", Attack_rate=", round(new_Attack_rate,3), "Total Healthcare Cost= $", round(total_cost,3), "Incremental cost effectiveness ratio = ", round(ICER,3)))
}
## [1] "Base case attack rate = 0.95 Base case cost = 4748966.747"
## [1] "Vaccine compliance = 0.4 Vaccine efficacy = 0.4 , Attack_rate= 0.764 Total Healthcare Cost= $ 3852833.06 Incremental cost effectiveness ratio =  482.763"
## [1] "Base case attack rate = 0.95 Base case cost = 4748966.747"
## [1] "Vaccine compliance = 0.4 Vaccine efficacy = 0.6 , Attack_rate= 0.666 Total Healthcare Cost= $ 3378568.928 Incremental cost effectiveness ratio =  483.081"

COMMENT

We register an attack rate of 0.764 (76.4%) when the vaccine efficacy is 40%. Total health care cost is $3,852,833.06 and ICER is $482.76 per case averted. When the vaccine efficacy is 60%, we register an attack rate of 0.666 (66.6%). In this case, the total health care cost is $3,378,568.92 and ICER is $483.08 per case averted.

vaccine_compliance = seq (0.5, 0.6, 0.1)   # Vaccine compliances
vaccine_efficacy = seq (0.4, 0.6, 0.2)     # Vaccine efficacies
cost_per_vaccine = 20                      # Average vaccination cost per person

par(mfrow=c(2,2))
for (i in vaccine_efficacy) 
for (ii in vaccine_compliance)  {
  vaccine_efficacy = i
  vaccine_compliance = ii 
  vac.total = vaccine_efficacy * vaccine_compliance
  initial_values = c(S=(X*(1-vac.total))/N,E=V/N,I=Y/N, R=(Z+X * vac.total)/N)
  output = lsoda (initial_values, timepoints, seir_model, parameter_list)
  
  # plot intervention
  plot (S~time, data=output, type="l", lwd=2, col=1, ylab="", ylim=c(0,1), main=paste("vaccine efficacy=",i, "compliance=",ii))
  lines (E~time, data=output, type="l", lwd=2, col=2, ylab="")
  lines (I~time, data=output, type="l", lwd=2, col=3, ylab="")
  lines (R~time, data=output, type="l", lwd=2, col=4, ylab="")
  legend ("right", legend=c("S","E","I","R"), col=c(1,2,3,4), lwd=2, cex=0.7)

   # Cost
   vaccine_cost = cost_per_vaccine * vac.total * X
   new_Attack_rate = sum(output[,4])/(1/gamma_value)
 
   new_treatment_cost = treatment_cost * new_Attack_rate * N 
   total_cost = vaccine_cost + new_treatment_cost
  
   # Incremental cost effectiveness ratio
   ICER = (total_cost - basecase_cost) / {(new_Attack_rate * N) - (Attack_rate * N)}
   
   #Summary
   print (paste("Basecase attack rate =", round(Attack_rate,3), "Basecase cost =", round(basecase_cost,3)))
   print (paste("Vaccine compliance=",ii, "Vaccine efficacy=",i, "Attack_rate=",round(new_Attack_rate,3), "Total Healthcare Cost= $", round(total_cost,3),"Incremental cost effectiveness ratio = ", round(ICER,3)))
}
## [1] "Basecase attack rate = 0.95 Basecase cost = 4748966.747"
## [1] "Vaccine compliance= 0.5 Vaccine efficacy= 0.4 Attack_rate= 0.716 Total Healthcare Cost= $ 3618629.858 Incremental cost effectiveness ratio =  482.913"
## [1] "Basecase attack rate = 0.95 Basecase cost = 4748966.747"
## [1] "Vaccine compliance= 0.6 Vaccine efficacy= 0.4 Attack_rate= 0.666 Total Healthcare Cost= $ 3378568.928 Incremental cost effectiveness ratio =  483.081"
## [1] "Basecase attack rate = 0.95 Basecase cost = 4748966.747"
## [1] "Vaccine compliance= 0.6 Vaccine efficacy= 0.6 Attack_rate= 0.504 Total Healthcare Cost= $ 2590156.109 Incremental cost effectiveness ratio =  483.864"

COMMENT

We register an attack rate of 0.716 (71,6%) when the vaccine efficacy is 40% and compliance 50%, with total health care cost of $3,618,629.85 and ICER of $482.91 per case averted. We register an attack rate of 0.666 (66,6%) when the vaccine efficacy is 40% and compliance 60%, with total health care cost of $3,378,568.92 and ICER of $483.08 per case averted

We register an attack rate of 0.504 (50.4%) when the vaccine efficacy is 60% and compliance 60%, with total health care cost of $2,590,156.1 and ICER is $483.86 per case averted.

#Summary
print(ICER <- structure (list('vaccine efficacy 40%' = c(482.76, 482.91, 483.08),'vaccine efficacy 60%' =  c(483.08, 483.38, 483.83)), .Names = c("vaccine efficacy 40%", "vaccine efficacy 60%"), row.names = c("compliance 40%", "compliance 50%", "compliance 60%"), class="data.frame"))
##                vaccine efficacy 40% vaccine efficacy 60%
## compliance 40%               482.76               483.08
## compliance 50%               482.91               483.38
## compliance 60%               483.08               483.83
print(Attack_rate <- structure (list('vaccine efficacy 40%' = c(0.764, 0.712, 0.666),'vaccine efficacy 60%' =  c(0.666, 0.588, 0.504)), .Names = c("vaccine efficacy 40%", "vaccine efficacy 60%"), row.names = c("compliance 40%", "compliance 50%", "compliance 60%"), class="data.frame"))
##                vaccine efficacy 40% vaccine efficacy 60%
## compliance 40%                0.764                0.666
## compliance 50%                0.712                0.588
## compliance 60%                0.666                0.504