Individual-based simulation of the spatial and temporal dynamics of macroinvertebrate functional groups provides insights into benthic community assembly mechanisms

Basic principles

Two versions of the model represent differences in the settlement success of benthic organisms between the subtidal and the intertidal zone. The model represents inter- and intra-group competitive interactions within macroinvertebrate communities, by recreating the life cycle of individuals, from their recruitment as juveniles, to their growth and death. These basic features are further modified through a basic representation of predation and biogenic habitat modification in the form of epibiosis (i.e., living on the surface of a living organism) and sediment engineering (i.e., modifying, maintaining and creating habitats by modulating sediment characteristics).

The supply of larvae and the settlement of juveniles are key factors shaping benthic communities, but their complexity often renders simplifications particularly difficult (Pineda, Reyns & Starczak, 2009). Recruitment processes should, at least partly, depend on the size of the spawner pool (i.e., the sum of adult individuals within a specific area) and the salinity preferences of each organism. They are also expected to be influenced by biological traits, such as fecundity, dispersal distance and early development mode, whose values often form combinations indicative of important life history trade-offs (Kupriyanova et al., 2001). The recruitment of juvenile organisms is also significantly impacted by post-settlement mortality, which shows high levels of both intra- and inter-specific variation (Hunt & Scheibling, 1997). Information at the former level indicates that small initial increments in juvenile body size lead to significant increases in survival rates in the face of predation (Gosselin & Qian, 1997).

The successional dynamics of benthic systems have been found to feature space limitation due to adult–juvenile interactions and exploitative competition for food, with the functional role of organisms largely defined by their size (Van Colen et al., 2008). Smaller organisms are expected to show higher growth rates and be better competitors for limited amounts of space, while larger organisms can occupy larger areas and should be competitively superior in the face of food limitation (Alexandridis et al., 2017b). Competitively equivalent individuals also compete for food and space, which should affect their growth rates (Côté, Himmelman & Claereboudt, 1994). Trophic distinctions, like the one between suspension- and deposit-feeders, can be invalid in view of the highly facultative feeding behavior of benthic organisms (Snelgrove & Butman, 1994). Trophic interactions could, instead, be dictated by expert knowledge and theoretically anticipated allometries, which should increase the realism and potential for stability of model representations (Brose, Williams & Martinez, 2006b).

The majority of predators appear to be larger than their prey and predator size tends to increase with prey size (Cohen et al., 1993). On the other hand, predator–prey body-size ratios are generally the lowest, just over two on average, for marine invertebrates, compared to other taxonomic groups and habitat types. It is possible that the energetic costs of prey capture and consumption set a limit to predator–prey size differences (Brose et al., 2006a). In the Rance estuary, fish predation of benthic organisms is restricted and predation by birds is highly seasonal and mostly limited to the intertidal zone. The mortality caused by predatory macroinvertebrates appears likewise to be limited in magnitude, as far as adult prey is concerned. This is due partly to the greater impact of predation on juveniles and partly to the partial ingestion of adults and the regenerative properties of many of them. Predation pressure should depend on the organisms’ defensive mechanisms and position in the sediment, along with their relative abundances (Desroy, 1998).

Benthic macroinvertebrates in the Rance estuary appear to differ with regard to their use of space, between those that are buried in the sediment and those that occupy its surface supported by hard substrates (Alexandridis et al., 2017b). The dominance of soft bottoms in the estuary indicates that the latter organisms partly occur due to the benefits provided by the phenomenon of epibiosis (Wahl, 1989). The settlement of basibiotic (i.e., providing substrate for other organisms) individuals should be facilitated by the habitat modification effect of other basibionts, while the settlement success of epibionts (i.e., organisms living on the surface of other organisms) should depend on their specific substrate preferences. Still, sediment engineering is expected to play a much more prominent role in soft bottom systems (Meadows, Meadows & Murray, 2012). Its impact on benthic communities can be summarized by the mobility-mode hypothesis, which groups organisms into sediment stabilizers and destabilizers (Posey, 1987). The dominance of the largest members of a group in an area should allow them to modify sediment characteristics in a way that facilitates the settlement of the group’s own members, while inhibiting members of the other group.

Entities, state variables and scales

