I was discussing this project with Avery Kruger from UBC (I think he is @avery-kruger on Github) at the CSEE meeting this week in Toronto, and he asked if it was possible to fit an early burst evolutionary model in MRFtools. I hadn't encountered the model before, but I think, after doing some reading on this (my main reference was 'Early Burst Models' at LibreTexts Biology by Luke Harmon), that it should be possible to implement it relatively easily both as a "tips only" model and as a model across a full tree (including internal nodes). I'm just using this as a placeholder for my notes until I can add this to the code.
Based on Harmon's notes, the model is basically a Brownian motion model with exponentially transformed of branch lengths. First define the transformation function $L$ as:
$L(t,b) = b^{-1}(1-e^{-bt})$
Here $b$ is assumed to be a positive constant $b\ge 0$; in Harmon's notes he states that b should be negative, but I find it easier to think of it as an exponential decline function where the rate of decline increases with b, so I re-arranged the signs of the function. It also makes it clear that this is the cumulative distribution function for an exponential distribution scaled by $b^{-1}$.
Note that, in the limit as $b\to 0$, this just reduces to $L(t) = t$, so standard branch lengths are the limiting behaviour for this model (I.e. it will reduce to a Brownian motion model).
This function tells us how long a branch is that starts at time 0 and goes until time t. However, if we have internal nodes, they will usually occur at some time after time 0. As such, we can define a branch length from some node $i$ to one of its descendants $j$ as:
$$l_{ij} = L(t_j,b) - L(t_i,b)$$
$$
= e^{-b(t_i)} L(t_j-t_i, b)
$$
Once you have the transformed lengths between all pairs of connected nodes (I.e. internal nodes + their direct descendants), then you can calculate the precision matrix $Q$ as:
$Q_{ij} = Q_{ji} = -1/l_{ij}$ if $i$ is an ancestor of $j$, 0 if not
$Q_{ii} = - \sum_{j \ne i} Q_{ij}$
I.e. the standard Brownian motion model of evolution across the new branch lengths.
I also think that this implies the possibility of other varying rate models that could be implemented, as any $L$ function that starts at 0 at time zero, and increases monotonically with time should be a viable model. I could see something like a hinge model could be used to model evolution that starts quickly by slows to a lower (but still positive) rate with increasing branch lengths.
I was discussing this project with Avery Kruger from UBC (I think he is @avery-kruger on Github) at the CSEE meeting this week in Toronto, and he asked if it was possible to fit an early burst evolutionary model in
MRFtools. I hadn't encountered the model before, but I think, after doing some reading on this (my main reference was 'Early Burst Models' at LibreTexts Biology by Luke Harmon), that it should be possible to implement it relatively easily both as a "tips only" model and as a model across a full tree (including internal nodes). I'm just using this as a placeholder for my notes until I can add this to the code.Based on Harmon's notes, the model is basically a Brownian motion model with exponentially transformed of branch lengths. First define the transformation function$L$ as:
Here$b$ is assumed to be a positive constant $b\ge 0$ ; in Harmon's notes he states that b should be negative, but I find it easier to think of it as an exponential decline function where the rate of decline increases with b, so I re-arranged the signs of the function. It also makes it clear that this is the cumulative distribution function for an exponential distribution scaled by $b^{-1}$ .
Note that, in the limit as$b\to 0$ , this just reduces to $L(t) = t$ , so standard branch lengths are the limiting behaviour for this model (I.e. it will reduce to a Brownian motion model).
This function tells us how long a branch is that starts at time 0 and goes until time t. However, if we have internal nodes, they will usually occur at some time after time 0. As such, we can define a branch length from some node$i$ to one of its descendants $j$ as:
Once you have the transformed lengths between all pairs of connected nodes (I.e. internal nodes + their direct descendants), then you can calculate the precision matrix$Q$ as:
I.e. the standard Brownian motion model of evolution across the new branch lengths.
I also think that this implies the possibility of other varying rate models that could be implemented, as any$L$ function that starts at 0 at time zero, and increases monotonically with time should be a viable model. I could see something like a hinge model could be used to model evolution that starts quickly by slows to a lower (but still positive) rate with increasing branch lengths.