A hybrid behavioural rule of adaptation and drift explains the emergent architecture of antagonistic networks

Abstract

1. Introduction

Antagonistic interactions, such as herbivory, parasitism and predation, are important to the provision of ecosystem function and service [1]. It represents the process of resource exploitation in ecological networks [2] and can divide species into clusters where consumers within a cluster are likely to share the same function and exploit similar resources [3–6]. This process of resource exploitation can also create a nested structure where species connected to specialists are embedded in the set of species connected to generalists [7]. Such clustering (namely, compartmentalization) and nested architectures can have profound effects on the stability of ecological communities [7–14]. Specifically, compartmentalization tends to stabilize ecological networks by containing the effect of perturbations within modules [3,11,15]. By contrast, although nested structure can foster high species richness [12] and enhance resilience against perturbations [9,10,16], it reduces species persistence in mutualistic networks [17] and destabilizes the community [14]. Despite their important role in securing ecosystem functions and services during perturbations, mechanisms that can account for the level of compartmentalization and nestedness close to those observed in real ecological networks remain poorly understood.

Both adaptive interaction switching [18] and evolutionary history [19,20] have been put forward as important factors in shaping the architecture of ecological networks. As the antagonistic interaction largely reflects the relationship between resources and consumers [2], it is the consumer's adaptive decision-making behaviour that decides which resources to exploit and, therefore, the structure of a bipartite ecological network. This adaptive nature of selecting and switching interaction partners (i.e. tuning target resources) is essential for the survival of consumers that are competing for available resources [21–23]. Adaptive interaction switches occur when the quantity and quality of available resources changes. Consumers prefer to select highly profitable resources rather than consuming all resources available to them as specified in optimal foraging theory [24]. Concurrently, they will exploit abundant resources over rare ones to avoid risk [25,26]. This adaptive behaviour of interaction switching could well be the predominant process that pushes an ecological network towards becoming compartmentalized or nested, forming a potentially stable and resilient complex community [18,27]. Furthermore, as evolutionary history can affect the traits of both consumers and resources, it often determines the availability and accessibility of resources on a physiological level [28,29]. Consequently, the architecture of an ecological network is also shaped by the imprint from its evolutionary history [20,30,31].

To date, although quite a few models have been put forward for explaining the emergence of function, compartmentalization and nestedness in antagonistic networks, those can be further improved and refined, especially regarding their predictive power, to incorporate all essential forces of evolution [32,33]. Evolution is largely driven by mutation, natural selection and random drift. Adaptive interaction switching implemented in most of these models only emphasizes adaptation from one or two particular processes, specifically by using the classical optimality methodology in evolutionary ecology [33] or optimization-based analytical treatment for adaptive behaviour [18,27]. To capture the essence of network evolution, the role of random drift needs to be appreciated, which is an important non-adaptive force to counterbalance the process of optimization. Real networks are often suboptimal, and ignoring random drift in optimality models can lead to the exaggeration of network architectures. To this end, a hybrid model that emphasizes both adaptation and drift could offer a more complete picture of network evolution.

Here we explore whether a model implementing both the processes of adaptive interaction switching and random drift can obtain a realistic level of compartmentalization, node-degree distribution and nestedness in real antagonistic networks. To do so, we first incorporated adaptive interaction switching and random drift in a modified Lotka–Volterra model, with the resource consumption depicted by the Holling type II functional response for multiple species. Consumers in the model are allowed to not only selectively eliminate the unfit resources from their diets based on the benefits and abundances of these resources (i.e. an adaptive process), but also randomly try out new resources (i.e. random drift). This hybrid behavioural rule depicts adaptation as Alfred Russel Wallace's natural selection via the elimination of the unfit and random drift as the innovation in foraging behaviours, making it distinctive from other rules of adaptive interaction switching. This hybrid rule, thus emphasizes the interplay between adaptation and drift in driving the emergence of realistic network architecture. We go further by examining the sensitivity of model outputs to a wide range of initial structures and parameter values, representing a diverse evolutionary history, and model performance by comparing model predictions with observed modularity, nestedness and node-degree distributions from 61 real antagonistic networks collated from literature. This model not only highlights the adaptive nature of ecological networks (i.e. the importance of density-dependent behavioural regulation) but also the role of random drift in explaining the emergence of network structures.

