## Abstract

Theoretical predictions for biodiversity patterns are typically derived under the assumption that ecological systems have reached a dynamic equilibrium. Yet, there is increasing evidence that various aspects of ecological systems, including (but not limited to) species richness, are not at equilibrium. Here, we use simulations to analyse how biodiversity patterns unfold through time. In particular, we focus on the relative time required for various biodiversity patterns (macroecological or phylogenetic) to reach equilibrium. We simulate spatially explicit metacommunities according to the Neutral Theory of Biodiversity (NTB) under three modes of speciation, which differ in how evenly a parent species is split between its two daughter species. We find that species richness stabilizes first, followed by species area relationships (SAR) and finally species abundance distributions (SAD). The difference in timing of equilibrium between these different macroecological patterns is the largest when the split of individuals between sibling species at speciation is the most uneven. Phylogenetic patterns of biodiversity take even longer to stabilize (tens to hundreds of times longer than species richness) so that equilibrium predictions from neutral theory for these patterns are unlikely to be relevant. Our results suggest that it may be unwise to assume that biodiversity patterns are at equilibrium and provide a first step in studying how these patterns unfold through time.

## 1. Introduction

The unified Neutral Theory of Biodiversity and biogeography, proposed by Hubbell [1] aims to provide a theoretical framework for explaining biodiversity patterns across all spatial and temporal scales, from the local and short-term scale traditionally investigated by ecologists to the regional and long-term scale more often contemplated by biogeographers and evolutionary biologists. This neutral theory (abbreviated as NTB hereafter) assumes classically [1] (i) that biotic communities are essentially governed by random population drift (i.e. demographic stochasticity hypothesis), (ii) that all individuals irrespective of species share the same *per capita* birth, death, migration and speciation rates (neutral hypothesis) and (iii) that the number of individuals in the system is constant through time (zero-sum hypothesis).

Although few biologists would accept these hypotheses as factually correct or sufficient to describe in detail how communities are organized and assembled and how they diversify over time (see [2–4] for a summary of the major criticisms), the neutral theory is surprisingly successful at fitting or mimicking a number of well-known biodiversity patterns, such as species abundance distributions (SAD) [5–7], species area relationships (SAR) [8–10], beta-diversity patterns [11,12], spatial distributions [10,13], levels of endemism on islands [14] and temporal variability of populations [15]. Whether or not the NTB accurately fits patterns of biodiversity and biogeography in real communities [16], Hubbell's neutral model is now recognized either as a useful approximation or as a null model [17] that helps to reveal the importance of other processes, such as environmental stochasticity [18–21], demographic trade-offs [22] or competitive interactions [23,24].

In its spatially implicit, classical version, the NTB operates at two loosely defined spatial scales, the local community—where birth, death and migration processes are solely responsible for the dynamics—and the much larger metacommunity—where speciation processes come into play. The local community receives, every generation through migration, a certain proportion of individuals from the metacommunity at large, which is essential to counterbalance the local extinctions that inevitably occur due to random drift. At the metacommunity scale, extinctions are no longer counterbalanced by migration, but instead by speciation events, which can be modelled according to one of several phenomenological modes (detailed below): point mutation [1], random fission [1,25] or peripheral isolate [26,27], the latter serving as a tuneable intermediate between the other two modes. More elaborate speciation modes also exist, such as protracted speciation [28] and those inspired by models of reproductive isolation and genetic differentiation [29–31].

At the heart of NTB is the concept of dynamic equilibrium. A dynamic equilibrium is expected to occur eventually at the metacommunity scale between speciations and extinctions, as both the size of the metacommunity and the *per capita* speciation rates are held constant through time. Virtually, all analytical predictions under NTB have been derived under the assumption of a dynamic equilibrium, by first writing a master equation that describes how the pattern of interest should vary through time, and then solving this equation at equilibrium. At the metacommunity scale, the resulting steady-state SAD is the log-series distribution with a point mutation mode of speciation [1,6] and the broken stick distribution with a random fission mode [25]. After setting the metacommunity SAD at equilibrium, it is possible to derive analytically the local community SAD at equilibrium for a given migration rate [6,32–34]. Analytical expressions for the SAR and patterns of beta-diversity [11,35] typically assume that the SAD, on which they rely, also has reached equilibrium.

Non-steady-state predictions of NTB have rarely been investigated, analytically or numerically. Time-dependent (out-of-equilibrium) predictions are available for the local community SAD, in general, from any starting point [36,37] and for situations near the stationary state [38]. Time-dependent predictions for the expected value and variance of the Simpson diversity index at the local and global scale have also been derived analytically for several spatial configurations [39,40]. The temporal behaviour for other biodiversity patterns has not received any analytical treatment to date and has only tangentially been documented through computer simulations [29,41–43].