The only entities are cells that make up a square grid of 60 × 60 cells. In order to avoid edge effects, the grid wraps horizontally and vertically into a torus. Cell size corresponds to the theoretical area exclusively occupied by an individual that belongs to one of the small groups. The dimensions of the grid represent an arbitrary sampling area of the real system. One time step corresponds to one year, starting right before spring larval dispersal.

A cell variable indicates whether each cell of the model is occupied and by a member of which functional group. There are twenty groups in the model (Table 2). Ten of them belong to the infauna (i.e., animals living inside the sediment) (FG1–8, FG10–11), five belong to the epifauna (i.e., animals living on the seabed’s surface) (FG12–16) and one can be part of either the infauna or the epifauna (FG9); all of these groups represent algae/detritus feeders. Four more groups represent predator/scavenger organisms (FGP1–4).

The infauna consists of four small groups (FG1–4), which can occupy one cell, five intermediate groups (FG5–9), which start with one cell and can occupy its eight immediate neighbors during their growth, and two large groups (FG10–11), which can, through the same procedure, occupy one cell and its twenty-four closest neighbors. The group that belongs to the infauna or the epifauna (FG9) represents basibiotic organisms of intermediate size. Epifauna, including individuals of the latter group but staying small in size, can settle in cells occupied by this group. In total, epifauna consists of four small groups (FG9, FG12–14), which can occupy one cell, and two large groups (FG15–16), which start with one cell and can occupy its eight immediate neighbors during their growth.

Process overview and scheduling

The first of the model’s actions represents the process of recruitment. Juveniles of the eleven infaunal groups settle randomly in empty cells, with probabilities defined by the contribution of each group to the infaunal spawner pool, along with the relative abundance of sediment stabilizing and destabilizing groups. Infaunal juveniles experience random post-settlement mortality and those that die are removed from the system. Juveniles of the six epifaunal groups settle randomly in cells that are occupied by infaunal adults of the basibiotic group and are empty of epibionts, with probabilities defined by the contribution of each group to the epifaunal spawner pool. Epifaunal juveniles experience random post-settlement mortality and those that die are removed from the system. During the first time step, all cells are considered empty, which corresponds with the initial conditions after the estuary was re-opened in 1966, and all groups of the infauna and the epifauna have equal contributions to the respective spawner pool.

The second action represents the growth in terms of occupied cells of juveniles that belong to groups with intermediate and large size and the process of inter- and intra-group competition for food and space that this entails. First, the juveniles of intermediate infaunal groups, then only those infaunal juveniles of the basibiotic group that settled near other group members and finally the juveniles of large infaunal groups grow in random order within each of these categories. The juveniles of the two large epifaunal groups are next to grow, also in random order within each category, first those of the group that is associated with hard substrate and then those of the group that is associated with gravel. During these growth processes, juveniles that grow earlier are allowed to grow over non-growing individuals and juveniles that grow later. Those of the latter that have not been grown over are allowed to grow over all other individuals. Juveniles are in all cases not allowed to grow over individuals within the same category.

The third action of the model represents the process of ageing by one year of all individuals that survived the previous time step.

During the fourth and final action, all individuals that could not grow to their full predefined size, were grown over, reached their lifespan during the current time step or were epibionts of deceased basibionts that could not be replaced, die and are removed from the system. Basibiotic individuals that die of ageing and have epibionts of the same group, take the age of their oldest epibiotic basibiont and retain the rest of their epibionts. Individuals of potential prey groups are then randomly selected to be removed due to predation. This is done in decreasing order of the predators’ size and starts with their most abundant prey group. If this is less abundant than each predator’s total prey, individuals from its next most abundant prey group are additionally removed, until the number of individuals removed is the closest possible to each predator’s total prey.

Subtidal settlement probabilities

