Some time ago, I was studying an article from Tamara Münkemüller et al. 2015: "Phylogenetic niche conservatism – common pitfalls and ways forward". It´s a nice study where the ability of current tools in order to identify the evolutionary model that drives the changes in a trait (in this case, an ecological suitability) is discussed.
They simulate the evolution of a trait under several evolutionary processes and then, try to apply the tools that are commonly used to test what kind of process is behind the evolution of the trait, allowing us to infer several patterns or processes as phylogenetic niche conservatism. Their results are really interesting, and I recommend you to take a look to the whole research!
However, I would like to focus in one of the main pictures of the article:
Looks really nice! As a beginner, it was very difficult to me to understand this awesome figure. From my biological background, to study all this mathematical stuff requires a lot of time that I should invert in laboratory or field work... but it´s also very important to know the underlying statistical background of all the analysis that I would like to apply to my data in the future! So I decided to try to replicate some of those boxes by myself using R language. Many times, it´s the best way to learn!
The figure shows trait evolution over time under Ornstein–Uhlenbeck 1 parameter in the left column, Ornstein–Uhlenbeck evolutionary at middle column, and Ornstein–Uhlenbeck extended (multiple optimum) at right column, with different selection strength in each row. This means that the first row is actually a Brownian Motion process, since alpha (selection strength) is 0.
The graphic is the result of applying the OU process equation:
x(t+dt)= x(t)+ alpha(theta(t)-x(t))dt + (sigma dt E)
where , x(t) is the trait value at time=t, alpha is the selection strength with which a specific trait value (optimum) is wanted, theta is the optimal value of the trait in time=t, sigma is Brownian motion rate (the amount of random trait variability per time) and E is a normal distribution(0,1).
Following this equation, I developed a simple and intuitive R function that simulate an OU evolutionary process for a hypothetical character in two species: eMotion (from "evolutionary motion")
You can provide each parameter of the OU equation for the two species, and then you should define the number of generations (time) and replicates (number of simulations per species).
Here you have the code to load the function into R console:
eMotion <- function (sigma, theta, alpha, sigma2, theta2, alpha2, generations, n) {
plot(seq(1:generations),seq(1:generations), ylim=c(-200,200), type="n",
xlab="generations", ylab="trait value", main=paste("Sp1- s=", sigma,
" th=", theta, " a=", alpha," Sp2- s=", sigma2, " th=", theta2,
" a=", alpha2, sep=""))
trait_a=numeric()
simulations= list(1:n)
for (j in 1:n){
a=0
for (i in 1:generations) {
a= a + (sigma * rnorm(1)) + (alpha * ( theta - a))
trait_a[i]=a
}
aa=data.frame(trait_a)
simulations[j]=aa
}
for(k in 1:n){
bb=data.frame(simulations[k])
lines(bb, col="red")
}
#######
trait_a=numeric()
simulations= list(1:n)
for (j in 1:n){
a=0
for (i in 1:generations) {
a= a + (sigma2 * rnorm(1)) + (alpha2 * ( theta2 - a))
trait_a[i]=a
}
aa=data.frame(trait_a)
simulations[j]=aa
}
for(k in 1:n){
bb=data.frame(simulations[k])
lines(bb, col="blue")
}
}
Once that you have loaded the function, you can use it to create simple plots to observe the effect of each parameter of OU equation in resultant graphic. For instance, here we describe an OU process of a trait in two species. The species 1 (red) follows an OU process of character evolution, with a sigma value (Brownian motion rate) of 2 and an optimum in 100, while species 2 (blue) follows an OU model with a sigma value of 0.5 and an optimum in -100. In both cases, alpha is 0.002.
eMotion(sigma, theta, alpha, sigma2, theta2, alpha2, generations, n)
eMotion(2, 100, 0.002, 0.5,-100, 0.002, 10000, 20) 
The function automatically generates this simple plot, where you can observe that species 1 (red), is looking for a trait value around 100 (optimum), while species 2 (blue) is trying to get an optimum near to -100. Although both species have the same alpha (selection strength), the red one has a Brownian rate (sigma) much higher than the blue one. That explains why species 1 (red) has a higher range of variation than species 2 (blue), where all values are always near to the optimum (-100).
You can use the function to explore other options in parameter values, which could be useful in order to teach about evolutionary processes or to understand the Ornstein–Uhlenbeck equation. Here you have one more example, enjoy!
eMotion(1, 0, 0, 1,0, 0.01, 10000, 20) 