# Simulating BM evolution. 
library(ape)

# First of all you need to randomly draw from a normal distribution with a mean of zero you can use the function rnorm

displace<-rnorm(100, mean=0, sd=1)
plot(displace)

# The sum of these displacements represents the trajectory of the trait through time

x<-cumsum(displace)
plot(x, type="l")

# Now create a plot with 5 separate runs of Brownian motion. You can imagine these as 5 independent species (sharing no common ancestry).

plot(cumsum(rnorm(100)), type="l", ylim=c(-30,30), ylab="Trait")
lines(cumsum(rnorm(100)), col="red")
lines(cumsum(rnorm(100)), col="blue")
lines(cumsum(rnorm(100)), col="green")
lines(cumsum(rnorm(100)), col="purple")

# If we do this 50 times we can see a very important fact about Brownian motion and that is that variance is proportional to time.
plot(cumsum(rnorm(100)), type="l", ylim=c(-30,30), ylab="Trait")
for(i in 1:50) lines(cumsum(rnorm(100)))  

# The standard deviation is proportional to the rate so if we increase the rates (sd) we increase the disparity generated.

# Obviously the above simulations were assuming the species were independent but we know that they are not, for example species within the same genus share evolutionary history and this affects the evolution of their traits.

# First lets make a simple 4 taxon tree

tree<-read.tree(text="(((species1:4.2, species2:4.2):3.1, species3:7.3):6.3,species4:13.5);")
plot(tree)
axisPhylo()

# Looking at this tree we see that species 1,2 & 3 share about 6 million years of history during which they all follow the same BM trajectory after which species 3 splits from 1 and 2 with both lineages starting on a new BM trajectory - we can simulate this. We need the branch length information: the longer the branch the greater the expected variance - with this we determine the standard deviation of the normal distribution.

# The total displacement of a character after time v is drawn from a normal distribtion with a mean of zero (i.e. net change in trait is 0) and a variance proportional to v (branch length). 

# Start with the first branch from the root node 5 to 6 and estimate the trait value at node 6. A root state of 1 and a BM rate of 0.5 (often called sigma-squared) use rnorm and the branch lengths to estimate the displacement.
x<-rnorm(1, mean=0, sd=sqrt(0.5*tree$edge.length[1]))

# You obviously need to convert this displacement to a trajectory, so if the root state is 1 then you add the displacement to the root state to get the value at the node.

node6est<-1+x

# You then repeat this process across the tree calculating the displacement with sqrt(0.5*tree$edge.length[i]) and adding it to the ancestral node. 

# Fortunately you don't need to do it by hand or work out how to automate it as there are several methods available in R: simchar in geiger and simulate in OUCH these two will simulate using your existing trait matrices to replicate their variance and covariance and for simulating traits with a given rate/variance you can use rTraitCont in the ape package or fastBM in the phytools package.

# Code for simulating Brownian motion with a tree.

speciation=rbinom(1000,1,0.01) #for each of the 1000 time steps determine if a speciation event has occured by randomly sampling from a binomial distribution.

species=1+cumsum(speciation) # total number of speciation events

bm.var=1 # rate of Brownian motion
time=(1:1000)

states=matrix(NA,1000,species[1000])
states[1,1]=0

plot(time[1],states[1,1],xlim=c(0,1000),ylim=c(-80,80),pch=".", cex = 2.5, xlab="", ylab="") # plots the point in the first time of the 1000 time bins  i.e. root state
mtext("Time", side=1, line=3, cex=2)
mtext("Trait Value", side=2, line=2, cex=2)

for (i in 1:999){
	if (speciation[i]==1){
		a=sample(species[i]-1,1)
		print(a)
		states[i,species[i]]=states[i,a]
		}
	states[1+i,1:species[i]]=states[i,1:species[i]]+((rbinom(species[i],1,0.5)-0.5)*bm.var)
	
	for (r in 1:species[i]){
		#points(time[i],states[i,r],pch=".",col=r, cex = 2.5)
		lines(c(time[i-1],time[i]),c(states[i-1,r],states[i,r]),pch=".",col=r, cex = 2.5)
		}
		#slow down simulation so you can see it happening
	for (z in 1:(10^4)){}
	}
	
# Now try it with a faster rate - so change bm.var >1