Competitive dominance in plant communities: Modeling approaches and theoretical predictions
CCompetitive dominance in plant communities: Modeling approachesand theoretical predictions
Jos´e A. Capit´an a,b, ∗ , Sara Cuenda c , David Alonso b a Complex systems group. Department of Applied Mathematics. Universidad Polit´ecnica de Madrid. Av. Juan de Herrera, 6.28040 Madrid, Spain b Theoretical and Computational Ecology Lab. Center for Advanced Studies, Blanes (CEAB-CSIC). C. Acc´es Cala St.Francesc 14, 17300 Blanes, Spain c Universidad Aut´onoma de Madrid. Facultad de Ciencias Econ´omicas y Empresariales. Depto. An´alisis Econ´omico:Econom´ıa Cuantitativa. C. Francisco Tom´as y Valiente 5, 28049 Madrid, Spain
Abstract
Quantitative predictions about the processes that promote species coexistence are a subject of activeresearch in ecology. In particular, competitive interactions are known to shape and maintain ecologicalcommunities, and situations where some species out-compete or dominate over some others are key todescribe natural ecosystems. Here we develop ecological theory using a stochastic, synthetic frameworkfor plant community assembly that leads to predictions amenable to empirical testing. We propose twostochastic continuous-time Markov models that incorporate competitive dominance through a hierarchyof species heights. The first model, which is spatially implicit, predicts both the expected number ofspecies that survive and the conditions under which heights are clustered in realized model communities.The second one allows spatially-explicit interactions of individuals and alternative mechanisms that canhelp shorter plants overcome height-driven competition, and it demonstrates that clustering patternsremain not only locally but also across increasing spatial scales. Moreover, although plants are actuallyheight-clustered in the spatially-explicit model, it allows for plant species abundances not necessarilyskewed to taller plants.
Keywords:
Continuous-time Markov processes, Birth-death-immigration processes, Hierarchicalcompetition, Spatially-explicit stochastic dynamics
1. Introduction
Classical coexistence theory (Gause, 1934; MacArthur and Levins, 1967) assumes that the moresimilar two species are in their niche requirements, the more strongly they will compete over shared re-sources, an idea that can be traced back to Darwin (1859). Ever since, competition-similarity hypotheseshave been at the forefront of theoretical explanations for species coexistence (MacArthur and Levins,1967; Abrams, 1983). This framework predicts that large species differences should be selected during ∗ Corresponding author.
Email addresses: [email protected] (Jos´e A. Capit´an), [email protected] (Sara Cuenda), [email protected] (David Alonso)
Preprint submitted to Journal of Theoretical Biology January 22, 2020 a r X i v : . [ q - b i o . P E ] J a n ommunity assembly to reduce competition. Therefore, trait and/or phylogenetic overdispersion haveoften been regarded as signatures of competitive interactions. However, progress in our understandingof how species differences influence the outcome of competitive interactions (Chesson, 2000; Mayfieldand Levine, 2010; HilleRisLambers et al., 2011) shows that this theoretical framework is too simplisticbecause it disregards the balance between stabilizing and equalizing species differences (Chesson, 2000).Stabilizing mechanisms are based on species differences that cause them to be limited more by their ownconspecifics than by their competitors, favoring species when they drop to low densities, which, in turn,promotes species coexistence. Fitness inequality, by contrast, promotes species dominance over po-tential competitors. In the absence of stabilizing species differences, superior competitors would driveother species to extinction through competitive exclusion. In communities controlled by fitness equal-izing mechanisms, species with similar trait values should be selected through competitive dominance,resulting in high levels of trait clustering even in the absence of environmental filtering (Mayfield andLevine, 2010; Kraft et al., 2015). This theoretical framework suggests that significant trait clusteringat local sites may be a fingerprint of biotic (competitive) interactions controlling the composition ofecological communities.Accurately separating the effect of biotic interactions from environmental filters as structuring agentsof community assembly is not trivial (van der Plas et al., 2015). Despite the undeniable success ofspecies distribution models (Peterson et al., 2011), there is an increasing recognition of the need forsimple, process-based models to make robust predictions that help understand species responses to en-vironmental change (Wisz et al., 2013). However, we still lack clear evidence for the role of bioticinteractions in shaping species assemblages. Studies based on species randomization models have at-tempted to separate the outcomes of competitive exclusion and environmental filtering by assuming thecompetition-similarity hypothesis as a given (Webb et al., 2002; Dini-Andreote et al., 2015). To re-late quantitatively clustering or over-dispersion patterns to competitive interactions, in this contributionwe defined and analyzed a toolbox of stochastic, process-based models designed to describe ecologicalcommunities where competitive dominance is the main driver of species interactions. Here we extendedearlier works on synthetic, stochastic approaches to community assembly (McKane et al., 2000; Haege-man and Loreau, 2011; Capit´an et al., 2015, 2017). Our models are trait-based and competition betweenheterospecifics is determined by a hierarchy in trait values. Our framework was designed to model plantgrowth in the presence of competition for light, although it could be extended to more general settings.We propose two stochastic models that explicitly incorporate trait hierarchies to model competitivedominance. The first one is spatially implicit and uses plant height as a proxy of competition for light.In situations where light is a scarce resource, plants tend to grow taller and out-compete shorter indi-viduals, yielding height-clustered communities formed by tall plants. This might not occur when lightis abundant. A similar effect is to be expected in situations where hydric resources are limited. Our2patially-implicit model can be used to mathematically describe the conditions under which height clus-tering is to be expected in situations where hierarchical competition drives community dynamics. Inaddition, our framework yields two predictions that can be confronted against data. First, the modelpredicts the expected fraction of species that survive, as in Serv´an et al. (2018) but using a stochasticframework of community dynamics. Second, when communities are affected by low dispersal rates, ourmodel yields communities significantly clustered about tall species.Our second model is a spatially-explicit extension of the previous one which also takes into accountthat competitive hierarchies can be traded-off by other alternative mechanisms (such as allelopathy)not related to height. Taller individuals are better competitors for light, whereas shorter individualsdevelop additional strategies to overcome competition, for example, they allocate more energy in alle-lochemics (Givnish, 1982). We explicitly incorporate this trade-off between potential growth and otheralternative mechanisms in the spatial model to show that the model predicts the dominance of eithertaller, mid-sized or shorter plants. A second, testable prediction of this model takes advantage of it beingspatially-explicit: the model predicts that clustering patterns persist for local communities of differentsizes. We finally discuss a number of practical implications of our work, and how the predictions ofour models can be actually tested against ecological community data. Such comparison with macro-ecological data of plant diversity is left to Capit´an et al. (2020).
2. Spatially-implicit plant competition model
Although community assembly results from the interplay between speciation, ecological drift, dis-persal, and selection acting across space and time (Vellend, 2010), our stochastic approach disregardsspeciation, does not account for environmental factors and assembles communities based only on dis-persal, ecological drift and asymmetric competition. Our first model is spatially implicit and considersspecies competitive dominance as measured by trait differences. In origin, this model was devised tostudy hierarchical plant competition for light, hence a functional trait that is usually related to compe-tition is species height. However, we present the model in general terms, without mentioning explicitlyplant heights, so that the theoretical framework can be applied in wider contexts.For all the species in a pool of richness S , a given quantitative, standardized species trait can alwaysbe sorted in ascending order, ≤ t ≤ t ≤ · · · ≤ t S ≤ . Our model assumes that this orderingdetermines a hierarchy of competitive dominance. In order to quantify if species i dominates over species j or vice versa , competitive interactions are chosen as signed trait differences, ρ ij = δ ij + ρ ( t j − t i ) , (1)where the scale factor ρ measures the ratio between inter- vs . intraspecific effects. Here δ ij = 1 if3 = j and otherwise. For this choice of species labeling, ρ ij are positive for i < j and negative for i > j . However, these strengths must be interpreted in competitive terms, i.e., the net effect of a signed ρ ij is opposite to the growth of species i . This means that species i is benefited from the presence ofspecies j if ρ ij < (i.e., if t i > t j ), and conversely, the population size n i decreases if j is present and ρ ij > ( t i < t j ). As a consequence, interactions are hierarchically ordered so that the first species isout-competed by the remaining species and the S -th one out-competes the rest of the community.As shown in Appendix A, the results reported in this contribution are robust to re-defining interactionstrengths as ρ ij = (1 − β ) δ ij + ρ ( t j − t i ) + β, (2)where adding a constant β to off-diagonal entries such that ρ ≤ β ≤ − ρ implies that ≤ ρ ij ≤ (recall that t i are standardized trait values, hence ≤ t i ≤ ). This choice, which makes all strengthsequally signed and positive, is in agreement with the intuitive idea that bi-directional strengths ρ ij and ρ ji must have the same sign to stand for competition. Note also that the effect of a positive ρ ij is to decreasethe growth rate of species i . It is important to notice that the use of signed (1) or unsigned (2) traitdifferences as proxies for competition does not matter as to the predictions derived in this contribution. Community dynamics is mathematically described as a birth-death-immigration Markov process incontinuous time (McKane et al., 2000; Haegeman and Loreau, 2011; Capit´an et al., 2015). Let n i be the i -th species abundance. At each time step, one of the following events affecting n i can take place: (1)immigrants of species i arrive from the pool at rate µ , (2) an individual of species i reproduces or dies atrates α + and α − , respectively, (3) two individuals of the same species i compete at a rate α/K , where α = α + − α − , (4) two individuals of distinct species i and j compete at a rate α | ρ ij | /K , resulting in anincrease of population size n i if ρ ij < and a decrease if ρ ij > .The continuous-time Markov process is completely described through the master equation, a sys-tem of coupled differential equations satisfied by the probability P ( n , t ) of observing a vector n =( n , . . . , n S ) formed by S species abundances. For a birth-death-immigration process, the master equa-tion can be expressed as ∂P ( n , t ) ∂t = S (cid:88) i =1 (cid:8) q + i ( n − e i ) P ( n − e i , t ) + q − i ( n + e i ) P ( n + e i , t ) − [ q + i ( n ) + q − i ( n )] P ( n , t ) (cid:9) , (3)the rate q + i ( q − i ) representing the overall birth (death) probability per unit time, and e i = ( δ ij ) Sj =1 beingthe i -th vector of the canonical basis of R S . According to the elementary processes described above, and4aking into account the sign of competitive interactions, the overall birth rate of species i is expressed as q + i ( n ) = µ + α + n i + αn i K (cid:88) ji | ρ ij | n j , (5)which implies that n i decreases by pure death events, by intraspecific competition or by interspecificcompetition with any dominant species j —which satisfy ρ ij > .The master equation balances two contributions in the rate of variation of P ( n , t ) : the probability ofvisiting state n at time t grows at rate q + i ( n − e i ) due to the elementary transition n − e i → n , and alsoincreases at rate q − i ( n + e i ) because of the transition n + e i → n . On the other hand, the probability P ( n , t ) decreases for all the transitions starting from state n to any other state.In the deterministic limit, we can expand the master equation in terms of species densities x i = n i / Ω ,where Ω is a measure of the volume occupied by the system and is used as a parameter for the vanKampen expansion of the master equation (van Kampen, 2011). Up to leading order in Ω , the stochasticmodel reduces to the Lotka-Volterra dynamics with immigration, d ˆ x i d ˆ t = ˆ x i (cid:18) − S (cid:88) j =1 ρ ij ˆ x j (cid:19) + λ, (6)where population abundances have been re-scaled by K , ˆ x i = x i /K , time is measured in units of α − ( ˆ t = αt ), the re-scaled immigration rate λ = µ/ ( αK ) is non-dimensional, and ρ ij is defined by eitherEq. (1) or Eq. (2).Gillespie’s stochastic simulation algorithm (Gillespie, 1977) is an exact simulation methodology thatcan be used to compute a sample trajectory of the stochastic process, based only on the knowledge ofthe elementary events driving the stochastic dynamics. Instead of directly applying Gillespie’s method,we simulated the stochastic process using an adaptive step (“tau-leaping”) algorithm to speed up simu-lations (Cao et al., 2006). Model dynamics can cause the extinction of some species in the steady-state regime. Hence wecan measure the observed local diversity relative to species pool richness (we refer to this ratio p c , thefraction of species that survive, as “coexistence probability”). Simulation results (see Appendix B for5 .1 1 10 Scaled competitive overlap, 〈ρ〉 S C oe x i s t en c e p r obab ili t y , p c S = 100S = 250S = 500 Figure 1. First prediction of the implicit model.
Average fraction p c of coexisting species (coexistence probability) incommunities as a function of their scaled competitive strength, (cid:104) ρ (cid:105) S . Simulation parameters are α + = 50 , α − = 0 . , µ = 5 and carrying capacity K = 50 (down triangles, inset) or K = 1000 (up triangles). There is a threshold in competitionover which the fraction of species that survive starts declining, and all curves show the same explicit form when averagecompetition is scaled with species richness as (cid:104) ρ (cid:105) S . details) are depicted in Figure 1, which shows full coexistence for values of the average competitivestrength (cid:104) ρ (cid:105) below an extinction threshold. Above that value, as interspecific competition increases,coexistence probability shows a power-law decay whose exponent is controlled by the immigration rate µ . Different pool sizes S yield different decay curves, but these curves “collapse” into a single curvewhen represented as a function of the competitive strength scaled by the species pool richness S , p c ∼ ( (cid:104) ρ (cid:105) S ) − γ . (7)Observe that the curve collapse eliminates the variability in S . This is important because, in practice,presence-absence plant data arising from different geographic regions will be available, and for each plotwe could measure the fraction of species that survive, p c , regarding each geographic region as a differentspecies pool. Therefore, empirical coexistence probabilities, which arise from different species poolsizes, can be fitted together. It is important to remark that, although we already reported the power-lawdecay of p c for non-hierarchical competition in Capit´an et al. (2015), here we show that the threshold incompetition still remains for hierarchical interactions, and the curve collapse is new to this paper.The exponent γ of the power-law decay is determined by immigration, as Fig. 2 shows. Increasingimmigration rate made the exponent γ decrease, in agreement with the fact that high immigration mustincrease the fraction of species that survive. For all the immigration rates shown in Fig. 2a, a power-lawfit to the range where p c < yielded γ < . In a more realistic scenario, where immigration rates weretaken proportionally to species abundances N i in the pool [ µ i = κN i , Hubbell (2001)], the effect ofimmigration still holds: large values of κ produce smaller values of the (unsigned) exponent γ (Fig. 2b).6 .1 1 10 Scaled competitive overlap, 〈ρ〉 S C oe x i s t en c e p r obab ili t y , p c µ = 250 µ = 50 µ = 5 µ = 0.1 κ = 50 κ = 5 κ = 0.1 ba Figure 2. Immigration rate controls the exponent of the power-law decay of coexistence probability. a , as µ increases,the exponent γ of the relation p c ∼ ( (cid:104) ρ (cid:105) S ) − γ decreases. Remaining simulation parameters are: S = 100 , α + = 50 , α − = 0 . , K = 1000 , and σ = 1 . Averages were taken over realizations of the spatially-implicit model. In panel b ,species immigration rates are proportional to theoretical abundances in the species pool, µ i = κN i . Abundances N i weredrawn from a log-series distribution with parameter α = 0 . (Hubbell, 2001). For the sake of comparison, the correspondingcurves for constant µ (cf. panel a, N i = 1 ) are reproduced as dashed lines. The slope is controlled by immigration rates aswell in this case. The stochastic model predicts the existence of an extinction threshold in competition above whichcommunity diversity is strictly smaller than richness’ pool, S [cf. Fig. 1 and Capit´an et al. (2015)].We have shown numerically that coexistence probability curves collapse, for different values of S , whenrepresented as a function of the scaled competitive overlap (cid:104) ρ (cid:105) S . Here we provide analytical calculationsthat support this scaling for coexistence curves in the limit of large pool sizes.We focus on a deterministic community model that illustrates the scaling and whose equilibriumabundances can be obtained analytically. To make equations solvable, we consider a uniform distributionof trait values, t i = i/S , i = 1 , , . . . , S (note that < t i ≤ for all i ). Competitive interactions ρ ij arecalculated as signed trait differences according to (1). Then it holds that ρ ij = δ ij + ρS ( j − i ) .Given a matrix ( ρ ij ) of pair-wise competitive interactions, our stochastic birth-death dynamics yields,in the deterministic limit, the Lotka-Volterra equations (6) for non-scaled species population densities x i , dx i dt = αx i (cid:32) − K S (cid:88) j =1 ρ ij x j (cid:33) + µ, i = 1 , , . . . , S. (8)We here assume, for the sake of simplicity, that the immigration rate µ is equal to zero. This assumptionis not expected to be determinant in the limit of low immigration rates (where most communities operate).7nterior equilibrium abundances of dynamics (8) with ρ ij = δ ij + ρS ( j − i ) satisfy the linear system x i + ρS S (cid:88) j =1 ( j − i ) x j = K, i = 1 , , . . . , S, (9)which can be written in matrix form, M x = K , with x = ( x i ) , = (1 , , . . . , T and M = r r · · · ( S − r − r r · · · ( S − r − r − r · · · ( S − r ... ... ... . . . ... − ( S − r − ( S − r − ( S − r · · · . (10)Here we have defined r := ρ/S .Because of competitive dominance, the hierarchy in traits introduces a hierarchy in competitionstrengths that induces an ordering of equilibrium abundances, x < x < · · · < x S , the more abundantthe species the larger the trait t i is. On the other hand, given that the rhs of (9) is a constant ( K ), densitieswill be proportional to the carrying capacity K . In addition, for r = 0 the system becomes trivial andthe solution of (9) is x i = K , i = 1 , , . . . , S . Therefore, we look for solutions of the form x i = Kf ( r ) [1 + y i g ( r )] , (11)where coefficients y i and functions f ( r ) , g ( r ) are to be determined and satisfy f (0) = 1 and g (0) = 0 to fulfill the requirement x i = K for r = 0 . Since the solution of the system (9) must be unique (thedeterminant of M is non-zero, see Appendix A), by finding non-trivial expressions for y i , f and g wewould have calculated the single, interior equilibrium point of the dynamics.The full calculation of equilibrium abundances is left to Appendix A. The result is x i = 6 K [2 + (2 i − S − rS ]12 + r S ( S − . (12)This expression yields a threshold in competition at which the lowest species density becomes equal tozero. This occurs when x = 0 , which implies the condition ρS = SS − . In the limit of large communitysizes, S → ∞ , the threshold in competition values above which the first species drops to zero abundancesatisfies ρS = 2 . (13)At this point, the first species goes extinct. The community is then formed by S − extant speciesand, after reaching a new steady state, according to (13), the condition for the second species to become8xtinct reduces to ρ ( S −
1) = 2 . Assuming that S − i (cid:29) , iteration of this argument implies that theextinction of the i -th species takes place when ρ = 2 / ( S − i ) . Once the i -th species have gone extinct,the fraction of extant species (relative to the species pool richness) is p c = 1 − i/S . Thus, eliminating i in the two latter relations yields the functional dependence between the probability of coexistence andthe scaled competitive overlap, p c = 23 ( (cid:104) ρ (cid:105) S ) − , (14)valid for (cid:104) ρ (cid:105) S ≥ / (note that, according to the interaction matrix of this particular case, (cid:104) ρ (cid:105) = ρ/ for S (cid:29) ). Coexistence probability decays as a power law (with exponent equal to − ) as a function of thescaled competitive overlap (cid:104) ρ (cid:105) S . Since no extinction occurs for (cid:104) ρ (cid:105) S < / [cf. Eq. (13)], we can write p c = , (cid:104) ρ (cid:105) S < , ( (cid:104) ρ (cid:105) S ) − , (cid:104) ρ (cid:105) S ≥ , (15)which is a continuous function at the extinction threshold (cid:104) ρ (cid:105) S = 2 / . Importantly, the deterministicmodel predicts a collapse of coexistence probability curves in terms of the product (cid:104) ρ (cid:105) S .As shown in the previous subsection by stochastic simulation, the collapse of power-law curves p c ∼ ( (cid:104) ρ (cid:105) S ) − γ is valid even when demographic stochasticity is included. This has an important practicalimplication: the collapse allows us to fit empirical fractions of surviving species arising from differently-sized species pools into a single function to adjust. The deterministic example discussed here predicta power-law exponent γ = 1 in the absence of immigration. Fig. 2 shows that this exponent has todecrease as immigration becomes increasingly important.Because the deterministic prediction (15) for the fraction of species that survive is independent of K and α , we expect a weak dependence on those parameters in the low immigration limit, which can beconsidered as a perturbation to the µ = 0 case, even if stochasticity is taken into account. Therefore,variations on those parameters will not significantly alter the curves depicted in Fig. 1 (the inset, inparticular, shows how the power-law decay for K = 50 is close to that of K = 1000 ). The thresholdin competition that limits coexistence cannot be significantly increased by, for instance, augmenting K .Accordingly, we have checked by stochastic simulation that the power-law curves do not significantlychange for K > in the limit of low immigration, as expected.
In this subsection we compare the species of the realized communities along the stochastic processwith a null model which assumes no species interactions (Triad´o-Margarit et al., 2019). Randomizationtests (see Appendix C for details) produce synthetic community assemblages that would be independentof any selection driven by environmental factors or biotic interactions. Confronting the empirical compe-9
Scaled competitive overlap, ⇥⇤ S S = 100S = 250S = 500 C oe x i s t en c e p r obab ili t y , p c ⌅ ~ . A egeanand W . T u r k e y f o r e s t s A l p s f o r e s t s A t l an t i c f o r e s t s B a l k an f o r e s t s B a l t i c f o r e s t s C an t ab r i an f o r e s t s C a r pa t h i an f o r e s t s C e l t i c f o r e s t s C en t r a l E u r opean f o r e s t s D i na r i c M oun t a i n s f o r e s t s E ng li s hLo w l and s f o r e s t s I be r i an f o r e s t s It a li an f o r e s t s N . A t l an t i c m o i s tf o r e s t s N . E . S pa i nand S . F r an c e f o r e s t s N . W .I be r i an m on t ane f o r e s t s P annon i an f o r e s t s P i ndu s M oun t a i n s f o r e s t s P o B a s i n f o r e s t s S a r m a t i c f o r e s t s S c and i na v i an f o r e s t andg r a ss l and s S c and i na v i anand R u ss i an t a i ga S . W .I be r i an M ed i t e rr anean f o r e s t s T y rr hen i an - A d r i a t i c f o r e s t s W . E u r opean f o r e s t s p - v a l ued i s t r i bu t i on Significant overdispersion Significant clustering
50 1000 p - v a l ued i s t r i bu t i on
50 1000 50 1000
Carrying capacity, K
50 1000 50 1000 µ = 0.1 µ = 1 µ = 5 µ = 10 µ = 25 b c Figure 3. Second prediction of the implicit model.
Model randomization tests for different immigration rates and twocarrying capacity values ( K = 50 and K = 1000 ). Here we chose (cid:104) ρ (cid:105) = 0 . and S = 100 , remaining parameters were α + = 50 and α − = 0 . . p -value distributions of test realizations are shown as Tukey boxplots. The closer the distributionis to , the larger is the fraction of cells where trait clustering is significant. For parameter values yielding low scaledimmigration rates [ µ/ ( αK ) (cid:28) ; this holds, for example, for µ (cid:46) and K = 1000 ], the model indicates a clear signature ofclustering. tition average (cid:104) ρ (cid:105) C for realized species assemblages with the distribution of the competition average (cid:104) ρ (cid:105) Q of the null model yields a probability p = Pr ( (cid:104) ρ (cid:105) Q ≤ (cid:104) ρ (cid:105) C ) . Probabilities significantly close to zero areindicative of trait clustering and probability values significantly close to one reflect trait overdispersion.Figure 3 shows the second prediction of our model. We observe that, in the limit of small (scaled)immigration rates [ λ = µ/ ( αK ) (cid:28) ], realized local communities are significantly more clustered thanexpected according to the null hypothesis (i.e., in the absence of interactions). The larger the carryingcapacity, the more significant is the signal of clustering. As immigration becomes more important, thesignal of clustering weakens and model communities show a broad distribution of p -values. At the high-est immigration rates, model communities are essentially random samples of the pool since immigrationoverrides competition in this regime (Etienne and Alonso, 2005). Therefore, as natural communities areexpected to operate in a low-immigration regime, our implicit model based on competitive dominancepredicts trait clustering.
3. An spatially-explicit extension of the model
It is known that taller individuals are better competitors for light (local shading depresses growth)and show higher colonization potential, while shorter individuals allocate more energy in allelopathiccompounds. Height hierarchies, however, as assumed in our spatially-implicit model, lead to the se-lection species that invest on potential growth (if t j > t i , then ρ ij = ρ ( t j − t i ) > and therefore theabundance of i decreases, hence taller plants are selected).10n this section we extend the model introduced in subsection 2.1 to incorporate spatially explicit localinteractions, as well as alternative mechanisms different than potential growth for plants to resist localheterospecific competitors. Thus the new model allows the study of properties, such as clustering, ata local level, and renders a wider variety of species abundance distributions, not necessarily skewed totaller plants. The system is structured on a hexagonal lattice of size N L where each site i can accommodate atmost a single plant individual of species s i ∈ { , , . . . , S } , with s i = 0 if the site is empty. Thus attime t the system is defined by the state vector s = ( s , s , . . . , s N ) . Species traits are ordered as in theimplicit model, where t s i < t s j if and only if s i < s j .Let d i be the transition probability rate of the death of an individual at site i if the site is occupied,and b ( σ ) i the transition probability rate of a birth event of species σ at site i if the site is empty: d i ( s ) = α − + αK e (cid:88) j ∈N i (cid:2) δ s i s j + ρ ( t s j − t s i )Θ s j − s i + ξ ( t s i − t s j )Θ s i − s j (cid:3) , (16) b ( σ ) i ( s ) = µ e + (cid:88) j ∈N i δ s j σ (cid:18) α + + αK e ρ (cid:88) k ∈N j ( t s j − t s k )Θ s j − s k (cid:19) . (17)Here N i is the set of neighbors of site i , and Θ n is the integer form of the Heaviside step function (definedhere as Θ n = 1 if n = 1 , , . . . and Θ n = 0 if n = 0 , − , − , . . . ). The model incorporates competitiondriven by alternative mechanisms (such as allelopathy), which is based on height differences as in theoriginal hierarchical competition term. The new term for competition is controlled by parameter ξ .Parameters α + , α − and ρ are the same as in the implicit model: α + and α − are the intrinsic probabilityrates of reproducing and dying, respectively, for each individual, while the terms with ρ (and ξ ) accountfor the intensity of pairwise interactions with respect to intraspecific competition.The meaning of K and µ in the implicit model (4) and (5), however, cannot be directly extrapo-lated to the spatially explicit model. In the implicit model, for each individual, its potential number ofcompetitors scales with KS (since each species abundance scales with K and all species can interactamong each other), whereas in the spatially explicit model the total number of potential competitors ofany species, K e , is fixed by the lattice ( K e = 6 in our case). Therefore, in order to use in the explicitmodel non-dimensional immigration rates comparable to the implicit ones, µ/ ( αK ) , we need values ofthe spatial model immigration rate µ e ∼ µK e / ( KS ) .Using expressions (16) and (17), the master equation that describes the continuous-time Markov11rocess can be written as ∂P ( s , t ) ∂t = N (cid:88) i =1 (cid:110) Θ s i (cid:104) b ( s i ) i ( s − s i e i ) P ( s − s i e i , t ) − d i ( s ) P ( s , t ) (cid:105) + (1 − Θ s i ) S (cid:88) σ =1 (cid:104) d i ( s + σ e i ) P ( s + σ e i , t ) − b ( σ ) i ( s ) P ( s , t ) (cid:105) (cid:111) . (18)For ξ = 0 , this model shows a threshold and a power-law decay of the coexistence probability p c as in Eq. (13), and it also reproduces the trait clustering described in section 2.4 (as we show inthe next subsection). The main drawback of this model with respect to the previous one is related tothe simulation of large communities. Whereas in the implicit model computation time grows with thenumber of species, S , in the spatially-explicit model grows with the lattice size, N L . This imposes hugelimits in the maximum number of individuals that the lattice can accommodate, and therefore in the sizeof the simulated communities. This is why we prefer to stick to the implicit model when not dealingwith local properties. The results that can be derived from the spatially-explicit extended model are all related to the clus-tering predicted by the implicit model. The first prediction, new to the explicit model, is related to thepersistence of trait clustering when species are aggregated at different spatial scales. This is importantbecause real individual plants interact at small spatial scales (1 to 1000ha), so local communities haveto comply with this spatial resolution in order for our battery of models to be able to capture signalsof competitive interactions. Our spatially-explicit model can help explain signals of height clustering atdifferent aggregation scales.Randomization tests like the ones described in Appendix C were conducted for different aggregationsizes on the simulated lattice. For that purpose, we divided the lattice into cells using different gridsizes and considered each cell as a local community. To keep sampling efforts comparable at differentcell sizes, we sampled a fixed number of lattice sites ( ) for all levels of resolution, and considered itas a distinct species assemblage. For each sample, we identified the number of distinct species presentand conducted randomization tests considering the whole lattice as the species pool from which speciescan arrive to communities. This way we obtained Fig. 4, which shows that this model was also able tocapture significant levels of clustering when species in the lattice were aggregated into cells of differentsizes.The second result extends the clustering to heights other than the highest ones. Besides height hier-archies, as assumed in our spatially-implicit model, the spatially-explicit stochastic model incorporatesalso alternative mechanisms that trade-off with growth. In this model, potentially taller plants are moreprone to reproduce and contribute to the death of neighboring shorter species, but these shorter species12 igure 4. First prediction of the explicit model. Significant local clustering across increasing spatial scales in modelrealizations. The simulation rectangle ( × sites organized in a hexagonal lattice) is divided into L × L cells ofdifferent sizes for L = 7 , , , , , . Each cell is regarded as a community for randomization tests. Cell size ismeasured as the square root of the number of sites in each cell. The shaded area represents the region where clustering issignificant ( p < . ). For most aggregation scales, the whole p -value distribution falls within the significance region. Modelparameters are α + = 50 , α − = 0 . , µ = 10 − , ρ = 0 . and S = 100 , with ξ = 0 in the left panel and ξ = 0 . in the rightpanel. The number of p -values in each boxplot was constant for the sake of comparison. Figure 5. Second prediction of the explicit model.
Example of species abundance distributions yielded by the explicitmodel, showing that lower, intermediate and higher species can predominate, depending on the relative values of ρ and ξ .Model parameters are: S = 200 , α + = 50 , α − = 0 . , µ = 10 − and ξ = 0 . . Here we used a simulation lattice, preservingthe hexagonal shape, formed by × sites. can also cause the death of taller individuals due to allelopathic effects (as an example of a non-size-related, alternative competition mechanism). Computer simulations show that the balance of these twomechanisms can end up selecting plant sizes clustered around an optimal potential height that can beeither shifted toward lower or higher values depending on the choice of model parameters, as shown inFig. 5. 13 . Discussion In this contribution we proposed a mathematical framework based on height hierarchies to modelplant community dynamics, which we analyzed in full detail to derive a number of theoretical predic-tions, namely: (i) when competition is only considered in terms of height hierarchy, there is a thresholdvalue of the average competitive overlap above which the expected fraction of extant species observedin species assemblages decays as a power-law whose exponent is essentially determined by immigrationrates; (ii) in the limit of low immigration and large carrying capacity, local communities are expectedto be clustered around similar height values; (iii) this clustering significantly remains in local communi-ties of different sizes; and (iv) when competition for light is traded off by other alternative mechanisms(such as the energy invested in allelopathic compounds), the abundances can be clustered around taller,middle-sized or smaller species in realized communities.Our theory represents a strong simplification of actual plant dynamics. Competitive hierarchies areseldom hard-wired. Real plant communities are obviously much more complex, but simple models canbe used to gain valuable insights into the functioning of complex plant communities. A careful descrip-tion of the heterogeneity and variability involved in the complex phenomena determining plant commu-nity assembly, particularly at larger scales, may require a considerable number of detailed variables andmore complex theoretical approaches. However, although attention to detail is essential to science, trueunderstanding of the causal relationships involved in the dynamics of a system is impossible withoutexamining models with only a handful of key aggregated variables that make model predictions andanalysis tractable.The first model, which only deals with species competition in a hierarchical, spatially implicit way,assumes that all individuals can interact with the rest. As apparent from Eqs. (4) and (5), interactions fa-vor the reproduction of individuals belonging to taller species and the out-competition of shorter speciesindividuals, and discourages large populations (compared to K ) of a single species, thus promoting di-versity. This leads to clustering around taller species. We have devised a second, spatially explicit model,which extends the implicit model to a lattice and includes non-hierarchical competition. Spatially-explicit competitive interactions helps us, on the one hand, analyze clustering at different spatial scalesand, on the other hand, unveil species abundance distribution of highly clustered species around small,medium or tall species. The clustering observed in this model at different community sizes is caused bynearest-neighbor interactions solely.All the predictions derived from our theoretical framework are amenable to empirical testing innatural communities. For that purpose, we just need plant diversity data as well as height estimates for allthe species under consideration. Empirical studies are commonly designed to recover species abundancesand to measure functional traits. Ideally, to calculate unambiguously the expected fraction of species thatsurvive (first prediction of our spatially-implicit model), diversity data must contain species presence-14bsence in different locations or regions. In Capit´an et al. (2020) we actually test the predictions of ourplant community models using available data for herbaceous plant communities realized across severalEuropean ecologically distinct regions. The application of our theoretical analysis and the confrontationof model predictions against plant community data shows that our simplified framework can be actuallyused to describe plant diversity across biogeographical scales, and also to unveil signals of competitivedominance in mid-latitude regions. We refer the reader to Capit´an et al. (2020) for further details aboutthis analysis. Code availability
Computer code to analyze data and run the models used to generate our results and support the claimsreported in the manuscript will be made available upon request with no restrictions.
Acknowledgments
The authors thank Mercedes Pascual for her useful comments, and are indebted to Rohan Arthur andHan Olff for their constructive criticism of earlier versions of the manuscript. This work was fundedby the Spanish ‘Ministerio de Econom´ıa y Competitividad’ under the projects CGL2012-39964 andCGL2015-69043-P (DA, JAC) and the Ram ´on y Cajal Fellowship program (RYC-2010-06545, DA).JAC acknowledges partial financial support from the Department of Applied Mathematics (UniversidadPolit´ecnica de Madrid).
Author contributions
JAC, SC and DA conceived the theory; JAC and SC conducted simulations; JAC, SC and DA ana-lyzed results; JAC, SC and DA wrote the paper.
AppendicesA. Equivalence of signed and unsigned trait differences
This Appendix provides a proof that Eq. (12) is the solution, for an arbitrary community size S , ofthe linear system that defines equilibrium densities for the deterministic dynamics (8). We focus here inthe case of positive strengths [ β ≥ , cf. Eq. (2)] and obtain the result (12) as a particular case.At equilibrium, the linear system to solve is (1 − β ) x i + r S (cid:88) j =1 ( j − i ) x j + β S (cid:88) j =1 x j = K (A.1)15or i = 1 , , . . . , S and S ≥ . It can be expressed in matrix form as M x = K , where stands for acolumn vector whose S entries are equal to one, and M = β + r β + 2 r · · · β + ( S − rβ − r β + r · · · β + ( S − rβ − r β − r · · · β + ( S − r ... ... ... . . . ... β − ( S − r β − ( S − r β − ( S − r · · · . (A.2)We first show that the determinant of the linear system matrix M can be computed explicitly forarbitrarily-sized matrices and verifies | M | = (1 − β ) S − [1 + β ( S − − ( S − β ) + r S ( S − / .Let us first introduce the column vectors u j = ( S − j − , j − S, e S − j − ) T , j = 1 , , . . . , S − , where e k = ( δ ik ) S − i =1 , j = 1 , , . . . , S − , are the (row) vectors of the canonical basis of R S − . Note that the( S − j + 1 )-th entry of u j is the one that is equal to 1.It is easy to check that u j are right-eigenvectors of M , all of them corresponding to the eigenvalue λ = 1 − β . The entries of matrix M = ( m jk ) satisfy m jk = (1 − β ) δ jk + r ( k − j ) + β . We use thevector notation m j = ( m jk ) Sk =1 for the j -th row of M . In particular, listing only the entries for columns1, 2 and S − j + 1 , we can write the first, second and ( S − j + 1 )-th rows as m = (1 , β + r, . . . , β + r ( S − j + 1 − , . . . ) , m = ( β − r, , . . . , β + r ( S − j + 1 − , . . . ) , m S − j +1 = ( β − ( S − j ) r, β − ( S − j − r, . . . , , . . . ) . (A.3)Then one can easily check that m u j = ( S − j −
1) + ( β + r )( j − S ) + β + r ( S − j ) = (1 − β )( S − j − , m u j = ( β − r )( S − j −
1) + ( j − S ) + β + r ( S − j −
1) = (1 − β )( j − S ) , m S − j +1 u j = [ β − r ( S − j )]( S − j −
1) + [ β − r ( S − j − j − S ) + 1 = 1 − β. (A.4)Now let k / ∈ { , , S − j + 1 } . If k > S − j + 1 it holds that m k = ( β − ( k − r, β − ( k − r, . . . , β − r ( k − S + j − , . . . ) , (A.5)whereas for k < S − j + 1 the k -th reads m k = ( β − ( k − r, β − ( k − r, . . . , β + r ( S − j + 1 − k ) , . . . ) . (A.6)16he entry at column S − j + 1 coincides in both cases, which allows us to write the remaining rowproducts as m k u j = [ β − ( k − r ]( S − j −
1) + [ β − ( k − r ]( j − S ) + r ( S − j + 1 − k ) + β = 0 (A.7)for k / ∈ { , , S − j +1 } . Hence we conclude that M u j = (1 − β ) u j for j = 1 , , . . . , S − , as stated. Theset of eigenvectors clearly forms a basis of the proper subspace associated to the eigenvalue λ = 1 − β ,which has dimension S − . Therefore, the calculation of the determinant | M | can be transformed into a × problem if the matrix is written in an appropriate R S basis.The matrix for the basis transformation we choose contains the S − aforementioned eigenvectors asthe first S − columns. The two remaining columns are set as convenient linearly independent columns.Using a block matrix notation, the matrix P for the basis transformation is P = A × ( S − U U S − ( S − × , (A.8)where subscripts denote matrix dimensions and sub-matrices are defined as A × ( S − = S − S − · · · − S − S · · · − − , (A.9) n × m is the n × m zero matrix and U n = ( δ i,n − j +1 ) denote a square, anti-diagonal matrix with size n andentries equal to one along the anti-diagonal. Since | P | = − , the columns of P form a basis of R S and,given that the S − first columns are linearly independent eigenvectors of M with eigenvalue λ = 1 − β ,the representation M (cid:48) of matrix M in the transformed basis is M (cid:48) = P − MP = (1 − β ) I S − C ( S − × × ( S − Q × , (A.10) I n being the n × n identity matrix. The sought determinant amounts to determining the sub-matrix Q × since | M | = | M (cid:48) | = (1 − β ) S − | Q × | .The inverse of the basis transformation matrix P can be written as P − = ( S − × U S − U B × ( S − (A.11)17here B × ( S − = · · · S − S − − − · · · − S − S . (A.12)In order to check that (A.11) is the inverse of P , we calculate the product PP − = U A × ( S − U S − + U B × ( S − ( S − × U S − . (A.13)Note that U n = I n , so it only remains to check that A × ( S − U S − + U B × ( S − = × ( S − . It can bedone in a straightforward way.To complete the calculation we simply take into account Eqs. (A.8) and (A.11) into the product P − MP and decompose M in four blocks such that M = M M M M , (A.14)where M is the × sub-matrix formed by the first two columns and rows, M is the corresponding × ( S − sub-matrix, M has dimensions ( S − × and M is the remaining ( S − × ( S − square sub-matrix. After block matrix multiplication, from the resulting product P − MP we identify thelower-right block as Q × = U M U + B × ( S − M U . (A.15)Recalling that, by definition, M = β + rβ − r and M = r − − − − ... ... − S − S − S − S + β ... ... , (A.16)we substitute matrices into (A.15) and finally arrive at Q × = − r S − (cid:88) k =1 k ( k + 1) − r S − (cid:88) k =1 k r + r S − (cid:88) k =1 k r S − (cid:88) k =1 k ( k + 1) + β S − (cid:88) k =2 k S − (cid:88) k =1 k − S − (cid:88) k =1 k − S − (cid:88) k =1 k . (A.17)18fter calculating explicitly of the sums, the determinant of M follows trivially, and reads | M | = (1 − β ) S − | Q × | = (1 − β ) S − (cid:20) β [ S − − ( S − β ] + 112 r S ( S − (cid:21) . (A.18)as anticipated.We now obtain a general expression for equilibrium densities for an arbitrary β ≥ . The solutionformally reduces to (12) when β = 0 . For β > we assume that r ( S − ≤ β ≤ − r ( S − in orderfor the maximum strength to verify β + r ( S − < and the smallest strength β − r ( S − to be positive.This way, intra-specific interactions dominate over inter-specific effects and the dynamical system willbe stable Chesson (2000). In particular, if β satisfies r ( S − ≤ β ≤ min { − r ( S − , − / ( S − } ,so that S − − ( S − β > , the determinant | M | will always be a positive function.Once we have computed the determinant of M , we look for solutions of (A.1) of the form x i = K (cid:48) f ( r, β ) [1 + y i g ( r, β )] (A.19)with the additional assumption (cid:80) Si =1 y i = 0 (to be checked later for consistency). Here K (cid:48) is the solutionof (A.1) in the case r = 0 , i.e., K (cid:48) = K/ (1 − β + βS ) . Compared to the case β = 0 , it can be interpretedas an effective carrying capacity.Observe that, according to Cramer’s rule, we can set f ( r, β ) = | M | . In particular, this non-vanishingexpression (even for β = 0 ) implies that a single solution to the system exists. In order to find g ( r, β ) and y i we sum up all the equations of the system (A.1). It holds K (cid:48) f ( r, β ) (cid:32) − β + βS − r S (cid:88) i =1 i (cid:33) + r S (cid:88) j =1 jx j = K, (A.20)where we have used that S (cid:88) i =1 x i = K (cid:48) Sf ( r, β ) , (A.21)which is a consequence of the assumption (cid:80) Si =1 y i = 0 . From (A.20) we get r S (cid:88) j =1 jx j = K + K (cid:48) f ( r, β ) (cid:20) rS ( S + 1)2 − (1 − β + βS ) (cid:21) . (A.22)Substitution of this expression into the system (A.1) yields x i = K (cid:48) f ( r, β ) (cid:20) rS − β (cid:18) i − S + 12 (cid:19)(cid:21) , (A.23)19hich, for β = 0 , is precisely (12) with g ( r ) = rS = ρ and y i = i − S + 12 . (A.24)As can be easily checked, (cid:80) Si =1 [ i − ( S + 1) /
2] = 0 , consistently with our previous assumption. Thiscompletes the calculation of the equilibrium densities yielded by the deterministic dynamics (8). Observethat equilibrium abundances are strictly increasing, x < x < · · · < x S , as expected (taller speciesdominate).Finally, note that the terms K (cid:48) /f ( r, β ) and rS/ (1 − β ) appearing in the solution (A.23) are alwayspositive within the range ρ (1 − /S ) ≤ β ≤ min { − ρ (1 − /S ) , − / ( S − } , yielding to positiveequilibrium densities as long as the expression in square brackets [Eq. (A.23)] remains positive. Moreimportantly, as for the β = 0 case, the hierarchical structure of equilibrium densities is maintained inthe β > situation. Hence the results regarding the coexistence probability curves for the deterministicmodel remain unchanged and, plausibly, will also be recovered by an stochastic version of the modelwith transitive (and equally signed) competitive interaction strengths. B. Computation of coexistence probability
The simulation methodology used to obtain Fig. 1 was the following. We first fixed the averagevalue (cid:104) ρ (cid:105) = ρ of the competition matrix ( ρ ij ) and drew its entries, in terms of trait values, as explainedin Appendix D.For a given realization of the competition matrix, we simulated a stochastic trajectory along a timespan ∆ t = 10 , which was divided into two sub-intervals. The first ∆ t/ time units were left to reach thesteady state, and no averages were taken. We have checked that, for the set of model parameters studied,the steady state was always reached after the first period, irrespective of the initial condition chosen.During the second time window we obtained uncorrelated measures of coexistence probability (es-timated as the fraction of extant species, s , relative to the species pool size S ), yielding a mean value p c = (cid:104) s (cid:105) /S along the trajectory.To get the curves for coexistence probability reported in Fig. 1 we also averaged p c over in-dependent samples of the competition matrix. Error bars were calculated as the standard deviation ofthe average over these matrix realizations. We computed the observed value of (cid:104) ρ (cid:105) in the simulation asan average over sampled matrices. Error bars for (cid:104) ρ (cid:105) were also calculated as standard deviations overrealizations. 20 . Randomization tests As a null model for plant community assembly, we considered that species in realized communitiesalong the stochastic process were randomly sampled from the species pool, so the null model effectivelyassumes no species interactions (Triad´o-Margarit et al., 2019). We generated a stochastic trajectory ofthe birth-death-immigration model for S = 100 species in the species pool and a value of the averagecompetitive overlap (cid:104) ρ (cid:105) . The competition matrix for the species pool was generated and, after the steady-state was reached ( ∆ t = 5 ), model communities were sampled after regular time intervals ( t s = 0 . ).We refer to each of these sampled communities as ‘local communities’. Given a community C observedalong the stochastic process, using the surviving species we computed the actual mean competitiveoverlap (cid:104) ρ (cid:105) C . Then we randomized it by sampling from the pool ‘synthetic’ communities with the samerichness s as the empirical one. For each synthetic community Q , we measured its average competitionstrength, (cid:104) ρ (cid:105) Q = 2 s ( s − s (cid:88) i =1 s (cid:88) j = i +1 | ρ Qij | , (C.1)where ( ρ Qij ) is the competition matrix restricted to species pairs present in the synthetically sampledcommunity. We took independent samples, which yielded distributions of (cid:104) ρ (cid:105) Q that were wellapproximated by Gaussian functions.For each community C , a randomization test was performed based on independent samples ofthe pool matrix to yield the corresponding p -value. This procedure was repeated times to take intoaccount a series of independent realizations of the pool competition matrix, each one yielding a list of p -values. The complete list of ×
50 = 1000 p -values was depicted as a boxplot in Fig. 3 for differentimmigration rates and carrying capacities. From those p -value distributions one can infer whether modeltrait values exhibit significant levels of clustering, overdispersion, or none of them. D. Competitive interactions generation in model simulations