Simulating trees from the fossilized birth-death process

sim.fbd functions

FossilSim can be used to simulate trees directly from a fossilized-birth death process using the sim.fbd family of functions, whose interface is similar to the sim.bd functions in TreeSim. These functions return a list of SAtree objects, since the simulated trees may contain sampled ancestors represented as tips with zero-length branches.

Simulating trees from a time-homogeneous fossilized-birth death process

As in TreeSim, trees can be simulated while conditioning either on the number of sampled extant taxa (using the sim.fbd.taxa function), or on the total age of the process (using the sim.fbd.age function). These functions require the fossil sampling rate parameter psi as an argument, but otherwise take the same arguments as the sim.bd.taxa and sim.bd.age functions in TreeSim.

trees <- sim.fbd.taxa(n = 10, numbsim = 10, lambda = 3, mu = 2, psi = 2)

plot(trees[[1]])

A stratigraphic range representation of the tree and fossil samples can be plotted using the rangeplot.asymmetric function, which assumes that all speciation events are asymmetric.

rangeplot.asymmetric(trees[[1]])

Note that when simulating a complete tree (complete = TRUE), tips in the resulting SAtree object represent extinction events, whereas they represent fossil sampling events in the tree without unsampled lineages (complete = FALSE). The complete field of resulting the SAtree object is set to the value of the complete function argument.

trees <- sim.fbd.taxa(n = 10, numbsim = 10, lambda = 3, mu = 2, psi = 2, complete = T)

print(trees[[1]]$complete)

## [1] TRUE

rangeplot.asymmetric(trees[[1]])

Simulating trees from an episodic fossilized-birth death process

Trees can also be simulated from a time-varying, piecewise constant fossilized birth-death process using the sim.fbd.rateshift.taxa function. Apart from a vector of interval specific psi parameters, this function also takes most of the same arguments as the corresponding sim.rateshift.taxa function in TreeSim. At the moment, this function assumes the sampling fraction of extant taxa (frac) is equal to 1.0 in all time intervals.

trees = sim.fbd.rateshift.taxa(n = 10, numbsim = 1, lambda = c(2,1), mu = c(0,0.3), psi = c(1,0.1), times = c(0,0.3))

rangeplot.asymmetric(trees[[1]])