How can we model the evolution of continuous traits?

This is not an easy task, but several approaches exist. We’ll focus on two models that have hugely influenced the field of evolutionary biology and remain important today: Brownian motion, and the Ornstein-Uhlenbeck process.

We will explore these related models on the basis of a tutorial by Sam Price, which I have slightly modified. All credit goes to Sam; you can find her original tutorial here, with updated functions at this location.

Simulation of Brownian motion

Now, what is Brownian motion (BM)? You may have heard about BM in the context of the random motion of particles suspended in liquid or gas. It’s basically a simple continuous stochastic (=probabilistic) process that describes the direction and magnitude of motion at any given point in time. Direction and magnitude of motion are completely independent from the previous history of motion.

We can apply this model to the evolution of a trait. For example, we can say that at each instant in time the trait value can increase or decrease, and the direction and the magnitude of change is not influenced by the previous history of the trait.

The total displacement of a trait after time t has elapsed is drawn from a normal distribtion with a mean of zero, so that the net change of the trait is 0). The variance in our trait is proportional to t, meaning the more time has passed the larger the accrued variance in observations.

Let’s investigate this by means of some simple, non-phylogenetic simulations.

We will first simulate trait displacement (=changes), using the rnorm() function to generate a normal distribution with a mean of zero for 1000 units of time.

changes <- rnorm(1000)
plot(changes, type="l", col="blue", ylab="Trait Change")

The cumulative sum of these changes represents the path of a trait through time — the resulting trajectory describes the trait history of a specific lineage through time.

sum.changes <- cumsum(changes)
plot(sum.changes, type="l", col="blue", xlab="Time", ylab="Trait Value")

Now proceed to generate a plot with 5 separate BM run — representative of the trait histories of 5 independent species sharing no common ancestry.

plot(cumsum(rnorm(1000)), type="l", ylim=c(-80,80), xlab="Time", ylab="Trait Value")
lines(cumsum(rnorm(1000)), col="red")
lines(cumsum(rnorm(1000)), col="blue")
lines(cumsum(rnorm(1000)), col="green")
lines(cumsum(rnorm(1000)), col="purple")

What patterns emerge? How does the overall variance change with time? Please describe in a short paragraph

What would happen if we were to change the parameters for the normal distribution and increased the standard deviation (sd)? The defualt is set to sd=1. Let’s change sd and find out.

#generate normal distributions
changes <- rnorm(1000) #as before with sd=1
heated.up <- rnorm(1000, sd=2) #new sd-value

#compare resulting BM simulations
plot(changes, type="l", col="blue", ylim=c(-6,6), ylab="Trait Change")
lines(heated.up, col="red")

plot(cumsum(changes), type="l", col="blue", ylim=c(-80,80), xlab="Time", ylab="Trait Value")
lines(cumsum(heated.up), col="red")

It appears the BM simulation run with the normal distribution created with larger SD shows larger displacements. The individual, instantaenous “jumps” are much bigger.

But what are the conseqences for the overall variance in the trait when we run these simulations for multiple lineages?

To explore this in more detail, let’s reduce the amount of typing and write a little function that summarizes our previous steps. Let’s start with a function that generates a normal distribution, calculates the cumulative sums, and then plots the results:

bm.sim0 <- function(sd){
  x <- rnorm(1000, sd=sd)
  sim <- cumsum(x)
  plot(sim, type="l", xlab="Time", ylab="Trait")

#BM simulation with sd=2