An animated GIF of Brownian evolution on a tree using phytools

Yesterday, I tweeted an animated GIF of Brownian motion evolution on a phylogeny.

Brownian motion on a phylogeny (a GIF), using #Rstats #phytools. pic.twitter.com/EYdvTGBzre

— Liam Revell (@phytools_liam) July 22, 2024

Here’s how this GIF was created (modified a bit so that the graphic stands vertically for this blog format, rather than horizontally).

## load phytools
library(phytools)
## set simulation parameters
N<-10 ## number of taxa
nlev<-200 ## number of steps
## simulate tree
tree<-pbtree(n=N,scale=1)
## simulate evolution on the tree 
mapped<-make.era.map(tree,
  limits=seq(0,1,length.out=nlev+1))
x<-fastBM(tt<-map.to.singleton(mapped),
  internal=TRUE)
## compute plot without plotting
plotTree(tt,plot=FALSE,ftype="off")
pp<-get("last_plot.phylo",envir=.PlotPhyloEnv)
## set tolerance
tol<-1e-6
## function to mix colors
mix_colors<-function(cols) colorRampPalette(cols)(3)[2]
## all colors
cols<-setNames(hcl.colors(n=Ntip(tree)),tree$tip.label)
## plotting loop
for(i in c(1:nlev,rep(nlev,20))){
  dev.hold()
  ## subdivide plotting area
  layout(matrix(c(1,2),2,1),heights=c(0.4,0.6))
  ## plot tree from 1 to i
  par(mar=c(1.1,3.1,2.1,1.1))
  plot(NA,xlim=c(0,1),ylim=c(range(pp$yy)),bty="n",
    axes=FALSE)
  ii<-intersect(which(pp$xx[pp$edge[,2]]>((i+tol-1)/nlev)),
    which(pp$xx[pp$edge[,2]]<=((i+tol)/nlev)))
  for(j in 1:nrow(pp$edge)){
    if(pp$xx[pp$edge[j,2]]<=(i+tol)/nlev){
      lines(pp$xx[pp$edge[j,]],pp$yy[pp$edge[j,]],type="s")
    }
  }
  ## add points at tips
  for(j in 1:length(ii)){
    daughters<-getDescendants(tt,pp$edge[ii[j],2])
    COL<-mix_colors(cols[daughters[daughters<=Ntip(tree)]])
    points(pp$xx[pp$edge[ii[j],2]],pp$yy[pp$edge[ii[j],2]],
      pch=16,col=COL,cex=1.5)
  }
  ## plot random walk from 1 to i
  par(mar=c(3.1,3.1,1.1,1.1))
  plot(NA,xlim=c(0,1),ylim=range(x),bty="n",las=1)
  for(j in 1:nrow(pp$edge)){
    if(pp$xx[pp$edge[j,2]]<=(i+tol)/nlev){
      lines(pp$xx[pp$edge[j,]],x[pp$edge[j,]])
    }
  }
  ## add points at tip
  for(j in 1:length(ii)){
    daughters<-getDescendants(tt,pp$edge[ii[j],2])
    COL<-mix_colors(cols[daughters[daughters<=Ntip(tree)]])
    points(pp$xx[pp$edge[ii[j],2]],x[pp$edge[ii[j],2]],
      pch=16,col=COL,cex=1.5)
  }
  dev.flush()
  Sys.sleep(0.05) ## sleep R
}

In an interactive R session, readers should run the above code with a plotting device that’s more taller than it is wide & a pretty cool animated visualization will come out!

To convert to a GIF, readers should have ImageMagick installed, and then can embed the main (outer) for loop of the code above within the code chunk below, and run it (adjust res, etc., as desired). Caution: it will normally be wise to run this within an otherwise empty directory because the last line cleans up by deleting all .png files in the working directory!!

png(file="bm-%03d.png",width=6,height=8,units="in",res=200)
## outer for loop here
dev.off()
system("magick convert -delay 10 -loop 0 *.png bm-anim.gif")
file.remove(list.files(pattern=".png"))

That’s it, folks!