Visualizing a discrete-character-dependent multi-regime OU process on a tree with phytools

Working on something else, I was recently reminded of the fact that phytools includes a function to simulate a multi- \(\theta\) (i.e., multi-“regime”) Ornstein-Uhlenbeck process called multiOU.

multiOU does a continuous-time simulation, but can (optionally) output node states, so I thought it would be fun to show how the function can be used (in combination with lots of non-branching nodes in our tree) to graphically illustrate multi-optimum OU.

This doesn’t serve any special purpose outside the visual, but perhaps someone currently preparing to present at Evolution 2026 in Cleveland in a couple of weeks will find this useful for their slideshow!

## load phytools
library(phytools)

We can start by simulating a tree. Here I’m going to scale the total length of this tree to 100 time units.

tree<-pbtree(n=40,scale=100)

Next I’d like to add a bunch of non-branching nodes along the length of my tree from the root to every tip. I’m going to evenly space these & the way I’m going to do it is by wrapping a function call of phytools::make.era.map with another map.to.singleton. This first generates “era” mapped regimes on the tree in "simmap" format and then converts these to nodes.

tt<-map.to.singleton(make.era.map(tree,
  limits=seq(0,100,length.out=1001)))

We can see that this worked as follows.

plotTree(tt,ftype="off",lwd=1)
get("last_plot.phylo",envir=.PlotPhyloEnv)->pp
points(pp$xx,pp$yy,pch=16,col=make.transparent("blue",0.2),
  cex=0.5)

(I hope it’s evident that the blur of blue dots along each edge of the tree & non-branching nodes!)

Next we’re going to simulate our regimes on the tree as a continuous-time Markov chain using phytools::sim.history. I’m just going to simulate two regimes, a and b, with a switching rate between regimes of \(q=0.02\). Keep in mind, we should do this on our tree with unbranching nodes, tt.

## our Q matrix
Q<-matrix(c(-0.02,0.02,0.02,-0.02),2,2,
  dimnames=list(letters[1:2],letters[1:2]))
## our regime tree
s.tt<-sim.history(tt,Q,anc="a")

## Done simulation(s).

To get a sense of the history we’ve simulated, let’s plot it. This can be done in the standard way because phytools::plot.simmap has no problem with unbranching nodes.

cols<-setNames(hcl.colors(n=2),letters[1:2])
plot(s.tt,cols,ftype="off")

So far, so good.

Now, we’ll set our simulation conditions for the multi-regime OU. The function allows us to simulate different \(\alpha\), \(\sigma^2\), and \(\theta\) for each of our regimes, but here we’ll hold the former two parameters, \(\alpha\) & \(\sigma^2\), constant between regimes and set \(\theta = [-1,1]\) as follows.

theta<-setNames(c(-1,1),letters[1:2])
sigsq<-setNames(rep(0.4,2),letters[1:2])
alpha<-setNames(rep(0.7,2),letters[1:2])
x<-multiOU(s.tt,alpha=alpha,sig2=rep(sigsq,2),
  theta=theta,a0=theta["a"],internal=TRUE)

Note that the vector, x, that I simulated contains not only the states for the tips but also the conditions of all the thousands of internal nodes of the tree – the overwhelming majority of which are unbranching.

head(x)

##          t4         t25         t26          t8          t9         t18 
## -1.47750858  0.01917008  1.05348083  1.71793727  0.71903373  1.17666809

length(x)

## [1] 14574

Now, for our plot!

For fun, I’m also going to add the stationary distribution at the tips based on our generating values of \(\theta=[-1,1]\), \(\alpha=0.7\), and \(\sigma^2=0.4\).

layout(matrix(c(1,2),2,1),heights=c(0.4,0.6))
plot(s.tt,cols,ftype="off",mar=c(1.1,4.1,2.1,2.1),
  xlim=c(0,110))
mtext("a)",adj=0,line=0)
par(mar=c(5.1,4.1,2.1,2.1))
phenogram(s.tt,x,ftype="off",
  colors=make.transparent(cols,0.5),
  xlim=c(0,110),cex.axis=0.8,las=1)
mtext("b)",adj=0,line=1.5)
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)*10
  polygon(x=stat.norm+100,y=seq(min(x),max(x),
    length.out=200),border=FALSE,
    col=make.transparent(cols[i],0.5))
}

That’s pretty much the idea. Cool, right?