simulateoutbreak: Simulate transmission and evolutionary dynamics
In seedy: Simulation of Evolutionary and Epidemiological Dynamics

Description Usage Arguments Details Value Examples

Simulate within-host evolutionary dynamics on top of an SIR transmission process and generate genomic samples.

simulateoutbreak(init.sus, inf.rate, rem.rate, mut.rate, nmat = NULL, 
                 equi.pop = 10000, shape=flat, init.inf = 1, inoc.size = 1, 
                 samples.per.time = 1, samp.schedule = "random", 
                 samp.freq = 500, full=FALSE, mincases = 1, 
                 feedback = 500, glen = 1e+05, ref.strain = NULL, ...)

`init.sus`	Initial number of susceptible individuals.
`inf.rate`	SIR rate of infection.
`rem.rate`	SIR rate of removal.
`mut.rate`	Mutation rate (per genome per generation).
`nmat`	Connectivity matrix. Entry [i,j] gives the relative rate at which person i may contact person j.
`equi.pop`	Equilibrium effective population size of pathogens within-host.
`shape`	Function describing the population growth dynamics. See Details.
`init.inf`	Initial number of infected individuals.
`inoc.size`	Size of pathogen inoculum at transmission.
`samples.per.time`	Number of samples taken at sampling times.
`samp.schedule`	How should sampling be conducted? Accepted values are: "calendar" - samples are taken from all current infectives every `samp.freq` generations; "individual" - samples are taken from each infective at intervals of `samp.freq` after infection; "random" - samples are taken at one time between infection and removal for each infective.
`samp.freq`	Number of generations between each sampling time (see `samp.schedule`).
`full`	Should ‘full’ genomic sampling be returned? That is, should a vector of genotypes and their respective frequencies be stored from each individual's sampling times?
`mincases`	Minimum final size of outbreak to output. If final size is less than this value, another outbreak is simulated.
`feedback`	Number of generations between simulation updates returned to R interface.
`glen`	Length of genome.
`ref.strain`	Initial sequence. By default, a random sequence of length `glen`.
`...`	Additional arguments to be passed to the `shape` function.

Population growth dynamics are defined by the function called by 'shape'. This function returns the expected population size at each time step, given the total simulation time. By default, the population is expected to grow exponentially until reaching an equilibrium level, specified by equi.pop (flat). Alternatively, the population can follow a sinusoidal growth curve, peaking at runtime/2 (hump). User-defined functions should be of the form function(time,span,equi.pop,...), where span is equal to the duration of infection in this setting.

Returns a list of outbreak data:

`epidata`	A matrix of epidemiological data with columns: person ID, infection time, removal time, source of infection.
`sampledata`	A matrix of genome samples with columns: person ID, sampling time, genome ID.
`libr`	A list with an entry for each unique genotype observed. Each entry is a vector of mutation positions relative to the reference genome.
`nuc`	A list with an entry for each unique genotype observed. Each entry is a vector of nucleotide types (integer between 1 and 4).
`librstrains`	A vector of unique genotype IDs corresponding to the `libr` object.
`endtime`	End time of the outbreak.

W <- simulateoutbreak(init.sus=10, inf.rate=0.002, rem.rate=0.001, mut.rate=0.0001, 
                      equi.pop=2000, inoc.size=1, samples.per.time=10, 
                      samp.schedule="calendar", samp.freq=500, mincases=3) 
sampledata <- W$sampledata
epidata <- W$epidata

# Calculate distance matrix for observed samples
distmat <- gd(sampledata[,3], W$libr, W$nuc, W$librstrains)

# Now pick colors for sampled isolates
colvec <- rainbow(1200)[1:1000] # Color palette
refnode <- 1 # Compare distance to which isolate?
colv <- NULL # Vector of colors for samples
maxD <- max(distmat[,refnode])

for (i in 1:nrow(sampledata)) {
  	colv <- c(colv, 
              colvec[floor((length(colvec)-1)*(distmat[refnode,i])/maxD)+1])
}


plotoutbreak(epidata, sampledata, col=colv, stack=TRUE, arr.len=0.1, 
             blockheight=0.5, hspace=500, label.pos="left", block.col="grey", 
             jitter=0.004, xlab="Time", pch=1)