This repository contains R code that is used to generate the figures displayed in the paper
Schreiner CL, Nuismer SL, Basinski AJ. When to vaccinate a fluctuating wildlife population: Is timing everything? Journal of Applied Ecology. 2020 Feb;57(2):307-19.
Most directories contain a simulation script to generate the graphed data, and a graphing script to produce the EPS figures. Some figures do not require much simulation (e.g. Figure 1), and in these cases the simulation and graphing are combined into a single script. Data that appears in figures 6 and 8 are simulated using C++ - see the notes section for details on how to simulate these. The Mathematica file contains derivations and expressions for average population size, pathogen R0, etc, for our model.
Parameters for the simulations are specified in a dataframe called parmat. Time units are days. The parameter names are
- d - rate of host death from natural causes
- gamv - rate at which hosts transition from class Sv to V
- Nv - number of hosts exposed to vaccine during each campaign
- tv - the timing of annual vaccination, measured as days into the year
- tb - the duration of the birthing season
- NPeak - the peak population size, achieved once per year
- Rp - basic reproduction number of the pathogen
- Bp - transmission coefficient of the pathogen
- gamp - rate at which infected hosts transition out of infected class
- pmu - probability an infected host dies from infection (rather than recovers)
Many of these simulations were run using parallel computation. In these cases, we show effective cpu minutes (number of cores used x total simulation time)
- Figure 1: < 1 minute
- Figure 2: 154 cpu minutes
- Figure 3: 316 cpu minutes (when 25 tb values used. Figure generated with 365 values)
- Figure 4: 137 core minutes
- Figure 5: < 1 minute
- Figure 6: 16 core minutes for c++ simulations; 20 minutes for Graph.r
- Figure 7: < 1 minute
- Figure 8: 16 core minutes for c++ simulations; 20 minutes for Graph.r
- Table: 80 cpu minutes
Supplemental Figures are included in the code for Figure 2 (Fig S1), Figure 4 (Fig S2), and within Fig_6_DENS (a density dependent version of Figure 6).
-The birth rate "b" parameter is specified by choosing parameter values for NPeak, d, and tb.
-Similar to the note above, we specify the transmission coefficient "Bp" by specifying the value of R0, and then solving for the value of Bp that R0 implies.
These figures are generated by scripts written in C++. The "Sim.cpp" script simulates the outcome of vaccination. The script initializes S, Ip, and P state variables at the values that are predicted by the ODE system in the absence of vaccination, when simulated to a stable limit cycle. The stable limit cycle is calculated in the Graph_JustPath.r script within the Initial_Condition directory. Graph_JustPath.r prints the values of the state variables at time t%%365 = 0. These values are used in the pars.h header file to intialize the simulations at the stable limit cycle.
When compiled, the Sim.cpp script simulates the host population in the absence of vaccination until time [StartVaccinationTime] is exceeded. These initial years of simulation are used as baseline to evaluate the fractional reduction in mean number of pathogen infected hosts, relative to when vaccination is applied annually.
These scripts can be run by compiling Sim.cpp, then running the compiled code. The Graph.r script (in the Fig_6 and Fig_8 directories) reads in the simulated data, and produces a 2x2 figure that shows the reduction in mean abundance of pathogen-infected hosts, as well as the probability of elimination, with respect to the timing of vaccination.
Instead of running each script separately, we've also included a makefile that
- compiles and runs Sim.cpp
- runs Graph.r to produce the figure The makefile is used in emacs by opening it, then typing M-x compile. Enter the command make -k buildrun, and the simulations and graphs will be automatically generated.
Because the datasets associated with Gillespie simulations (Fig_6, Fig_8) were large, we compressed them into tar.xz files. You must first uncompress this data, or resimulate it, before running the graphing script. This is also true for the density-dependent version of Figure 6 (Fig_6_DENS).