Scale-invariant topology and bursty branching of evolutionary trees emerge from niche construction

Abstract

Phylogenetic trees describe both the evolutionary process and community diversity. Recent work has established that they exhibit scale-invariant topology, which quantifies the fact that their branching lies in between the two extreme cases of balanced binary trees and maximally unbalanced ones. In addition, the backbones of phylogenetic trees exhibit bursts of diversification on all timescales. Here, we present a simple, coarse-grained statistical model of niche construction coupled to speciation. Finite-size scaling analysis of the dynamics shows that the resultant phylogenetic tree topology is scale-invariant due to a singularity arising from large niche construction fluctuations that follow extinction events. The same model recapitulates the bursty pattern of diversification in time. These results show how dynamical scaling laws of phylogenetic trees on long timescales can reflect the indelible imprint of the interplay between ecological and evolutionary processes.

Keywords: evolution; molecular phylogeny; niche construction; scaling laws.

Fig. 1.

(A) A balanced tree. All nodes have exactly two descendants. The topological relation $C\left(A\right)\sim A\ln A$ at large $A$. (B) A maximally unbalanced tree. For any node, only one of the two children continues branching. $C\sim A^2$ at large $A$. Actual phylogenetic trees have topology and scaling behavior in between the two extreme cases, as studied in ref. .

Fig. 2.

(A) Averaged $C\left(A\right)$ calculated for a typical tree generated by the Niche Inheritance Model, with ${\sigma }_n=0$. The dots scatter along a straight line in the linear-logarithmic scale, indicating $\overline{C}/A\sim \ln A$. (B) Averaged typical $C\left(A\right)$ calculated with ${\sigma }_n=2$. The scale is double logarithmic. Fitting the well-averaged region, $A<200$, to a power function $\overline{C}\sim A^{\eta }$ gives an exponent of $\eta =1.501$, with the 95% CI being $\left[1.496,1.505\right]$. The red line stands for the fitted function. (C) Dependence of averaged $C\left(A\right)$ on ${\sigma }_n$, with $r_{ϵ}=0$. As ${\sigma }_n$ increases, the apparent power-law region of $\overline{C}\left(A\right)$ also stretches. Since subtrees with $A>10^4$ are undersampled, those data points are not shown in the main plot. The full data are shown in C, Inset. All of the following $C\left(A\right)$ plots are handled in a similar fashion. Other parameters for all sets of simulations are $r_{ϵ}=0$, ${\mu }_n=0$, $R_0=10$, and $n_0=1$ for the root node.

Fig. 3.

(A) Dependence of averaged $C\left(A\right)$ on $r_{ϵ}$, with ${\sigma }_n=2.5$. As $r_{ϵ}$ approaches zero, $\overline{C}\left(A\right)$ expands the power-law range. Other parameters are ${\mu }_n=0$, $R_0=10$, and $n_0=1$ for the root node. (B) The same data as in A plotted in the linear-logarithmic scale. The scaling turns $\overline{C}\left(A\right)\sim A\ln A$, as $r_{ϵ}$ becomes large. Insets in A and B show the full data range, respectively.

Fig. 4.

Critical scaling of $\overline{C}\left(A\right)$ as $r_{ϵ}$ decreases, indicated by the data collapse. By tuning $\eta$ and $b$, we reach a data collapse of the eight $\overline{C}\left(A\right)$ datasets corresponding to eight values of $r_{ϵ}$ ranging from 0 to 0.1 across four orders of magnitude. $b=0.36$ and $\eta =1.53$ gives the best result. The tail matches the expected $x^{1-\eta }\ln x$ behavior, which is indicated by the straight reference line in the double-logarithmic scale. Other parameters for all datasets are ${\sigma }_n=2.5$, ${\mu }_n=0$, $R_0=10$, and $n_0=1$ for the root node.

Fig. 5.

Dependence of EAD on ${\sigma }_n$, with $r_{ϵ}=0$. Data are shifted for clarity. Increasing ${\sigma }_n$ changes the power-law exponent approximately from $\alpha =2$ to $\alpha =1$. Since large subtrees are not well sampled, data beyond $k>10^3$ become noisy and, thus, are not informative. Inset shows the EADs plotted with full data. Other parameters used are ${\mu }_n=0$, $R_0=10$, and $n_0=1$.

Fig. 6.

(A) Averaged $C\left(A\right)$ for the modified model with a constant waiting time. Theoretical predictions from mean field analysis in SI Appendix for small niche construction strengths are also plotted. None shows power-law behavior. (B) Averaged $C\left(A\right)$ for the modified model with a constant bifurcation rate. Mean-field analysis no longer applies to this case, and power-law behavior is still not recovered. For both plots, other parameters used are $r_{ϵ}=0$, ${\mu }_n=0$, $R_0=100$, and $n_0=1$. The reference line is the power-law function $C=A^{1.51}$. The simulated tree has $10^6$ nodes, and the data are averaged over 10 repetitions. The Insets in A and B show the full data range, respectively.