Diversity and coevolutionary dynamics in high-dimensional phenotype spaces

TL;DR

This paper investigates long-term coevolutionary dynamics in high-dimensional phenotype spaces, revealing that diversity increases exponentially with dimensionality and dynamics tend to stabilize at saturation levels.

Contribution

It combines microevolutionary competition models with long-term coevolutionary analysis in high-dimensional spaces, revealing new patterns of diversity and stability.

Findings

Diversity at saturation increases exponentially with phenotype space dimension.

Coevolutionary dynamics are fast and non-stationary at intermediate diversity levels.

Community diversity tends to stabilize as it reaches saturation.

Abstract

We study macroevolutionary dynamics by extending microevolutionary competition models to long time scales. It has been shown that for a general class of competition models, gradual evolutionary change in continuous phenotypes (evolutionary dynamics) can be non-stationary and even chaotic when the dimension of the phenotype space in which the evolutionary dynamics unfold is high. It has also been shown that evolutionary diversification can occur along non-equilibrium trajectories in phenotype space. We combine these lines of thinking by studying long-term coevolutionary dynamics of emerging lineages in multi-dimensional phenotype spaces. We use a statistical approach to investigate the evolutionary dynamics of many different systems. We find: 1) for a given dimension of phenotype space, the coevolutionary dynamics tends to be fast and non-stationary for an intermediate number of coexisting lineages, but tends to stabilize as the evolving communities reach a saturation level of diversity; and 2) the amount of diversity at the saturation level increases rapidly (exponentially) with the dimension of phenotype space. These results have implications for theoretical perspectives on major macroevolutionary patterns such as adaptive radiation, long-term temporal patterns of phenotypic changes, and the evolution of diversity.

Figures

a). This type of complicated dynamics is characterized by average evolutionary speeds $v>10^{-2}$. In the same system, but now with a value of $m_{max}$ that lies above the saturation level $M_{\sigma,d}$, the diversity evolves to the saturation level, at which the community consists of ca. 30 coexisting phenotypic clusters (Fig. 2b). In this saturated state, the average evolutionary speed is much lower than $10^{-2}$, and the community exhibits much more stationary coevolutionary dynamics (that would eventually converge to a coevolutionary equilibrium). Moreover, the saturated community exhib

A. Over time the evolving community becomes more diverse due to adaptive diversification, and as a consequence the nature of the coevolutionary dynamics of the community changes. In the example shown in Figure 5A, the coevolutionary dynamics are fast for low to intermediate levels diversity, and then slow down as the community acquires more and more species, until eventually the community reaches a coevolutionary equilibrium at the diversity saturation level. Again, the slowdown of the evolutionary speed during an adaptive radiation appears to occur exponentially with an increase in diversity.

Equations

$$f(\mathbf{x},\mathbf{y})=1-\frac{\alpha (\mathbf{x},\mathbf{y})K(\mathbf{x})}{K(\mathbf{y})}.$$

$$\displaystyle s_{i}(\mathbf{x})\equiv\frac{\partial f(\mathbf{x},\mathbf{y})}{\partial y_{i}}{\big{|}_{\mathbf{y}=\mathbf{x}}}=-\frac{\partial\alpha(\mathbf{x},\mathbf{y})}{\partial y_{i}}{\big{|}_{\mathbf{y}=\mathbf{x}}}+\frac{\partial K(\mathbf{x})}{\partial x_{i}}\frac{1}{K(\mathbf{x})},$$

$$\frac{d\mathbf{x}}{dt}=\mathbf{M}(\mathbf{x})\cdot \mathbf{s}(\mathbf{x}).$$

$$\alpha (\mathbf{x},\mathbf{y})=\exp \left[\sum _{i,j=1}^d{b}_{ij}(x_i-y_i)x_j-\sum _{i=1}^d\frac{(x_i-y_i{)}^2}{2{\sigma }_i^2}\right].$$

$$K(\mathbf{x})=\exp \left(-\frac{{\sum }_i^d{x}_i^4}{4}\right).$$

$$\frac{d{x}_i}{dt}=\sum _{i=1}^d{b}_{ij}{x}_j-x_i^3,i=1,...,d.$$

$$\frac{d{N}_r(t)}{dt}=N_r(t)\left(\frac{1-{\sum }_{s=1}^m\alpha ({\mathbf{x}}_s,{\mathbf{x}}_r)N_s(t)}{K({\mathbf{x}}_r)}\right),r=1,...,m.$$

$$f({\mathbf{x}}_1,...,{\mathbf{x}}_m,{\mathbf{x}}_r^{\prime })=1-\frac{{\sum }_{s=1}^m\alpha ({\mathbf{x}}_s,{\mathbf{x}}_r^{\prime })N_s^{\ast }}{K({\mathbf{x}}_r^{\prime })}.$$

$$s_{ri}=\sum _s{N}_s^{\ast }\left(-\frac{1}{K({\mathbf{x}}_r)}{\frac{\partial \alpha ({\mathbf{x}}_s,{\mathbf{x}}_r^{\prime })}{\partial x_{ri}^{\prime }}|}_{{\mathbf{x}}_r^{\prime }={\mathbf{x}}_r}+\frac{\alpha ({\mathbf{x}}_s,{\mathbf{x}}_r)}{K^2({\mathbf{x}}_r)}\frac{\partial K({\mathbf{x}}_r)}{\partial x_{ri}}\right),i=1,...,d.$$

$$\frac{d{x}_{ri}}{dt}=N_r^{\ast }{s}_{ri},i=1,\dots ,d,r=1,\dots ,m.$$

$$v=\sum _{r=1}^m\frac{N_r\sqrt{{\sum }_{i=1}^d(d{x}_{ri}/dt{)}^2}}{{\sum }_{r=1}^m{N}_r}$$

$$C=\frac{⟨[Pg-⟨Pg⟩][X-⟨X⟩]⟩}{{\sigma }_{Pg}{\sigma }_X},$$

$$\frac{\partial N(x,t)}{\partial t}=N(x,t)\left(\frac{1-\int \alpha (y,x)N(y,t)dy}{K(x)}\right)+D\sum _{i=1}^d\frac{{\partial }^{2}N(x,t)}{\partial x_i^2},$$

$$M_{{\sigma }_a,d}\approx C_{\sigma }\frac{L^d}{{\sigma }^d}.$$

$$M_{{\sigma }_a,d}=M_{{\sigma }_b,d}{\left(\frac{{\sigma }_b}{{\sigma }_a}\right)}^d\text{and }{M}_{\sigma ,d_1}=M_{\sigma ,d_2}^{d_1/d_2},$$

$$\left(\begin{array}{ccc}0.407& 0.498& 0.287\\ -0.199& -1.102& -0.305\\ 1.387& -0.896& 0.341\end{array}\right)$$

$$\left(\begin{array}{cccc}-1.289& 0.682& 0.217& -0.093\\ -0.223& -0.035& 0.697& -0.117\\ -0.563& 0.434& -0.953& -0.198\\ 0.119& 0.398& 0.183& 0.530\end{array}\right)$$

Full text

Diversity and coevolutionary dynamics in high-dimensional phenotype spaces

Michael Doebeli∗ & Iaroslav Ispolatov*∗∗*

∗ Departments of Zoology and Mathematics, University of British Columbia,

6270 University Boulevard, Vancouver B.C. Canada, V6T 1Z4; doebeli$@$zoology.ubc.ca

∗∗ Departamento de Fisica, Universidad de Santiago de Chile

Casilla 302, Correo 2, Santiago, Chile; jaros007$@$gmail.com

Supporting Material: Appendix A

RH: Diversity and coevolutionary dynamics

Corresponding author: Michael Doebeli, Department of Zoology and Department of Mathematics, University of British Columbia, 6270 University Boulevard, Vancouver B.C. Canada, V6T 1Z4, Email: doebeli@zoology.ubc.ca.

Abstract We study macroevolutionary dynamics by extending microevolutionary competition models to long time scales. It has been shown that for a general class of competition models, gradual evolutionary change in continuous phenotypes (evolutionary dynamics) can be non-stationary and even chaotic when the dimension of the phenotype space in which the evolutionary dynamics unfold is high. It has also been shown that evolutionary diversification can occur along non-equilibrium trajectories in phenotype space. We combine these lines of thinking by studying long-term coevolutionary dynamics of emerging lineages in multi-dimensional phenotype spaces. We use a statistical approach to investigate the evolutionary dynamics of many different systems. We find: 1) for a given dimension of phenotype space, the coevolutionary dynamics tends to be fast and non-stationary for an intermediate number of coexisting lineages, but tends to stabilize as the evolving communities reach a saturation level of diversity; and 2) the amount of diversity at the saturation level increases rapidly (exponentially) with the dimension of phenotype space. These results have implications for theoretical perspectives on major macroevolutionary patterns such as adaptive radiation, long-term temporal patterns of phenotypic changes, and the evolution of diversity.