The settlement probabilities of all infaunal and epifaunal groups that are employed during the model’s initialization are defined by each group’s early development mode (T2), dispersal distance (T3) and fecundity (T4). Specifically, groups with brooded early development mode (i.e., development of eggs or juveniles under parental care) (FG4, FG5, FG12, FG14) or long dispersal distance and high fecundity (FG7, FG8, FG9, FG10, FG13, FG15), along with groups with any of the last two trait levels and laid early development mode (i.e., development of eggs actively deposited on the substrate) (FG11), have settlement probabilities that are three times the settlement probabilities of groups with none of the aforementioned trait levels (FG2). The latter groups have settlement probabilities that are half the settlement probabilities of groups with planktonic early development mode (i.e., development in the water column without the ability to swim against a current) and either long dispersal distance (FG1) or high fecundity (FG6, FG16), along with groups with short dispersal distance, low fecundity and laid early development mode (FG3).

The settlement probabilities of infaunal and epifaunal groups that are employed during the model’s action of recruitment use the settlement probabilities employed during the model’s initialization, multiplied by a factor proportional to each group’s contribution to the respective spawner pool. The settlement probabilities of infaunal groups are additionally determined by each group’s position in the sediment (T9) and role in sediment engineering (T14), along with the relative abundance of sediment stabilizing and destabilizing groups. This represents the impact of sediment engineering and is applied as follows. First, the settlement probabilities of mobile stabilizers and groups that live deep in the sediment (FG2, FG3, FG6, FG8) are divided by 2. The rest of the rules apply only to the remaining infaunal groups. If the effective sediment stabilizers (intermediate and large sessile stabilizers FG5 and FG10) outnumber the effective destabilizers (intermediate and large destabilizers FG7 and FG11), the settlement probabilities of the stabilizers (FG5, FG10) and the basibiotic group FG9 are multiplied by 2 and those of the destabilizers (FG4, FG7, FG11) are divided by 2. Otherwise, the settlement probabilities of the stabilizers (FG5, FG10) and the basibiotic group FG9 are divided by 2 and those of the destabilizers (FG4, FG7, FG11) are multiplied by 2. The settlement probability of the infaunal group that is neutral with respect to epibiosis and sediment engineering (FG1) is not modified during this process.

Intertidal settlement probabilities

The settlement probabilities of all infaunal and epifaunal groups that are employed during the model’s initialization are partly defined by each group’s early development mode (T2), dispersal distance (T3), and fecundity (T4). The same multiplication factors are therefore used as those employed during the initialization of the subtidal version of the model. These factors are additionally modified based on each group’s tolerance for low salinity and tidal exposure (T5). Specifically, the settlement probabilities of all stenohaline (i.e., not tolerating long tidal exposure and salinities that differ greatly from those of the open sea) groups (FG1, FG2, FG4, FG5, FG6, FG9, FG12, FG14, FG15, FG16) are divided by 2.

The settlement probabilities of infaunal and epifaunal groups that are employed during the model’s action of recruitment use the settlement probabilities employed during the model’s initialization, multiplied by a factor proportional to each group’s contribution to the respective spawner pool. The settlement probabilities of infaunal groups are additionally determined by each group’s role in sediment engineering (T14) and the relative abundance of sediment stabilizing and destabilizing groups. This represents the impact of sediment engineering and is applied similarly to the subtidal version of the model; only this time the settlement probabilities of euryhaline (i.e., tolerating salinities that differ greatly from those of the open sea) and emersed (i.e., tolerating long tidal exposure) groups (FG3, FG7, FG8, FG10, FG11) are modified. If the effective sediment stabilizers (intermediate and large, euryhaline and emersed stabilizers FG8 and FG10) outnumber the effective destabilizers (intermediate and large, euryhaline and emersed destabilizers FG7 and FG11), the settlement probabilities of the stabilizers (FG8, FG10) are multiplied by 2 and those of the destabilizers (FG3, FG7, FG11) are divided by 2. Otherwise, the settlement probabilities of the stabilizers (FG8, FG10) are divided by 2 and those of the destabilizers (FG3, FG7, FG11) are multiplied by 2.

Post-settlement mortality

The mortality rates experienced by juveniles of all groups after their settlement on the seabed are defined by each group’s body size (T7). The selected values represent extreme and intermediate levels observed in nature. A fraction of the juveniles’ abundance equivalent to 90% for the small groups, 50% for the intermediate groups and 10% for the large groups are removed from the system following their settlement.

Predation

