In jameshay218/driftSim: Short-sighted evolution of influenza cellular receptor binding avidity shapes influenza epidemic dynamics
library(driftSim)
library(ggplot2)
attach(exampleParameters)
knitr::opts_chunk$set(echo = FALSE, fig.width=7,fig.height=6,
message=FALSE, results=FALSE,warning=FALSE)
p = viruspars["p"] #' parameter to control degree by which changes in binding avidity affect probability of escape from immune response
r = viruspars["r"] #' parameter to control degree by which previous exposure reduce probability of immune escape
b = viruspars["b"] #' parameter to control the shape of the relationship between probability of successful replication and changes in binding avidity
a =viruspars["a"] #' controls rate of changes of relationship between probability of successful replication and change in binding avidity.
n = viruspars["n"] #' average number of virus copies: n is number of offspring per virus replication
v = viruspars["v"] #' v is number of virions initialyl transmitted
nv = n*v
q = viruspars["q"] #' parameter to control the shape of the relationship between binding avidity and immune escape (shift on the x-axis)
V = seq(0, 3, by = 0.005) #' Binding avidity
N_reinfect = as.numeric(25)
max_reinfect = N_reinfect
delta = as.numeric(5)
#' Get a colour scale gradient blue to red
f_array <- NULL #' Survival prob
df_array <- NULL #' Derivative of survival prob
for(k in 0:(N_reinfect-1)){
probTarget = exp(-p*(V+q)) #' probability of being targetted by immune system. As binding avidity increases, this probability decreases
probEscape = 1-probTarget #' probability of escape
immJ = r*(k- delta) #' Strength of host immune respone. As k increases, virus must escape more antibodies. As delta increases, this effect diminishes as infecting virus is further from host immunity.
if(immJ < 0) immJ = 0
f = (probEscape)^(immJ) #' probability of escape from immunity
if(k >= 1) f_dash= immJ*p*((1-probTarget)^(immJ-1))*(probTarget) #' derivative of this. ie. rate of change of relationship between binding avidity and probability of immune escape
else f_dash= rep(0, length(V))
f_array[[k+1]] <- f
df_array[[k+1]] <- f_dash
}
probInf <- exp(-a*(V^b)) #' probability of successful infection within a host. as binding avidity increases in a naive host, chance of successfully replicating decreases
probInf_dash= -a*b*(V^(b-1))*exp(-a*(V^b)) #' rate of change of this relationship
rho_array <- NULL
dV_array <- NULL
probRep_array <- NULL
for(i in 1:length(df_array)){
R0 = f_array[[i]]*probInf*n
rho = 1 - R0^-v
rho[rho < 0] = 0
rho_array[[i]] <- rho
dV = df_array[[i]]*probInf + f_array[[i]]*probInf_dash
dV_array[[i]] <- dV
probReplication = f_array[[i]]*probInf
probRep_array[[i]] = probReplication
}
rho0 = max(rho_array[[1]])
rho1 = max(rho_array[[2]])
rho2 = max(rho_array[[3]])
probRep_data <- NULL
for(i in 1:length(probRep_array)){
data <- data.frame(x=V,y=probRep_array[[i]],z=as.character(i))
probRep_data <- rbind(probRep_data,data)
}
colours <- NULL
blue <- rev(seq(0,1,by=1/N_reinfect))
red <- seq(0,1,by=1/N_reinfect)
for(i in 1:N_reinfect){
colours[i] <- rgb((i-1)/(max_reinfect-1),0,(N_reinfect-i)/(max_reinfect-1))
}
A <- ggplot() + geom_line(data=probRep_data,aes(x=x,y=y,colour=z)) + scale_color_manual(values=colours) +
theme(plot.title=element_text(hjust=0.5),
panel.grid.major=element_blank(),
panel.grid.minor=element_blank(),
panel.background=element_blank(),
axis.text.x=element_text(colour="gray20",size=12),
axis.text.y = element_text(colour="gray20",size=12),
text = element_text(colour="gray20",size=14),
axis.line=element_line(colour="gray20"),
axis.line.x =element_line(colour="gray20"),
axis.line.y=element_line(colour="gray20"),
legend.position = "none",
axis.title=element_text(size=12)
) +
scale_x_continuous(expand=c(0,0)) +
scale_y_continuous(limits=c(0,1),expand=c(0,0)) +
ylab("Probability of Successful\n Replication Within a Host (g(V))") +
xlab("Binding Avidity") + ggtitle("Plot A")
rho_data<- NULL
for(i in 1:length(rho_array)){
data <- data.frame(x=V,y=rho_array[[i]],z=as.character(i))
rho_data<- rbind(rho_data,data)
}
B <- ggplot() + geom_line(data=rho_data,aes(x=x,y=y,colour=z)) + scale_color_manual(values=colours) +
theme(plot.title=element_text(hjust=0.5),
panel.grid.major=element_blank(),
panel.grid.minor=element_blank(),
panel.background=element_blank(),
axis.text.x=element_text(colour="gray20",size=12),
axis.text.y = element_text(colour="gray20",size=12),
text = element_text(colour="gray20",size=14),
axis.line=element_line(colour="gray20"),
axis.line.x =element_line(colour="gray20"),
axis.line.y=element_line(colour="gray20"),
legend.position = "none",
axis.title=element_text(size=12)
) +
scale_x_continuous(expand=c(0,0)) +
scale_y_continuous(limits=c(0,1),expand=c(0,0)) +
ylab(expression(Probability~of~Infection~Between~Hosts~rho)) +
xlab("Binding Avidity") +
ggtitle("Plot B") +
geom_hline(yintercept = rho0,linetype="longdash",colour="dodgerblue4")+
geom_hline(yintercept = rho1,linetype="longdash",colour="dodgerblue4")+
geom_hline(yintercept = rho2,linetype="longdash",colour="dodgerblue4")
df_data<- NULL
for(i in 1:length(df_array)){
data <- data.frame(x=V,y=df_array[[i]],z=as.character(i))
df_data<- rbind(df_data,data)
}
#' Derivative of f
C <- ggplot() + geom_line(data=df_data,aes(x=x,y=y,colour=z)) + scale_color_manual(values=colours) +
theme(plot.title=element_text(hjust=0.5),
panel.grid.major=element_blank(),
panel.grid.minor=element_blank(),
panel.background=element_blank(),
axis.text.x=element_text(colour="gray20",size=12),
axis.text.y = element_text(colour="gray20",size=12),
text = element_text(colour="gray20",size=14),
axis.line=element_line(colour="gray20"),
axis.line.x =element_line(colour="gray20"),
axis.line.y=element_line(colour="gray20"),
legend.position = "none",
axis.title=element_text(size=12)
) +
scale_x_continuous(expand=c(0,0)) +
scale_y_continuous(expand=c(0,0)) +
ylab(expression(d*f/d*V)) +
xlab("Binding Avidity") + ggtitle("Plot C")
dV_data<- NULL
for(i in 1:length(dV_array)){
data <- data.frame(x=V,y=dV_array[[i]],z=as.character(i))
dV_data<- rbind(dV_data,data)
}
#' Infection
D <- ggplot() + geom_line(data=dV_data,aes(x=x,y=y,colour=z)) + scale_color_manual(values=colours) +
theme(plot.title=element_text(hjust=0.5),
panel.grid.major=element_blank(),
panel.grid.minor=element_blank(),
panel.background=element_blank(),
axis.text.x=element_text(colour="gray20",size=12),
axis.text.y = element_text(colour="gray20",size=12),
text = element_text(colour="gray20",size=12),
axis.line=element_line(colour="gray20"),
axis.line.x =element_line(colour="gray20"),
axis.line.y=element_line(colour="gray20"),
legend.position = "none",
axis.title=element_text(size=12)
) +
scale_x_continuous(expand=c(0,0)) +
scale_y_continuous(expand=c(0,0)) +
ylab(expression(d*beta/d*V)) +
xlab("Binding Avidity") +
geom_hline(yintercept=0,linetype='longdash',colour="gray20") + ggtitle("Plot D")
f_data<- NULL
for(i in 1:length(f_array)){
data <- data.frame(x=V,y=f_array[[i]],z=as.character(i))
f_data <- rbind(f_data,data)
}
E <- ggplot() + geom_line(data=f_data,aes(x=x,y=y,colour=z)) + scale_color_manual(values=colours) +
theme(plot.title=element_text(hjust=0.5),
panel.grid.major=element_blank(),
panel.grid.minor=element_blank(),
panel.background=element_blank(),
axis.text.x=element_text(colour="gray20",size=12),
axis.text.y = element_text(colour="gray20",size=12),
text = element_text(colour="gray20",size=14),
axis.line=element_line(colour="gray20"),
axis.line.x =element_line(colour="gray20"),
axis.line.y=element_line(colour="gray20"),
legend.position = "none",
axis.title=element_text(size=12)
) +
scale_x_continuous(expand=c(0,0)) +
scale_y_continuous(expand=c(0,0)) +
ylab("Probability of Evading Immune System, f(j,V)") +
xlab("Binding Avidity") + ggtitle("Plot E")
probRep_dat <- data.frame(x=V,y=probInf)
F <- ggplot() + geom_line(data=probRep_dat,aes(x=x,y=y,colour="red")) +
theme(plot.title=element_text(hjust=0.5),
panel.grid.major=element_blank(),
panel.grid.minor=element_blank(),
panel.background=element_blank(),
axis.text.x=element_text(colour="gray20",size=12),
axis.text.y = element_text(colour="gray20",size=12),
text = element_text(colour="gray20",size=14),
axis.line=element_line(colour="gray20"),
axis.line.x =element_line(colour="gray20"),
axis.line.y=element_line(colour="gray20"),
legend.position = "none",
axis.title=element_text(size=12)
) +
scale_x_continuous(expand=c(0,0)) +
scale_y_continuous(expand=c(0,0)) +
ylab("Probability of Successful Replication, g(V)") +
xlab("Binding Avidity") + ggtitle("Plot F")
#' Plot antigenic distance against j at a given binding avidity.
d <- seq(0,N_reinfect,by=0.1)
fixedV <- 0.5
delta_array <- NULL
for(k in 0:(N_reinfect-1)){
probT = exp(-p*(fixedV+q))
probE = 1- probT
immJ1 = r*(k-d)
immJ1[immJ1 < 0] <- 0
probSurvival1 = 1 - (n*exp(-a*(fixedV^b))*(probE^immJ1))^-v
probSurvival1[probSurvival1 < 0] <- 0
delta_array[[k+1]] <- probSurvival1
}
delta_data<- NULL
for(i in 1:length(delta_array)){
data <- data.frame(x=d,y=delta_array[[i]],z=as.character(i))
delta_data <- rbind(delta_data,data)
}
G <- ggplot() + geom_line(data=delta_data,aes(x=x,y=y,colour=z)) + scale_color_manual(values=colours) +
theme(
plot.title=element_text(hjust=0.5),
panel.grid.major=element_blank(),
panel.grid.minor=element_blank(),
panel.background=element_blank(),
axis.text.x=element_text(colour="gray20",size=12),
axis.text.y = element_text(colour="gray20",size=12),
text = element_text(colour="gray20",size=14),
axis.line=element_line(colour="gray20"),
axis.line.x =element_line(colour="gray20"),
axis.line.y=element_line(colour="gray20"),
legend.position = "none",
axis.title=element_text(size=12)
) +
scale_x_continuous(expand=c(0,0)) +
scale_y_continuous(limits=c(0,1),expand=c(0,0)) +
ylab("Probability of Infection Between Hosts") +
xlab("Antigenic Distance to Host Immunity") + ggtitle("Plot G")
immJ_array <- NULL
immJ2 <- seq(0,r*N_reinfect,by=1)
for(k in 0:(N_reinfect-1)){
probT1 = exp(-p*(fixedV+q))
probE1 = 1- probT1
probSurvival2 = 1 - (n*exp(-a*(fixedV^b))*(probE1^immJ2))^-v
probSurvival2[probSurvival2 < 0] <- 0
immJ_array[[k+1]] <- probSurvival2
}
immJ_data<- NULL
for(i in 1:length(immJ_array)){
data <- data.frame(x=immJ2/r,y=immJ_array[[i]],z=as.character(i))
immJ_data <- rbind(immJ_data,data)
}
H <- ggplot() + geom_line(data=immJ_data,aes(x=x,y=y,colour=z)) + scale_color_manual(values=colours) +
theme(plot.title=element_text(hjust=0.5),
panel.grid.major=element_blank(),
panel.grid.minor=element_blank(),
panel.background=element_blank(),
axis.text.x=element_text(colour="gray20",size=12),
axis.text.y = element_text(colour="gray20",size=12),
text = element_text(colour="gray20",size=14),
axis.line=element_line(colour="gray20"),
axis.line.x =element_line(colour="gray20"),
axis.line.y=element_line(colour="gray20"),
legend.position = "none",
axis.title=element_text(size=12)
) +
scale_x_continuous(expand=c(0,0)) +
scale_y_continuous(limits=c(0,1),expand=c(0,0)) +
ylab("Probability of Infection Between Hosts") +
xlab("immJ") + ggtitle("Plot H")
Foreword
The purpose of this vignette is to describe the theory behind the model of the driftSim package. Here, we will describe the model that governs the within-host dynamics of influenza binding avidity adaptation. Another way of understanding these equations is to check out the second (parameters) panel of the driftSim shiny app, which can be run with driftSim::runSimulationApp(). A useful piece of background reading can be found in a paper by Yuan et al.
Note that this document does not describe the transmission model itself, which is described here.
Model equations
The plots in the shiny app are governed by the equations shown below. These equations make up a model to represent how changes in binding avidity impact the various stages from host infection through to reinfection. The model also considers the immune history of the host, $\textbf{h}$, and the level of immunity, $j$, a host has, which is represented by the blue-red lines in the plots (blue $j=0$, red $j=j_{max}$). For the purposes of these plots, a higher value of $j$ simply represents higher levels of antibody-mediated immunity.
1. Probability of evading immune system
After entering a host, the virus must first evade the immune response (Equation 1). Here, the probability of escaping the immune response increases as binding avidity increases as shown in plot A. The virus must also evade the host's antibody-mediated immunity conferred from all previous infections ($\textbf{h}$), adjusted by the antigenic distance between the infecting virus and the viruses that elicited the host's previous antibody response.
E + ggtitle("Plot A")
As hosts will have encountered a number of infections in their lifetime, the infection history $\textbf{h}$ for a given host is the vector of all viruses that the host has seen. This part of the model is described in greater detail in Section 3.3 of the accompanying transmission model vignette.
In summary, the level of immunity that a host has against an infecting virus is the sum of the host's antibody repetoire less the minimum antigenic distance between the infecting virus and all past infections in the host's immune history, represented by the term $r(j - \delta_{il})$ below. $i$ represents the infecting strain and $l$ represents the most antigenically similar strain in the host's infection history. The definition of antigenic distance is analagous to the concepts described by Smith et al. 2004
The rate of change of $f$ with respect to binding avidity is shown in plot B.
C + ggtitle("Plot B")
\begin{equation}
f(k,\delta_{il},V_i) = [1-e^{-p(V_i+q)}]^{r(j - \delta_{il})}
\end{equation}
Where $V_i$ is the binding avidity of virus $i$; $j$ is the total antibody titre of the host; $\delta_{il}$ is the antigenic distance between strain $i$ and the host immunity virus, $l$; and the remaining parameters are defined at the bottom of this document (or with ?exampleParameters from the package).
2. Probability of successful replication within host
Binding avidity also affects how well a virus is able to replicate within the host, as described by Equation 2. This relationship is shown in plot C. The naive case (ie. $j=0$) is shown in plot D.
\begin{equation}
g(V_i) = e^{-aV_i^b}
\end{equation}
cowplot::plot_grid(A + ggtitle("Plot C"), F + ggtitle("Plot D"),ncol=2)
3. Probability of successful within host infection
To successfully infect a host, the virus must escape the immune system (Equation 1) and successfully replicate (Equation 2). The probability of a successful within-host infection is given by:
\begin{equation}
\phi(H_j,V_i) = f(j,\delta_{il},V_i) \cdot g(V_i) = [1-e^{-p(V_i+q)}]^{r(j - \delta_{il})} \cdot e^{-aV_i^b}
\end{equation}
Where $H_j$ represents the host with immunity $j$ and a known infection history such that $H_j : {j, \delta_{il}}$.
4. Within host reproductive number
The within-host reproductive number is given by the product of the probability of successful within host infection, and the number of offspring virions produced per event:
\begin{equation}
R_{in} = n \cdot \phi(H_j,V_i)
\end{equation}
5. Infectiousness
The infectiousness of a particular virus is therefore related to the within-host reproductive number and the number of initially infecting virions as follows:
\begin{equation}
\rho = 1 - (\frac{1}{R_{in}})^{-v} = 1-(\frac{1}{\phi(H_j, V_i)})^{-nv}
\end{equation}
Plot E below shows the probability of infection between hosts as a function of binding avidity.
B + ggtitle("Plot E")
6. Transmission rate
The transmission rate between hosts, $\beta$, is therefore given by the product of the infectiousness of that virus and the contact rate between hosts. This relationship is shown in plot E. The rate of change of $\beta$ with respect to binding avidity is shown in plot F.
\begin{equation}
\beta = c \cdot \rho
\end{equation}
D + ggtitle("Plot F")
7. Derivatives
The rate of binding avidity adaptation of a virus within a given host is assumed to be proportional to the derivate of $R_{in}$ with respect to binding avidity:
\begin{equation}
\frac{dV}{dt} = \frac{dR_{in}}{dV} = g(V_i)f'(J,V_i) + g'(V_i)f(J,V_i)
\end{equation}
Where $J = j - \delta_{il}$, and the derivaties for $g(V)$ and $f(J, V)$ are given by:
\begin{equation}
g'(V) = -abV^{b-1}e^{-a(V^b)}\
f'(J, V) = prJ(1-e^{-p(V+q)})^{rJ-1}e^{-p(V + q)}
\end{equation}
Parameter descriptions
- $p$: parameter to control degree by which changes in binding avidity affect probability of escape from immune response
- $r$: parameter to control degree by which previous exposure reduce probability of immune escape
- $b$: parameter to control the shape of the relationship between probability of successful replication and changes in binding avidity
- $a$: controls rate of changes of relationship between probability of successful replication and change in binding avidity
- $c$: per day contact rate
- $n$: number of offspring per virus replication event
- $v$: number of virions initially transmitted
- $q$: parameter to control the shape of the relationship between binding avidity and immune escape (shift on the x-axis)
- $\delta$: antigenic distance between two viruses
jameshay218/driftSim documentation built on Jan. 24, 2020, 6:14 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
library(driftSim) library(ggplot2) attach(exampleParameters) knitr::opts_chunk$set(echo = FALSE, fig.width=7,fig.height=6, message=FALSE, results=FALSE,warning=FALSE) p = viruspars["p"] #' parameter to control degree by which changes in binding avidity affect probability of escape from immune response r = viruspars["r"] #' parameter to control degree by which previous exposure reduce probability of immune escape b = viruspars["b"] #' parameter to control the shape of the relationship between probability of successful replication and changes in binding avidity a =viruspars["a"] #' controls rate of changes of relationship between probability of successful replication and change in binding avidity. n = viruspars["n"] #' average number of virus copies: n is number of offspring per virus replication v = viruspars["v"] #' v is number of virions initialyl transmitted nv = n*v q = viruspars["q"] #' parameter to control the shape of the relationship between binding avidity and immune escape (shift on the x-axis) V = seq(0, 3, by = 0.005) #' Binding avidity N_reinfect = as.numeric(25) max_reinfect = N_reinfect delta = as.numeric(5) #' Get a colour scale gradient blue to red f_array <- NULL #' Survival prob df_array <- NULL #' Derivative of survival prob for(k in 0:(N_reinfect-1)){ probTarget = exp(-p*(V+q)) #' probability of being targetted by immune system. As binding avidity increases, this probability decreases probEscape = 1-probTarget #' probability of escape immJ = r*(k- delta) #' Strength of host immune respone. As k increases, virus must escape more antibodies. As delta increases, this effect diminishes as infecting virus is further from host immunity. if(immJ < 0) immJ = 0 f = (probEscape)^(immJ) #' probability of escape from immunity if(k >= 1) f_dash= immJ*p*((1-probTarget)^(immJ-1))*(probTarget) #' derivative of this. ie. rate of change of relationship between binding avidity and probability of immune escape else f_dash= rep(0, length(V)) f_array[[k+1]] <- f df_array[[k+1]] <- f_dash } probInf <- exp(-a*(V^b)) #' probability of successful infection within a host. as binding avidity increases in a naive host, chance of successfully replicating decreases probInf_dash= -a*b*(V^(b-1))*exp(-a*(V^b)) #' rate of change of this relationship rho_array <- NULL dV_array <- NULL probRep_array <- NULL for(i in 1:length(df_array)){ R0 = f_array[[i]]*probInf*n rho = 1 - R0^-v rho[rho < 0] = 0 rho_array[[i]] <- rho dV = df_array[[i]]*probInf + f_array[[i]]*probInf_dash dV_array[[i]] <- dV probReplication = f_array[[i]]*probInf probRep_array[[i]] = probReplication } rho0 = max(rho_array[[1]]) rho1 = max(rho_array[[2]]) rho2 = max(rho_array[[3]]) probRep_data <- NULL for(i in 1:length(probRep_array)){ data <- data.frame(x=V,y=probRep_array[[i]],z=as.character(i)) probRep_data <- rbind(probRep_data,data) } colours <- NULL blue <- rev(seq(0,1,by=1/N_reinfect)) red <- seq(0,1,by=1/N_reinfect) for(i in 1:N_reinfect){ colours[i] <- rgb((i-1)/(max_reinfect-1),0,(N_reinfect-i)/(max_reinfect-1)) } A <- ggplot() + geom_line(data=probRep_data,aes(x=x,y=y,colour=z)) + scale_color_manual(values=colours) + theme(plot.title=element_text(hjust=0.5), panel.grid.major=element_blank(), panel.grid.minor=element_blank(), panel.background=element_blank(), axis.text.x=element_text(colour="gray20",size=12), axis.text.y = element_text(colour="gray20",size=12), text = element_text(colour="gray20",size=14), axis.line=element_line(colour="gray20"), axis.line.x =element_line(colour="gray20"), axis.line.y=element_line(colour="gray20"), legend.position = "none", axis.title=element_text(size=12) ) + scale_x_continuous(expand=c(0,0)) + scale_y_continuous(limits=c(0,1),expand=c(0,0)) + ylab("Probability of Successful\n Replication Within a Host (g(V))") + xlab("Binding Avidity") + ggtitle("Plot A") rho_data<- NULL for(i in 1:length(rho_array)){ data <- data.frame(x=V,y=rho_array[[i]],z=as.character(i)) rho_data<- rbind(rho_data,data) } B <- ggplot() + geom_line(data=rho_data,aes(x=x,y=y,colour=z)) + scale_color_manual(values=colours) + theme(plot.title=element_text(hjust=0.5), panel.grid.major=element_blank(), panel.grid.minor=element_blank(), panel.background=element_blank(), axis.text.x=element_text(colour="gray20",size=12), axis.text.y = element_text(colour="gray20",size=12), text = element_text(colour="gray20",size=14), axis.line=element_line(colour="gray20"), axis.line.x =element_line(colour="gray20"), axis.line.y=element_line(colour="gray20"), legend.position = "none", axis.title=element_text(size=12) ) + scale_x_continuous(expand=c(0,0)) + scale_y_continuous(limits=c(0,1),expand=c(0,0)) + ylab(expression(Probability~of~Infection~Between~Hosts~rho)) + xlab("Binding Avidity") + ggtitle("Plot B") + geom_hline(yintercept = rho0,linetype="longdash",colour="dodgerblue4")+ geom_hline(yintercept = rho1,linetype="longdash",colour="dodgerblue4")+ geom_hline(yintercept = rho2,linetype="longdash",colour="dodgerblue4") df_data<- NULL for(i in 1:length(df_array)){ data <- data.frame(x=V,y=df_array[[i]],z=as.character(i)) df_data<- rbind(df_data,data) } #' Derivative of f C <- ggplot() + geom_line(data=df_data,aes(x=x,y=y,colour=z)) + scale_color_manual(values=colours) + theme(plot.title=element_text(hjust=0.5), panel.grid.major=element_blank(), panel.grid.minor=element_blank(), panel.background=element_blank(), axis.text.x=element_text(colour="gray20",size=12), axis.text.y = element_text(colour="gray20",size=12), text = element_text(colour="gray20",size=14), axis.line=element_line(colour="gray20"), axis.line.x =element_line(colour="gray20"), axis.line.y=element_line(colour="gray20"), legend.position = "none", axis.title=element_text(size=12) ) + scale_x_continuous(expand=c(0,0)) + scale_y_continuous(expand=c(0,0)) + ylab(expression(d*f/d*V)) + xlab("Binding Avidity") + ggtitle("Plot C") dV_data<- NULL for(i in 1:length(dV_array)){ data <- data.frame(x=V,y=dV_array[[i]],z=as.character(i)) dV_data<- rbind(dV_data,data) } #' Infection D <- ggplot() + geom_line(data=dV_data,aes(x=x,y=y,colour=z)) + scale_color_manual(values=colours) + theme(plot.title=element_text(hjust=0.5), panel.grid.major=element_blank(), panel.grid.minor=element_blank(), panel.background=element_blank(), axis.text.x=element_text(colour="gray20",size=12), axis.text.y = element_text(colour="gray20",size=12), text = element_text(colour="gray20",size=12), axis.line=element_line(colour="gray20"), axis.line.x =element_line(colour="gray20"), axis.line.y=element_line(colour="gray20"), legend.position = "none", axis.title=element_text(size=12) ) + scale_x_continuous(expand=c(0,0)) + scale_y_continuous(expand=c(0,0)) + ylab(expression(d*beta/d*V)) + xlab("Binding Avidity") + geom_hline(yintercept=0,linetype='longdash',colour="gray20") + ggtitle("Plot D") f_data<- NULL for(i in 1:length(f_array)){ data <- data.frame(x=V,y=f_array[[i]],z=as.character(i)) f_data <- rbind(f_data,data) } E <- ggplot() + geom_line(data=f_data,aes(x=x,y=y,colour=z)) + scale_color_manual(values=colours) + theme(plot.title=element_text(hjust=0.5), panel.grid.major=element_blank(), panel.grid.minor=element_blank(), panel.background=element_blank(), axis.text.x=element_text(colour="gray20",size=12), axis.text.y = element_text(colour="gray20",size=12), text = element_text(colour="gray20",size=14), axis.line=element_line(colour="gray20"), axis.line.x =element_line(colour="gray20"), axis.line.y=element_line(colour="gray20"), legend.position = "none", axis.title=element_text(size=12) ) + scale_x_continuous(expand=c(0,0)) + scale_y_continuous(expand=c(0,0)) + ylab("Probability of Evading Immune System, f(j,V)") + xlab("Binding Avidity") + ggtitle("Plot E") probRep_dat <- data.frame(x=V,y=probInf) F <- ggplot() + geom_line(data=probRep_dat,aes(x=x,y=y,colour="red")) + theme(plot.title=element_text(hjust=0.5), panel.grid.major=element_blank(), panel.grid.minor=element_blank(), panel.background=element_blank(), axis.text.x=element_text(colour="gray20",size=12), axis.text.y = element_text(colour="gray20",size=12), text = element_text(colour="gray20",size=14), axis.line=element_line(colour="gray20"), axis.line.x =element_line(colour="gray20"), axis.line.y=element_line(colour="gray20"), legend.position = "none", axis.title=element_text(size=12) ) + scale_x_continuous(expand=c(0,0)) + scale_y_continuous(expand=c(0,0)) + ylab("Probability of Successful Replication, g(V)") + xlab("Binding Avidity") + ggtitle("Plot F") #' Plot antigenic distance against j at a given binding avidity. d <- seq(0,N_reinfect,by=0.1) fixedV <- 0.5 delta_array <- NULL for(k in 0:(N_reinfect-1)){ probT = exp(-p*(fixedV+q)) probE = 1- probT immJ1 = r*(k-d) immJ1[immJ1 < 0] <- 0 probSurvival1 = 1 - (n*exp(-a*(fixedV^b))*(probE^immJ1))^-v probSurvival1[probSurvival1 < 0] <- 0 delta_array[[k+1]] <- probSurvival1 } delta_data<- NULL for(i in 1:length(delta_array)){ data <- data.frame(x=d,y=delta_array[[i]],z=as.character(i)) delta_data <- rbind(delta_data,data) } G <- ggplot() + geom_line(data=delta_data,aes(x=x,y=y,colour=z)) + scale_color_manual(values=colours) + theme( plot.title=element_text(hjust=0.5), panel.grid.major=element_blank(), panel.grid.minor=element_blank(), panel.background=element_blank(), axis.text.x=element_text(colour="gray20",size=12), axis.text.y = element_text(colour="gray20",size=12), text = element_text(colour="gray20",size=14), axis.line=element_line(colour="gray20"), axis.line.x =element_line(colour="gray20"), axis.line.y=element_line(colour="gray20"), legend.position = "none", axis.title=element_text(size=12) ) + scale_x_continuous(expand=c(0,0)) + scale_y_continuous(limits=c(0,1),expand=c(0,0)) + ylab("Probability of Infection Between Hosts") + xlab("Antigenic Distance to Host Immunity") + ggtitle("Plot G") immJ_array <- NULL immJ2 <- seq(0,r*N_reinfect,by=1) for(k in 0:(N_reinfect-1)){ probT1 = exp(-p*(fixedV+q)) probE1 = 1- probT1 probSurvival2 = 1 - (n*exp(-a*(fixedV^b))*(probE1^immJ2))^-v probSurvival2[probSurvival2 < 0] <- 0 immJ_array[[k+1]] <- probSurvival2 } immJ_data<- NULL for(i in 1:length(immJ_array)){ data <- data.frame(x=immJ2/r,y=immJ_array[[i]],z=as.character(i)) immJ_data <- rbind(immJ_data,data) } H <- ggplot() + geom_line(data=immJ_data,aes(x=x,y=y,colour=z)) + scale_color_manual(values=colours) + theme(plot.title=element_text(hjust=0.5), panel.grid.major=element_blank(), panel.grid.minor=element_blank(), panel.background=element_blank(), axis.text.x=element_text(colour="gray20",size=12), axis.text.y = element_text(colour="gray20",size=12), text = element_text(colour="gray20",size=14), axis.line=element_line(colour="gray20"), axis.line.x =element_line(colour="gray20"), axis.line.y=element_line(colour="gray20"), legend.position = "none", axis.title=element_text(size=12) ) + scale_x_continuous(expand=c(0,0)) + scale_y_continuous(limits=c(0,1),expand=c(0,0)) + ylab("Probability of Infection Between Hosts") + xlab("immJ") + ggtitle("Plot H")
Foreword
The purpose of this vignette is to describe the theory behind the model of the driftSim package. Here, we will describe the model that governs the within-host dynamics of influenza binding avidity adaptation. Another way of understanding these equations is to check out the second (parameters) panel of the driftSim shiny app, which can be run with driftSim::runSimulationApp(). A useful piece of background reading can be found in a paper by Yuan et al.
Note that this document does not describe the transmission model itself, which is described here.
Model equations
The plots in the shiny app are governed by the equations shown below. These equations make up a model to represent how changes in binding avidity impact the various stages from host infection through to reinfection. The model also considers the immune history of the host, $\textbf{h}$, and the level of immunity, $j$, a host has, which is represented by the blue-red lines in the plots (blue $j=0$, red $j=j_{max}$). For the purposes of these plots, a higher value of $j$ simply represents higher levels of antibody-mediated immunity.
1. Probability of evading immune system
After entering a host, the virus must first evade the immune response (Equation 1). Here, the probability of escaping the immune response increases as binding avidity increases as shown in plot A. The virus must also evade the host's antibody-mediated immunity conferred from all previous infections ($\textbf{h}$), adjusted by the antigenic distance between the infecting virus and the viruses that elicited the host's previous antibody response.
E + ggtitle("Plot A")
As hosts will have encountered a number of infections in their lifetime, the infection history $\textbf{h}$ for a given host is the vector of all viruses that the host has seen. This part of the model is described in greater detail in Section 3.3 of the accompanying transmission model vignette.
In summary, the level of immunity that a host has against an infecting virus is the sum of the host's antibody repetoire less the minimum antigenic distance between the infecting virus and all past infections in the host's immune history, represented by the term $r(j - \delta_{il})$ below. $i$ represents the infecting strain and $l$ represents the most antigenically similar strain in the host's infection history. The definition of antigenic distance is analagous to the concepts described by Smith et al. 2004
The rate of change of $f$ with respect to binding avidity is shown in plot B.
C + ggtitle("Plot B")
\begin{equation} f(k,\delta_{il},V_i) = [1-e^{-p(V_i+q)}]^{r(j - \delta_{il})} \end{equation}
Where $V_i$ is the binding avidity of virus $i$; $j$ is the total antibody titre of the host; $\delta_{il}$ is the antigenic distance between strain $i$ and the host immunity virus, $l$; and the remaining parameters are defined at the bottom of this document (or with ?exampleParameters from the package).
2. Probability of successful replication within host
Binding avidity also affects how well a virus is able to replicate within the host, as described by Equation 2. This relationship is shown in plot C. The naive case (ie. $j=0$) is shown in plot D.
\begin{equation} g(V_i) = e^{-aV_i^b} \end{equation}
cowplot::plot_grid(A + ggtitle("Plot C"), F + ggtitle("Plot D"),ncol=2)
3. Probability of successful within host infection
To successfully infect a host, the virus must escape the immune system (Equation 1) and successfully replicate (Equation 2). The probability of a successful within-host infection is given by: \begin{equation} \phi(H_j,V_i) = f(j,\delta_{il},V_i) \cdot g(V_i) = [1-e^{-p(V_i+q)}]^{r(j - \delta_{il})} \cdot e^{-aV_i^b} \end{equation}
Where $H_j$ represents the host with immunity $j$ and a known infection history such that $H_j : {j, \delta_{il}}$.
4. Within host reproductive number
The within-host reproductive number is given by the product of the probability of successful within host infection, and the number of offspring virions produced per event: \begin{equation} R_{in} = n \cdot \phi(H_j,V_i) \end{equation}
5. Infectiousness
The infectiousness of a particular virus is therefore related to the within-host reproductive number and the number of initially infecting virions as follows: \begin{equation} \rho = 1 - (\frac{1}{R_{in}})^{-v} = 1-(\frac{1}{\phi(H_j, V_i)})^{-nv} \end{equation}
Plot E below shows the probability of infection between hosts as a function of binding avidity.
B + ggtitle("Plot E")
6. Transmission rate
The transmission rate between hosts, $\beta$, is therefore given by the product of the infectiousness of that virus and the contact rate between hosts. This relationship is shown in plot E. The rate of change of $\beta$ with respect to binding avidity is shown in plot F. \begin{equation} \beta = c \cdot \rho \end{equation}
D + ggtitle("Plot F")
7. Derivatives
The rate of binding avidity adaptation of a virus within a given host is assumed to be proportional to the derivate of $R_{in}$ with respect to binding avidity:
\begin{equation} \frac{dV}{dt} = \frac{dR_{in}}{dV} = g(V_i)f'(J,V_i) + g'(V_i)f(J,V_i) \end{equation}
Where $J = j - \delta_{il}$, and the derivaties for $g(V)$ and $f(J, V)$ are given by:
\begin{equation} g'(V) = -abV^{b-1}e^{-a(V^b)}\ f'(J, V) = prJ(1-e^{-p(V+q)})^{rJ-1}e^{-p(V + q)} \end{equation}
Parameter descriptions
- $p$: parameter to control degree by which changes in binding avidity affect probability of escape from immune response
- $r$: parameter to control degree by which previous exposure reduce probability of immune escape
- $b$: parameter to control the shape of the relationship between probability of successful replication and changes in binding avidity
- $a$: controls rate of changes of relationship between probability of successful replication and change in binding avidity
- $c$: per day contact rate
- $n$: number of offspring per virus replication event
- $v$: number of virions initially transmitted
- $q$: parameter to control the shape of the relationship between binding avidity and immune escape (shift on the x-axis)
- $\delta$: antigenic distance between two viruses
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.