Stable Coexistence in a Field-Calibrated Individual-Based Model of Mangrove Forest Dynamics Caused by Inter-Speciﬁc Crown Plasticity

: (1,2) In this theoretical study, we apply MesoFON, a ﬁeld-calibrated individual-based model of mangrove forest dynamics, and its Lotka–Volterra interpretations to address two questions: (a) Do the dynamics of two identical red mangrove species that compete for light resources and avoid inter-speciﬁc competition by lateral crown displacement follow the predictions of classical competition theory or resource competition theory? (b) Which mechanisms drive the dynamics in the presence of inter-speciﬁc crown plasticity when local competition is combined with global or with localized seed dispersal? (3) In qualitative support of classical competition theory, the two species can stably coexist within MesoFON. However, the total standing stock at equilibrium matched the carrying capacity of the single species. Therefore, a “non-overyielding” Lotka–Volterra model rather than the classic one approximated best the observed behavior. Mechanistically, inter-speciﬁc crown plasticity moved heterospeciﬁc trees apart and pushed conspeciﬁcs together. Despite local competition, the community exhibited mean-ﬁeld dynamics with global dispersal. In comparison, localized dispersal slowed down the dynamics by diminishing the strength of intra-/inter-speciﬁc competition and their difference due to a restriction in the competitive race to the mean-ﬁeld that prevails between conspeciﬁc clusters. (4) As the outcome in ﬁeld-calibrated IBMs is mediated by the competition for resources, we conclude that classical competition mechanisms can override those of resource competition, and more species are likely to successfully coexist within communities.


Introduction
Individual-based models (IBMs) simulate the complex life-cycle of individual organisms, including their usage of dynamic resources and their variation among each other within and between life-stages [1]. Frequently, the dynamics of field-calibrated individualbased models of plant communities have deviated from the behavior of a mean field [2], where "individual organisms encounter one another in proportion to their average abundance across space" [3] or, more stringently formulated, "where the interactions are global, that is mortality depends on global density and/or dispersal is infinite" [4]. This has been demonstrated by comparing the dynamics of two models: the original IBM in which interactions among plants and/or the dispersal of seeds were local (in line with natural conditions) and a model in which interactions/dispersal were artificially extended and made global. Pacala and Deutschman (1995) [5], for instance, developed a mean-field model in which the horizontal spatial heterogeneity inherent in the individual-based forest growth model SORTIE was removed and only the vertical heterogeneity among trees was retained. Due to an absence of gaps, their mean-field model maintained only approximately half of the standing crop and it lost diversity in the succession about twice as fast as the full spatial model [5]. In general terms, spatial patterns assumed to indiscriminately arise in the IBMs from local competition and localized seed dispersal have been made responsible for the deviations [2].
This has led to the belief that real plant communities cannot be simulated by meanfield models and that IBMs are necessary to accomplish this task. On the other hand, the persistence of this belief has prevented the contribution of individual-based modeling to the ultimate question of plant ecology: Why do so many plant species coexist in plant communities?

Theories of Plant Competition
For more than 90 years, coexistence and biodiversity theorists have sought to address this question largely by way of mean-field models. In Appendix A, we give an overview of the opposing theories on the coexistence of competing species: (1) the classical competition theory founded on the Lotka-Volterra (LV) model [6,7], (2) the resource competition theory, including (2.1) the resource ratio hypothesis [8] and (2.2) the resource co-limitation theory. Recently, mounting evidence of nitrogen/phosphorus co-limitation (N/P co-limitation), thus supporting theory 2.2, has been found in terrestrial, freshwater, and marine ecosystems from meta-analyses of large numbers of fertilization experiments [9,10].
Within this study, we point out the major difference among the two main theories in terms of the competition for a single resource: In the classical theory, stable coexistence at an over-yielding equilibrium-one that exceeds the carrying capacities of the single speciesis enabled when intra-specific competition is stronger than inter-specific competition. The coexistence is defined here by biotic interactions alone and is, thus, valid for any number of resources. In resource competition theory, two species deplete a common resource until a resource level R* is finally reached, at which only the best competitive species can survive [8]. Hence, in the competition for one resource, coexistence no longer exists. There is, in fact, only one best competitor. We assume here that resource co-limitation theory (2.2) condenses to the R* theory (2.1) in the single resource case (compare [11]).

Recent Developments of Individual-Based Modeling
The work of Adams et al. (2011) [12] has cast doubts on the belief that IBMs and meanfield models are distinct: In their study on Scots Pine, the effect of space (local interaction) on tree density and basal area was minor, the trajectories of density for the mean-field and spatial models were almost indiscernible, and the basal area at equilibrium was only 10% lower in the mean-field model. The authors stated that such small discrepancies between spatial IBMs and mean-field models are commonly observed in temperate forests (Deutschman et al., 1999;Busing and Mailly, 2004, cited in [12]). Larger discrepancies only developed when either the tree density or the neighborhood size was artificially lowered or when the interaction strength was unnaturally raised. After all, the work of Adams et al. [12] opens the way to a contribution of the IBM to the competition theory.
Localized seed dispersal is the second important process that might induce deviations from mean-field behavior in IBMs, because it leads to conspecific clustering [13]. The results of Adams et al. [12] were less clear in this regard. Local dispersal is generally considered an equalizing coexistence mechanism [14][15][16][17] that, by definition, slows the community dynamics but does not alter the competitive outcome among the species. However, not much is known about the mechanisms by which local dispersal prompts this slowing of dynamics.
During the past two decades, individual-based modelers, in particular, the developers of the SORTIE model, have made considerable efforts to replace the old approximation of a tree crown as a rigid cylinder by a new approximation in which a crown is able to change its shape flexibly (cylinder, ellipsoid, and cone) and move laterally [18]. While the most critical biological shortcoming of the SORTIE model could be overcome [19], the modification brought about another dilemma. SORTIE, already complicated prior to the modification, turned into a mathematically intractable model. As a solution, the authors translated all spatial stochastic processes of SORTIE to a macroscopic partial differential equation, the von Foerster equation, using a perfect plasticity approximation (PPA) [18].
Intrigued by the work of Strigul et al., Grueters et al. (2014) [20] aimed to gain insight into the role crown plasticity could play in mangrove ecosystems. The authors developed MesoFON, a new individual-based and field-calibrated model of mangrove forest dynamics that advances beyond current models by describing the crown plasticity of mangrove trees. They applied a simpler approach to crown plasticity than the SORTIE developers. Their routines took advantage of the Fields-Of-Neighborhood (FON) approach, a general-purpose, position-dependent competition index ( [21], see [22] for the concepts of competition indices), and they accounted for the trunk bending and the differential side branch growth mechanism as Strigul et al. did [18], but, in contrast to those researchers, they ignored the flexibility of the crown shape, because this is not as important in mangroves as in temperate forests, which may consist of coniferous species possessing conical crowns and deciduous species owning ellipsoidal crowns.
More recently, efforts have been made to integrate MesoFON into an ecosystem model to be able to simulate the impact of large-scale threats, such as eutrophication and altered sedimentation/hydrology [23], on the wide range of valuable ecosystem goods and services provided by mangroves (compare [24], Appendix A in particular).
As all field-calibrated IBMs, MesoFON is implicitly resource-based. "Implicitly resource-based" here means that the resource (co-)limitation is accounted for in the IBM via mathematical functions, but the uptake of mineral resources from the soil into plants and their return to the soil by the decay of litter is not simulated explicitly mechanistically therein. The competition for above-and below-ground resources is dealt with separately in the MesoFON model assuming resource co-limitation. In the first application of the model, Grueters et al. [20] exposed two plant functional types of the red mangrove (Rhizophora mangle L.), one with plastic and the other with rigid crowns, to two disturbance regimes, namely (1) without disturbance and (2) with hurricane impacts returning every 5 years. The major finding of this simulation study was that disturbances strongly promoted the beneficial effects of crown plasticity. Without disturbances, the overall high competitive strength constrained the effects of plasticity in the dense stands, whereas in the treatment with hurricane impacts, crown shifts were particularly advantageous because of their contribution to gap closure.

