3.6: Simulating Brownian motion on trees

To simulate Brownian motion evolution on trees, we use the three properties of the model described above. For each branch on the tree, we can draw from a normal distribution (for a single trait) or a multivariate normal distribution (for more than one trait) to determine the evolution that occurs on that branch. We can then add these evolutionary changes together to obtain character states at every node and tip of the tree.

I will illustrate one such simulation for the simple tree depicted in figure 3.4b. We first set the ancestral character state to be $\bar{z}(0)$, which will then be the expected value for all the nodes and tips in the tree. This tree has three branches, so we draw three values from normal distributions. These normal distributions have mean zero and variances that are given by the rate of evolution and the branch length of the tree, as stated in equation 3.1. Note that we are modeling changes on these branches, so even if $\bar{z}(0) \neq 0$ the values for changes on branches are drawn from a distribution with a mean of zero. In the case of the tree in Figure 3.1, x1 ∼ N(0, σ2t1). Similarly, x2 ∼ N(0, σ2t2) and x3 ∼ N(0, σ2t3). If I set σ2 = 1 for the purposes of this example, I might obtain x1 = −1.6, x2 = 0.1, and x3 = −0.3. These values represent the evolutionary changes that occur along branches in the simulation. To calculate trait values for species, we add: xa = θ + x1 + x2 = 0 − 1.6 + 0.1 = −1.5, and xb = θ + x1 + x3 = 0 − 1.6 + −0.3 = −1.9.

This simulation algorithm works fine but is actually more complicated than it needs to be, especially for large trees. We already know that xa and xb come from a multivariate normal distribution with known mean vector and variance-covariance matrix. So, as a simple alternative, we can simply draw a vector from this distribution, and our tip values will have exactly the same statistical properties as if they were simulated on a phylogenetic tree. These two methods for simulating character evolution on trees are exactly equivalent to one another. In our small example, this alternative is not too much simpler than just adding the branches - but in some circumstances it is much easier to draw from a multivariate normal distribution.