While it is analytically convenient to assume that the metacommunity has reached equilibrium, is it biologically reasonable? This question has rarely been investigated by ecologists, but is fundamental for evolutionary biologists; the existence of an equilibrium (‘limit’) for species richness over evolutionary time, and whether this equilibrium has been reached, is still actively debated and questioned [31,44–51]. A dynamical study of how biodiversity patterns emerge through deep time and converge towards equilibrium at the metacommunity scale is thus timely, both in ecology and evolutionary biology. Ecologists are used to talking loosely about ‘the’ metacommunity equilibrium (in the singular), possibly assuming tacitly that all macroecological patterns (e.g. species richness, SAR and SAD) converge simultaneously to equilibrium. Yet, we do not actually know if these patterns converge to equilibrium at the same time, and if not, which ones converge first. We can anticipate that macroevolutionary (i.e. phylogenetic) patterns take longer than species richness to stabilize, for instance because early events in the diversification history gradually disappear from reconstructed phylogenies, the longer the system spends at equilibrium species richness [52]. However, how much longer do phylogenetic patterns take to reach equilibrium compared with macroecological patterns is not well known. Anecdotal evidence exists to suggest that the gamma statistic, an index that reflects whether the branching pattern is concentrated early/late in a clade's history relative to what is expected under a pure birth diversification model [53], may take up to five times longer to reach equilibrium than species richness [52]. By contrast, other studies have suggested that both gamma and measures of phylogenetic imbalance converge to equilibrium ‘shortly’ after Simpson's diversity index (which summarizes the SAD) [54].

The aim of this study is to measure how long it takes for biodiversity patterns at the metacommunity (whole system) scale to reach their respective equilibrium as the system diversifies according to neutral dynamics. This is important to interpret more aptly the biodiversity patterns seen in real clades/metacommunities, and potentially to identify new clues to determine whether the number of species in a clade/metacommunity has reached its equilibrium. We consider four types of biodiversity patterns: (i) the number of species in the metacommunity, (ii) the SAD at the global level, (iii) the entire SAR and (iv) measures of phylogenetic structure (such as branching tempo (gamma), imbalance (beta) and phylogenetic diversity). We investigate by computer simulation (using individual-based, spatially explicit cellular automata) a range of NTB scenarios, varying the speciation rate (high, medium and low), mode of dispersal (global panmictic or local) and mode of speciation (random fission, fixed fission and point mutation).

## 2. Material and methods

### (a) Computer simulations

The metacommunity dynamics were simulated as individual-based spatially explicit cellular automata (coded in C++, available upon request from the corresponding author), on a square grid of 1000 cells × 1000 cells, acting as a torus (i.e. every side is connected to its opposite side) to avoid edge effects. Every cell contained only one individual at a time and was initially populated (at the start of a simulation) with the same ancestor species. Apart from treating space explicitly, our simulations adopted the classical view of neutral theory where all three underlying hypotheses (see Introduction above) were adhered to. We used grids containing a million individuals as this was shown in preliminary tests to be large enough to make the biodiversity patterns considered here behave fairly consistently from one replicate simulation to another (figure 1) while being computationally tractable (each simulation taking between 12 and 36 h to complete). Although much larger grid sizes (numbering billions of cells not merely a million) would be desirable for realism sake, this would have been impractical computationally given that increasing the grid size by an order of magnitude leads to an approximate increase in computational length by two orders of magnitude (if one is interested to track the approach to equilibrium).

Birth, death and dispersal events in the metacommunity were implemented, in a slightly different manner to that proposed by Hubbell [1] but with identical dynamical behaviour, as follows. First, a cell was picked at random as the parent of a single offspring. Second, this offspring dispersed away from the parent according to a given dispersal kernel, either global or local (see below). Third, the owner of the destination cell was then killed off and replaced by the offspring. This triplet of events, birth–dispersal–death, was considered to take place simultaneously and guaranteed that the zero-sum dynamics was maintained at all times. Half a million of those triplets, occurring one after the other, made up one generation (as this is the expected number of triplets that will occur on average between the birth of an individual and the birth of its offspring in a grid of 1 million cells).

With global, panmictic, dispersal, the destination cell could be any grid cell with equal probability. With local dispersal, a random direction was chosen and the distance travelled (away from the parent) followed a Gaussian distribution with a mean of zero and a variance of 5.102. The continuous distance travelled along the *x*- and *y*-axes was then discretized (rounded to nearest integer) to identify the destination cell relative to the parent cell. If the destination cell happened to be the same as the parent cell, the algorithm (whether global or local) was repeated until they were different. Using the local kernel procedure, the average distance travelled by an offspring turned out to be 4.46 cells away from the parent.