Objectives of This Study
By contrast, in this study, we aimed to discern which of the above theories explains the competitive dynamics in the MesoFON IBM. We adopted an explorative case study approach [25][26][27] and designed a critical case [25] toward this end.
The chosen case consists of two equivalent species that are identical in all respects and inherit their entire behavior from Rhizophora mangle. Resource competition theory predicts that a community of two equivalent (or identical) species is only neutrally stable [28] and that it either fluctuates forever [28] or one of the species goes extinct in a random walk (own judgement). Given that MesoFON is implicitly resource-based, we expect competitive dynamics in MesoFON to follow this prediction. As detailed in the Materials and Methods (Section 2.1), the growth, competition, and dispersal attributes of Rhizophora mangle make this a typical case of a pioneer tree species with a fast initial, but later average, growth, shortdistance competitive interaction (due to the medium size of mature trees), and moderate dispersal distance, while we ignore the long tail of the species' dispersal kernel in this study (compare Section 2.1.1). It is emphasized that the growth parameters under the optimal environmental conditions applied here are not affected by the physiological, morphological, and ecological adaptations (salt tolerance, aerial roots, vivipary, and hydrochory) that allow the red mangrove to survive in intertidal tropical regions.
Furthermore, we compare community dynamics when local competition is applied alone (with global dispersal) or when it is combined with localized seed dispersal. We aim to verify that IBM dynamics solely with local competitive interactions, as observed in the field, are in accordance with the conclusions of Adams et al. (2011) [12] and resemble those of a mean-field. We expect that the combination of localized dispersal combined with local interaction would slow competitive dynamics and lead to some form of conspecific clustering. If so, we pay special attention to examine the mechanisms by which this behavior occurs.
In the end, we transform the case study into a critical one by injecting inter-specific crown plasticity, as we call it, into the scene. Toward this end, we establish exclusively a (hard-coded) avoidance of inter-specific neighbors by lateral crown displacement in MesoFON. This should raise intra-specific competition over inter-specific competition. If classical competition theory holds, the two equivalent (or identical) species should coexist at a stable equilibrium [6,7], whereas the community of equivalent species should only be neutrally stable if resource competition theory holds ( [28], see above). The setup of the cases study is graphically illustrated as part of the workflow diagram in Appendix B, Figure A1.
The typical case is already simplistic, even minimalistic, serving to prevent the potential danger of misidentifying mechanisms. Of course, the critical case was ultimately artificial. Even though the transformation into the critical case was justified by scientific evidence from field surveys that crown plasticity is stronger in mixed than in pure stands [29], the artificiality of the critical case requires us to evaluate the implications of our results for field surveys. We provide such a resolution in the final discussion section.

MesoFON-An Individual-Based Model of Mangrove Forest Dynamics
The MesoFON model simulates mangrove forest dynamics as they emerge from the recruitment, establishment, growth, and death of individual mangrove trees, and their spatially explicit above-/below-ground interaction.
In light of the objectives and the community under study, we chose to start each simulation with an empty plot of two hectares. Patch sizes of two hectares and beyond are common among Rhizophora mangle forests, in the Southwest of Florida, for instance [30]. Initially, a number of saplings was randomly placed on the plot, thereby mimicking a severe disturbance event at the beginning. We let the community pass through an initial re-colonization phase until it achieved a steady-state (after 500 annual time intervals, set to 1000 years to be on the safe side), and we maintained it in the steady-state for a further 7000 simulation years. This procedure allowed us to compare the state in which most mangrove communities reside, namely the re-colonization, with the state primarily considered in theoretical models, i.e., the steady-state.
The state of an individual tree is specified in mesoFON by its age, diameter at 1.37 m height (d 1.37 ), diameter growth during the preceding 5 years, x,y-position of its stem base, and the possibly shifted x,y-position of its crown center. An important new state variable is the intra-/inter-specific competition that a tree experiences from its neighbors. Other variables, such as tree height, basal area, and stem volume, are derived from the state variables.
The following series of processes was scheduled for individual trees, to be executed at the conclusion of each annual time interval: (1) tree recruitment, including propagule production, and random global or natural propagule dispersal by water localized around the parental trees, (2) tree growth, (3) growth reduction due to above-ground competition and below-ground competition, (4) lateral crown displacement, and (5) tree death. The processes were implemented as follows:

Tree Recruitment
Recruitment includes propagule production and random global or natural localized propagule dispersal. Under optimal conditions, the density of propagules D (m −3 ) per crown surface area A (m 3 ) is assumed to remain constant during the lifetime of a tree, so the propagule production rate increases with the size of the tree. The density parameter D is set to resemble sapling numbers established per hectare per year, as measured by Chen and Twilley (1998) [31]. Overall, seed production scales less with tree size using Meso-FON, relative to other models [20]. Due to competition, propagule production deviates downward, compared with this optimal behavior.
Irrespective of the dispersal regime (global or localized), the placement of saplings is the sole stochastic process in the model. The localized propagule dispersal by water around parental trees represents the single extension to the previous model [20]; thus, details of corresponding dispersal routines are described in the following discussion. We assume that flowers are primarily produced by terminal branches. Thus, propagules are released exclusively from random x,y-coordinates on the crown margins. After landing on the water surface at the x,y-position of release, a propagule is transported forward by water in steps of 1 m in the direction of the water current that is randomly chosen from 0 • to 360 • for each propagule. The specific behavior is defined by Equation (1): where Pos(x,y) n+1 is the x,y-position of a propagule at step n+1, x is a random number taken from a uniform distribution defined over the closed interval [0,1], P spread = 0.9 represents the spread probability, and → u is a unit vector pointing in the direction of the water current. The dispersal pattern established by this routine mimicks a negative exponential kernel [32]. In order to reduce the number of simulated individuals, trees are established as saplings rather than as seedlings in the model.
With the parameters defined above, the simulation creates a simple geometric distribution of propagule distances, all defined relative to a starting point at the crown margin of a parental tree. The mean and maximum dispersal distances of 100 propagules were computed as 9.09 and 39.91 m, respectively.
In the following, we aimed to assess the realism of the simulated dispersal kernel for Rhizophora mangle and its representativeness for larger parts of the tree flora. In principle, the simulated dispersal distances are in good qualitative agreement with recent findings on the Rhizophora genus. The majority of Rhizophora mucronata Lamk. propagules remained in a 20 m vicinity of the parental tree, and only few were carried more than 65 m away (Chan and Hursin 1985, cited in [30]). In release experiments carried out by Van der Stocken et al. [33], more than 50% of the R. mucuronata propagules were recovered at distances less than 20 m, and only 1% were recovered at distances of more than 90 m. In a further study, Sousa et al. [34] recovered all painted R. mangle propagules within a distance of 8.0 m four weeks after release. Self-planting is the first mechanism likely to be responsible for the observed short dispersal distances [35][36][37]. With the given spread probability, we assumed a self-planting probability of 10% here. By this, we took an intermediary position between researchers that claim a preponderance of self-planting [35] and those who state it as a rare event [38]. Further mechanisms that were made responsible for the short-distance dispersal of Rhizophora mucoronata propagules were the dense aerial root network acting as a physical barrier and the elongated propagule shape, whereas wave action and water flow velocity were identified as important antagonists in the flume tank experiments by Van der Stocken et al. [33].
Certainly, propagules of parental trees at the seaward-side or, more generally, propagules that could escape the aerial root network are prone to long-distance dispersal (LDD). Sengupta et al. [30] demonstrated that, under the influence of oceanic currents, dispersal distances of up to 5 km are most relevant for the supply-side ecology of R. mangle. However, it is not only oceanic currents that affect LDD. Flume tank experiments revealed that, depending on anatomical and morphological characteristics of the propagules, wind also has an important influence on LDD [39][40][41][42]. However, simulating the long tail of the dispersal kernel realistically was not among our objectives here.
On the other hand, the chosen dispersal kernel is representative of larger parts of the tree flora. Seed dispersal by wind (anemochory) and by animals in their coat (epizoochory) are often the most widespread natural dispersal modes [43]. While epizoochorous tree species used to exhibit much longer mean dispersal distances [44], those of anemochorous species are only slightly longer than the one chosen here. The 50% dispersal distance of those species was reported to be at 40 m (Table 1, Dispersal Type 4, [45]). Some anemochorous tree species even have mean dispersal distances comparable to those chosen for R. mangle