The potential prey of the predator groups consists of groups that are smaller or similar in size but no smaller than 1/3 of their own size (T7), are not buried deep in the sediment (T9) and their representative species are not protected by plates, shells or tubes. Specifically, FGP1 has groups FG4 and FG7 as its potential prey, while the potential prey of FGP2 consists of FG7 and FG11. The sessile predator group FGP4 is additionally limited to mobile epifaunal organisms (T10) and can, therefore, feed on groups FGP1 and FGP3. The availability of precise information on the diet of the representative species of predator group FGP3 in the Rance estuary allows the assignment of groups FG4 and FG5 as its potential prey. The abundance of predator groups at each time step is equal to 1/10 of the total abundance of their potential prey in the subtidal and 1/10 of the abundance of the most abundant group among their potential prey in the intertidal version of the model, both observed at the previous time step. The 1/10 value was derived by assuming 10% energy transfer efficiency between the two trophic levels, as well as similar values of biomass-relative productivity and individual biomass. The less and more conservative levels of total and most abundant potential prey were used in the subtidal and the intertidal version of the model, respectively. The total number of individuals of the potential prey groups that explicitly die due to predation and are removed from the system at each time step is equal to the abundance of their respective predators.

Basic principles

The spatial scales at which pre- and post-settlement processes take place are not expected to overlap considerably, as pre-settlement processes usually operate at much coarser scales than post-settlement ones (Fraschetti et al., 2002). Exchanges across spatial scales are mainly the result of larval dispersal; immigration and emigration of adults can in most cases be considered as trivial for the population dynamics of benthic organisms (Eckman, 1996).

Initialization

The grid of the coarse-scale model represents the bottom of the Rance estuary, with cell size corresponding to an area of approximately 0.2 km². A sediment type is attributed to each cell, based on a sedimentary map of the Rance estuary from 1994 (Bonnot-Courtois, 1997). Each cell is then assigned to the subtidal or the intertidal zone, based on its sediment type. Areas covered by gravel, coarse sand, intermediate/coarse sand, fine/intermediate sand, muddy sand and sandy mud are assigned to the subtidal zone. Areas covered by silty mud, mud, pure mud and salt marshes are assigned to the intertidal zone.

The subtidal version of the fine-scale model is the first to be loaded and the cells of the subtidal zone, in random order, ask it to initialize and export the generated model instances in files named after their own x and y coordinates. At the same time, each cell of the subtidal zone is attributed with the generated group abundances, which are printed out, as well as a color on the grid, which indicates whether these abundances are dominated by sediment stabilizers or destabilizers. The same procedure is then repeated for the cells of the intertidal zone and the intertidal version of the fine-scale model.

Process overview and scheduling

One time step corresponds to one year, starting right before spring dispersal. The cells of the subtidal zone, in random order, ask the subtidal version of the fine-scale model to import the model instances that were generated for them during the previous time step and set each infaunal and epifaunal group’s contribution to the respective spawner pool equal to the median abundance of each group within the cells themselves and their eight immediate neighbors that are part of the system. Within the same procedure, the cells ask the fine-scale model to move one step forward and export the generated model instances in files named after their x and y coordinates. At the same time, each cell is attributed with the generated group abundances, which are printed out, as well as a color on the grid, which indicates whether these abundances are dominated by sediment stabilizers or destabilizers. The same procedure is then repeated for the cells of the intertidal zone and the intertidal version of the fine-scale model.

Model simulation

A 10-year simulation of the model was replicated three times. The number of time steps was limited to 10, because the possibility of all sites in the Rance estuary to stay undisturbed and evolve concurrently should decrease significantly as the number of years increases. The choice of three replicate simulations was made for practical reasons. The simulation that produced the highest level of β-diversity (see section Sensitivity analysis) was singled out for a more detailed analysis of the model output (hereafter called the “benchmark simulation”).

Sensitivity analysis

Four modifications were independently applied to different elements of the model, in order to examine the role of the respective processes. The role of dispersal distance was examined first. The contribution of each infaunal and epifaunal group to the respective spawner pool in each cell of the coarse-scale model was derived from group abundances in all of the model’s cells, instead of just each cell’s immediate neighbors. The role of post-settlement mortality was examined next, by completely removing it from the fine-scale model. The role of sediment engineering was examined by eliminating its effect on the settlement probabilities used in the fine-scale model. Finally, the role of predation was examined by removing explicit predation mortality from the fine-scale model. Each of these model configurations was simulated for 10 years and replicated three times.