Keywords: Long-term evolution — Diversity and stability — Adaptive radiation

Introduction

One of the fundamental problems in evolutionary biology is to understand how microevolutionary processes generate macroevolutionary patterns. In particular, the emergence of macroevolutionary changes in the speed of evolution (Simpson, 1944, Gould and Eldredge, 1977), and of macroevolutionary changes in patterns of species diversity (Rosenzweig, 1995, Schluter, 2000) have long been of great interest. For example, Uyeda et al. (2011) have recently proposed that over macroevolutionary time scales, relatively short intermittent bursts of high rates of evolutionary change should alternate with long periods of bounded phenotypic fluctuations. Also, there is much discussion about whether species diversity saturates over evolutionary time in a given environment (Rosenzweig, 1995, Harmon and Harrison, 2015, Rabosky and Hurlbert, 2015). Phylogenetic analysis has been used to shed light on these questions (Pennell and Harmon, 2013), but mechanistic models in which short-term ecological interactions are extrapolated to yield long-term patterns of diversity and evolutionary change have only recently been developed. Most of these models have been used to study the long-term evolution of diversity by analyzing processes of community assembly emerging from short-term ecological dynamics (Loeuille and Loreau, 2005, 2009, Ito and Dieckmann, 2007, Allhoff et al., 2015, Rosindell et al., 2015, Gascuel et al., 2015). In particular, these papers have mainly focussed on how diversity changes over time, but not on how the nature of the coevolutionary dynamics of a given set of coexisting species changes as the diversity changes. In fact, in all these models, the evolutionary dynamics for a fixed amount of diversity, i.e., for a given set of species, converge to an equilibrium. However, if one wants to understand macroevolutionary changes in the “tempo and mode” (Simpson, 1944) of evolution, one not only needs to consider how diversity changes over evolutionary time, but also how such changes in diversity affect the nature of evolutionary dynamics (Pennell and Harmon, 2013). Indeed, there is evidence from evolution experiments with microbes that evolutionary dynamics in more diverse communities are qualitatively different from the evolutionary dynamics in less diverse communities (Lawrence et al., 2012). Here we present a theoretical investigation of the questions of how diversity affects the complexity of coevolutionary dynamics.

In general, the number of different phenotypes that affect ecological and evolutionary processes is an important quantity. For example, determining the dimensionality of niche space in ecological food webs is a classical problem (Cohen, 1977, Eklöf et al., 2013), and it has recently been shown that including more phenotypic dimensions in models for community assembly has a strong effect on the structure of the emerging food webs (Allhoff et al., 2015). Implicitly, the importance of the dimension of phenotype space is also acknowledged in phylogenetic research through the notion of “adaptive zones” (Simpson, 1944, Uyeda et al., 2011). In particular, it is thought that much of the extant diversity has evolved as a consequence of lineages entering new adaptive zones, which can be interpreted from the phenotypic perspective as an increase in the dimension of phenotype space. In general, given the large number of phenotypic properties that determine an individual’s life history and ecology in almost any species, one would expect that ecological interactions are generally determined by many phenotypic properties, and that selection pressures emerging from ecological interactions in turn affect many phenotypes simultaneously.

For example, comprehensive modelling of the metabolic network in E. coli cells comprises more than 2000 reactions (Yoon et al., 2012). These reactions are in turn controlled by thousands of genes in a complicated interaction network whose exact workings are largely unknown. Nevertheless, many of the genes contributing to this network of metabolic reactions will be under selection in any given environmental setting, and as a consequence, a large number of phenotypic properties have the potential to undergo evolutionary change. It is generally not known how exactly these phenotypic properties impinge on birth and death rates of individual organisms, and hence what exactly the ecological selection pressures are on these properties. Nevertheless, it seems clear that in general, many phenotypes will evolve at the same time, i.e., that evolution generally takes place in high-dimensional phenotype spaces.

We have recently argued that if evolution takes place in high-dimensional phenotype spaces, then the evolutionary dynamics, that is, the phenotypic change over evolutionary time, can be very complicated, i.e., non-stationary and often chaotic (Doebeli and Ispolatov, 2014, Ispolatov et al., 2015). In low-dimensional phenotype spaces, non-equilibrium evolutionary dynamics are less likely. However, if a species evolving on a simple attractor gives rise to diversification, the effective dimensionality of the evolving system increases, as the species that emerge from diversification coevolve, driven by both intra- and interspecific ecological interactions. Thus the total dimensionality of the resulting dynamical system describing multispecies coevolution is the number of species times the dimensionality of the phenotype space in which each species evolves. Based on our earlier results (Doebeli and Ispolatov, 2014, Ispolatov et al., 2015), one could then expect that due to the increase in dimensionality, diversification leads to more complicated evolutionary dynamics in each of the coevolving species. On the other hand, as a multispecies community becomes more diverse and evolves towards saturation, the available niches tend to get filled, and hence evolutionary change has to become highly coordinated between interacting species and thus constrained, potentially leading to simplified evolutionary dynamics. It is thus unclear how the nature of the evolutionary dynamics changes as the pattern of diversity changes during community assembly.

We investigate these issues by applying the framework of adaptive dynamics (Geritz et al., 1998) to a general class of competition models. The main question we address is, how does the complexity of long-term coevolutionary dynamics depend on the diversity of the coevolving community? We show that in low-dimensional phenotype spaces, there is a humped-shaped relationship between diversity and the complexity of evolutionary dynamics: in communities with low diversity, coevolutionary dynamics are often simple, i.e., stationary in the long-time limit; for intermediate degrees of diversity, non-stationary (complex) coevolutionary dynamics are common, and each of the species in the community evolves on a complicated trajectory in phenotype space; and for high amounts of diversity, coevolutionary dynamics become simple again, i.e., stationary. In particular, as communities reach diversity saturation, e.g. through adaptive diversification (Doebeli, 2011), coevolutionary dynamics change from complex to simple.

Our results are relevant for a number of issues concerning patterns of macroevolution. For example, the results suggest that during processes of adaptive radiation (Schluter, 2000, Gavrilets and Losos, 2009, Seehausen, 2015), evolutionary dynamics are more complicated early in the radiation than late in the radiation, a pattern that corresponds to the “early-burst” perspective of macroevolution that has attracted much attention in recent years (Gavrilets and Losos, 2009, Harmon et al., 2010, Slater and Pennell, 2013). Our results also show that the level at which diversity saturates depends on the dimensionality of phenotype space, with higher dimensions allowing for more diversity. This observation is in accordance with data from radiations in fishes (Seehausen, 2015) and points to the possibility of a microevolutionary mechanism for the “blunderbass theory” of temporal patterns of macroevolutionary changes and diversification (Uyeda et al., 2011): if evolution operates on the dimension of phenotype space on a very slow time scale, then on shorter time scales diversity may saturate and thereby generate relatively stationary evolutionary dynamics, whereas on longer time scales the dimension of phenotype space may increase, e.g. due to gene duplications, thus generating a new burst of non-equilibrium (co-)evolutionary dynamics until the diversity reaches a new saturation level. Such patterns of intermittent bursts have recently been found in the phylogenies of birds and echinoids (Brusatte et al., 2014, Hopkins and Smith, 2015), and the bursts have been attributed to the evolution of flight capabilities and of novel feeding techniques, respectively, both of which can be interpreted as an increase in the dimensionality of the relevant phenotype space. This perspective may also shed light on the question of whether diversity saturates or not (Harmon and Harrison, 2015, Rabosky and Hurlbert, 2015): diversity may saturate for a given dimension of phenotype space, but evolutionary innovation in the form of new phenotypic dimensions may intermittently generate room for additional bouts of evolutionary diversification.

Methods

Single-cluster adaptive dynamics