Tree Growth
Tree growth is simulated by using the Shugart growth function, which defines the annual increment in the diameter at a height of 1.37 m ∆d/∆t as follows [46,47]: where d is the diameter at a height of 1.37 m (d1.37) (cm) of a focal tree, d max is the maximum attainable d 1.37 in (cm), h is the tree height (cm), h max is the maximum achievable height (cm), G, b 2 , and b 3 are three auxiliary variables, (∆d) max denotes the maximum initial d 1.37 growth rate of a sapling, while f red represents the reduction in d 1.37 increment due to the competition for (co-)limiting resources in the environment. Equation (2) can be classified as a saturating function that, in the optimum case, asymptotically approaches d max and h max . It couples the annual diameter and the height increment to generate the overall growth for the year. The auxiliary variables b2 and b3 ensure the asymptotic behavior. The parameters given by Chen and Twilley [31], namely b 2 = 77.26 and b 3 = 0.396, do not respect this definition. Therefore, with the help of the reported d max = 100 cm and h max = 3000 cm [31], they were corrected here to values of 57.26 for b 2 and 0.2863 for b 3 . The parameter G = 267 was chosen in accordance with the value given by Chen and Twilley [31].
Because of the many morphological and physiological specialties common to true mangroves [38], one might consider the red mangrove to be a special case, like an outlier among tree species. However, in the subsequent paragraph, we aim to convince this kind of reader that the opposite is true when it comes to Rhizophora mangle growth. To this end, we compared its growth parameters with those of 279 tree species referred to by Shugart [47] in Tables 4.1-4.7. As shown in the box-and-whisker plots of Figure A2 in Appendix C, d max and h max of R. mangle were slightly above the medians of the other tree species, whereas b 2 and b 3 were slightly below the respective medians. Only the growth parameter G of the red mangrove lied outside the inter-quartile range of the set of tree species. It exceeds the upper 25% quartile, but it is still far from being an outlier. In conclusion, R. mangle is initially fast growing, but it is an intermediate species in terms of growth later on.

Growth Reduction Due to Competition
The "Field of Neighborhood" approach (FON, [21]) is used in MesoFON to simulate the competition for resources. In principle, the FON represents a position-dependent competition index that, properly parametrized, might resemble the behavior of other such indices (see Pretzsch [22] for an overview of commonly applied competition indices). The pursuing outline of the FON concepts follows that given by Grueters et al. [20].
The competitive strength exerted by a tree is presumed to decline exponentially from the trunk toward the finite margin of the circular FON. At radius r (m), it is given by where r 1.37 is the radius at a height of 1.37 m (m), R FON is the radius (m) of the FON, I max is the maximum intensity of competition (inside the trunk), and I min is the minimum intensity of competition at RFON, defined as a fraction of I max .
The rise in the FON radius (m) with increasing tree size is governed by an allometric relationship that (in the case of R. mangle) resembles crown radius allometry: Moreover, the local competitive strength at an x,y-location is found by adding the intensities of all n FONs that have nonzero intensity at x,y and overlap with the coordinate.
After all, the competitive influence F k A of n neighboring trees (n = k) on a focal tree k is calculated by summing the area-based integrals of all n FON overlaps and dividing each by the integral of the focal FON to normalize it: where A k is the FON area of the focal tree and da is the area variable over which we are integrating, and F k A is cut to a range of [0,1] in case the value calculated by Equation (6) falls outside this range.
The growth reduction factor then consists of the following components: In this equation, F k A, above encompasses the growth reduction due to the competition for light above-ground, while F k A, below covers negative effects of an over-investment in root biomass under the influence of below-ground competition ( [48], p. 105). The competition routines, in fact, are implemented via a duplicated Field-Of-Neighborhood approach. In order to simulate the behavior of size-asymmetric above-ground competition for light [49], we assumed that the competitive intensity in an above-ground FON decreases strongly with radius. In contrast, the commonly believed size-symmetric below-ground competition [20,49,50] is approximated by choosing a constant competitive intensity over the below-ground FON. With the chosen parameters, we assigned equal weightings to abovevs. below-ground competition. However, current parametrization efforts on Rhizophora apiculata suggest that the chosen weighting overestimates below-ground competition. Equation (7) further contains environmental controls, such as soil phosphorus (f red, Nut ) and soil pore salinity (f red,Salt ). The influence of salt and phosphorus levels on tree development was calculated, as described in Chen and Twilley [31]. We chose salinity and phosphorus levels that ensured a maximum possible red mangrove growth, which is characteristic for the subtidal, lower, and middle intertidal zones along shorelines and estuarine mouths in the Southwest of Florida [31]. As competition for the below-ground resource phosphorus was bypassed in this way, light remained the single resource for which trees compete in this study. At this point, it should be recalled that light has also been treated in the resource ratio theory as a limiting resource by defining the light intensity L* out measured below the canopy (in a monoculture) as an analog to R* [51,52] (compare Appendix A).
In conclusion of Sections 2.1.2 and 3, our species-and environment-related choices enabled us to study the dynamics of a typical pioneer tree species under conditions where it realizes its fast growth potential.

(Inter-Specific) Crown Plasticity
In the MesoFON model, trees can jointly shift their crown and above-ground FON positions away from the most severe competitive pressure produced by neighboring trees [20]. In the model version applied to this study, we have hard-coded the neighborhood avoidance to occur only among heterospecific trees (see Figure 1). Thus, trees express only inter-specific crown plasticity. This decision was made to reduce inter-specific competition and to transform this explorative case study into a critical case in which either classical competition or resource limitation theory holds (compare Introduction).

Figure 1.
Graphical illustration of inter-specific crown displacement between Rhizophora spec. 1 (red tree pair on the left) and Rhizophora spec. 2 (green tree pair on the right) due to trunk bending. Above-ground FONs are depicted as transparent blue cones on top of colored discs, while below-ground FONs are depicted as colored discs at the bottom. Two mechanisms responsible for crown displacement, trunk bending and differential side branch growth, are included within the simulation. Each is controlled by a scaling coefficient. Larger scaling coefficients reduce the costs of lateral movement, and they enable a species to escape more rapidly from the negative effects of inter-specific competitive pressure. Both coefficients were chosen to generate displacement patterns in accordance with those reported for species similar to R. mangle [20].

Tree Death
In the model, tree death occurs when tree growth is severely affected by competition or when it slows at the end of a tree's life (natural death). For tree death to occur, the average diameter growth must fall below the species-specific threshold of 0.05 cm yr −1 , over a period of five consecutive years.
A more thorough description of the model following the ODD protocol (ODD = overview, design concepts, detail; [53,54]) is provided by Grueters et al. [20].