Changes in the model’s behavior were investigated by depicting the evolution of modelled β-diversity, compared to the value that was observed at the level of functional groups in the Rance estuary in 1995 (data can be found in section S3 in the Supplemental Material). β-diversity was in all cases quantified as the variance of a Hellinger-transformed table of group abundances in different sites (Legendre & De Cáceres, 2013). In the case of the observations, the sites correspond to 113 stations (71 subtidal, 42 intertidal) that were sampled in the Rance estuary in 1995. In the case of the output of the five different model configurations, the sites correspond to 113 of the 230 cells of the coarse-scale model (71 subtidal, 42 intertidal; see Fig. 2 for the exact position of the cells in the model). These cells were used with the goal of approximating the location in the estuary of the sampled stations.

Correspondence analysis

Correspondence analysis (Legendre & Legendre, 1998) was first performed on the table of functional group abundances in 113 stations that were sampled in the Rance estuary in 1995. It was then performed on tables of group abundances in 113 cells selected out of the 230 cells of the coarse-scale model (see section Sensitivity analysis) in the 1st, 2nd, 3rd and 10th year of the benchmark simulation. These years were chosen to represent the qualitatively greatest changes in model output. All tables of group abundances were Hellinger-transformed before analysis. The goal was to compare patterns in the relative frequencies of functional groups across stations or cells, so scaling 2 was selected for the projection of groups on the first two axes of the reduced multivariate space (Borcard, Gillet & Legendre, 2011).

Spatial correlation

Mantel correlograms were used to quantify spatial correlation in the multivariate domain of functional group abundances. The technique is based on calculation of the normalized Mantel statistic between pairs of site dissimilarity matrices. One matrix in each pair quantifies differences in multivariate community composition and the other is derived by attributing the value 0 to pairs of sites that belong to the same distance class and the value 1 to all other pairs of sites. The process is repeated for each distance class and values of the Mantel statistic are tested by permutations. Mantel correlograms were produced for the 113 stations that were sampled in the Rance estuary in 1995 and the output of the 1st, 2nd, 3rd and 10th year of the benchmark simulation in 113 cells selected out of the 230 cells of the coarse-scale model (see section Sensitivity analysis). Distances in the model were measured between the cells’ centers of symmetry, by assuming cell dimensions of 450 m × 450 m. The tables of group abundances were Hellinger-transformed and Holm’s correction for multiple testing was applied to the permutation tests. The number of distance classes was in each case based on Sturge’s rule and the correlograms were restricted to distances that included all sites (Borcard, Gillet & Legendre, 2011).

Software

Simulations of both the fine- and the coarse-scale model were implemented in the multi-agent modeling environment NetLogo version 5.3.1 (Wilensky, 1999). Interactions between the two scales were realized through the NetLogo extension LevelSpace (Hjorth, Head & Wilensky, 2015). The source code of the subtidal and the intertidal version of the fine-scale model and the coarse-scale model can be found in section S4 in the Supplemental Material. The NetLogo model files and the GIS data of substrate types in the Rance estuary are also available in the Supplemental Material. All model analyses were performed using the statistical software R version 3.2.2 (R Core Team, 2015) with the packages vegan (Oksanen et al., 2015) and raster (Hijmans, 2016) and the function beta.div (Legendre & De Cáceres, 2013, Appendix S4).

Model simulation

A common pattern in all simulations of the standard model configuration was the initial presence of all functional groups in most cells of the coarse-scale model and the gradual dominance of different groups in different sets of cells (see section S5 in the Supplemental Material). During this transition, three obligate epibiotic groups (FG12–14) initially became rare and were eventually eliminated from the system. The median number of groups per cell of the coarse-scale model declined from 20 at initialization to 11 after 10 years of simulation in all three replicates. The maximum and minimum number of groups per cell in the 10th year was also the same in all replicate simulations, equal to 16 and four, respectively. The general trend was for the whole system to be covered by stabilizer-dominated cells. Pockets of resistance to this trend were, however, formed by destabilizer-dominated cells of the subtidal and the intertidal zones, within which the opposite trend could be observed.

