VE calculation when vaccine is antigenically different to circulating infection

vaccination_antigenic_mismatch.R

# Analysis of VE during epidemic when vaccine protection < infection protection -------------------------------------------------------------
# Adam Kucharski (2024)


# Define assumptions about exposure over time in epidemic -----------------

xx <- seq(0,6,1/30) # time points for model
peak_time <- 3 # peak in months post-vaccine campaign
peak_width <- 0.8 # variance in epidemic width (i.e. sd of distribution)
exposure_prop <- 0.6 # proportion of population exposure during wave

cumulative_prop <- exposure_prop*pnorm(xx,mean=peak_time,sd=peak_width)
incidence_exposure <- exposure_prop*dnorm(xx,mean=peak_time,sd=peak_width)
incidence_exposure <- incidence_exposure*exposure_prop/sum(incidence_exposure) # normalise so cumulative = exposure_prop

# VE estimation -----------------------------------------------------------
# Assume leaky protection: https://www.nature.com/articles/s41467-023-40750-8

vaccine_protect <- 0.5 # vaccine protection against infection
infection_protect <- 0.95 # homologous infection protection against infection

# Simulate incidence in the two groups
# Leaky vaccination model

#Define ICs
unexposed_vaccinated <- 1
unexposed_unvaccinated <- 1

incidence_store_vaccinated <- NULL
incidence_store_unvaccinated <- NULL

for(ii in 1:length(xx)){

  exposure_ii <- incidence_exposure[ii]

  # Vaccinated group
  incidence_vaccinated <- exposure_ii*(unexposed_vaccinated*(1-vaccine_protect) + (1-unexposed_vaccinated)*(1-vaccine_protect)*(1-infection_protect))
  unexposed_vaccinated <- unexposed_vaccinated - incidence_vaccinated
  #unexposed_vaccinated_infected <- unexposed_vaccinated - incidence_vaccinated

  incidence_store_vaccinated <- c(incidence_store_vaccinated,incidence_vaccinated)

  # Unvaccinated group
  incidence_unvaccinated <- exposure_ii*(unexposed_unvaccinated + (1-unexposed_unvaccinated)*(1-infection_protect)) 
  unexposed_unvaccinated <- unexposed_unvaccinated - incidence_unvaccinated

  incidence_store_unvaccinated <- c(incidence_store_unvaccinated,incidence_unvaccinated)

}

# Calculate observed VE over time
ve_estimate <- 1-incidence_store_vaccinated/incidence_store_unvaccinated

# Plot comparative incidence -------------------------------------------------------------
# Shows infection incidence in vaccinated and unvaccinated, plus overall exposure incidence

par(mfrow=c(1,1),mgp=c(2.5,0.7,0),mar = c(3.5,3.5,1,3.5))

# Plot distribution of initial titres
plot(xx,incidence_exposure,xlab="months",ylab="probability",
     xaxs="i",yaxs="i",bty="l",type="l",lwd=2,col="darkred",main="",xlim=c(0,7))
lines(xx,incidence_store_vaccinated,col="darkgreen")
lines(xx,incidence_store_unvaccinated,col="darkorange")

# Plot curves -------------------------------------------------------------

par(mfrow=c(1,1),mgp=c(2.5,0.7,0),mar = c(3.5,3.5,1,3.5))

# Plot VE
plot(xx,ve_estimate,xlab="months",ylab="",
     xaxs="i",yaxs="i",yaxt="n",bty="l",type="l",lwd=2,col="darkred",main="",xlim=c(-0.5,max(xx)),ylim=c(0,0.6))
lines(c(-2,max(xx)),c(0,0),lty=2)
axis(2,col="dark red",col.axis="dark red",at=seq(-1,1,0.1),labels = sapply(seq(-100,100,10),function(x){paste(x,"%",sep="")}))
mtext("estimated vaccine effectiveness", side=2, line=2,col="dark red",cex=1) # Label for 2nd axis

graphics::arrows(x0=0,y0=0.55,x1=0,y1=-0.2,length=0.1,lwd=2)
graphics::text(x=-0.2,y=0.58,labels="vaccine campaign",adj=0)

lines(c(-2,max(xx)),c(vaccine_protect,vaccine_protect),lty=3,col="black")
graphics::text(x=3,y=0.52,labels="true VE",adj=0.5,col="black")

# Add epidemic dynamics

par(new=TRUE)
plot(xx,incidence_exposure,col="dark orange",xaxt="n",
     yaxt="n",xlab="",ylab="",bty="l",ylim=c(0,0.02),type="l",lwd=2,xlim=c(0,max(xx)),xaxs="i",yaxs="i")

axis(4,col="dark orange",col.axis="dark orange",at=seq(0,1,0.01),labels = sapply(seq(0,100,1),function(x){paste(x,"%",sep="")}))
mtext("daily exposure risk", side=4, line=2,col="dark orange",cex=1) # Label for 2nd axis

dev.copy(png,paste0("ve_waning.png"),units="cm",width=15,height=14,res=150)
dev.off()