Conducting MesoFON Simulation Experiments
In this study, we simulated the behavior of two hypothetical equivalent red mangrove species, Rhizophora spec. 1 and Rhizophora spec. 2, both of which behave like Rhizophora mangle (L.). In one type of simulation run, we invoked inter-specific crown plasticity alone (combined with random global propagule dispersal). In this treatment, trees avoided neighbors belonging to the other species and moved their crowns away from inter-specific competition. In the second type of simulation run, we combined this behavior with localized propagule dispersal.
In order to differentiate between asymptotic stability (classical competition) and neutral stability (resource competition), we followed the guidelines of Grimm and Wissel [55] and investigated two stability properties on the community scale: "staying essentially unchanged" (constancy [55], reliability [56], or "temporal stability" [56,57]) and "returning to the reference state after a temporary perturbation" (resilience [55,56,58]) of stem volumes for Rhizophora spec. 1 and 2. As the two species are identical, we assumed that the reference state, i.e., the potential equilibrium, occurs at equal shares or equal stem volumes of the two species. Both kinds of simulation runs were initialized with equal sapling numbers of their respective species. Then, constancy, the lack of change in spite of minor perturbations [56] (here, due to demographic stochasticity), was measured by the mean/standard deviation of the stem volumes (as in [57]), both computed over the steady state phase of the simulation runs.
Following Wilson et al. [56], we defined the term "resilience" as the rate by which a system returns to an equilibrium following a pulse perturbation. In our case, the system is the community consisting of the two identical species, the equilibrium is the potential reference state referred to above, and assigning a fraction of Rhizophora spec. 1 trees to Rhizophora spec. 2 (or vice versa) in a community at the reference state comprises a pulse perturbation of the community's species composition.
This kind of perturbation is notably different from the more practicable disturbancerelated approach in the form of manipulative removal experiments commonly applied in field ecology, but it is readily available in IBMs and is ideally suited to examine the stability of the species composition in a community according to the definition given by [56]. For the sake of simplicity, we considered unequal initial sapling numbers (735/15, see below) as a displacement from the reference state. The recovery time began after the relatively short time interval to reach canopy closure and overcome initial oscillations of standing stocks and ends when the stem volumes of both species returned to the reference state. The resilience RL (m 3 yr −1 ) was calculated according to [59] using the following equation: where V e,o denotes the stand-based stem volume of Rhizophora spec. 1 or 2 in (m 3  The sole stochastic process in the model runs is the placement of the saplings. All model runs were initialized with 750 randomly placed saplings per hectare (ha) on 2 ha plots and with a constant homogeneous environment. Each type of run was repeated twice. This completes the description of the ecological situation [55] for which the stability properties are valid.

Parameterization of the Lotka-Volterra (LV) Model from Simulation Experiments
The design of simple individual-level models that correspond, in the long-range limit, to the classical Lotka-Volterra model of competition has been the most common approach to provide evidence for the importance of short-range interactions and dispersal [2,4,60].
For example, Neuhauser and Pacala [4] applied a continuous-time Markov process to an integer lattice that is occupied either by "species 1" or "species 2" on all lattice points and that has no vacancies (i.e., the high-density limit). The community dynamics on the lattice were governed by (1) a mortality process that is dependent on the conspecific density (multiplied by an intra-specific competition coefficient) and the heterospecific density (multiplied by an inter-specific competition coefficient), and (2) a reproductive replacement of death that is proportional to the densities of the two species weighted by their specific fecundities. A shift from local to global processes was achieved in the model by including more and more distant individuals as neighbors of a focal tree in respective density calculations. Overall, this approach is elegant, because the operation on the high-density limit avoids complications related to intrinsic growth rates (by ignoring them) and carrying capacities (by simplifying them to the total number of lattice points for both species).
The comparison of field-calibrated individual-based models, which simulate birth, growth, and death processes, and their likewise field-calibrated counterparts in the longrange limit (i.e., mean-field models) served the same objective as above [2,5,12,61]. However, with the exception of Pacala and Silander Jr. [61] (cited in 2), the calibrated mean-field models did not correspond to the Lotka-Volterra equations. Gathering LV parameters from the field, on the other hand, has turned out to be difficult (e.g., [62]) and is generally considered to be tedious, because single species, as well as neighbor removal experiments, in natural communities or mutual invasion experiments in spatially designed communities are required for this [63,64].
We avoided these difficulties here by following Pacala and Levin [2] and we adopted the approach to estimate intra-/inter-specific competition coefficients and other LV parameters from (manipulative) experiments in the IBM. Carrying capacities and intrinsic growth rates are readily obtained from single-species simulation experiments. Competition coefficients for the LV-model are readily available from community simulation experiments conducted with IBMs, because individual plants record competition from intra-and inter-specific neighbors at each time interval in these models.
Basically, the settings selected for the individual-based simulations should ensure that we achieved the simplification of environmental assumptions of the LV model in the long run (spatial and temporal homogeneity, no larger trophic structure, [28]). The biotic assumption that organisms are identical through time [28] is likely satisfied while the populations are in the steady-state and exhibit a constant size and age structure. However, the second biotic assumption of homogeneous unstructured populations, the so-called mean-field assumption, has to be carefully tested in the IBM simulation runs. Additionally, the approach as a whole requires us to critically evaluate the matching/scaling of individual-based processes contained in MesoFON and those contained in the LV model. We provide such an assessment in the following discussion.
Single-species experiments: We conducted single-species simulation experiments with the MesoFON model to parameterize the monospecific logistic growth curve of stem volume integrated over time t: where V i is the stem volume of species i at time t, V 0 is the initial stem volume for that species, K i is the species-specific carrying capacity, and r i is the stem volume growth rate. All former entities are defined in units of m 3 ha −1 , while the last is given in m 3 ha −1 yr −1 .
In principle, because of its sigmoidal form, the Shugart function that governs the stem diameter growth of an individual tree in the MesoFON Equation (2) might scale up to a logistic population growth. Indeed, stand volume trajectories of the JABOVA model (which also relies on the Shugart function) has repeatedly shown compliance with logistic growth, but other forms of population growth have also been reported [47].
Two separate simulation runs were conducted using the single species: (1) one underlying random global and (2) one underlying natural localized propagule dispersal. Equation (9) was fitted to the initial colonization phase where saplings grew almost in synchrony. We used the probabilistic optimization routine of the simulated-annealing method SANN within the optim-function of the statistical software R to fit the data [65]. We adjusted two of the three SANN control parameters following recommendations given by Cortez [66]. The initial temperature T was raised to 1000 to ensure that SANN operates as a global parameter search method. The maximum number of iterations was increased to 100,000 in order to enable convergence in the localized seed dispersal treatment. The resulting models were evaluated visually because SANN does not produce a goodness-of-fit measure (other than the residual sum of squares RSS).
Additionally, we recorded long-term mean values of the stem volume (K i ), the population size (K i,indiv ), and the average individual-based intra-specific competition (α ii,indiv ) from the MesoFON simulation runs. These annual recordings commenced in simulation year 1000 and are representative of the steady-state phase of population dynamics (with constant tree size structure). The multiplication of α ii,indiv by K i,indiv and the division by K i produces the volume-specific α ii for the monoculture, which can then be compared with α ii , values achieved in the community.
Community experiments: Additionally, we carried out simulation experiments with a community consisting of both species. The purpose of the experiments was to verify that the two identical species can stably coexist in the individual-based model. The goal was to parameterize the intra-and inter-specific competition coefficients using annual competition recordings of the individual trees in the IBM. Separate MesoFON simulation runs were conducted such that the underlying community was formed by random global or natural localized propagule dispersal. In order to cover the full range of species compositions, without risking extinction of one species by chance, we initialized the runs with 15/735 and 735/15 saplings per hectare with Rhizophora spec. 1 and 2. We restricted data collection to the steady-state phase of the community dynamics (>500 yrs) in those experiments.
Before we outline the methodology used for data analysis, we introduce the set of differential equations that comprises the competitive LV-model, commonly based on population sizes [67][68][69]. We avoid difficulties related to size differences between species by defining the LV-model in terms of the stem volume per hectare as a proxy for the most suitable "common currency," i.e., biomass per hectare. As Lotka (1925) [6] stated, users are free to switch units in the equations. The logistic growth equations for the interaction of species 1 and 2 are as follows: where: 1. r 1 and r 2 denote the intrinsic volume-specific growth rates for species 1 and 2 in m 3 ha −1 yr −1 , respectively; 2.
V 1 and V 2 represent the respective stem volumes in m 3 ha −1 ; 3. K 1 and K 2 are the respective carrying capacities, also given in m 3 ha −1 ; 4.
α 2←1 is the inter-specific competition coefficient of species 1 on species 2; 6. α 11 and α 22 are the intra-specific competition coefficients of species 1 and 2, respectively.
Similar to the simple IBM of Neuhauser and Pacala [4], all processes included in our model under the chosen scenario conditions, namely fecundity, consisting of propagule production and sapling establishment, tree growth, and death from competition (except for natural death due to tree age), are density-dependent in a simple manner, and, thus, MesoFON is supposed to correspond in the long-range limit to the Lotka-Volterra model. The question of whether the distances of competitive interaction and/or propagule dispersal effective in the field-calibrated model are sufficiently long-ranged is addressed in this study.
The derivation of competition coefficients from simulations is described as follows: MesoFON provides the intra-and inter-specific competition that an average individual tree experiences as a derived system state variable. This normalized competition reflects the fractional growth reduction due to competition exhibited by the average tree. We multiply these individual-based competition coefficients by the stem volume of the respective species to produce the population-wide competition coefficients. Each of these coefficients was subsequently plotted against the stem volume of the species that affects that coefficient, and we then applied linear regression techniques to determine each volume-specific competition coefficient used in Equations (10) and (11). The linear regression allowed us to validate whether the mean field assumption that underlies the LV-model was satisfied by the IBM. The mean field assumption states "that individual organisms encounter one another in proportion to their average abundance across space" [3]. We verified the assumption here by confirming that those population-wide competition coefficients (e.g., α ii · V i ) were proportional to the stem volumes (e.g., V i ) along the respective isoclines. Due to the large samples we took in the simulation runs (≈2 × 7000 simulation years, ≈600 trees), we refused to apply statistical tests on the parameters of the linear regressions in the verification (compare Lin et al. [70] for statistical consequences of the large sample problem). Instead, we compared the robust regression having a nonzero intercept with a robust regression through the origin, and we assumed the mean-field assumption to hold when differences among the regressions were hardly visible. The robust regression was chosen to remove the influence of several outliers recorded outside the isocline range. Additional normal linear regression was used to obtain an R 2 goodness-of-fit measure.
The four possible competitive outcomes of the two species case are commonly determined graphically by the cardinal points (intersections with the axes) of the two speciesspecific isoclines, in the phase-plane diagram graphing V 2 vs. V 1 . The zero growth isoclines are identified by setting Equations (10) and (11) to zero.
The "stable coexistence at equilibrium" occurs when (K 1 α 11 /α 1←2 ) > K 2 and (K 2 α 22 /α 2←1 ) > K 1 . As we deal with identical species in this study and as we can assume that K 1 = K 2 , these inequalities can be simplified to: α 11 > α 1←2 and α 22 > α 2←1 ; these relations are interpreted as the intra-specific competition must be greater than the inter-specific competition for the stable equilibrium to exist.
This scenario arises when inter-specific crown plasticity lowers the inter-specific competition sufficiently. On the other hand, if the species are identical, compete only for a single light resource, and all LV assumptions hold, we are on the nonspatial single-resource hyperbola in the LV parameter space [60] at the point where all competitive coefficients equal 1, and the scenario is, thus, only neutrally stable [28].
We present all results of the LV-analysis by way of specially adapted and highly illustrative phase-plane diagrams (Figures 5 and 7). The R-script used to draw these diagrams was modified after Poisot (2014) [71]. The black arrows located at grid points of the diagrams point in the direction of the local community shift, while the color scheme ranging from yellow to red indicates ascending rates of community change. Isoclines for the respective species are added to the diagram in the traditional way. Moreover, original data points are depicted as a smoothed scatter in the background. If the shape of the smoothed scatter suggested the existence at a classic equilibrium point, we fitted a multivariate adaptive regression spline (MARS) to the original data because MARS supports a bending in the data. The adaptive regression spline was retained if the value of its adjusted R 2 was superior. Otherwise, a major axis regression was fitted. Multivariate adaptive regression splines were constructed using the earth-package in R [72], whereas major axis regression models were constructed using the smatr-package in R [73]. In any case, graphs of the fitted model were added to the phase-plane diagram. When we found evidence for stable coexistence, we also plotted the means and standard deviations of the stem volumes at equilibrium (derived from a run with initial sapling numbers of 325/325 per hectare).
A summarizing workflow diagram of this study is provided in Appendix B, Figure A1.

Community Dynamics in the Individual-Based Model
The two identical Rhizophora species that were designed to avoid each other's competition were able to stably coexist in the individual-based model regardless of the dispersal regime. This is shown clearly in the stem volume trajectories of the runs (compare Figure 2).
In runs initialized with equal sapling numbers, the computed stem volumes were (1) V 1 = 309. 36  As expected, localized dispersal did reduce the speed of community dynamics. Apparently, much more time was required after pulse perturbation to return to the equilibrium in the localized dispersal regime (≈2460 yrs), relative to the random dispersal treatment (≈540 yrs) (compare Figure 2). Likewise, the resilience has been estimated to about 1.83 × 10 −3 m 3 yr −1 in the global dispersal treatment compared to about 0.41 × 10 −3 m 3 yr −1 with local dispersal. Of course, the question arises regarding which mechanisms caused this difference. We hope to clarify this question with the Lotka-Volterra parameterizations in the next sections.

Lotka-Volterra Parameters from Single-Species Experiments
This type of experiment was used to parameterize V 0 , r i , K i , and the intra-specific competition coefficient α ii . Table 1a presents the results of the respective model fittings. A stem volume trajectory is displayed in Figure 3a. It is characterized by an initial recolonization and an overshoot in both, of which the stem volume comprises even-aged trees, followed by a damped oscillation wherein the age similarity fades away and a steady state develops after 500 or 1000 simulation years, at the latest, in which tree size structure is constant and oscillation is absent except for some demographic stochasticity. We subsumed by definition the overshoot and the damped oscillation under the term re-colonization, and we differentiated this from the steady state.  (2)   Diagrams of the fitted models are shown in Figure 3b,c. The growth rate values r i , 5.119 and 5.201 m 3 ha −1 yr −1 , of the two dispersal regimes were very similar. In contrast, the Rhizophora species carrying capacity with localized dispersal was 16.7% larger than the carrying capacity with global dispersal. The visual appearance of Figure 3b and the low residual sums of squares (RSS = 6.76) indicate that Equation (9) fitted the data in the random dispersal treatment very well. Based on the visual appearance of Figure 3c, the large value of V 0 = 31.93 m 3 ha −1 , and the associated large RSS = 46.76, we can conclude that the logistic curve overestimated the low stem volumes and fitted the growth in the localized dispersal (initially slow, and then rapidly accelerating) in a less accurate manner.
The carrying capacities, recorded as long-term averages and representative of steadystate population dynamics, were much lower than those attained in the colonization phase, as described above (compare Figure 3a). The reductions in the random and the localized dispersal case were 17.4% and 29.7%, respectively (see Table 1b). As the competitive outcome unfolds in the steady-state phase, we restricted our LV-analysis to the carrying capacities realized in that phase. The only parameter that we retained from the initial colonization was the growth rate r i , which does not influence competitive outcomes.

Lotka-Volterra Parameters from Community Experiments with Global Dispersal
Community experiments were conducted to determine values of the missing parameters α 1<−2 and α 2<−1 , to achieve a full calibration of the Lotka-Volterra model. Even though these parameters were previously computed in the single-species experiments (compare Table 1b), we extended our analysis here to the intra-specific competition to confirm that the parameter values did not change relative to the values previously determined. As outlined in the "Materials and Methods" section, we obtained a volume-specific competition coefficient from the slope of a robust linear regression of the plot "populationwide competition vs. respective stem volume". Table 1c  With values of R 2 ranging between 0.916 and 0.937, the linear models fitted the data very well. The intercepts of the linear regressions were very small, with the exception of α 1←2 . As a result, differences between the robust regression with nonzero intercept, and the robust regression running through the origin, were hardly visible (compare Figure 4a,c,d), again with the exception of α 1←2 where the deviation was moderate but clearly visible (see Figure 4b). Hence, the visual inspection revealed that the mean field assumption was approximately achieved in the individual-based model. Therefore, we ignored the minor intercepts and assumed that the slopes represent the competition coefficients. As slopes of the intra-specific models were much greater than their inter-specific counterparts (+39% and +51%, respectively), we did obtain a stable equilibrium.
Subsequently, we constructed a phase-plane diagram calculated from the slopes. Figure 5a shows the results of this classic but explorative LV-analysis.
The stable equilibrium of the LV-model is in qualitative agreement with the IBM, but the location of the equilibrium is distinct among the models. The LV-model predicted the equilibrium to be located atV 1 = 398.94 m 3 ha −1 andV 2 = 337.98 m 3 ha −1 , whereas the equilibrium observed in the IBM was at V 1 = 309.36 ± 17.17 m 3 ha −1 and at V 2 = 314.79 ± 17.17 m 3 ha −1 (mean ± standard deviation). Consequently, with a value of 736.92 m 3 ha −1 , the classic LV-model predicted the total stem volume at equilibrium to exceed or over-yield the carrying capacity of the single species at 621.72 by 115.20 m 3 ha −1 , while, using the original data of the IBM, the equilibrium was at 624.15 m 3 ha −1 and, thus, very close to the capacity of the single species. This is confirmed by the major axis regression that fitted the original data very well (R 2 = 0.959). With an intercept of 622.13 ± 2.88 m 3 ha −1 and a slope of −0.9942 ± 0.0092 (mean ± 95% confidence interval), the graph connected the carrying capacities perfectly.
To overcome the defect of a falsely overyielding LV-model in further explorative analysis, we applied a non-overyielding Lotka-Volterra model that is defined by the following equation:

Lotka-Volterra Parameters from Community Experiments with Localized Dispersal
According to the stem volume trajectories of the MesoFON simulation runs, the equilibrium was obtained much later with the localized dispersal, and, therefore, the community was likely less resilient (see Figure 2). The derivation of the LV-parameters from the community experiments with local dispersal should clarify the underlying mechanisms.
Once again, we fitted robust linear regression models to the plots "population-wide competition vs. respective stem volume" to obtain volume-specific competition coefficients. Table 1c illustrates the results of the linear regression; Figure 6a-d show corresponding diagrams with the fitted straight lines.
The values of R 2 ranged from 0.874 to 0.938, indicating a good fit of the normal linear regression. However, in the case of localized dispersal, the slopes of both the intraand inter-specific cases were much lower than the corresponding slopes for the random dispersal cases, but intra-specific slopes were still larger than the inter-specific slopes. The intercepts of the robust regression in the inter-specific cases were again close to zero. Similarly, differences between the robust regression with a nonzero intercept and the robust regression passing through the origin were hardly visible in the inter-specific diagrams of Figure 6b,d. In contrast, robust regressions in the intra-specific cases had large intercepts. Such an intercept represents a competition that is inversely proportional to stem volume and converges to zero at infinity. We factored in the intercept by averaging the corresponding volume-specific competition above 100 m 3 ha −1 , adding this value to the constant intraspecific competition value of the slope. Even though this is a conservative procedure, it raised intra-specific over inter-specific competition so dramatically that the resulting stable equilibrium in the exploratively applied classic LV-model (diagram not shown) greatly overyielded the carrying capacity of the species at 617.42 m 3 ha −1 , and the equilibrium observed in the IBM specified by V 1 = 318.50 ± 17.35 m 3 ha −1 and V 2 = 294.43 ± 17.96 m 3 ha −1 (mean ± standard deviation), or at 612.93 m 3 ha −1 , in total. The location of the observed equilibrium is confirmed by the major axis regression that fitted the original data very well (R 2 = 0.983). With an intercept of 610.56 ± 0.88 m 3 ha −1 and a slope of −1.0044 ± 0.0028 (mean ± 95% confidence interval), the graph connected the carrying capacities almost perfectly. In contrast, MARS did not provide meaningful results.
If we exploratively applied the non-overyielding Lotka-Volterra model instead, we were again confronted with consequences of the large difference between the intra-and the inter-specific competition; in this case, the equilibrium exhibited a strong under-yield of the carrying capacities, with a large distance between the axis intersections of the isoclines (see Figure 7a). This did not correspond to observed community dynamics.
To address the issue in further exploration, we made a major assumption: We separated intra-specific competition into an effective mean-field part (i.e., the competition coefficient determined by the slope alone) and ignored the part that we believe reflects the inner-cluster competitive processes and does not influence the competitive outcome (i.e., the corresponding intercept). As Figure 7b shows, this approximation resulted in the convergence to the observed equilibrium at a rate that was close to the observed rate.

Discussion
The two identical Rhizophora species that were designed to avoid each other's competition were able to stably coexist in the individual-based model regardless of the dispersal regime. In our study, we avoided artificial methods to stabilize the community, such as external seed input or large plot dimensions [14]. Hence, the stable coexistence observed here is genuine.
The carrying capacities, recorded as long-term averages and representative of steadystate population dynamics, in the global and localized dispersal treatment were 17.4% and 29.7% lower than those attained in the early colonization phase. During initial colonization, the even-aged monoculture stands drove along the self-thinning line [74], therefore attaining the highest possible stem volumes, and they passed through an overshoot. An initial overshoot of a forest stand's biomass followed by a damped oscillation toward a lower steady state has been reported as a potential outcome of forest gap models, such as JABOWA, which have the Shugart growth function in common with MesoFON; the phenomenon is known to be a consequence of lagged mortality and/or lagged regeneration (Borrmann and Likens 1979a,b, Peet 1981 cited in [47]). In accordance with that, a tree dies in MesoFON only when slow diameter growth persists for five consecutive years (compare Section 2.1.5) and generative reproduction does not start until a certain tree age is reached (i.e., when the crown surface area A is sufficiently large that the no. of propagules = f red × D × A ≥ 1, compare Section 2.1.1). Surpassing standing stocks after (re-)colonization are acknowledged behaviors of mangrove forests (Twilley, personal communication). As several assumptions of the LV model are unlikely to be met during the overshoot and oscillatory phases (temporal homogeneity, mean-field assumption: outliers in Figures 4  and 6 likely originated from these phases), the validity of the Lotka-Volterra model is constrained to the steady-state phase that commences of the order of 500 years after, what we consider, a complete canopy removal by a large-scale disturbance event at the beginning of the simulation (compare Figure 2a). Given the typical frequencies of large-scale mangrove disturbance by hurricanes, this leads to the question of whether the Lotka-Volterra equations are applicable to mangrove forests at all. Rhizophora mangle, even though considered as a sensitive species, is severely damaged only when sustained hurricane wind speeds exceed 177 km/h (Category 3 to 5 hurricanes, [75]). For Florida, Costanza et al. (2008) [76] report return periods of 6, 26, and 77 years for Category 3, 4, and 5 hurricanes to hit coastal wetlands, respectively. Based on these values, it seems appropriate to answer the above question with "No." However, one should bear in mind that not all the mangrove area being hit by a hurricane would be damaged severely. Even the large Category 5 hurricane Andrew damaged >90% of the Everglades mangrove vegetation on a coastal transect only in the 40 km eyewall zone and >70% in the eyewall zone, including a 40 km wide band north of it [77].
Nevertheless, as we characterized R. mangle as a typical pioneer tree species, we may also look at forest disturbance more broadly. In their global review on this topic, Folking et al. (2009) [78] list the recurrence of major large-scale forest disturbance types: fire (spatial scale: <1 to >10 4 km 2 ): recurrence <10-1000 years; hurricanes (10 3 -10 5 km 2 ): 15-200 years, flooding (10-10 4 km 2 ): 50-100 years. These figures suggest that the LV model might be applicable to less fire-prone forests and to forests rarely impacted by hurricanes if we take into account the above argument regarding the eyewall zone.
As the mean field assumption was approximately achieved in the individual-based model, we consider the effects of crown plasticity combined with random dispersal as an archetype for localized behavior. This is outlined in the following discussion. The meanfield assumption is generally considered valid when the environment is homogeneous, when organisms are well mixed, and when their interactions are over long distances [3]. In our simulations, the environment was homogeneous, with random dispersal, and the community was well mixed even though the trees were sessile. However, the application of the local interaction, present in all the MesoFON runs, should have resulted in a fine-scale spatial structure and a "plant's eye view" that differs from large-scale spatial averages [79]. However, as long as the alteration of spatial structure does not scale up (see below) and a sufficiently large number of trees take a broad sample of the density, the composition of the average neighborhood is a "miniaturized" image of the global density average (as in configuration-field models, [80]). As the mean field was approximately achieved, we ignored the minor intercepts and assumed that the slopes represent the competition coefficients. As slopes of the intra-specific models were much greater than their interspecific counterparts (+39% and +51%, respectively), we did obtain a stable equilibrium. The suggested mechanism is: when displaced, the community is returned to the stable equilibrium by the high mortality of one species due to the higher intra-specific competition, just as a rolling ball is returned to the bottom of a deep bowl by the moment of the gravitational force (compare Figure 8a). When displaced, the community returns to the stable equilibrium by the high mortality of one species (red and blue arrows) due to the higher intra-specific competition, just as a rolling ball returns to the bottom of a deep bowl by the moment of the gravitational force. (b) Local dispersal: When displaced, the community slowly returns to the stable equilibrium by the slightly higher mortality of one species (red and blue arrows) due to the higher apparent intra-specific competition (at the margin of intra-specific clusters), just as a rolling ball slowly returns to the bottom of a shallow bowl.
The non-overyielding LV model matched the individual-based equilibrium in the global dispersal treatment far better than the classic model and did under-yield only slightly. The latter result may be explained by the fact that the accuracy of the equation is restricted for ratios of intra-vs. inter-specific competition sufficiently close to 1. A mechanistic explanation for this "goal-oriented" LV-model is suggested by the comparison of competition coefficients attained in single-species and community experiments. In accordance with our expectations, inter-specific coefficients in the community (1.93 and 2.08 × 10 −4 ) are far below those of single-species experiments (2.62 × 10 −4 ), but concurrently, intra-specific coefficients (both 2.92 × 10 −4 ) are above their single-species counterparts. Both alterations have certainly contributed to the strong stabilizing effects of crown plasticity. While the former change represents a decrease of approximately 27%, the latter represents an 11% increase in the competition. We thus conclude that it is necessary to estimate intra-specific competition coefficients in both community and single-species experiments, because they are not invariant and their comparison contributes to the mechanistic understanding of the underlying processes. Yet, the mechanistic interpretation of the pattern is: the implemented inter-specific neighborhood avoidance moved heterospecific trees apart from each other and pushed conspecifics together. Based on this result, the linear trend of the original data and the occurrence of a stable equilibrium at carrying capacity seem reasonable.
However, even though the constructed non-overyielding LV model provides a good prediction of the equilibrium and the speed approaching the equilibrium (determined by the width of the triangles between isoclines), it implies the occurrence of bending at the equilibrium and a reduction in the apparent carrying capacities. Both implications seemed to be absent from the original data, because MARS did not provide reasonable results. Moreover, the latter implication conflicts with the carrying capacity measured in the single-species experiments. Yet, we cannot rule out that the properties of the LV-model alter when one of the species becomes rare (data not shown). According to the discussion in Loreau (2010) [64], a non-overyielding equilibrium occurs only when isoclines are nonlinear and convex functions. Apparently, this was not the case here. From the perspective of an LV-modeler, a linear trend, as observed in the original data of the individual-based simulations, is not feasible.
In summary, we found that a non-overyielding model, rather than the classical Lotka-Volterra model, approximated most of the competitive behavior observed in the MesoFON model, even though the LV-model was not fully compatible with the IBM. The inter-specific crown plasticity acted as a strong stabilizing mechanism. Given the mean field assumption was approximately satisfied, it represented largely a localized process.
The approximate validity of the mean-field assumption was in line with mathematical descriptions given by Neuhauser and Pacala [4], but more importantly, it is in accordance with the work of Adams et al. [12], who stated that competitive interactions in real temperate forests are less local than those reported to cause species separation in theoretical studies conducted with spatially explicit model types other than field-calibrated IBMs (e.g., continuous-time Markov process [4]; interacting particle systems, stochastic point processes, noncontiguous patch models [60]).
Localized seed dispersal is generally considered an equalizing coexistence mechanism [14][15][16][17] that, by definition, slows the community dynamics but does not alter the competitive outcome among the species. In this study, we have analyzed two community stability properties, namely constancy and resilience. Both were affected by localized dispersal in the same direction, but to a different extent. While local dispersal lowered the constancy only slightly by about 4%, it reduced the resilience about 4.5-fold. The decline in resilience is consistent with the definition of an equalizing mechanism, even though this does not necessarily mean that the community is less stable with local dispersal. We expect a similar drop in the antagonist of resilience, the resistance [56], because this would lead to delayed responses, as described in the above definition. Resistance, however, is almost impossible to measure in the field [56], Unfortunately, this prevented us from including resistance in our analysis. Based on the gained understanding, we recommend defining resistance henceforth as the rate of a press perturbation that does not yet lead to a significant change in the community at the reference state, and determining resistance in future IBM studies by assigning increasing numbers of species 1 trees to species 2 until the threshold that leads to a change is passed. Nonetheless, the lower resilience and the presumably delayed community responses contradict the expectation that stabilizing and equalizing mechanisms act in the same direction and that both contribute to the coexistence of species.
However, in this study, plotting "population-wide competition vs. stem volume" provided evidence of the following mechanism: Individuals in a large portion of the plot area formed conspecific clusters and left the competitive race. The strength of intra-and inter-specific competition, as well as their difference, diminished, concurrently, in the inter-cluster mean field. When displaced, the locally dispersed community slowly returned to the stable equilibrium by the higher mortality of one species due to the slightly higher apparent intra-specific competition, just as a rolling ball is slower when returned to the bottom of a shallow bowl compared to the bottom of a deep bowl by the moment of the gravitational force (compare Figure 8b). In contrast to crown plasticity, local dispersal is a process that scales up to larger spatial dimensions.
Implications of the results for field surveys: Our case study revealed that an excess of inter-specific over intra-specific crown displacement facilitates the coexistence of two typical pioneer tree species ( [81,82] and Material and Methods section). However, as we pointed out in the introduction, the artificiality of our critical case (identical species exhibiting solely inter-specific crown plasticity) requires us to evaluate the implications of the obtained results for field surveys. Hence, in the following, we first address the question: Have any field studies provided evidence of excessive inter-specific crown plasticity and positive implications on species coexistence?
Plants, in general, are known to respond pro-actively and differently to the shade of different neighboring species via their red/far-red (R/FR) spectral perception system [83]. In accordance with that, a wide range of alterations in crown architecture at various spatial scales was found between mixed and pure stands in terrestrial laser scanning studies [84][85][86]. In addition, conspecific neighbors were often found to be stronger competitors than heterospecific neighbors [87][88][89].
The most direct evidence of stronger inter-specific plasticity, however, is offered by the survey on the variability of crown projection area in inter-vs. intra-specific environments conducted by Pretzsch et al. (2017) [29] on 80-year-old permanent plots containing mixed stands of European beech (Fagus sylvatica L.) x Sessile oak (Quercus petraea (Matt.) Liebl.) and pure stands of those two species. With 2326 beech and 1959 oak trees per hectare (4285 ha −1 in total), mixed stands comprised a much larger tree density than pure stands of Fagus sylvatica (3173 ha −1 ) and of Quercus petraea (2888 ha −1 ), respectively. This occurred partly at the expense of tree size [29]. Crowns of individual oak trees responded mostly with a wider projection area per stand area (cpa/sa, the latter derived from Voronoi tessellation) to the crowding in the mixed stands (0.90 in pure to 1.39 in mixed stands) and, to a lesser extent, with a rise in normalized crown displacement (3.10 in pure to 3.30 in mixed stands). In contrast, beech trees responded to it primarily with an exceeding surge in normalized crown displacement (4.0 in pure to 7.4 in mixed stands) and less so with an increase in cpa/sa (1.32 in pure to 1.48 in mixed stands). It might be due to the overall high packing density that the late-successional and shade-tolerant species was more plastic than the mid-successional light-demanding species in this study, whereas in most other studies, the opposite ranking was observed [90][91][92]. However, because crown displacement had been divided by d 1.37 and, thus, normalized to equal tree size (29), it is highly likely that European beech and, even though less likely, that Sessile oak expressed stronger interspecific crown plasticity values. As our study suggests, a stronger inter-specific crown displacement opens a window for stable coexistence, but its realization will likely depend on the difference between carrying capacities of the species and the excess of inter-specific over intra-specific crown plasticity for both species-irrespective of whether early-or late-successional species are more plastic. The carrying capacity of European beech, in particular, seems to exceed that of Sessile oak by far, as intensive forest management or frequent disturbances are needed to prevent Sessile oak from extinction [93]. If, in turn, the carrying capacities of early-and late-successional tree species resemble one another more closely, excessive inter-specific crown plasticity could represent an alternative coexistence mechanism beyond the competition-colonization tradeoff in the disturbance-mediated Mosaic cycle [28,46,64,94]. To the best of our knowledge, such a mechanism has not been studied yet. Its substantiation will require the careful interaction of field research and modeling efforts.
The second question we would like to address is the specification of the first question to mangroves. This is not an easy task, because field research on mangrove crown plasticity is still in its infancy. Vovides et al. (2018) [95] were the first researchers and, so far, the only ones who investigated the crown displacement of mangrove species. They measured stem positions and crown projection areas in eight cardinal directions for canopy and sub-canopy trees of Neotropical mangrove species (R. mangle, Laguncularia racemosa (L.) C.F.Gaertn., A. germinans) in Mexican mangrove forests along a salinity gradient. Wind direction/intensity exerted a strong influence on the crown displacement of canopy trees at all salinities, but for sub-canopy trees, a change in driving forces of crown displacement along the salinity gradient was detected: In the dense canopies at low salinity, neighborhood avoidance prevailed as a driver, whereas in sparse canopies at high salinities, crown overlap and neighborhood avoidance subsided-presumably due to stronger belowground competition [95]. Among the species, the late-successional Avicennia germinans did exhibit the strongest crown plasticity-reminding us of the beech-oak case study described by Pretzsch et al. [29] (see above). Examining differences between inter-and intra-specific crown displacement and/or their implications for mangrove species coexistence was not among the objectives of Vovides et al. [95]. However, as these data are freely available on the worldwide web (compare the link in [95]), this could be readily addressed by comparing response strengths of crowns in heterospecific and monospecific neighborhoods or by allowing crown plasticity parameters in MesoFON to differ between species-specific neighborhoods and retrieve optimum values for them via Genetic Algorithms (compare Grueters et al., 2019 [23]).

Conclusions
In qualitative support of the LV-model, two identical Rhizophora species that avoid competition with one another by lateral crown displacement can stably coexist in the "small world" of a field-calibrated individual-based model with a homogeneous environment.
The mean-field assumption was approximately satisfied in the random propagule dispersal treatment. In spite of that, the classic Lotka-Volterra model that was calibrated from the IBM simulation runs was invalid. Rather, the "non-overyielding" LV-model provided a more accurate approximation of the competitive behavior observed in the IBM. The underlying mechanism possessed the following characteristics: The implemented inter-specific neighborhood avoidance moved heterospecific trees apart from each other, and it pushed conspecifics together. Inter-specific crown plasticity represents the archetype of a local process that exerts a strong stabilizing effect when it is combined with random dispersal.
Community dynamics were slowed by localized propagule dispersal. This delayed the time to reach equilibrium and, thus, the recovery time of the community. The underlying mechanism for this behavior was: Individuals in a large portion of the plot area formed conspecific clusters and left the competitive race. The competitive race took place in the inter-cluster mean field, where the strength of intra-and inter-specific competition, as well as their difference, diminished. Our results suggest that the community dynamics of field-calibrated individual-based models deviate from the behavior of a mean field if local competition and localized seed dispersal combine, but no deviation occurs when local competition operates in isolation.
In this study, we found that the dynamics of a field-calibrated IBM followed predictions of classical competition theory, rather than resource competition theory. As fieldcalibrated IBMs are resource-based and, thus, the competitive outcome in such models is implicitly mediated by the competition for resources, we conclude that classical competition mechanisms can override those of resource competition theory. This implies that more species are able to coexist in a plant community than predicted by resource competition theory.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. The Three Theories of Plant Competition
Appendix A.1. Classical Competition Theory: The Lotka-Volterra (LV) Model The first rigorous and systematic effort to model the coexistence of competing species was independently made by Alfred J. Lotka and Vito Volterra in the 1920s [6,7]. Their mean-field competition model is based on a set of differential equations for species-specific logistic growth and it describes observed community dynamics in terms of species-specific growth rates, carrying capacities, and intra-vs. inter-specific competition coefficients. The main message of the Lotka-Volterra model in terms of species coexistence is: Stable coexistence at an over-yielding equilibrium-one that exceeds the carrying capacities of the single species-is enabled when intra-specific competition is stronger than inter-specific competition. Hence, coexistence in the model is defined by biotic interactions alone and is based on niche separation. Initially, the theory was applied in animal ecology [6] where species, having separate niches due to distinct food sources or diurnal activity periods, compete therefore more strongly intra-specifically and manage to coexist at an over-yielding equilibrium. In contrast, plant species are similar in their resource usage, and, thus, mechanistic resource-based evidence for niche separation at local sites is difficult to provide [96]. A stable coexistence at an over-yielding equilibrium was found among certain crop species [64]. However, because the LV-model is descriptive or, say, phenomenological [64], this is merely circumstantial evidence. The Lotka-Volterra theory remained under criticism for its phenomenological nature [62,64], but it was not until Tilman published his book on "Resource Competition and Community Structure" [8] that a major change toward mechanistic models took place. In his Monod-type model that was originally proposed for aquatic organisms, two species deplete a common resource until a resource level R* is finally reached, at which only the best competitive species can survive [8]. While this model was initially postulated for a (below-ground) mineral resource, light was later also treated as a limiting resource by defining the light intensity L* out measured below the canopy (in a monoculture) as an analog to R* [51,52]. In any case, the bad message inherent to this model was: In the competition for one resource, coexistence no longer existed. There is, in fact, only one best competitor.
However, in the extension toward multiple resources, it is still possible that species coexist at an equilibrium.
The conditions under which this happens are best explained using Tilman's famous graphical model representation (compare Figure 24, Case 3, p. 73 in [8]) in which two species compete for two resources and plot axes represent the supply rates of the two resources [69]. The behavior of each species is defined (1) by its zero growth isoclines that-under the assumption of Liebig's law of the minimum-have a rectangular shape and run parallel to the axes at low supply rates, and (2) by its consumption vector whose direction depends on the required resource ratio. An equilibrium exists only if the zero growth isoclines of the species cross, meaning that each species is the best competitor for the other resource. This equilibrium is then located at the intersection point of the two species-specific isoclines. Furthermore, this equilibrium can only be achieved if each species consumes its most limiting resource at a faster rate than the other one and when the ratio of the resources in the initial supply is balanced and intermediary between the different demands of the species (the consumption vectors, [97]).
On the other hand, if either resource is supplied in excess, the other resource becomes limiting, and-as in the single resource case-the best competitor for that resource survives [97]. In this situation, spatial and temporal heterogeneity of the environment might explain plant diversity and coexistence in an area [14,17,98]. However, irrespective of whether the environment is homogeneous or heterogeneous at an equilibrium state, the species richness cannot exceed the number of limiting resources [8,99]. Certainly, this prediction is at odds with numerous observations given the enormous diversity of plant species present in many areas [100], despite there being merely a handful of resources for which plants compete. This contradiction has been called the conundrum of biodiversity research [101].
As a way out of the conundrum, avoidance mechanisms have been introduced. In the spatial domain, the competition-colonization trade-off [14,102,103], and in the temporal domain, the storage effect [14,104] are thought to contribute to coexistence. From modeling studies on algal species [51], we can infer even more opportunities for coexistence, as these models yield a periodic state where more species than limiting complementary resources can coexist in a homogeneous environment. However, tests of the resource ratio hypothesis carried out at the Cedar Creek experimental facility revealed conflicting results [28]: While the responses of three out of four species to fertilization was consistent with the resource ratio hypothesis, successional changes in species composition along a soil nitrogen gradient were opposite to predictions. Beyond that, the R* theory has been criticized because it (falsely) focuses on the average concentration of a nutrient in a soil solution rather than on nutrient acquisition by the (species-specific) root length density and nutrient supply preemption, as soil scientists do [48].

Appendix A.2.2. The Resource Co-Limitation Theory
In response to this criticism, efforts have been made over the past two decades to develop an alternative resource limitation theory. Based on economic principles, natural selection should favor the co-limitation of multiple resources. Given the long evolutionary history of plants, resource co-limitation is expected to occur frequently. In accordance with this expectation, there is mounting evidence of N/P co-limitation in terrestrial, freshwater, and marine ecosystems from meta-analyses of large numbers of fertilization experiments [9,10]. As first outlined by Gleason and Tilman [105], co-limitation requires plants to flexibly adjust carbon and/or protein allocation to balance trade-offs between the uptake of different resources, store resources in excess to current demand for later usage, or, overall, to be plastic in their response to the supply of multiple resources [11,48]. Basically, two approaches were taken to formulate resource co-limitation: (1) a sum rule that assumes that resource cycles advance independently in series, and (2) a product rule that assumes an intimate bound between resource cycles such that a delay in one cycle constrains biosynthesis processes in another cycle [11]. The model of Wirtz and Kerimoglu [11] integrates these two approaches, as well as the Liebig law, in a seamless manner. The transition from the Liebig law to the sum and, finally, product rule goes hand in hand with a decline in curvature and a change in shape from a true rectangle to a rectangular hyperbola in the Tilman-like graphical representation (Figure 1 of [11]). The change in isocline shape, however, is modest. Irrespective of the internal (biosynthetic) co-limitation scheme, optimal eco-physiological regulation might restrict the apparent co-limitation to the sum rule that deviates even less from Liebig's law ( Figure 6 of [11]). Additionally, co-limitation modelslike resource ratio models-did yield a periodic state where more species than limiting complementary resources can coexist in a homogeneous environment [106]. Due to these arguments, co-limitation theory is unlikely to cause substantial changes to the conundrum of biodiversity research. Less clear are the effects that the high flexibility in stoichiometric ratios exerted by most plants has on biodiversity issues. Certainly, the thinking in (static) demand vectors should be given up and replaced by the thinking in demand functions being close together for various species at high resource supply rates, but drifting apart at low supply rates. It is likely that the resulting behavior will be in accordance with the results of fertilization experiments that the loss of plant species at high resource supplies is caused by a reduced niche dimension [107]. This is the ongoing debate from which IBM stayed away, as explained above.