As in Doebeli and Ispolatov (2014), we study a general class of models for frequency-dependent competition in which ecological interactions are determined by $d$-dimensional phenotypes, where $d\geq 1$. For simplicity, we consider homogeneous systems, so no spatial coordinates are included. The ecological interactions are described by a competition kernel $\alpha(\mathbf{x},\mathbf{y})$ and by a carrying capacity $K(\mathbf{x})$, where $\mathbf{x},\mathbf{y}\in\mathbb{R}^{d}$ are the $d$-dimensional continuous phenotypes of competing individuals. The competition kernel $\alpha$ measures the competitive impact that an individual of phenotype $\mathbf{x}$ has on an individual of phenotype $\mathbf{y}$, and we assume that $\alpha(\mathbf{x},\mathbf{x})=1$ for all $\mathbf{x}$. Assuming logistic ecological dynamics, $K(\mathbf{x})$ is then the equilibrium density of a population that is monomorphic for phenotype $\mathbf{x}$. The adaptive dynamics of the phenotype $\mathbf{x}$ is a system of differential equations for $d\mathbf{x}/dt$. To derive the adaptive dynamics, one defines the invasion fitness $f(\mathbf{x},\mathbf{y})$ as the per capita growth rate of a rare mutant phenotype $\mathbf{y}$ in the monomorphic resident $\mathbf{x}$ population that is at its ecological equilibrium $K(\mathbf{x})$:

[TABLE]

The expression for the invasion fitness reflects the fact that the growth rate of the mutant $\mathbf{y}$ is negatively affected by the effective density experienced by the mutant, $\alpha(\mathbf{x},\mathbf{y})K(\mathbf{x})$, discounted by the carrying capacity $K(\mathbf{y})$ of the mutant (see (Doebeli, 2011, Doebeli and Ispolatov, 2014) for more details). Note that $f(\mathbf{x},\mathbf{x})=0$ for all $\mathbf{x}$. The invasion fitness $f(\mathbf{x},\mathbf{y})$ gives rise to the selection gradients in the $i=1,...,d$ phenotypic components:

[TABLE]

The selection gradients in turn define the adaptive dynamics as a system of differential equations on phenotype space $\mathbb{R}^{d}$, which is given by

[TABLE]

Here $\mathbf{s}(\mathbf{x})$ is the column vector $(s_{1}(\mathbf{x}),...,s_{d}(\mathbf{x}))$, and $\mathbf{M}(\mathbf{x})$ is the mutational variance-covariance matrix. In this matrix, the diagonal elements contain information about the size and rate of mutations in each of the phenotypic dimensions, whereas the off-diagonal elements contain information about the covariance between mutations in two different phenotypic dimensions. This matrix essentially captures “evolvability” of a population and generally depends on the current resident phenotype $\mathbf{x}$, and influences the speed and direction of evolution. For simplicity, we assume here that this matrix is the identity matrix. For more details on the derivation of the adaptive dynamics (3) we refer to a large body of primary literature (e.g. (Dieckmann and Law, 1996, Geritz et al., 1998, Diekmann, 2004, Leimar, 2009, Doebeli, 2011)). We note that the adaptive dynamics (3) can be derived analytically as a large-population limit of an underlying stochastic, individual-based model that is again defined based on the competition kernel $\alpha(\mathbf{x},\mathbf{y})$ and the carrying capacity $K(\mathbf{x})$ (Dieckmann and Law, 1996, Champagnat et al., 2006, 2008).

Specifically, here we consider a class of systems that are defined by competition kernels of the form

[TABLE]

Here the coefficients $b_{ij}$ in the first sum on the right hand side are arbitrary and correspond to the simplest form of a generic, non-symmetric competition kernel that can generate non-stationary evolutionary dynamics. It can be interpreted as the lowest-order (non-trivial) term from a Taylor expansion of an unknown non-symmetric competition function. Adaptive dynamics of asymmetric competition has been studied quite extensively (e.g. Law et al. (1997), Kisdi (1999), Doebeli (2011)), and is necessary to generate single-species non-equilibrium dynamics in high-dimensional phenotype spaces (Doebeli and Ispolatov, 2013, 2014). The second sum on the right hand side represents “Gaussian competition”, according to which the competitive impact between individuals increases with phenotypic similarity between the competing individuals. The parameters $\sigma_{i}$ measure how fast the effect of competition declines as phenotypic distance in the $i$-component increases. For the carrying capacity we assume

[TABLE]

This implies that the carrying capacity imposes an element of stabilizing selection for the phenotype $\mathbf{x}=0$, at which the carrying capacity is maximal. Thus, the frequency-dependent component of selection is generated by the competition kernel, whereas the frequency-independent component of selection is due to the carrying capacity. With these assumptions, the adaptive dynamics (3) become

[TABLE]

We note that the terms $-x_{i}^{3}$ in (6) are due to the carrying capacity and serve to contain the trajectories of (6) in a bounded domain of phenotype space. Also, the Gaussian part of the competition kernel does not affect the adaptive dynamics of monomorphic populations, i.e., the $\sigma_{i}$ do not appear in (6), because the Gaussian part always has a maximum at the current resident, and hence the corresponding first derivative in the selection gradient (2) is 0.

The system of ODEs (6) describes the trajectory of an evolving monomorphic population in phenotype space $\mathbb{R}^{d}$. In Doebeli and Ispolatov (2014) we have shown that for general competition kernels $\alpha$ such trajectories can be very complicated, particularly when the dimension $d$ is large. With complex evolutionary dynamics, trajectories can be quasi-periodic or chaotic, and typically visit many different regions of phenotype space over evolutionary time. When $d$ is low the dynamics tend to be simpler, and often converge to an equilibrium attractor. We can assess the likelihood of equilibrium dynamics for a given dimension $d$ by choosing the $d^{2}$ coefficients $b_{ij}$ in (6) randomly and independently, e.g. from a normal distribution with mean 0 and variance 1, solving the resulting adaptive dynamics (6) and checking whether it converges to an equilibrium. If this is done repeatedly, we can approximate the probability of equilibrium dynamics as the fraction of runs that converged to an equilibrium. For $d=1$ the probability of equilibrium dynamics is of course 1, and for $d=2,3,4$, the resulting probabilities of equilibrium dynamics are approximately $85\%$, $81\%$ and $74\%$, respectively. These are the dimensions that we will primarily use in the analysis presented below, but we note that the probability of equilibrium dynamics goes to 0 for large $d$ (Edelman, 1997, Doebeli and Ispolatov, 2014).

Multi-cluster adaptive dynamics