After each birth–dispersal–death triplet, a speciation event was triggered with a fixed probability, set by the speciation rate (exact rates given below). Three speciation modes were investigated: point mutation, random fission and fixed fission. In the point mutation mode of speciation, a single individual, the offspring that just dispersed, was instantaneously changed into a new species. In the random fission mode, a random proportion of individuals (uniformly picked between 0 and 50%) from the parent species was allocated to the new species. If the proportion of individuals picked at random happened to be very low, so that not even the most common species in the grid would have enough individuals to allocate at least one individual to both of its daughter species, a new random proportion was picked until this was achieved. With global dispersal, allocation of individuals to the new species was done in a random manner spatially, whereas with local dispersal, the allocation was done so that all individuals were confined to the smallest possible square region of the grid, centred around the offspring cell. In the fixed fission mode of speciation, a new mode of speciation that we introduce here, the process happened in a similar fashion to random fission, except that the proportion of individuals allocated to the new species was fixed to 5%. When, occasionally, the parent species did not have enough individuals to undergo speciation at the required level of fission (random or fixed), another grid cell was picked at random until the species it belonged to fulfilled the requirement. The fixed fission mode of speciation is close to another intermediate mode of speciation, the peripheral isolate mode, where a fixed number of individuals (not a fixed proportion) is allocated to new species [26,27]. We preferred to work with a fixed proportion in order to better evaluate the role of unequal partition of individuals at speciation. Three *per capita* speciation rates—high, medium and low—were studied for each mode. For random fission and fixed fission, these were set at 1 × 10^{−6}, 2 × 10^{−7} and 4 × 10^{−8}, which translated, respectively, into 0.5, 0.1 and 0.02 speciation events occurring on average per generation (in the system as a whole). For point mutation, higher *per capita* speciation rates had to be used (set at 9 × 10^{−5}, 3 × 10^{−5} and 1 × 10^{−5}, which translated, respectively, into 45, 15 and 5 speciation events occurring on average per generation) to make sure that the biodiversity patterns converged to their stationary state in a reasonable amount of time, comparable with that of the other speciation modes. We have assumed, as in classical neutral theory, that all speciation events, once triggered, happened instantly and completely (no incipient stage), which facilitates phylogeny reconstruction (see below). Given that the number of speciation events per generation remains constant system-wide over time, the per lineage speciation rate will tend to decrease over time as diversification proceeds, until it matches the per lineage extinction rate (values are provided in the electronic supplementary material), at which point the dynamic equilibrium in species richness is reached.

Species go extinct in our simulations when the last individual of a species is replaced by the offspring of another species (i.e. through demographic stochasticity). This outcome is quite frequent right from the start with a point mutation mode of speciation because every new species starts its existence with just one individual. With the other two modes of speciation (random and fixed fission), the rate of extinction is initially nil and only starts to kick-in when sufficient speciations have occurred for the abundance of some species to dwindle down to just a few individuals.

A total of 18 different neutral scenarios (3 speciation rates × 3 speciation modes × 2 dispersal kernels) were investigated, and for each scenario 20 replicate simulations were run for 2 million generations (10^{12} triplets) each. The fate of every species in the grid (speciation or extinction) was recorded continuously as it occurred in order to be able to reconstruct the phylogeny at any point in time. Snapshots of the entire grid were also recorded on file at regular intervals (every 100 generations for the first 60 000 generations and every 2000 generations thereafter). This enabled us to measure all macroecological biodiversity patterns considered here in a temporal manner (total of 1570 possible data points) as the system approached equilibrium. All analyses were subsequently carried out in R (v. 3.1.2) [55] using specific packages where necessary (see below).

