Visualizing a multi-regime OU process on the tree using phytools

In a recent tweet about phytools’ new discrete character dependent multi-optimum OU model, I included an illustration of this process that was not derived from the original blog post.

More on a discrete character dependent multi-optimum OU model in #phytools: https://t.co/LpLHskJoo6. pic.twitter.com/JrmAlmQpbh

— Liam Revell (@phytools_liam) July 13, 2026

Just in case it might become useful later, I decided to post the code for that cool illustration here. In contrast to my prior use, here I decided to set different values for \(\alpha\) & \(\sigma^2\), as well as for \(\theta\), between the three regimes.

## load packages
library(phytools)

In this first chunk I’m going to generate the data and objects I need for the plot.

## simulate a tree
tree<-pbtree(n=100,scale=10)
## add singleton (unbranching) nodes
tt<-map.to.singleton(make.era.map(tree,
  limits=seq(0,10,length.out=1001)))
## set the transition matrix for the discrete
## trait
q<-0.2
Q<-matrix(c(
  -2*q,q,q,
  q,-2*q,q,
  q,q,-2*q),3,3,
  dimnames=list(letters[1:3],letters[1:3]))
## generate a character history of the discrete
## trait
s.tt<-sim.history(tt,Q,anc="a")

## Done simulation(s).

## set the parameters of the multi-regime OU
## process
theta<-setNames(c(-0.5,0.9,2.2),letters[1:3])
sigsq<-setNames(c(0.12,0.05,0.08),letters[1:3])
alpha<-setNames(c(0.7,0.4,0.8),letters[1:3])
## generate data for tips (and nodes) under the
## multi-regime OU process
x<-multiOU(s.tt,alpha=alpha,sig2=sigsq,
  theta=theta,a0=theta["a"],internal=TRUE)

Now for our plot.

## set colors
cols<-setNames(hcl.colors(n=3),letters[1:3])
## set plot layout
layout(matrix(c(1,2),2,1),heights=c(0.4,0.6))
## plot discrete character history on the tree
plot(s.tt,cols,ftype="off",mar=c(1.1,4.1,2.1,2.1),
  xlim=c(0,11))
mtext("a)",adj=0,line=0)
## set margins for second subplot
par(mar=c(5.1,4.1,2.1,2.1))
## plot continuous character history
phenogram(s.tt,x,ftype="off",
  colors=make.transparent(cols,0.5),
  xlim=c(0,11),cex.axis=0.8,las=1)
mtext("b)",adj=0,line=1.5)
## add the predicted equilibrium distribution
## at the tips
for(i in 1:length(theta)){
  stat_theta<-dnorm(seq(min(x),max(x),
    length.out=200),mean=theta[i],
    sd=sqrt(sigsq[i]/(2*alpha[i])))
  stat.norm<-stat_theta/max(stat_theta)
  polygon(x=c(0,stat.norm,0)+10,
    y=c(min(x),seq(min(x),max(x),
      length.out=200),max(x)),border=FALSE,
    col=make.transparent(cols[i],0.5))
}

Nothing to it!