Sensitivity analysis

The evolution of β-diversity during the 10-year simulations of the standard model configuration shows a clear increasing trend during the first nine years, after which it levels off (Fig. 3A). The levels reached in the 9th year were, at least in one of the three simulations (0.35), close to those observed in the Rance estuary (0.38). The replacement of local dispersal by its global counterpart drastically changed β-diversity in the modeled system (Fig. 3B), as no clusters of cells that were dominated by sediment stabilizers or destabilizers were formed. The very low levels of β-diversity that were reached in this case were the result of the distinction between the two tidal zones, along with stochasticity regarding the relative abundances of habitat modification groups in the fine-scale model’s initialization. Without sediment engineering, this stochastic effect was removed and all simulations had almost identical output (Fig. 3D). β-diversity increased at a very slow rate, as differences in community composition were gradually amplified through local dispersal. The removal of post-settlement mortality allowed small groups to overwhelmingly dominate the fine-scale model, resulting in unrealistic abundance levels and low β-diversity (Fig. 3C). The removal of predation mortality had a minimal impact on β-diversity (Fig. 3E). Still, the values that were observed in the 10th year were slightly higher than those of the standard model configuration.

Correspondence analysis

The relative position of groups as they are projected on the first two axes of the reduced multivariate space illustrates their similarity with regard to their relative frequencies across sites or cells (Fig. 4). The first two axes of the correspondence analysis that was performed on observed group abundances (Fig. 4A) represent only about 32% of the total variation, whereas the same value for the output of the 1st, 2nd, 3rd and 10th year of the benchmark simulation (Fig. 4B–4E) equals 72%, 84%, 79% and 73%, respectively, all above 70%.

Some of the patterns that can be seen in the observed group associations are also evident in the model output. The observed association of groups FG1 and FG12, both of which are neutral with regard to epibiosis and sediment engineering, can be seen after the first simulation year. The same holds true for the association of the intertidal groups FG3, FG8 and FG13. Groups FG12 and FG13 gradually get eliminated and groups FG1 and FG8 converge toward the majority of the groups, but the separation of group FG3 remains a constant feature of the model output. The basibiotic group FG9 is positioned near the gravel-associated groups FG6 and FG15 in both the observations and the initial simulation years, before epibiosis is largely eliminated from the modeled system. The association of subtidal prey group FG4 and, to a lesser extent, intertidal prey groups FG7 and FG11 with predatory groups FGP2, FGP3 and FGP4 is also a pattern shared by the observations and the entire simulation period after the initialization, as is the association of mud-associated, stabilizer groups FG2, FG5 and FG10.

Spatial correlation

The Mantel correlogram of the observed group abundances demonstrates significantly similar community composition within distances of 0.5 and 1.5 km and significant compositional dissimilarity at distances of 3.5 and 7.5 km (Fig. 5A). Spatial correlation in the 1st year of the benchmark simulation was demonstrated as significantly dissimilar community composition at distances around 5 and 6 km (Fig. 5B). Similar levels of compositional dissimilarity were observed at slightly shorter distances in the subsequent simulation years (Figs. 5C and 5D), accompanied by significantly similar community composition within a distance of 1 km. The latter feature was retained in the 10th year (Fig. 5E), along with significant compositional dissimilarity at distances of 3 km and around 7 and 8 km, in a pattern that is similar to the one observed in the empirical data.

Mechanistic approach

This study aimed to explore the potential to reproduce the dynamics of benthic macroinvertebrate communities through a mechanistic representation of the system, based primarily on theoretical expectations about community assembly mechanisms. The first step was taken with the reduction of the system’s components through a systematic and testable procedure that retained sufficient information on the organisms’ functional role (Alexandridis et al., 2017a). The verification of the basic assumptions of the emergent group hypothesis (Hérault, 2007) affirmed the removal of functionally equivalent variability, allowing the subsequently built model to only reproduce the stochastic variation in a community’s composition that is relevant to its functioning.

