Threshold character evolution visualization along the branches of a tree

I'll be presenting on ancestral character estimation under the threshold model at the upcoming symposium entitled "Mathematics for an evolving biodiversity" in Montreal, Quebec in a couple of weeks. Last night, I thought it might be neat to take my existing function, bmPlot (which does discrete-time simulation & visualization under the threshold model), and modify it so that we could also visualize evolution of a threshold character up the branches of the tree.

The way I did this was by tracking the evolution of liability, translating liability into the threshold trait, and then aggregating adjacent generations with the same state into a stochastic-map style tree. Code for the new version of bmPlot is here.

One of the predictions of the threshold model is that the character should change back & forth across the threshold many times very rapidly when near the threshold. (Actually, infinitely many times in theory, under continuous time; remember, ours will be a discrete-time simulation.) That means that (compared to a 'nucleotide-model' type simulation) there some periods when the character changes very rapidly many times; and other long periods when the character doesn't change at all. Let's see how that bears out in simulation.

> tree<-pbtree(b=b,n=50,t=1000,type="discrete")
simulating with both taxa-stop (n) and time-stop (t) is
performed via rejection sampling & may be slow

51 trees rejected before finding a tree

> par(mfrow=c(2,1))
> par(mar=c(4.1,4.1,1.1,1.1))
> X<-bmPlot(tree,type="threshold",thresholds=c(0,0.5,2), anc=1,sig2=1/1000,ngen=max(nodeHeights(tree)), return.tree=TRUE)
> plotSimmap(X$tree,lwd=3,ftype="off",colors=X$colors)

So our prediction appears to be exactly borne out - large parts of the tree, including very long branches, have few to no changes; however in areas when a lineage is close to one of the thresholds, it tends to change back & forth rapidly. Cool.