Here we are interested in the question of how diversification and subsequent coexistence of species (also called phenotypic clusters or simply clusters through the text) affects the evolutionary dynamics. While the Gaussian term in the competition kernel (4) does not affect the adaptive dynamics of single monomorphic populations, this term is crucial for determining whether evolutionary diversification occurs. For one-dimensional phenotype spaces ($d=1$) this is very well known and is encapsulated in the concept of evolutionary branching (Metz et al., 1992, Geritz et al., 1998, Doebeli, 2011). An evolutionary branching point is an equilibrium point of (6) that is both an attractor for the adaptive dynamics and a fitness minimum. The reason that such points exist in the competition models considered here is precisely that the Gaussian term does not affect the adaptive dynamics, but does affect the curvature of the fitness landscape, i.e., the second derivative of the invasion fitness (1). In particular, small enough $\sigma_{i}$’s in the Gaussian term will make any equilibrium point a fitness minimum, and hence will give rise to evolutionary diversification. Evolutionary branching in scalar traits has been described in a plethora of different models (for an overview we refer to Eva Kisdi’s website at the Department of Mathematics and Statistics at the University of Helsinki, http://www.mv.helsinki.fi/home/kisdi/addyn.htm). In high-dimensional phenotype spaces, equilibrium points of (6) can also be fitness minima along some directions in phenotype space. For this to happen the Hessian matrix of second derivatives of the invasion fitness (1), evaluated at the equilibrium, must have positive eigenvalues. Indeed, in higher dimensional phenotype spaces the conditions for the existence of positive eigenvalues of this Hessian matrix, and hence for diversification, generally become less stringent (Doebeli and Ispolatov, 2010, Débarre et al., 2014, Svardal et al., 2014).

Importantly, evolutionary diversification can also occur from non-equilibrium adaptive dynamics trajectories (Ito and Dieckmann, 2014, Ispolatov et al., 2016). If the adaptive dynamics (6) exhibit non-equilibrium dynamics, the crucial quantity determining whether diversification occurs is again the Hessian matrix of second derivatives of the invasion fitness (1), but now restricted to the subspace of phenotype space that is orthogonal to the selection gradient (Ispolatov et al., 2016). Essentially, diversification can occur in orthogonal directions in which this Hessian has positive curvature, and hence in which the invasion fitness has a minimum. Because the population is still evolving along the selection gradient, elucidating the exact conditions for diversification requires a careful analysis (Ito and Dieckmann, 2014). In the present context, the implication of these results is that, just as with equilibrium adaptive dynamics, diversification can occur along non-equilibrium trajectories of (6) if the $\sigma_{i}$ in the Gaussian term of the competition kernel are small enough, i.e., if the frequency dependence generated by Gaussian competition is strong enough (Ispolatov et al., 2016).

To investigate the process of diversification and the subsequent coevolutionary dynamics, we extend the adaptive dynamics (6) to several coexisting phenotypic clusters as follows. We assume that an evolving community consists of $m$ monomorphic populations, each given by a phenotype $\mathbf{x}_{r}$, $r=1,...,m$, with phenotypic components $x_{ri}$, $i=1,...d$ (where $d$ is the dimension of phenotype space). Let $N_{r}$ be the population density of cluster $\mathbf{x}_{r}$. Then the ecological dynamics of the $m$ clusters are given by the system of logistic differential equations

[TABLE]

Let $N_{r}^{*}$, $r=1,...,m$ denote the equilibrium of system (7) (more generally, for the purposes of deriving the adaptive dynamics, the quantities $N_{r}^{*}$ are suitable time averages of population densities over the ecological attractor of (7); however, our extensive numerical simulations indicated that (7) always converges to an equilibrium). Making the traditional adaptive dynamics assumption that ecological dynamics occur on a faster time scale than evolutionary dynamics, we calculate the invasion fitness function in cluster $r$ based on the densities $N_{r}^{*}$ of the various clusters:

[TABLE]

Here $\mathbf{x}_{1},...,\mathbf{x}_{m}$ describe the phenotypic state of the resident population, and $\mathbf{x}_{r}^{\prime}$ denotes the mutant trait in cluster $r$, $r=1,...,m$.

Taking the derivative of (8) with respect to $\mathbf{x}_{r}^{\prime}$ and evaluating it at the resident, $\mathbf{x}_{r}^{\prime}=\mathbf{x}_{r}$, yields the components of the selection gradient $\mathbf{s}_{r}$ for the cluster $r$ as:

[TABLE]

For coevolutionary adaptive dynamics, one has to take into account that the rate of mutations in each evolving phenotypic cluster is proportional to the current population size of that cluster (Dieckmann and Law, 1996), and hence that the speed of evolution is influenced by the population size. In the single-cluster system such consideration only rescales time without affecting the geometry of the trajectory and thus is usually ignored. However, in the multi-cluster system, instead of assuming that the mutational process is described by the identity matrix as in (3), we now assume that in each cluster $r$, the mutational variance-covariance matrix $M_{r}$ is a diagonal matrix with entries $N_{r}^{*}$. This generates the following $d\cdot m$ differential equations describing the adaptive dynamics in the coevolving community:

[TABLE]

For the multicluster adaptive dynamics, the equation (8,9,10) replace their single-cluster analogs (1,2,3). It is important to note that the Gaussian part of the competition kernel $\alpha$ not only affects whether diversification occurs, but in contrast to the adaptive dynamics of single monomorphic populations, the Gaussian term will indeed affect the coevolutionary adaptive dynamics (10) of the phenotypic clusters that coexist after diversification has occurred, because it affects both the ecological dynamics (7) and the selection gradient (9).

Numerical procedure

To study diversification and subsequent multi-cluster adaptive dynamics, we implemented the following iterative numerical scenario:

Step 1: Each simulation run is initiated with a randomly generated $d\times d$ matrix of the coefficients $b_{ij}$ for the competition kernel (4). The coefficients are drawn from a Gaussian distribution with zero mean and $d^{-1/2}$ variance. As explained in (Doebeli and Ispolatov, 2014), this is done to keep the sum of the $d$ terms $\sum_{j=1}^{d}b_{ij}x_{j}$ in (6) of order $x_{i}$, i.e. independent of $d$. Then a certain number of clusters, given by a parameter $m_{0}$, each with population size of order $1$, are randomly placed near the phenotype [math], i.e., near the maximum of the carrying capacity.

Step 2: For a given set of phenotypic clusters, the population dynamics of all clusters is solved using the ecological dynamics (7). The system of differential equations is integrated using a 4th-order Runge-Kutta algorithm for $\sim 10^{3}$ time steps of duration $dt\sim 10^{-2}$ to ensure convergence to the equilibrium (or, in case there is no such convergence, to ensure a correct calculation of the time average of the various population densities). If the population density of a given cluster falls below the threshold $N_{min}\sim 10^{-8}$, the cluster is eliminated from the system. During the ecological dynamics the evolutionary dynamics is frozen and evolutionary time does not advance.

Step 3: After calculating the $N_{r}^{*}$, $r=1,...,m$ (where $m$ is the current number of clusters), the adaptive dynamics of the phenotypes of the clusters is advanced via (9,10) using a 4th-order Runge-Kutta algorithm with a typical time-step $d\tau\sim 10^{-2}$, by which the evolutionary time is advanced as well. After this evolutionary time step, the ecological dynamics are recalculated, potentially preceded by the following step 4, which is only performed if the corresponding time condition is satisfied.

Step 4: The level of diversity, i.e., the number of clusters in the system, is controlled as follows. Each $\tau_{c}$ time units the distances between clusters are assessed. If the distance between two or more clusters is below a threshold $\Delta x\sim 10^{-3}$, these clusters are merged, preserving the total population size of the merged clusters and the position of their centre of mass. Immediately after this comparison step, the total number of clusters is compared to the target number of clusters, which is given by a system parameter $m_{max}$. If the current number of clusters is below $m_{max}$, a new cluster is created by randomly picking an existing cluster, splitting it in half and separating the two new clusters in a random direction in phenotype space by the distance of the merging threshold, $\Delta x$.

Step 5: In our simulations, we take measurements at regular time intervals (ranging from $\tau_{m}\sim 1-10$ time units). One of the main quantities of interest is the average per capita evolutionary speed $v$ in the evolving community, which is the average of the norms of the vectors of trait variation (evolution) rates in each cluster, weighted by the cluster population size, computed as

[TABLE]

This quantity is a strong indicator of the nature of the evolutionary dynamics of the coevolving system. In particular, our very extensive numerical simulations indicate that when the average speed falls below $10^{-2}$, then the system eventually exhibits equilibrium evolutionary dynamics. In contrast, when the average evolutionary speed remains high, the coevolving system tends to exhibit complicated, non-equilibrium dynamics, with the majority of the clusters exhibiting large fluctuations in phenotype space over evolutionary time. An example of such non-equilibrium coevolution is given in the next section. Other measurements include the position and population size of all clusters in the system, and the number of “distinct” clusters separated by a “visible” distance $\Delta X=0.1$. These measurements can also be averaged over time.

For any given simulation run initiated by step 1 above, steps 2-5 were repeated iteratively until a specified final simulation time is reached, or until evolution comes to a halt, which by our definition occurs when the average evolutionary speed falls below a threshold, $v<10^{-4}$. Our general approach consisted of simulating many different systems according to the above scheme, and then computing statistical characteristics such as the fraction of runs that result in non-equilibrium dynamics, or the average evolutionary speed as a function of the level of diversity (see Results section).

One crucial feature of our algorithm is the periodic generation of new clusters in step 4, which mimics diversification events, i.e., evolutionary branching. Diversification is thus modeled by simply adding new phenotypic clusters at certain points in time and close to existing clusters. This mimics the sympatric split of an ancestral lineage. Sympatric diversification is a theoretically robust phenomenon (Doebeli, 2011) and our procedure represents a shortcut for this phenomenon necessitated by computational feasibility. If such splitting is not feasible given the current ecological circumstance, the new cluster will not diverge phenotypically from the ancestor, and hence will be merged again with the ancestor (see below). Alternatively, newly generated clusters may go extinct ecologically. In either case, speciation was not successful. Thus, in our models it is the ecological circumstances that determine whether speciation can occur or not, but the process of speciation itself (i.e., the splitting) is performed in a simplified manner. If speciation is successful and the newly generated clusters diverge and persist ecologically, then diversity has increased (unless other clusters go extinct). We note that by construction, the maximal level of diversity in a given simulation run, i.e., the number of different clusters, cannot exceed the parameter $m_{max}$. Therefore, this parameter allows us to control the level of diversity in a given simulation.

There are in principle other, less artificial ways to model diversification. In particular, stochastic, individual-based based models and partial differential equation models (Champagnat et al., 2006, 2008) have been used to describe the evolutionary dynamics of phenotype distributions. In such models, diversification is an emergent property that is reflected in the formation of new modes in the evolving phenotype distributions. While these techniques are very useful in general, they are currently not computationally feasible for the statistical approach that we employed here, which requires systematic simulation of many different systems. Also, they would not allow for control of the level of diversity, as the number of phenotypic modes would simply be an emergent property of the evolving system. Nevertheless, we have used these alternative techniques to illustrate the robustness of salient results using particular examples. A more detailed description of these techniques is given in the Appendix. Another alternative would be to assume that new clusters (species) are assigned phenotypes that are chosen randomly in phenotype space, rather than close to an existing cluster. This could correspond to immigration of new species into an existing community. However, we would not expect this to affect our main results, because with complicated evolutionary dynamics, the initial phenotypic position of a given cluster becomes irrelevant after some time.

Finally, we note that the merging of clusters (species) is done solely for computational reasons and has no biological meaning (apart from designating organisms that are closely related and phenotypically very close as belonging to the same species). Merging of clusters only occurs shortly after a new cluster is seeded close to an existing one, and only if the new cluster does not diverge from the existing one (i.e., only if the ecological conditions for diversification are not satisfied). If divergence is successful, the clusters will never again get close enough to other clusters to be merged because of the repelling force of frequency-dependent competition. Thus, the only function of merging is to prevent the number of clusters from artificially becoming very large.

Results

The parameter that controls the level of diversity in our simulations is $m_{max}$, which is the maximal number of different phenotypic clusters allowed to be present at any point in time in an evolving community (see step 4 in the Methods section). Our first result is obtained by allowing this parameter to be very large, so that we can estimate the number of clusters that eventually coexist by simply running the simulations for a long time and recording the number of clusters at which the diversity equilibrates. We denote by $M_{\sigma,d}$ the equilibrium number of clusters for a given phenotypic dimension $d$ and strength of the Gaussian component $\sigma$ in the competition kernel (4). We found that such equilibrium level of diversity increases exponentially with the dimension $d$ of phenotype space, and decreases with the strength $\sigma$ (Figure 1). Here and below we assume for simplicity that the $\sigma_{i}$ are the same in all phenotypic directions, $\sigma_{i}=\sigma$ for $i=1,...,d$. In the Appendix we indicate scaling relationships that hold for $M_{\sigma,d}$ as functions of the parameters $\sigma$ and $d$. In general, diversity is only maintained if $\sigma\lesssim 1$, which is roughly the scale of the phenotypic range set by the carrying capacity (5). Only if $\sigma\lesssim 1$, the equilibrium level of diversity increases exponentially with increasing dimension of phenotype space, Figure 1.

Our main results are now obtained based on the observation that by fixing the parameter $m_{max}$ at a value $\leq M_{\sigma,d}$ for a given $d$ and $\sigma$, the community will typically evolve to a diversity level $m_{cluster}$ of approximately $m_{max}$. That is, if the diversity is constrained to be below the maximal level of diversity possible for a given set of parameters, then the diversity will typically evolve to the value set by the constraint. Note that this is an “average” statement about many simulations runs, i.e., many different choices of the coefficients $b_{ij}$ and stochastic realizations of cluster splitting. While some simulation runs will result in a diversity that is lower than $m_{max}$ (which may reflect an intrinsic state of the system for the given set of coefficients, or a long-living metastable state which has not yet reached its full diversity), most runs will evolve to the level of diversity that is prescribed by this parameter.

This allows us to then assess, for a given level of diversity, the nature of the coevolutionary dynamics that unfolds in communities with that level of diversity. Two paradigmatic examples are shown in Figure 2. We first set the level of diversity $m_{max}=12$, which is far below the saturation level $M_{\sigma,d}$ for the given system. Starting from very few clusters the diversity quickly evolves to the level set by $m_{max}$, and the coexisting clusters then exhibit complicated, non-stationary evolutionary dynamics, with all clusters undergoing sustained and irregular fluctuations in phenotype space (Fig. 2a). This type of complicated dynamics is characterized by average evolutionary speeds $v>10^{-2}$. In the same system, but now with a value of $m_{max}$ that lies above the saturation level $M_{\sigma,d}$, the diversity evolves to the saturation level, at which the community consists of ca. 30 coexisting phenotypic clusters (Fig. 2b). In this saturated state, the average evolutionary speed is much lower than $10^{-2}$, and the community exhibits much more stationary coevolutionary dynamics (that would eventually converge to a coevolutionary equilibrium). Moreover, the saturated community exhibits a characteristic pattern of over-dispersion in phenotype space due to competitive repulsion caused by the Gaussian component of the competition kernel (see also Fig. A1 in the Appendix).

To obtain a more systematic characterization of the coevolutionary dynamics as a function of the diversity of the evolving community, we ran, for a given dimension of phenotype space $d$ and strength of competition $\sigma$, 100 simulations with randomly chosen coefficients $b_{ij}$ for each $m_{max}=1,...,M$, where $M$ is some number that is larger than the saturation level of diversity $M_{\sigma,d}$. For each run, we recorded the average per capita evolutionary speed $v$ and the number of phenotypic clusters, i.e., the level of diversity, present at the end of $1000$ evolutionary time units (averaged over the last 4 time units). We classified the dynamics into equilibrium dynamics if the average speed $v$ was $<10^{-2}$, and non-equilibrium dynamics otherwise. As mentioned earlier, this was based on individual inspection of many simulation that ran longer than $1000$ time units, which showed that the threshold $10^{-2}$ is a very good indicator of whether the coevolutionary system eventually equilibrates.

Our main results are shown in Figures 3 and 4. The general pattern is that the probability of non-equilibrium dynamics increases as diversity increases from single-cluster communities to communities with a few clusters (Figure 3). For intermediate diversity, the fraction of non-equilibrium dynamics remains high. For communities with high diversity, the fraction of non-equilibrium dynamics starts to decrease, and almost almost all communities with a diversity close to the saturation level $M_{\sigma,d}$ exhibit equilibrium coevolutionary dynamics. To illustrates these trends, we give a more detailed account of the average velocities $v$ defined in (11) in the coevolving communities (Fig. 4). It shows that there is an exponential decrease in the average speed as the diversity increases, and that there is a substantial fraction of low-diversity communities that exhibit equilibrium dynamics. The exact shape of these patterns depends on $d$ and $\sigma$ (Figures 3 and 4), but whenever diversification is possible, the overall trend is that non-equilibrium dynamics are most likely at intermediate levels of diversity, and that high levels of diversity tend to generate equilibrium coevolutionary dynamics.

The patterns shown in Figs. 3 and 4 are based on many different simulated communities with different levels of diversity. However, similar patterns can be observed in simulations of single communities as they evolve from low to high diversity, i.e., as they undergo an adaptive radiation. Such a radiation, starting from a single phenotypic cluster, is shown in Fig. 5A. Over time the evolving community becomes more diverse due to adaptive diversification, and as a consequence the nature of the coevolutionary dynamics of the community changes. In the example shown in Figure 5A, the coevolutionary dynamics are fast for low to intermediate levels diversity, and then slow down as the community acquires more and more species, until eventually the community reaches a coevolutionary equilibrium at the diversity saturation level. Again, the slowdown of the evolutionary speed during an adaptive radiation appears to occur exponentially with an increase in diversity. This can also be seen by running a given community defined by a given set of coefficients $b_{ij}$ for different values of the parameter $m_{max}$, determining the level of diversity possible in the evolving community. The evolutionary speed exponentially decreases with the diversity given by $m_{max}$ (Fig. 5B). We currently do not have a mechanistic explanation for the exponential decay in evolutionary rates with increasing diversity. It is informative to watch the process of diversification and subsequent evolutionary slowdown unfold dynamically. To verify that the observed dynamical pattern is not an artifact of the adaptive dynamics approximation, we performed the individual-based and partial differential equation simulations of the same system. The movies in Videos in the Appendix, corresponding to the scenario used for Figures 2B and 5A, confirm that all three methods produce qualitatively similar evolutionary pictures. The detailed descriptions of the individual-based and partial differential equation methods are given in the Appendix.

Another interesting, although perhaps not so surprising observation for single adaptive radiations concerns the rate of accumulation of new species in the evolving ecosystem. Figure 5C shows the number of species as a function of time during the adaptive radiation scenario used for Figure 5A, illustrating that the rate of diversification is highest at the beginning of the radiation, and then slows down as the community evolves towards the diversity saturation level. The details of these dynamics depend on system parameters, and in particular on the rate at which new species are introduced into the system, but the qualitative behaviour of diversification rates, which are initially high and then slow down, is common to all adaptive radiations generated by our models.

Discussion

We investigated the expected long-term evolutionary dynamics resulting from competition for resources in models for gradual evolution in high-dimensional phenotype spaces. In reality, most organisms have many different phenotypic properties that impinge on their ecological interactions in generally complicated ways, and here we assumed that multi-dimensional phenotypes determine logistic ecological dynamics through the competition kernel and the carrying capacity. We then used a coevolutionary adaptive dynamics algorithm to extend the ecological dynamics to macroevolutionary time scales, and we used a statistical approach to capture general properties of the ensuing evolutionary dynamics.

If the negative frequency-dependence generated by the competition kernel is strong enough, competition results in repeated adaptive diversification, and hence in communities of coevolving phenotypic species. By randomly choosing many different competition kernels, we showed that the complexity of the coevolutionary dynamics in such communities is expected to be highest for intermediate levels of phenotypic diversity. In particular, as the evolving communities increase in diversity towards the saturation level, i.e., the maximal number of different species that can coexist, the evolutionary dynamics becomes simpler, and communities at the saturation level are expected to exhibit a coevolutionary equilibrium. We also showed that the diversity saturation level increases exponentially with the dimension of phenotype space.

We have used a statistical approach to determine the expected long-term evolutionary dynamics resulting from competition for resources. We have assumed that multi-dimensional phenotypes determine logistic ecological dynamics through the competition kernel and the carrying capacity, and we then used a coevolutionary adaptive dynamics algorithm to extend the ecological dynamics to macroevolutionary time scales. If the negative frequency-dependence generated by the competition kernel is strong enough, competition results in repeated adaptive diversification, and hence in communities of coevolving phenotypic species. By randomly choosing many different competition kernels, we showed that the complexity of the coevolutionary dynamics in such communities is expected to be highest for intermediate levels of phenotypic diversity. In particular, as the evolving communities increase in diversity towards the saturation level, i.e., the maximal number of different species that can coexist, the evolutionary dynamics becomes simpler, and communities at the saturation level are expected to exhibit a coevolutionary equilibrium. We also showed that the diversity saturation level increases exponentially with the dimension of phenotype space.

Our interpretation of these findings is that in low-dimensional phenotype spaces such as the ones considered here, evolutionary dynamics of single species are expected to converge to an equilibrium (Doebeli and Ispolatov, 2014). However, as diversity increases, the different phenotypic clusters will “push” each other around evolutionarily due to frequency-dependent competition. This occurs mostly due to the repulsive nature of pairwise interaction induced by the Gaussian term in the competition kernel (4): clusters that move further apart decrease competition felt from each other. For example, a splitting of a cluster stuck in an attractive fixed point of the adaptive dynamics creates two offspring which may become moving again if the repulsion between clusters is stronger than the attraction of the fixed point. As long as diversity is not very high, i.e., as long as there is enough available niche or unoccupied phenotype space, this typically results in non-equilibrium coevolutionary dynamics, thus leading to an increase in evolutionary complexity with phenotypic diversity. As the diversity keeps increasing towards saturation levels, which for each phenotypic dimension is determined roughly by the ratio of the widths of the carrying capacity and the competition kernel (see Video 2), the available carrying capacity niche gets filled, so that the evolving clusters “have nowhere to go” evolutionarily. An analogy with gas-liquid-solid phase transitions may illustrate this in the following way: As in the dynamics of molecules, the adaptive dynamics of phenotypic clusters contains a pairwise-repulsive term, which originates from the Gaussian term in the competition kernel. A few-cluster regime qualitatively corresponds to the gas phase, when the range of the repulsive interaction is significantly less than the typical distance between clusters. As the number and thus density of clusters increases, the repulsive interaction becomes more relevant, constraining the individual motion of clusters and resulting in a liquid-like behaviour, where clusters are predominantly localized and occasionally hop to a new location. Finally, the maximum cluster density creates a crystal-like structure, albeit not necessarily entirely symmetric due to the randomly generated $b_{ij}$ terms in the adaptive dynamics. The motion of individual clusters is heavily constrained by its neighbours via mutual repulsion, while the collective motion of an ensemble of clusters is limited by the carrying capacity function. Thus, phenotypic saturation leads to a state in which the coevolving clusters are strongly constrained evolutionarily by the other clusters in the community, and hence to coevolutionary equilibrium dynamics.

Some empirical support for an initial increase in the complexity of evolutionary dynamics with the number of species in an ecosystem comes from the laboratory evolution experiments of Lawrence et al. (2012), who showed that the speed of adaptation to novel environments is higher in bacterial species that are part of microbial communities with a small number of competitors than when evolving in monoculture. However, our results are seemingly in contrast to previous theoretical results about the effect of diversity on evolutionary dynamics (Johansson, 2008, de Mazancourt et al., 2008). These authors essentially argued that while a single species is free to evolve in response to changes in the environment, evolution of the same species is more constrained in a community of competitors, in which other species are more likely to evolutionarily occupy new niches. Hence diversity is expected to slow down evolution. However, these models only describe evolution in 1-dimensional phenotypes, and may thus miss the complexity arising in higher-dimensional spaces. Moreover, even in higher-dimensional spaces, the arguments for evolutionary slowdown presented in (Johansson, 2008, de Mazancourt et al., 2008) essentially correspond to our observation of a slow-down when diversity reaches saturation, at which point evolutionary change in each species is indeed constrained due to competing species occupying all available niches. Our approach also needs to be distinguished from approaches based primarily on ecological dynamics, as in Shtilerman et al. (2015). In these approaches, emerging ecological communities are also modelled by periodically adding new species, but there is no underlying phenotype space that would determine competitive interactions. Instead, every time a new species added, its interaction coefficients with the already existing species are chosen according to a specific, randomized procedure. This leads to interesting results, such as saturating levels of diversity after initially fast and fluctuating increases from low levels of diversity. However, since there is no underlying phenotype space, this approach does not reveal the evolutionary dynamics of continuous phenotypes, and in particular, it does not yield any information about the effects of the dimension of phenotype space on the evolutionary dynamics or on the amount of diversity at saturation.

There has been much interest in recent years in determining the effects of phylogenetic relationships on the functioning of ecosystems (e.g. Ives and Godfray (2006), Cadotte et al. (2013), Nuismer and Harmon (2015)). The intuitive notion is that phylogenetic information has predictive power for ecological interactions if recently diverged species are more likely to interact than those that diverged long ago. More specifically, Nuismer and Harmon (2015) have argued that phylogenetic information is most likely to be relevant for ecosystem dynamics if ecological interactions are based on phenotypic matching, so that species with more similar trait values are more likely to interact strongly. Our models have a component of phenotypic matching due to the Gaussian part of the competition kernel, but they also have a strong component of different types of interactions due to the “random” part of the competition kernel given by the coefficients $b_{ij}$. As we have shown, it is this non-Gaussian part of the competition kernel that causes the complicated coevolutionary dynamics, and it is this complexity in turn that makes phylogenetic signal largely irrelevant in our models.

A full phylogenetic analysis of the macroevolutionary dynamics generated by our models is beyond the scope of this work, but we can provide some basic insights based on the complicated evolutionary dynamics in phenotype space that the different phenotypic clusters (species) perform when there is an intermediate number of clusters in the coevolving community. An example of this is shown in the movie in Figure A1A. Here, after an initial phase of diversification, the community contains 12 coevolving clusters. These clusters move on a complicated evolutionary trajectory, with each cluster undergoing large evolutionary changes without further diversification. No matter what the phylogenetic relationship between these clusters (as given by their emergence from the single initial cluster), it is clear that because of the large evolutionary fluctuations in phenotype space of each cluster (species), there will be no consistent correlation between phylogenetic relationship and phenotypic distance. Even if there were such a correlation (positive or negative) at a particular point in time, it would change over time due to the large evolutionary fluctuations of each cluster over time. This is illustrated in Figure A1B, which shows that no persistent correlation pattern between phylogenetic and phenotypic distance should be expected in communities with an intermediate amount of diversity. In particular, recently diverged species are not more likely to interact than those diverged less recently, because the evolving community has a short “phenotypic memory” due to complicated evolutionary dynamics.

However, when further diversification is allowed, so that the system reaches its saturation level of diversity, the coevolving community not only becomes more diverse, but the evolutionary dynamics slows down, leading to ever smaller phenotypic fluctuations. In particular, new clusters emerging towards the end of the assembly of the evolutionarily stable community will stay phenotypically closer to their phylogenetically most closely related clusters, i.e., to their parent or sister species. Therefore, in the last phase of community assembly a positive correlation between phylogenetic and phenotypic distance can be expected to build up at least to some extent. This is illustrated in Figure A1B. Thus, weak phylogenetic signals are expected to develop towards the end of community assembly.

Regarding adaptive radiations, two observations emerge from our models. The first concerns the classical notion that rates of diversification should decline over the course of a radiation (Schluter, 2000, Gavrilets and Losos, 2009), a pattern that seems to have good empirical support (Schluter, 2000, McPeek, 2008, Rabosky and Lovette, 2008, Gavrilets and Losos, 2009). Our models confirm this pattern of declining rates of diversification (Figure 5). The second observation is that rates of evolution should generally slow down with an increase in diversity. This should not only be true when different ecosystems are compared (Figures 3,4), but also during an adaptive radiation in a single evolving community (Figure 5). Thus, we would expect the evolutionary dynamics to be faster and more complicated early in an adaptive radiation, and to slow down and eventually equilibrate late in the radiation. This corresponds to the so-called “early-burst” model of macroevolution (Gavrilets and Losos, 2009, Harmon et al., 2010) in the context of adaptive radiations. This model predicts that when lineages enter novel “adaptive zones” (Simpson, 1944), such as novel ecological niches, evolutionary rates in the lineage should be fast initially and then slow down as the adaptive zone gets filled with diverse phenotypes. Harmon et al. (2010) found little evidence for the early-burst model when analyzing a large set of data from many different clades. Nevertheless, these authors noted that younger clades have higher rates of evolution than older clades, which points to the fact that evolutionary rates may slow down with clade age. Moreover, few clades in their data set correspond to the type of very fast adaptive radiation envisaged and observed in our models, and they did not consider high-dimensional phenotypes. Finally, Harmon et al. (2010) note that groups with a larger proportion of sympatric species early in their history would be more likely to exhibit an early-burst pattern. In our models, adaptive radiations occur in complete sympatry and indeed produce the early burst pattern.

According to Slater and Pennell (2013), the jury on early-burst models is still out, and in fact substantial evidence for this model has accumulated in recent years. For example, Cooper and Purvis (2010) reported an early burst in body size evolution in mammals, Weir and Mursleen (2013) observed an early-burst pattern in the evolution of bill shape during adaptive radiation in seabirds, Gonzalez-Voyer and Kolm (2011) and Arbour and López-Fernández (2013) reported early-burst patterns in morphological and functional evolution in cichlids, and Benson et al. (2014) described patterns of early bursts in the evolution of dinosaur morphology.

Uyeda et al. (2011) have incorporated the early-burst concept into a macroevolutionary perspective in which over very long evolutionary time scales, rare but substantial phenotypic bursts alternate with more stationary periods of bounded phenotypic fluctuations, somewhat reminiscent of the concept of punctuated equilibrium (Gould and Eldredge, 1977) when applied to rates of phenotypic evolution (Pennell et al., 2014). We think that the models presented here could provide a microevolutionary basis for such a perspective if they are extended by considering evolutionary change in the dimension of the phenotype space that determines ecological interactions. Such an extended theory would have three time scales: a short, ecological time scale, an intermediate time scale at which co-evolution and single diversifications take place in a given phenotype space, and a long time scale at which the number of phenotypic components increases (or decreases). Our hypothesis would then be that in such systems, periods of bounded evolutionary fluctuations near diversity saturation levels for a given dimension of phenotype space would alternate with bursts of rapid evolutionary change, brought about by an evolutionary increase in phenotypic dimensions and the subsequent increase in diversity and acceleration in evolutionary rates until a new saturation level is reached. The resulting long-term evolutionary dynamics would thus show periods of relative phenotypic stasis alternating with periods of fast evolution. This picture would fit very well with the “blunderbass” pattern envisaged in Uyeda et al. (2011). These authors proposed that the intermittent bursts in evolutionary rates are caused by lineages encountering novel “adaptive zones” (Simpson, 1944). Novel adaptive zones would correspond to the opening up of new habitats or new resources, which would in turn correspond to new phenotypes that determine use of the novel adaptive zone. Alternatively, novel adaptive zones could also be generated by the emergence of novel sets of regulatory mechanisms allowing novel uses of already existing habitats and resources (as e.g. when a trade-off constraint is overcome through gene duplication). In either case, novel adaptive zones would correspond to an increase in the dimensionality of ecologically important phenotypes.

It is interesting to note that such intermittent burst patterns have in fact been observed in phylogenetic data, and that they seem to be connected to novel, ecologically important phenotypes. Hopkins and Smith (2015) have shown that evolutionary rates in echinoids reveal at least two instances of rapidly accelerating and subsequently declining evolutionary rates, i.e., two intermittent bursts. Moreover, these bursts appear to be associated with the evolution of novel feeding strategies (Hopkins and Smith, 2015, Slater, 2015). Also, Brusatte et al. (2014) have shown that an evolutionary burst occurs in the dinosaur-bird transition, and it is tempting to conjecture that this burst was caused by the increase in phenotype dimensionality due to the proliferation of flight capabilities.

There is also good empirical support for our finding that the level at which diversity saturates increases with the dimension of phenotype space. Seehausen (2015) has argued that essentially, the high number of different ecologically relevant traits is the basis for the spectacular radiations of cichlids in African lakes. In conjunction with ecological opportunity, genetic and phenotypic flexibility, which appears to be at least in part due to gene duplications, has allowed this group of fish to reach a much higher diversity than other groups, such as cichlids in rivers or whitefish in arctic lakes, in which fewer phenotypes appear to be ecologically relevant (Hudson et al., 2011, Seehausen, 2015). In this context, we note that incorporating the evolution of the dimension of phenotype space may also shed light on the ongoing debate about whether diversity saturates over evolutionary time or not (Harmon and Harrison, 2015, Rabosky and Hurlbert, 2015). It seems that the answer could be “yes and no”: diversity saturates in the intermediate term for a given dimension of phenotype space, but does not saturate in the long term if the dimension of phenotype space increases over long evolutionary time scales, thus generating recurrent increases in saturation levels.

Our study has a number of limitations that should be addressed in future research. It is currently impractical to perform the statistical analysis presented here for phenotype spaces with dimensions higher than $4$ due to computational limitations. Our results indicate that the diversity saturation level, i.e., the maximal number of coexisting phenotypic clusters, increases rapidly with the dimension $d$ of phenotype space, which makes simulations of communities at saturation levels unfeasible. Nevertheless, we expect the salient result that coevolutionary dynamics slow down as communities reach the saturation level to be true in any dimension as long as the Gaussian component of competition in (4) affects all phenotypic directions. Also, in our approach we have assumed that the phenotypes determining competitive interactions are the same for intra- and inter-specific competition. This may be a fair assumption for closely related species, such as those coevolving in an adaptive radiation. However, for competition in more general ecosystems it may also be relevant to assume that from a total set of $d$ phenotypes, different subsets determine competition within a species and competition with various other species. In addition, to describe general ecosystems and food webs, it will be important to include not just competitive interactions, but also predator-prey and mutualistic interactions, each again determined by potentially high-dimensional phenotypes. Also, throughout we have assumed a simple unimodal form of the carrying capacity to represent the external environment. More complicated forms of the carrying capacity, and hence of the external fitness landscape will likely generate even richer patterns of coevolutionary dynamics and diversification. Finally, we have assumed throughout that evolving populations are well-mixed, and it will be interesting so see how the results generalize to spatially structured ecosystems. All these extensions remain to be developed.

We are of course aware of the fact that we did not include genetic mixing due to sexual reproduction in our models, and our method of describing diversification by simply adding new phenotypic clusters, although fairly standard, does not take into account the actual process of speciation. In sexual populations, adaptive diversification due to disruptive selection, as envisioned here, requires assortative mating, and the conditions for the evolution of various types of assortative mating, as well as for the likelihood of speciation once assortment is present, have been studied extensively (e.g. (Doebeli, 2011)). A general, if crude conclusion from this work is that when there is enough disruptive selection for diversification to occur in asexual models, then it is likely that adaptive speciation also occurs in the corresponding sexual models, although factors such as the strength of assortment, population size and linkage disequilibrium may become important. It would in principle be possible to incorporate sexual reproduction into the models presented here, e.g. along the lines of Gascuel et al. (2015). Our previous results (Doebeli and Ispolatov, 2010, Ispolatov et al., 2016) indicate that adaptive diversification is generally more likely in high-dimensional phenotype spaces, and we think that the present models serve well as a first approximation to study adaptive diversification and coevolutionary dynamics in evolving communities.

Ultimately, the applicability and relevance of our models for understanding macroevolutionary patterns in nature depends in part on being able to determine evolutionary rates of high-dimensional phenotypes from phylogenetic data, which appears to be a difficult problem (Denton and Adams, 2015, Adams, 2013, 2014). Nevertheless, overall we think that our approach of incorporating microevolutionary processes based on ecological interactions in high-dimensional phenotype spaces into statistical models for macroevolutionary dynamics has the potential to shed new light on a number of fundamental conceptual questions in evolutionary biology.

Acknowledgments

M. D. was supported by NSERC (Canada). I. I. was supported by FONDECYT grant 1151524 (Chile). Both authors contributed equally to this work.

Correlation between phylogenetic and phenotypic distance

For each pair of clusters (species) in an evolving community we define the phylogenic distance between them, $Pg$, as the number of links in the path between them on the phylogenic tree. To measure this distance, we add the following scheme to our evolutionary algorithm:

• •

The system is initialized with a single cluster.

• •

Each cluster splitting event produces two offspring separated by the distance 2. The distance between an offspring and all its existing neighbours is incremented by one.

• •

When two recently split cluster that failed to diverge are merged, the distance between the newly produced common cluster and each of its neighbours is calculated as the minimum of the distances of the two merged clusters minus one. This reflects the observation that merging events only happen with newly split clusters.

As a result, at any given time we know phylogenic distances between all pairs of clusters currently present in the system. To quantify the relation between the phenotypic and phylogenic similarity, we compute the correlation $C$ between phylogenetic and phenotypic distance as follows:

[TABLE]

where $Ph$ and $X$ are phylogenic and phenotypic distances between clusters, $\langle\ldots\rangle$ define the average over all pairs of clusters present in the system and $\sigma_{Pg}$ and $\sigma_{X}$ are the standard deviations of distances.

The above scheme allows us to track the correlation between phylogenetic and phenotypic distance over time, as illustrated in Figure A1. Fig. A1A shows the time dependence of $C$ for the simulation shown in Video 1, and in Fig. A1B shows the time dependence of $C$ for the simulation shown in Video 2. During the early phase of community assembly the correlation $C$ rapidly decays due to complicated coevolutionary dynamics of the emerging clusters. When the diversity of the coevolving community is kept intermediate (by setting the parameter $m_{C}$ to intermediate values, as in Video 1), the correlation between phylogenetic and phenotypic distance itself undergoes fluctuations around 0 (Fig. A1A). This is because the clusters in the community with intermediate diversity undergo large phenotypic fluctuations while their phylogenetic relationship is constant, because no further diversification (or extinction) occurs. However, when the diversity is allowed to reach saturation levels (by setting $m_{C}$ to a large value, as in Video 2), a positive correlation between phylogenetic and phenotypic distance develops in the final stages of community assembly, i.e., as the coevolving community reaches the saturation diversity and hence undergoes much smaller phenotypic fluctuations (Fig. A1B). Note that the correlation is still close to 0 during the early stages of community assembly, but some correlation remains at the end due to clusters emerging in the last phase of community assembly, which tend to stay phenotypically closer to their sister species because evolutionary dynamics become slow and stable.

Individual-based simulations

Individual-based realizations of the model were based on the Gillespie algorithm (Gillespie, 1976) and consisted of the following steps:

The system is initialized by creating a set of $K_{0}\sim 10^{3}-10^{4}$ individuals with phenotypes $\operatorname{\mathbf{x}}_{k}\in\mathbf{R}^{d}$ localized around the initial position $\operatorname{\mathbf{x}}_{0}$ with a small random spread $|\operatorname{\mathbf{x}}_{k}-\operatorname{\mathbf{x}}_{0}|\sim 10^{-3}$. 2. 2.

Each individual $k$ has a constant reproduction rate $\rho_{k}=1$ and a death rate

$\delta_{k}=\sum_{l\neq k}A(\operatorname{\mathbf{x}}_{l},\operatorname{\mathbf{x}}{{}_{k}})/[K_{0}K(\operatorname{\mathbf{x}}_{k})]$, as defined by the logistic ecological dynamics. 3. 3.

The total update rate is given by the sum of all individual rates, $U=\sum_{k}(\rho_{k}+\delta_{k})$. 4. 4.

The running time $t$ is incremented by a random number $\Delta t$ drawn from the exponential distribution $P(\Delta t)=U\exp(-\Delta tU)$. 5. 5.

A particular birth or death event is randomly chosen with probability equal to the rate of this event divided by the total update rate $U$. If a reproduction event is chosen, the phenotype of an offspring is offset from the parental phenotype by a small mutation randomly drawn from a uniform distribution with amplitude $\mu=10^{-3}-10^{-2}$. 6. 6.

The individual death rates $\delta_{k}$ and the total update rate $U$ are updated to take into account the addition or removal of an individual. 7. 7.

Steps 4-6 are repeated until $t$ reaches a specified end time.

The movie in Video A2 shows the dynamics of the individual-based model corresponding to the adaptive dynamics simulation shown in Video A1, which is the same as the scenario used for Video 2 in the main text (note that the movie in Video A1 runs for $t=1200$ time units, whereas the movie in Video 2 runs for $t=400$ time units).

Partial differential equation models

A deterministic large-population limit of the individual-based model is obtained as the partial differential equation (PDE)

[TABLE]

where $N(\operatorname{\mathbf{x}},t)$ is the population distribution at time $t$ (Champagnat et al., 2006). The second term of the right hand side is a diffusion term that describes mutations, with the diffusion coefficient typically set to $D\sim 10^{-4}-10^{-3}$. Local maxima of the solution $N(x,t)$ can be interpreted as positions of the centers of the phenotypic clusters. Their dynamics are shown in Video A3. For any given scenario, the corresponding adaptive dynamics solution can be used to determine the single- or few-cluster trajectory, and hence to approximately determine the region occupied by the system in phenotype space over time. Note that the deterministic PDE model is invariant with regard to the coordinate change $\operatorname{\mathbf{x}}\rightarrow-\operatorname{\mathbf{x}}$, and hence its solutions must be symmetric with regard to simultaneous reflection on all coordinate axes. To numerically solve the PDE model (13) we chose a lattice noticeably larger than the corresponding adaptive dynamics attractor. The number of bins $B$ in each dimension of this lattice is strongly constrained by memory limitations: An efficient implementation requires computing and storing an array of $B^{d}\times B^{d}$ values of the competition kernel $\alpha(\operatorname{\mathbf{y}}_{i},\operatorname{\mathbf{x}}_{j})$ for the pairwise interactions between all pairs of sites $i$ and $j$. With $B=25-30$ to achieve a reasonable spatial resolution, the memory constraint makes the PDE implementation feasible only for $d=2,3$.

The movie in Video A3 shows the dynamics of the partial differential equation model corresponding to the scenarios shown in Videos A1 and A2.

Scaling relationship for the diversity at saturation

The number of clusters at the diversity saturation level, $M_{\sigma,d}$, can be estimated to be proportional to the volume of the available phenotype space with the linear dimension $L$, divided by the volume occupied by each cluster, which has a typical linear size $\sigma$:

[TABLE]

Hence, the following scaling relationships hold:

[TABLE]

where $\sigma_{a}$ and $\sigma_{b}$ denote different strengths of competition, and $C_{\sigma}$ is a constant of order 1 that takes into account the “imperfect packing” occurring when $\sigma$ and $L$ have similar magnitude. Based on this, the equilibrium level of diversity is expected to increase exponentially with increasing dimension of phenotype space (as illustrated Figure 1), and with increasing frequency-dependence (i.e., decreasing $\sigma$). In general, diversity is only maintained if $\sigma\lesssim 1$, which is roughly the scale of the phenotypic range set by the carrying capacity given by eq. (5) in the main text.

1 Specific sets of coefficients used

The following set of coefficients $b_{ij}$ determining the competition kernel were used for Figures 5A in the main text and for the movies.

[TABLE]

The following set of coefficients $b_{ij}$ determining the competition kernel were used for Figure 5B in the main text:

[TABLE]