The second step involved the study of associations of biological traits with environmental variables and with each other, seeking support for ecological theories that would allow the definition of functional relationships in the marine benthos (Alexandridis et al., 2017b). Each of the community assembly mechanisms that are represented by these relationships encompasses a variety of processes that could potentially be modeled in more detail. The level of representation was dictated by the available trait, environmental and theoretical knowledge. Hence, biological traits were used as proxies for the role of functional groups in a set of theoretically anticipated community assembly mechanisms.

Ideally, the rules of interaction in IBMs are formulated in terms of fitness maximization, providing a representation that is more general and reliable than empirical formulations (Stillman et al., 2015). In cases where entire communities and their complex webs of interactions need to be represented, IBMs have had to settle for more implicit representations of this first principle. Our rules of interaction are partly phenomenological, but they represent well-established ecological theories that use fitness maximization as their basis. Furthermore, the algorithmic representation of interactions allowed the incorporation of expert knowledge, which is widely available but often difficult to formulate mathematically.

Inter-scale modeling

The transfer of knowledge from the level of individual organisms, where most observations and experiments are performed, to the level at which biodiversity patterns are typically observed is one of the central problems in ecology (Denny & Benedetti-Cecchi, 2012). The nonparametric up-scaling approach of Cipriotti et al. (2015) is an example of how this issue can be addressed in IBMs. The important state variables of the fine-scale model define a state space, which is divided into a finite number of discrete states. Simulation runs of this model, covering a range of pertinent initial conditions, states and drivers, define transition matrices that are used by the coarse-scale model. The resulting up-scaling is not dynamic and it is restricted to the range of the simulations. Mumby & Dytham (2006) opted, instead, for the development of two dynamically linked, spatially structured models. Computational limitations imposed a mismatch between the two scales, which was bridged by using a toroidal lattice with no boundaries and a mean representation of some processes.

The representation of the dynamic link between processes that operate at different spatial scales is important for the reproduction of the main consequence of cross-scale interactions, namely nonlinear dynamics with threshold values (Peters, Bestelmeyer & Turner, 2007). This study identified transfer processes and spatial heterogeneity at intermediate scales as important components of the link between fine- and coarse-scale patterns and processes. Similarly to this conclusion, the model presented here consists of two spatially distinct representations that are linked by larval dispersal and cell variability. The cells correspond to large ecosystem areas and each of them is represented by a community of much smaller surface. The link between community and ecosystem is implemented through aggregate community measures fed into the ecosystem model, whose output in turn influences community assembly parameters.

Emerging patterns

The output of the fine-scale model can be considered to represent α-diversity. Diversity at this scale is maintained as a result of inter-group competitive trade-offs and intra-group inhibition. Pre-emptive competition for space is represented explicitly, while exploitative competition for food is implied in the overgrowth competition. The trait of body size is central to the definition of each group’s role in competition, additionally controlling predatory interactions and post-settlement mortality. The role of epibiosis and sediment engineering is represented through the modification of settlement probabilities, which are initially defined by a combination of reproduction-related traits. The generated abundance patterns were characterized by the gradual dominance of a few groups, through a process that represents competitive exclusion.

The turnover of group abundances in the coarse-scale model can be considered to represent β-diversity. Abundance patterns at this level were shaped by diversity among the cells and the process of local dispersal. Inter-cell diversity was driven by the distinction between the subtidal and the intertidal zone and the random initial dominance of stabilizers or destabilizers in each cell. The few areas of destabilizer-dominated cells that were left after 10 years of simulation did not seem able to resist the general trend of cells becoming dominated by sediment stabilizers. This trend could be the result of allowing the system to evolve undisturbed and/or due to overestimating the effect of sediment engineering, particularly sediment stabilization, which might be contingent on groups reaching certain density levels (Posey, 1987).

The model was able to generate levels of β-diversity that were not far from those observed in the Rance estuary in 1995. Exploratory longer simulations of the model indicate that it would probably fail at sustaining high diversity levels, mostly due to the eventual dominance of sediment stabilizers. Local dispersal and sediment engineering were indicated as major drivers of β-diversity. They played this role by facilitating the creation of clusters of cells that were dominated by sediment stabilizers or destabilizers. Early post-settlement mortality also contributed to the maintenance of high levels of β-diversity, but its impact on community composition appeared to be more profound. It prevented small functional groups from dominating and restricted total group abundances to realistic levels. It can probably be largely attributed to adult–juvenile interactions and predation (Desroy, Retière & Thiébaut, 1998).