2. Model and methods

Let us consider an antagonistic network, consisting of m resource species and n consumer species. The population dynamics of resource i is controlled by its own density-dependent recruitment minus the loss due to feeding by consumers, whereas the population dynamics of consumer j is governed by the increase rate due to exploiting resources (depicted by Holling's type II functional response) minus its mortality. This yields the following Lotka–Volterra resource–consumer model:

2.1

where R_i and N_j are the population size of resource i and consumer j, respectively; r_i and c_i the intrinsic growth rate and the density-dependent coefficient of resource i; d_j the mortality of consumer j. The last term depicts the functional response of resources to the exploitation of consumers [34–36], satisfying all rules for prey switching (equivalent to the last functional response listed in the table of Morozov & Petrovskii [37]). Specifically, the binary interaction matrix {a_ij} indicates whether resource i is exploited by consumer j (a_ij = 1) or not (a_ij = 0); the preference matrix {v_ij} depicts the probability of whether consumer j decides to exploit resource i once encountered; the benefit matrix {b_ji} represents the benefit received by consumer j from consuming an individual of resource i, also known as the conversion efficiency; h denotes a consumer's handling time of a resource individual and is assumed to be equal for all consumers on all resources (h = 0.1). To keep the model simple, direct competition within the same trophic level is ignored as its impact on population dynamics is often much weaker than antagonistic interactions of resource exploitation [6], thus emphasizing indirect resource competition.

We define the hybrid rule of adaptive interaction switch and random drift as follows. At each time step, a consumer stops wasting its limited foraging time on the resource that contributes the least to its fitness gain, and a randomly selected consumer also starts to exploit a new resource. This hybrid rule of interaction switch depicts the process of natural selection as the elimination of the unfit and the behavioural innovation of trying out new resources. Specifically, the above model (equation (2.1)) was numerically solved with a time step of 0.01; at each time step, a randomly selected consumer species stops exploiting the resource that contributes the least to its per capita population growth rate (i.e. the resource with the minimum non-zero numerator in the consumer's functional response, b_jia_ijv_ijR_i in equation (2.1)); at the same time step, a consumer starts to exploit a randomly selected new resource. This hybrid behavioural rule not only emphasizes the adaptive process that the consumer gradually improves its resource utilization efficiency by retaining highly beneficial and abundant resources and eliminating less beneficial and rare ones [36], but also allows for behavioural innovation that new resources can be exploited by consumers via the random drift of interactions.

In the simulation, the entries of binary interaction matrix were initially randomly assigned 0 or 1, with the number of interactions being equal to the observation from the real networks and also ensuring no isolated species in the network. Moreover, multiple values of parameters and initial interaction matrices were used to ensure the robustness of model predictions. When running the model we randomly assigned values to the model parameters and chose only to report the results for parameters that ensure the persistence of all species in the network (see the electronic supplementary material). The performance of the model was evaluated by comparing predictions with the observed modularity (i.e. the extent to which nodes cluster into compartments [38]), nestedness (i.e. specialist species interact with a subset of the partners of generalist species [7]) and node-degree distribution (depicting the asymmetry of network topology [39–41]) of 61 real networks (33 host–parasitoid and 28 plant–herbivore), collated from published materials. For each input of observed network size and interaction numbers, the model was run from t = 0 to 300 (each time unit of t includes 100 time steps).

The modularity of interaction matrix was calculated by using the software Netcarto based on simulated annealing [42], with the statistical significance tested by the null model F that has fixed row and column sums equal to the observations (similar to the SIM9 in [43]). Nestedness was measured by Aninhado v. 3.0 [44], with the significance tested by null models ER (random matrices with equal connectance to the observations) and CE (probability of an interaction equals (Sr_i/N_c + Sc_j/N_r)/2, where Sr_i and Sc_j are sums of row i and column j, N_c and N_r are the number of columns and rows, respectively). The predicted modularity and nestedness were the average of 200 matrices for t = 101, 102 … 300 after the dynamics of the network has reached its stable equilibrium, and the node-degree distribution calculated for the interaction matrix at t = 300. Reduced major axis regression was used to compare observed with predicted modularity and nestendess, and the Kolmogorov–Smirnov test was used to compare the observed with predicted node-degree distributions. We further used a general linear (statistical) model to partition the variance of observed modularity into proportions explained by network size and connectance, with/without considering the hybrid rule of interaction switch.

3. Results

Through the interaction switch, the modularity of a network gradually converged to a stable equilibrium similar to the observed modularity of the real network, and as illustrated in figure 1 this equilibrium is also independent of the initial interaction matrix and assigned model parameters (see also figure 2; for details see the electronic supplementary material, appendix S1). The predicted levels of modularity for the 61 real networks were not significantly different from their observed values as demonstrated for two sets of model parameters using the reduced major axis regression (figure 2a: slope = 0.95, t-test for y = x, t = −1.15, p = 0.26; figure 2b: slope = 0.97, t = −1.19, p = 0.24; see the electronic supplementary material, table S1 for the prediction and observation for each real network), with more than 90% variation of observed modularity explained by the model (r² > 0.9). Surprisingly, the trajectory of network structure was not necessarily evolving towards a higher level of modularity; rather, the predicted networks showed a significantly lower level of modularity than the modularity of initial random networks (t = 18.066, p < 0.001; t = 15.48, p < 0.001; figure 2c,d).

Open in a separate window

Figure 1.

The effect of initial network structure and model parameters on the dynamics of modularity over time. Simulations of (a,b) were run for the network from Thompson & Townsend [45] and Labunets [46] (PH22 and HP20 in electronic supplementary material, table S1), respectively, with parameters assigned as in Case II of the electronic supplementary material, figure S1. Simulations of (c) were run for the network PH22 using various demographic and benefit parameter ranges, while simulations of (d) were run for PH22 using different values of handling time. Dashed lines represent the observed modularity of the real network.

Open in a separate window

Figure 2.

Observed, initial and predicted modularity for 28 plant–herbivore and 33 host–parasite networks. (a,b) Observations versus predictions from different model parameters (Cases I and II in the electronic supplementary material, figure S1). (c,d) The comparison between these two modularity predictions and the initial random matrices at the beginning of each simulation. Predictions are the averages of the last 200 time units after the dynamics stabilizes. Error bars indicate double standard deviation. Solid lines are results from the reduced major axis regression ((Case I): slope = 0.95, r² = 0.897, t = −1.15, p = 0.26; (Case II): slope = 0.97, r² = 0.91, t = −1.19, p = 0.24) and the dashed lines indicate perfect agreement between predictions and observations (y = x).

Predictions of nestedness also agreed well with the observations. Of the 61 real networks, 55 were significantly nested when using the null model ER, and 49 were significantly nested when using the null model CE (electronic supplementary material, table S2). In addition, the predicted NODFs of the 61 networks were significantly higher than the nestedness of their initial random networks (t = −19.36, p < 0.001; t = −19.033, p < 0.001; figure 3c,d).

Open in a separate window

Figure 3.

Observed, initial and predicted nestedness for 28 plant–herbivore and 33 host–parasite networks. (a,b) Observations versus predictions from different model parameters (Cases I and II in the electronic supplementary material, figure S1). (c,d) The comparison between these two nestedness predictions and the initial random matrices at the beginning of each simulation. Predictions are the averages of the last 200 time units after the dynamics stabilizes. Error bars indicate double standard deviation. Solid lines are results from the reduced major axis regression (Case I: slope = 0.95, r² = 0.58, t = −0.977, p = 0.333; Case II: slope = 0.94, r² = 0.64, t = −1.458, p = 0.15; and the dashed lines indicate perfect agreement between predictions and observations (y = x).

Our model also produced the node-degree distribution resembling the observations (figure 4). The Kolmogorov–Smirnov test revealed no significant difference between the observed and predicted node-degree distribution for more than 98% of the real networks (inset of figure 4), with only one distribution of consumer degrees differing significantly from the observation (see the electronic supplementary material, table S1). Evidently, the model has successfully predicted, with admirable precision, the observed levels of compartmentalization, nestedness and node-degree distributions, suggesting that the model has captured the essential process of structural emergence in antagonistic networks, namely the adaptive interaction switch plus random drift.

Open in a separate window

Figure 4.

Comparison of the predicted and observed node-degree distributions of consumer species (a) and resource species (b) for the network from Townsend et al. [47] (PH19 in the electronic supplementary material, table S1). ‘Random’ is the degree distribution of the initial random network at the beginning of the simulation. Inset shows the distribution of the p-value of the Kolmogorov–Smirnov (ks) test for 61 real networks (see the electronic supplementary material, table S1). Predictions are made from Case II parameters assignment (electronic supplementary material, figure S1).

The general linear model that explains the observed level of modularity in these 61 networks by the number of resource species, the number of consumer species and the number of interactions in the network can also explain a substantial amount (52%) of modularity variation (F_3,57 = 22.5, p < 0.01; adjusted r² = 0.518). This still suggests an extremely high predictive power of our model, adding nearly 40% variation explained over the general linear model. Network size does not play a significant role in determining network architecture, here specifically modularity (less than 7% variance explained, F_2,58 = 3.102, p = 0.052; figure 5). Instead, the role of network size on network architecture is indirectly realized through affecting network complexity (i.e. explaining the number of interactions in a network, 63% variance explained, F_2,58 = 52.32, p < 0.001; figure 5). Network complexity alone is an important factor explaining compartmentalization (29% variance explained, F_1,59 = 25.86, p < 0.001) and when together with network size can explain 52% variance of modularity (F_3,57 = 22.49, p < 0.001; figure 5). With network size and complexity as the input, our model that incorporated the hybrid rule of interaction switch can explain 90% variance of modularity (F_1,59 = 514.09, p < 0.001).

Open in a separate window

Figure 5.

A concept framework of the emergence of network architecture, based on general linear models of explaining the level of compartmentalization (modularity) by network size, complexity (the number of interactions; electronic supplementary material, table S1) and interaction switch. Dashed lines indicate factors not included in the model; double-lined arrows indicate model inputs; dotted lines indicate either parameter insensitivity (figure 1 and electronic supplementary material, figure S1) or insignificant variation explained. Percentage indicates variation explained (measured by the adjusted r² in the general linear model).

4. Discussion

Consumer–resource interaction is the mainstay of ecosystem function in food webs, both in antagonistic (e.g. parasitic and predation) and mutualistic (e.g. pollination and seed dispersal) networks [48,49]. The optimal foraging theory predicts that a consumer can maximize its energy intake rate by optimally allocating searching and handling time; in doing so, the consumer will only target highly beneficial resources and discard wasting time on low-benefit resources [24]. A recent update on the optimal foraging theory suggests that imprint from past experience and hunger aversion can also make consumers prefer abundant resources to those rare ones even if the abundant resources are less profitable than the rare ones [25,26]. The component of adaptive interaction switch in the model reflects this updated optimal foraging process of both profit seeking and hunger aversion as the decision of a switch depends on both the benefit (b_jia_ijv_ij) and the abundance (R_j) of exploited resources.

Although such a behavioural rule of adaptive interaction switching that can rapidly change interacting partners has been suggested a key force of shaping the structure of ecological networks [18,32,33,50], the optimization via adaptive interaction switching alone would lead to unrealistic extreme structures in ecological networks, such as perfectly nested structures in mutualistic networks or a higher level of compartmentalization than null model predictions in antagonistic networks, exaggerating the level of asymmetry in network topology. This highlights the counterbalancing role of random drift in pushing network structures back from the extremes to a realistic level. The hybrid behavioural rule not only included the adaptive process of eliminating the unfit interacting partners but also the random drift of interactions between species. Importantly, the simple model implementing both adaptive rewiring and random drift of interactions successfully predicted the observed level of compartmentalization and nestedness in real antagonistic networks (figures ​(figures11​1–3), with the proportion of variance explained drastically improved from the general linear model based on network size and connectivity. The hybrid behavioural rule of both adaptive rewiring and random drift could well be the key processes for structure emergence in real ecological networks.

Results from this model, together with progress made in literature, further allow us to elucidate key drivers behind the emergence of network architecture (figure 5). First, evolutionary history, community assembling process and environmental characteristics (e.g. productivity and heterogeneity) could largely determine the number of species (i.e. network size) and the composition/turnover of species and traits (i.e. model parameters) that a community can hold [4,5,30,51–54]. Second, extracting the contribution of network size and complexity, the hybrid behavioural rule alone can explain nearly 40% variation of observed modularity. Such a hybrid behaviour rule of adaptation and drift can occur rapidly at a pace even faster than the typical ecological time scale (e.g. host switch in parasites [21]) and is the most important determinant of network architecture (figure 5). The hybrid behavioural rule, together with network size and complexity, can account for more than 90% variation of network architectures. We suggest that this model almost perfectly explains fundamental network architecture as only 10% variance is unaccounted for—which could be due to many other stochastic factors or sampling artefacts.

The conceptual framework proposed in figure 5 also pinpoints two future research directions: (i) factors determining the number of species that a community can hold (i.e. network size); and (ii) mechanisms of how network size indirectly affects network architectures through directly dictating network complexity. To this end, the interplay of environment, evolutionary process (e.g. adaptation limits) and community assembly rule (e.g. species packing of generalist/specialist) could together determine the ceiling of species richness and the interaction complexity within a community. Studies on species packing in local and regional communities (e.g. on Darwin's naturalization hypothesis and the integration of alien species in native species assemblages) could shed light on these research directions [4].

For a long time, ecologists have sought a solution to the diversity–stability debate [55] with comparisons such as binary versus weighted interaction matrices, antagonistic versus mutualistic interactions, and random versus non-random interactions [3,6,11,14,22,33,55,56]. However, most early models in this debate depict a rigid network with a constant interaction matrix, with Lyapunov stability a proxy of ecological stability, portraying whether the perturbation caused by small changes in population sizes amplifies or dampens [14]. This is inconsistent with the dynamic and adaptive nature of ecological systems. The interaction matrix in recent models of adaptive interaction switch is constantly changing, reflecting an adaptive network and ecosystem converging or responding to changes [22,27]. By allowing consumers to readjust their exploited resources via updating the interaction matrix, our model successfully captured the essence of the structural emergence and the process of stabilization in antagonistic networks (see the electronic supplementary material, figure S7), suggesting that a network model implementing all relevant evolutionary processes is a better proxy than a rigid network for ecological communities and that the adaptive interaction switch is important for forecasting the response of ecosystems to environmental changes and perturbations [18,57], without ignoring the behavioural innovation from random drift that can further broaden the evolutionary trajectory. To this end, we should also consider structural stability [58] or more generally evolutionary stability of ecological networks. Since Lyapunov and evolutionary stability reflect different aspects of interaction networks, they have different implications for understanding network resilience. This is highlighted in evolutionary invasion analysis (e.g. adaptive dynamics) where ecological and evolutionary stability are clearly distinguished and handled separately (e.g. [59]). The complexity–stability debate will be better resolved once these two stability concepts are differentiated for the unique value that each captures in describing interaction networks and the resilience thereof.

The adaptive interaction switch allows the abundance of species to fluctuate without necessarily leading to extinction when facing perturbations [21] and is a vital adaptive behaviour that can rebalance the network back to its equilibrium (figure 1) and buffer ecosystems against perturbations [22]. It is, thus, a strong structuralizing and stabilizing force that compartmentalizes antagonistic networks, which can also form nested structures in mutualistic networks [33,36]. Species do select and adjust which other species to interact with, in response to changes in ambient environment and resource availability [21–23]. Fast ecological processes (e.g. density- and benefit-dependent interaction switching, selective foraging and importantly randomly drift), together with the slow evolutionary processes behind network size and complexity (e.g. assemblage history and the coevolution of trait complementarity), are dominant forces that give rise to the realized network architecture. The lack of either one will seriously hamper our understanding and predictive power of how ecosystems function. Our results highlight the need to incorporate random drift in current models based on adaptive interaction switch for maximizing the predictive power. Lessons should be learned from population genetics that considers evolution as driven by adaptation (natural selection plus mutation), genetic drift and gene flow. Future network models could further consider the process of colonization, extinction and speciation to initiate a shift in modelling from closed to open and adaptive systems.

Supplementary Material

Acknowledgements

References