10.4: Simulating Birth-Death Trees

We can use the statistical properties of birth-death models to simulate phylogenetic trees through time. We could begin with a single lineage at time 0. However, phylogenetic tree often start with the first speciation event in the clade, so one can also begin the simulation with two lineages at time 0 (this difference relates to the distinction between crown and stem ages of clades; see also Chapter 11).

To simulate our tree, we need to draw waiting times between speciation and extinction events, connect new lineages to the tree, and prune lineages when they go extinct. We also need a stopping criterion, which can have to do with a particular number of taxa or a fixed time interval. We will consider the latter, and leave growing trees to a fixed number of taxa as an exercise for the reader. Our simulation algorithm is as follows.

simulation algorithm

Assume that we have a certain number of “living” lineages in our tree (1 or 2 initially), a current time (tc = 0 initially), and a stopping time tstop.

1. Draw a waiting time ti to the next speciation or extinction event. Waiting times are drawn from an exponential distribution with rate parameter Nalive * (λ + μ) where Nalive is the current number of living lineages in the tree.
2. Check to see if the simulation ends before the next event. That is, if tc + ti > tstop, end the simulation.
3. Decide whether the next event is a speciation event [with probability λ/(λ + μ)] or an extinction event [with probability μ/(λ + μ)]. This can be done by drawing a uniform random number ui from the interval [0, 1] and assigning speciation to the event if ui < λ/(λ + μ) and extinction otherwise.
4. If (3) is a speciation event, then choose a random living lineage in the tree. Attach a new branch to the tree at this point, and add one new living lineage to the simulation. Return to step 1.
5. If (3) is an extinction event, choose a random living lineage in the tree. That lineage is now dead. As long as there is still at least one living lineage in the tree, return to (1); otherwise, your whole clade has gone extinct, and you can stop the simulation.

This procedure returns a phylogenetic tree that includes both living and dead lineages. One can prune out any extinct taxa to return a birth-death tree of survivors, which is more in line with what we typically study using extant species. It is also worth noting that entire clades can – and often do – go extinct under this protocol before one reaches time tstop. Note also that there is a much more efficient way to simulate trees (Stadler 2011).

We can think about phylogenetic predictions of birth-death models in two ways: by considering tree topology, and by considering tree branch lengths. I will consider each of these two aspects of trees below.