The gradual elimination of the three epibiotic groups could indicate that epibiosis cannot by itself sustain observed levels of epifaunal diversity, calling for the inclusion of hard substrate types. However, the resulting rareness of groups is also a feature of the observations, which disappears from the model output only when these groups are eliminated. This discrepancy between observations and the final years of the simulation could be due to the concurrent development of cells in the model. Assemblages that are at different successional stages should sustain group rareness and be a more realistic representation of a system that is subject to frequent perturbations, such as wave action or extreme cold (Desroy, 1998). Marine benthic systems are known to demonstrate such transient dynamics (Beukema, Essink & Dekker, 2000).

Some associations of group distribution in the observations could also be seen in the model output, indicating the adequate representation of certain community assembly mechanisms, such as environmental filtering caused by tidal zonation and trophic interactions. It is not clear whether the similarity of observations with patterns produced by the model due to biogenic habitat modification reflects a good depiction of this phenomenon or just mimicry of environmental filtering caused by substrate type. The same holds true for patterns of spatial correlation. Compositional similarity at small distances and dissimilarity at intermediate distances emerged early in the simulation, due to tidal zonation. Compositional dissimilarity at large distances, on the other hand, emerged only toward the end of the simulation, probably due to the formation of clusters of synchronized cells. Similar patterns have been observed in other estuarine macroinvertebrate systems, where they were also attributed to synchronous community development, driven by local dispersal (Ysebaert & Herman, 2002). However, the role of substrate, especially in interaction with biogenic habitat modification, is not clear.

Future improvements

It appears that the potential of the model to explore the assembly mechanisms of benthic communities could be drastically increased through a few extensions. The combination of sediment engineering with a distinction between different substrate types could allow the model to sustain high diversity levels. This task would be facilitated by the addition of a realistic disturbance regime, which would keep the model away from equilibrium and promote the persistence of rare groups. The explicit representation of these processes would help disentangle the role of substrate type, habitat modification and physical disturbance in shaping marine benthos. The model might also benefit from an improved formulation of trophic interactions. A more detailed description of trophic strategies could be combined with known links between community structure and organic fluxes (Herman et al., 1999), toward the definition of individual-based interactions. Information on primary production would help define organic fluxes, while its combination with hydrodynamics could extend competition to larger distances (Fréchette, Butman & Geyer, 1989). The latter extension could additionally allow for a more explicit representation of larval dispersal, based on the groups’ reproductive trait values.

Any understanding of community assembly mechanisms that is gained by this and future versions of the model should eventually allow its transition to a predictive modeling tool, which could improve projections of benthic biodiversity responses to environmental change. The impact of change could be readily predicted by the model for certain of its aspects, while others might require modifications to the model’s formulation. For instance, the response of benthic communities to changes in the tidal zonation of the Rance estuary due to alternative operation schemes of the tidal power station could be modeled thanks to the inclusion of the environmental filtering effect of tidal zonation in the functional groups’ traits. The same holds true for the response of benthic communities to expected additions of invasive species with certain combinations of the available traits. The responses of benthic biodiversity to increased temperatures or ocean acidification would, on the other hand, require an extension of the model’s functional scope, in order to incorporate these aspects of environmental change. Ultimately, the use of the model for predictive purposes would require its precise parameterization and validation through context-oriented approaches (Kubicek et al., 2015).

Modeling community dynamics largely owes its value to the well-established but poorly understood link between biodiversity and ecosystem functioning (Srivastava & Vellend, 2005). The study of this link requires biodiversity patterns to be explicitly associated with ecosystem functions, such as energy and elemental cycling, the provision of habitat or the modification of physical properties of the system (Frid et al., 2008). These functions have long been the subject of statistical analyses based on many of the biological traits that were here assumed to represent the organisms’ role in community assembly mechanisms (Bremner, Paramor & Frid, 2006). The role of these traits can be extended to that of indicators of ecosystem functioning, in order to associate the latter with modeled biodiversity. This task can be assisted by theoretical (Brown et al., 2004) or empirical (Brey, 2010) links between traits and ecosystem functions.