Stabilizing selection of seasonal influenza receptor binding in populations with partial immunity

1. Foreword

This vignette describes the potential output files of the driftSim package. When running the function run_simulation, there is an option to specify 5 output files. The contents of these, in the order that they are specified in the vector, are given below. Furthermore, any included functions that use these output files are also described.

## Basic example of producing these host files
attach(exampleParameters)

sim_duration <- 365 ## Duration of simulation in days
version <- 4
flags <- as.logical(c(1,1,1,0,1,0,0,1))
scenario_descriptions(4) ## Which version of the simulation to run (1-4)
viruspars["kc"] <- 0.2
output <- run_simulation(flags=flags, 
                         hostpars=hostpars, 
                         viruspars=viruspars, 
                         deltaVMat=deltaVMat,
                         iniKs=NULL,
                         start=0,end=365,
                         input_files="",
                         output_files=c("SIR.csv","voutput_characteristics.csv","voutput_distances.csv","hosts.csv","hostKs.csv"),                                   VERBOSE=TRUE,scenario=4)

2. SIR dynamic outputs

• Default filename: “SIR.csv”
• Description: A data frame with 4 columns giving counts for the number of susceptible, infected, recovered individuals and incidence at each time step. Note that the data frame has no column names.
• File type: .csv file
• Usage: plot_SIR()

sir <- read.table("SIR.csv",header=FALSE,sep=",")
plot_SIR(sir)

3. Virus characteristics

• Default filename: “voutput1.csv”
• Description: A data frame with 15 columns containing the following information:

1. vid: the unique virus identifier
2. birth: the time of birth in the simulation (ie. time of infection)
3. death: the time of virus death (ie. time of recovery)
4. parentid: the unique virus identifier of the parent
5. bindingAvidityIni: the binding avidity at the time of initial infection
6. bindingAvidityFinal: the binding avidity at the time of death
7. infectionNo: if the infected host has had n infections previously, then this number is n+1. Note that this number is fixed at the time of infection
8. distToParent: antigenic distance to the parent virus
9. hostImmunity: the level of unmediated immunity of the host, \(k\) at the time of infection
10. hostInfections: the number of infections that the host has had at the end of the simulation
11. immJ: the effective immunity of the host against this virus at the time of infection
12. distHost: the smallest antigenic distance between this virus and the viruses of the host’s infection history at the time of infection
13. distRoot: the antigenic distance to the root virus at the end of the virus’s life
14. changeFromV: the amount of antigenic change elicited as a byproduct of binding avidity adaptation over the virus’s entire life
15. changeFromR: the amount of antigenic change elicited as from antigenic drift over the virus’s entire life

• File type: .csv file
• Usage: various. For example:

## Read in virus output data
tmp <- data.table::fread("voutput_characteristics.csv",data.table=TRUE)
dt <- data.table::data.table(tmp)

## Get mean and 95% CI of distance to root over time
means <- data.frame(dt[,list(mean=mean(distRoot)),by=birth])
lower <- data.frame(dt[,list(lower=quantile(distRoot,c(0.025))),by=birth])
upper <- data.frame(dt[,list(upper=quantile(distRoot,c(0.975))),by=birth])

## Also take a random virus that was extant at the
## end of the simulation, and use this to look at
## an example antigenic change profile
virus_data <- data.table::fread("voutput_characteristics.csv",data.table=FALSE)
## Note that extant viruses are given a death time of -1
samps <- virus_data[(virus_data$death == -1 & virus_data$birth > 0),]
i <- sample(nrow(samps),1)
startVid <- virus_data[virus_data$vid==samps[i,"vid"],"vid"]
index <- 1
nextVid <- 1

## Trace back through this virus lineage until we reach the
## root virus.
x <- y <- z <- NULL
while(nextVid != 0){
  nextVid <- virus_data[virus_data$vid==startVid,"parentid"]
  x[index] <- virus_data[virus_data$vid==startVid,"birth"]
  y[index] <- virus_data[virus_data$vid==startVid,"bindingAvidityFinal"]
  z[index] <- virus_data[virus_data$vid==startVid,"distRoot"]
  index <- index + 1
  startVid <- nextVid
}
dat <- data.frame(time = x, V=y ,delta=z )
dat1 <- dat[with(dat,order(time,delta)),]
dat1 <- dat1[dat1$time > 0,]

par(mfrow=c(2,1))
plot(means,type='l',lwd=1,col="blue", ylim=c(0,40),
     xlab="Time (days)",ylab="Average Antigenic Distance to Root",
     main="Mean and 95% CI of antigenic distances to root over time")
lines(lower,col="blue",lwd=0.5)
lines(upper,col="blue",lwd=0.5)

plot(dat1$delta~dat1$time,type='l',col="red",
     xlab="Time (days)", ylab="Antigenic distance from root",
     main="Example antigenic distance over time of an extant lineage")

4. Pairwise antigenic distances of all viruses

• Default filename: “voutput2.csv”
• Description: An n x n matrix giving the pairwise antigenic distance between any two viruses from the simulation. Saves only a random sample of 1000 viruses. It is possible to change the size of this sample in the source code, but it is currently hard coded. The matrix is in chronological order (ie. top left is most recent pair of viruses)
• File type: .csv file
• Usage: For (crude) example:

voutput2 <- data.table::fread("voutput_distances.csv",data.table=FALSE)
## An uninformative plot, but just to show an example use
image(as.matrix(voutput2))

5. Host population characteristics

• Default filename: “hosts.csv”
• Description: An data frame with four columns containing the following information on the entire living host population from the simulation:
  1. state: Infection status at the end of the simulation (0 = Susceptible; 1 = Infected; 2 = Recovered)
  2. last_vid: virus identifier of the last virus to infect the host. 0 if not applicable
  3. cur_inf: if the host was infected at the end of the simulation, this is the virus ID of the currently infecting virus
  4. hostK: the level of unmediated host immunity, \(j\) (referred to as \(k\) in the code), at the end of the simulation.
• File type: .csv file
• Usage: A key use is to generate a distribution of \(j\) values (referred to as \(k\) in the code) from which a new host population can be created and a new simulation seeded.

hostKDist <- generateHostKDist("hosts.csv",100000)
hist(hostKDist)

6. Host K distribution over time

• Default filename: “hostKs.csv”
• Description: A matrix containing the number of individuals in each \(j\) value (columns) (referred to as \(k\) in the code) over the course of the simulation. The frequency at which these are saved is specified by saveFreq in the hostpars argument for run_simulation. See ?exampleParameters.
• File type: .csv file
• Usage: Can be used to look at the distribution of unmediated host immunity, \(j\), over the course of the simulation.

hostKs <- read.table("hostKs.csv",sep=",")
day <- 25*hostpars["saveFreq"]

barplot(as.numeric(hostKs[25,1:30]),
        xlab="j",ylab="Frequency",
        main=paste0("Frequency of host j values at day ",day))