To validate our simulation approach, we compared our estimates of species richness at equilibrium with those expected for point mutation and random fission, given the size of our metacommunity (*J*_{M} = 1 million individuals) and *per capita* speciation rates (*ν*). For the point mutation mode of speciation, the number of species expected at equilibrium is given by *S*_{eq} = 1 + *θ* ln(1 + [*J*_{M} − 1]/*θ*) [1], where *θ* is the fundamental biodiversity number. Under the Moran model with overlapping generation, our case, *θ* = *J*_{M} *ν*/(1 − *ν*), or more exactly (*J*_{M} − 1)*ν*/(1 − *ν*), to account for the fact that an offspring could not take over the parental cell. Importantly, we do not calculate *θ* as Hubbell [1] originally did (*θ* = 2*J*_{M} *ν*), which only applies to the Wright–Fisher model with non-overlapping generations [56,57]. For the random fission mode, the number of species expected at equilibrium is given by *S*_{eq} = *J*_{M} *ν*^{1/2} [58,59]. Very close agreement (less than 1% discrepancy) was found between the number of species estimated at equilibrium in our simulations with global dispersal and the above theoretical predictions for the random fission and point mutation modes of speciation, with little variation between replicate runs (electronic supplementary material, table S2).

### (b) Measuring biodiversity patterns

Species richness is the only biodiversity pattern that is unambiguously measured by a single metric. All other biodiversity patterns investigated here (SAD, SAR and phylogenetic structure) have a complex shape that is harder to summarize. Whenever possible we attempted to use synthetic metrics that would capture whole biodiversity patterns and thus reliably measure how quickly they approach their stationary shape. We also decided to use several metrics that are pertinent to more specific aspects of a pattern (e.g. SAR behaviour at large scale, abundance of commonest species in SAD) to discern which of these aspects stabilizes first and whether this behaviour is consistent across neutral scenarios. We used the last metric to stabilize, whichever it turned out to be, for each scenario separately, to characterize how quickly the whole pattern of relevance approached equilibrium. We also looked at whether the same metric (for a given biodiversity pattern) was consistently the last to stabilize.

The SAD was studied through rank-abundance plots [60], constructed with logged raw abundances and linear ranks. A direct link exists between rank-abundance plots and empirical cumulative distribution functions used in statistics [61]. To make this link more rigorous, data points sharing the same abundance were collapsed into single data points (with a rank equal to the minimum rank represented and a combined weight equal to the number of collapsed data points). Several metrics were then extracted from these rank-abundance plots to summarize their most salient features: three of these represented the average abundance (after log transformation) at the beginning, at the end and in the middle of the plot; another three described the slope of the plot (how quickly abundance decreased) for the same three plot sections. Two statistics summarizing the whole SAD were also measured: the exponentiated Shannon index (*e*^{H}) [62] and the inverse Simpson index (1/*D*) [63]. These indices (in units of species richness) belong to the Hill series of ecological diversity measures [64] and both express how long and flat the rank-abundance plot is. They differ in how severely rare species are down-weighted compared with abundant species. Their associated measures of evenness *E*_{H} (*e*^{H}/*S*) [65] and *E*_{D} ([1/*D*]/*S*) [66] were also measured.

SARs were calculated using fully nested sets of ‘square’ areas, spanning the entire spectrum of possible areas, from 1 cell up to 1 million cells, evenly spread on a log scale. Four hundred starting points, located regularly across the grid every 50 grid cells in each direction were used for these nested sets (proceeding in a southeasterly direction each time), in order to produce a smooth SAR pattern (in log–log space). The average SAR pattern across the 400 nested sequences was smoothed further, using a cubic spline smoothing approach (with 7 d.f., package stats) combined with an adaptive filtering algorithm [67] with five overlapping windows of 19 data points each, to obtain more accurate slope values. Several metrics were then extracted: three of these described the SAR slope at the beginning (small scale), the end (large scale) and the middle of the relationship. Another metric measured the average species richness observed (geom. mean across the 400 replicates) in the middle of the SAR (‘square’ area of 2000 cells). A synthetic metric that applies to both local and global scenarios was measured that estimates the whole area under the log–log SAR. In addition, two metrics that only applied to local dispersal scenarios were measured: the slope value at the point of inflection in the triphasic pattern of increase and a relative index of spatial structure, which compared the actual SAR (in log–log space) with the one expected if individuals were randomly drawn from the entire grid (instead of being spatially contiguous) (relative index = [area under the rarefied SAR/area under actual SAR] − 1).

Finally, patterns of phylogenetic structure were summarized as follows. For each time point considered for the other biodiversity patterns, the reconstructed phylogeny [68] for the extant species at the time was inferred (using our own R scripts, available upon request, and the file of all extinctions, speciations and ancestral relationships recorded during a simulation run) and the following measures calculated using available R packages: phylogenetic imbalance (beta-splitting parameter, package apTreeshape 1.4-5) [69], the relative position of phylogenetic nodes, aka branching tempo (gamma statistic, package ape 3.2) [70] and phylogenetic diversity (Faith's value, total evolutionary time in generations along all tree branches, except the stem branch that leads to the first branching event in the reconstructed phylogeny, package picante 1.6-2) [71].

### (c) Timing the approach to equilibrium

The generalized Weibull function, a very flexible monotonic growth function with six parameters [72], was fitted to the median timeline of every biodiversity metric (described above) by minimizing the sum of squared residuals. The median timeline simply joined together the median data points across the 20 replicate simulations conducted for every scenario with the same parameter values. For species richness, we also fitted the timeline for each simulation run separately. Achieving the optimal fit with such a flexible function is difficult. A genetic algorithm using derivatives [73] (package rgenoud 5.7-12) was, therefore, used to explore the parameter space thoroughly to get close to an optimal fit. Obtaining a precise value for the time a given metric actually reaches equilibrium is hard to achieve, given the dynamic nature of the equilibrium. Instead, we have measured the time it took to approach within a certain relative ‘distance’ of the equilibrium. After initial testing, we decided to measure the time taken (abbreviated as *t*_{90} hereafter) to cover 90% of the path to equilibrium (i.e. 90% of the differential between the minimum/maximum value close to the start and the equilibrium value, which is estimated through curve fitting). For monotonic trendlines, the differential is simply the difference between the starting value (typically a minimum, in rarer cases a maximum) and the equilibrium value. For non-monotonic trendlines (e.g. gamma statistic), the differential is no longer calculated from the initial value, but from the actual minimum (alternatively maximum) to the equilibrium value. Choosing a higher proportion of the path to equilibrium (e.g. 95%) would land in a flatter, more horizontal part of the growth function, and thus inherently generate less precise estimates. A target of 90%, therefore, was a good compromise, high enough to be close to the actual equilibrium, but not so high as to produce results too imprecise or inconsistent to measure. Times to approach equilibrium are either expressed in absolute units (number of generations since the start) or relative terms (dividing the absolute *t*_{90} for a given biodiversity metric with the absolute *t*_{90} for species richness). Relative values of *t*_{90} make it easier and more direct to compare speeds of approaching equilibrium from one biodiversity pattern to another, and whether the lags observed are consistent across the 18 neutral scenarios studied here. For instance, a relative *t*_{90} of 1 means that the given metric and species richness both approach their equilibrium at the same time (no lag) and a relative *t*_{90} of 2 means that the given metric takes twice as long to approach equilibrium compared with species richness (lag of 100%).

## 3. Results

As expected, the trajectory of species richness through time, the number of species at equilibrium, and the time needed to reach this equilibrium depended on the mode of speciation and speciation rates (figure 1). Examination of the macroecological patterns (SAD and SAR) showed that the shape of these patterns varied through time, sometimes widely, and that they did not necessarily reach equilibrium at the same time as species richness (figure 2). Similarly, phylogenetic patterns varied through time, and took much longer to stabilize than the other patterns (figure 3). Below, we focus on detailing the time it took for these respective patterns to approach equilibrium.

For each speciation mode, increasing the speciation rate increased the number *S*_{eq} of species present at equilibrium in the metacommunity, as expected (figure 4*a*). Under the random fission mode of speciation, this trend is predicted to be a power law with exponent ½ (see Material and methods and [58,59]), which was recovered in our simulations (figure 4*a*). Under the point mutation mode of speciation, the relationship was also well approximated by a power law, although analytical predictions show that it is not strictly a power law (figure 4*a*). The impact of changing the speciation rate was much more important for the point mutation mode of speciation (power-law exponent = 0.9) than for the other two modes of speciation (power-law exponent = 0.5 for both). The negative relationship between speciation rate and the time to approach equilibrium, *t*_{90}, was also conveniently approximated by a power law (figure 4*b*). Here, too, the impact of changing the speciation rate was much more important for the point mutation mode of speciation (power-law exponent = −0.9) than for the other two modes of speciation (power-law exponent = −0.55 for fixed fission and −0.5 for random fission). Changing the speciation rate, thus, approximately had the same effect in amplitude (but opposite in sign) on species richness at equilibrium on the one hand, and on the time it takes to approach equilibrium on the other, irrespective of speciation mode. The number of species present at equilibrium was found to be slightly higher (significantly so based on just 20 replicates) with local dispersal compared with global dispersal (Kruskal–Wallis tests not shown; note the negligible overlap in equilibrial values in the electronic supplementary material, table S2). The time to approach equilibrium in terms of species richness, however, was not significantly different between the two dispersal kernels.

The SAR often took longer to approach equilibrium than species richness and this trend varied considerably between speciation modes. The relative time, *t*_{90}, for the SAR to approach its stationary shape (as evidenced by the last SAR measure to stabilize) compared with the time it took for species richness to stabilize (figure 5) was slightly less than 1 (i.e. no delay), for the random fission mode (0.91, geometric mean of the highest relative *t*_{90} across the three speciation rates), slightly more than 1 for the fixed fission mode (1.07), and 3.8 for the point mutation mode (took almost 4 times longer), all with local dispersal. Scenarios with global dispersal took less time to reach a stationary SAR than the same scenarios with local dispersal, except for random fission (0.89 compared with 0.91 for random fission, 0.96 compared with 1.07 for fixed fission and 2.7 compared with 3.8 for point mutation).

It took slightly longer for the SAD to approach equilibrium than did the SAR, irrespective of speciation mode. The relative time for the SAD to approach its stationary shape (as evidenced by the last SAD metric to stabilize) compared with the time it took for species richness to stabilize (figure 5) is approx. equal to 1 for the random fission mode (i.e. still no delay), 2.25 for the fixed fission mode (took more than twice as long compared with 0–7% longer for SAR) and 5.7 for the point mutation mode (compared with 2.7–3.8 for SAR). However, for the SAD (unlike for the SAR), there is no marked difference in timings between the two dispersal kernels (local and global).

Overall, these results (on SAR and SAD) suggest that the delay that tended to exist between the time it took for macroecological patterns of biodiversity to reach their equilibrium and for species richness to do so got larger the more uneven the split was between sibling species at speciation (i.e. from random fission to fixed fission to point mutation). However, it is worth noting that it was not always the same metric that converged last to equilibrium, and this was contingent on the mode of speciation. For SARs, the number of species halfway through the SAR was often the last metric to stabilize, except for random fission with local dispersal (for which the slope at largest scale converged more slowly) and for point mutation with global dispersal (for which the slope at the smallest scale converged more slowly). For SADs, the metrics that most often converged last were either the Simpson index of diversity or the evenness index based on it. Note also that the metric converging the fastest towards its equilibrium also tended to be contingent on the speciation mode, in some cases it converged even more quickly than species richness itself (electronic supplementary material, tables S3 and S4).

The majority of phylogenetic patterns of biodiversity typically took longer than macroecological patterns to converge towards their equilibrium (electronic supplementary material, table S5), the only exception being the measure of imbalance for the random fission mode that took hardly any time to approach equilibrium (relative *t*_{90} much lower than 1). Of the three measures of phylogenetic structure, imbalance converged the fastest, followed much later by gamma (figure 6). By contrast, phylogenetic diversity did not show any sign of converging, even at the end of the simulations (figure 3), and consequently, its *t*_{90} could not be estimated with any confidence (leaving open the question of whether this pattern would ever reach a stationary value). The relative time it took for phylogenetic patterns to approach equilibrium was also contingent on the speciation mode, but did not follow the same trend shown for macroecological patterns (where random fission stabilized first, followed by fixed fission and eventually point mutation). For imbalance, random fission stabilized first, well before species richness did (rel. *t*_{90} of 0.21 and 0.29 with global and local dispersal, respectively), followed by point mutation (rel. *t*_{90} of 6.43 for global and 6.85 for local) stabilizing shortly after the SAD did (rel. *t*_{90} of 5.70) and eventually fixed fission (rel. *t*_{90} of 26.1 for global and 26.2 for local) stabilizing more than 10 times later than SADs (rel. *t*_{90} of 2.25) (figure 6). For gamma, which took considerably more time to approach equilibrium than macroecological patterns (figure 6), the trend was opposite to the macroecological trend, with point mutation taking comparatively the least amount of time to approach equilibrium (rel. *t*_{90} of 28.8 for global and 62.5 for local), followed by fixed fission (rel. *t*_{90} of 102.1 for global and 96.1 for local) and random fission (rel. *t*_{90} of 114.2 for global and 101.4 for local).

## 4. Discussion

We showed that biodiversity patterns tend to follow a consistent pattern of succession, with species richness converging to equilibrium first, then SAR, then SAD and eventually phylogenetic patterns, often after some considerable time. We showed further that macroecological patterns of biodiversity tended to reach their equilibrium at the same time solely under the random fission mode of speciation and to diverge increasingly in their approach to equilibrium as the speciation mode (first fixed fission then point mutation) partitioned individuals more unevenly between sibling species at speciation. Phylogenetic patterns did not follow the same trend, and for gamma in particular, the trend was opposite to that observed for macroecological patterns: it took longer for this statistic to approach equilibrium, the more even the partition of individuals was between sibling species at speciation. For phylogenetic diversity, no signs of convergence could be detected at all, cautioning against taking any PD value as anything but transient.

### (a) Number of species at equilibrium

Under neutral theory, it is well understood that increasing dispersal will tend to increase the number of species existing at the local scale [1]. Here, we have shown that limited dispersal, while decreasing local scale diversity, tends to increase global scale diversity, as previously shown [10], though in our case the effect is slight enough (3–4% difference in equilibrium species richness, electronic supplementary material, table S2) for the predictions assuming panmictic dispersal at the metacommunity scale to remain a good approximation in terms of species richness. This slight increase in global species richness may be explained by a tendency for the stochastic dynamic of individual species in the system to slow down, as the dispersal kernel contracts. Indeed, the probability for a species to increase (or decrease) in abundance by one individual after each birth–dispersal–death event would tend to get lower as the number of potential ‘competitors’ for a given grid cell decreases (due to local dispersal). These lower transition probabilities in turn would tend to increase the average longevity of species, which coupled with a constant *per capita* speciation rate, would lead more species to coexist.

The link between speciation rate and metacommunity species richness is well known under neutral theory. However, the implications for macroevolutionary diversity dynamics have not been discussed as much. Under neutral theory, systems with higher speciation rates reach higher diversity at equilibrium. Although this is a classical result for neutral theory, in macroevolution the two quantities (species richness at equilibrium and speciation rate) are usually considered independent. This results in two hypotheses to explain species richness [47,48,74]: either being controlled by ‘ecological limits’ (equilibrium hypothesis) or being driven jointly by the speciation/extinction rates and the age of a system (diversification hypothesis). Neutral theory provides an interesting alternative hypothesis, whereby these classical hypotheses are no longer opposed but play together. Here, the upper ‘ecological limit’ set ultimately by the total number of individuals and the speciation rate, is combined with the system's age to predict species richness at the current point in a clade's history.

### (b) Time to approach the species richness equilibrium

Under neutral theory, the speciation rate affects not only equilibrium species richness, but also the speed at which this equilibrium is reached, with faster dynamics associated with higher speciation rates. Our results suggest a useful rule of thumb: for a given speciation mode, if changing the speciation rate leads to an increase in species richness at equilibrium by a certain factor, the time to approach this equilibrium is expected to be reduced by the same factor (i.e. if *S*_{eq} doubles, the time to equilibrium is halved). This rule of thumb works very well for the point mutation and random fission modes, slightly less well (but still useful as an approximation) for the fixed fission mode. This rule of thumb has some interesting implications for our understanding of equilibrium in nature. Under neutral theory, if we assume that ecological systems (metacommunities) of equivalent areas are at equilibrium and closed (negligible immigration from other metacommunities), the magnitude of their differences in terms of species richness, would parallel the magnitude of their differences in terms of time taken to reach equilibrium. However, species richness typically differs more than ecosystem age in nature, suggesting that part of the difference in species richness may be linked to different positions on the trajectory towards equilibrium. Taking boreal and tropical forests as an exemplar, they cover a similar proportion of land, yet their tree species richness (*ca* 161 spp. for boreal, *ca* 44 000 spp. for tropical forests) [75] differ by more than 2 orders of magnitude. Given that the boreal forests originated at most 10 Ma [76], we would have to conclude using our rule of thumb that tropical forests reached their equilibrium number of species in less than 100 000 years (2 orders of magnitude less than 10 Myr) and remained at equilibrium ever since. This unrealistic prediction therefore means *ad absurdum* that boreal forests are far from having reached their equilibrium, as predicted by classical neutral theory.

### (c) Macroecological patterns of biodiversity

When extinctions eventually balance speciations as predicted by NTB, species richness will have reached its dynamic equilibrium. It was previously unknown whether other biodiversity patterns would have reached their equilibrium at the same time. Here we show that they typically do not reach equilibrium simultaneously. The three major macroecological patterns, species richness, SAD and SAR, only coincide in their approach to equilibrium for random fission, a speciation mode shown to fit empirical data less well than other modes of speciation [25,26,58]. Furthermore, we show that the lag occurring between macroecological patterns follows a consistent successional pattern, whereby species richness stabilizes first, SAR second and SAD third. If one assumes that the SAD has reached equilibrium (for instance to fit the local SAD), it is safe to imply that the SAR and species richness would have also reached equilibrium and this can be taken advantage of to infer the equilibrium SAR from the metacommunity SAD at equilibrium for instance [35,77]. On the other hand, only having evidence that species richness has reached equilibrium, is insufficient to tell how close SAR and SAD actually are to reaching equilibrium; they may have reached equilibrium already or still have a long way to go. Future simulation studies on NTB will, therefore, need to pay more attention to the length of their simulations and which equilibrium (species richness, SAR, SAD), if any, their simulations are expected to have reached.

At first it may seem counterintuitive that the SAR may take longer than species richness to stabilize and reach equilibrium (as we have shown for point mutation), even when dispersal is panmictic; after all SAR only describes how species richness decreases with the extent of the area considered. However, it is well known that the shape of the SAR is affected by the pattern of commonness and rarity (i.e. the SAD) in the system [77]. Under global panmictic dispersal (the equilibrium SAR is then concave throughout), the slope at large scale will be affected by the proportion of extremely rare species (such as singletons) and the overall pattern of increase in the SAR (how quickly its slope decreases) will be affected by the evenness of species abundances. Under local dispersal, the ‘triphasic’ shape of the SAR will also be influenced by the SAD, and in addition by the non-random spatial distribution of species across the grid. The time needed for this spatial distribution to reach equilibrium may explain why SARs with local dispersal took longer to approach their equilibrium compared with SARs with global dispersal.

### (d) Phylogenetic structure

Previous work on phylogenetic patterns under NTB has focused on predictions at equilibrium, or after an arbitrary long simulation time [1,54,78]. This has resulted in the widely admitted conclusion that NTB can explain the imbalance (negative beta) but not the branching times (negative gamma) typically observed in real phylogenies [54]. However, here we have shown that phylogenetic patterns take so long to reach equilibrium (at least 20 times longer than species richness in the case of gamma), that the patterns calculated at equilibrium or any arbitrary long time (e.g. [41]) are irrelevant for assessing whether the neutral theory can explain phylogenetic structure of real assemblages. By examining how phylogenetic patterns vary through time, we find here that NTB can explain phylogenetic patterns, at least under some specific speciation modes. The random fission and point mutation modes are unable to produce realistic patterns of phylogenetic structure, with random fission converging very early towards nil beta values (indicative of null, ‘pure-birth’, diversification), and point mutation producing positive gamma values consistently through time (figure 3). Our fixed fission mode, however, resulted in negative values of phylogenetic imbalance and gamma, similar to those typically observed in real phylogenies, around the time the system was reaching its species richness equilibrium (figure 3). Other intermediate speciation modes could also produce such realistic phylogenetic patterns. The peripheral mode, for example, unlike random fission and point mutation, produces realistic species longevities [26,27,58,59] and could potentially also produce realistic phylogenies. The recently proposed mode of speciation by genetic differentiation [31] has also been shown to produce realistic trees under conditions of non-equilibrium species richness. Even non-neutral models including competitive interactions suggest that realistic values of beta and gamma are attained early in the diversification process, before species richness equilibrium is reached [79]. Lastly, the Unified Theory of Ecology and Macroevolution [80] that adds selection and thus fitness differences between individuals to the NTB, can also generate realistic diversification patterns, in the form of decelerating Lineage-Through-Time plots towards the present, provided protracted speciation also occurs.

### (e) Neutral Theory of Biodiversity and empirical data

Our results have implications for the fit of NTB to empirical data. Fitting NTB to empirical data has been useful not only for testing whether the neutrality assumption is reasonable (e.g. [18]), but also for understanding the role of environmental heterogeneity versus dispersal limitation in explaining the spatial distribution of species [81–83], the role of speciation versus dispersal limitation in driving the richness of local assemblages [84], the role of negative density-dependence in local communities [12] and for estimating the size of regional species pools that are compatible with neutral theory and local phylogenies [85].

Such approaches have invariably used equilibrium predictions, such as those for beta-diversity patterns [81], local SADs [84] and phylogenies [85], with few notable exceptions [31]. Given how much these patterns vary along their approach to equilibrium and the repeated evidence that natural systems are not at equilibrium, we argue that previous conclusions obtained by fitting neutral patterns at equilibrium to empirical data may not be tenable. For example, under point mutation, beta is higher at equilibrium than prior to it, and also higher with a higher speciation rate; thus, assuming that the system is at equilibrium, when it is not, results in a significant underestimation of speciation rates (and thus of the fundamental biodiversity number) in ABC approaches that use beta as a summary statistic [85]. Another example: the overwhelming majority of marine ecosystems across the globe had SADs better fitted by a Poisson-lognormal distribution than a Poisson-gamma distribution (the latter being expected under various neutral scenarios) suggesting that non-neutral processes were at play [86]. However, one would need to check that the expected shape of the SAD is Poisson-gamma before as well as at equilibrium SAD, before drawing definitely such conclusions.

## 5. Conclusion

Solid evidence that species richness may have reached equilibrium at the metacommunity scale in real taxa is hard to find whether it is from the fossil record [48] or from analysing extant species phylogenies [46,87,88]. Consequently, evolutionary biologists still argue for and against the existence of a limit (i.e. equilibrium) to species richness and whether it is ever going to be reached [49,50]. Given that real metacommunities may not have reached equilibrium in terms of species richness, it would be unwise for users of NTB to continue assuming that other biodiversity patterns, taking even longer to converge to equilibrium, are themselves at equilibrium. This is most emphatically the case for the metacommunity SAD, the macroecological pattern of biodiversity that takes the longest to stabilize compared with species richness (even assuming an intermediate mode of speciation like fixed fission) and for phylogenetic patterns. This change in paradigm—acknowledging that biodiversity patterns are not necessarily, and most often not probably, at equilibrium—will undoubtedly complicate how neutral theory is fitted to empirical data, but will also ultimately contribute to improving its predictive power.

## Data accessibility

Simulation data are too large (several Terabytes in archive format) to be uploaded. Estimated timings of approaching equilibrium for all metrics and scenarios studied here have been uploaded as part of the electronic supplementary material.

## Authors' contributions

O.M., C.D. and H.M. devised the study. O.M. carried out the research. O.M. coded the simulations with contributions from C.D. O.M. performed the analyses and wrote the first draft of the manuscript. O.M., C.D. and H.M. revised critically the manuscript.

## Competing interests

We have no competing interests

## Funding

O.M. was funded by an EU Marie Curie Intra-European Fellowship (project MEDIATEMP).

## Acknowledgement

We would like to thank Tiago Quental, James Rosindell and an anonymous reviewer for useful comments.

## Footnotes

One contribution of 11 to a theme issue ‘The regulators of biodiversity in deep time’.

- Accepted January 7, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.