seedy (version 1.3)

simulateoutbreak: Simulate transmission and evolutionary dynamics

Description

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

Usage

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, ...)

Arguments

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 i

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.

Value

• 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.

Details

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.

Examples

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)