Adaptive foraging fosters food web stability in environments disturbed by bioaccumulative pollutants
Authors: Constanza Vega-Olivares, Rodrigo Ramos-Jiliberto, Pablo Moisset de Espanés, Andreu Rico
Abstract
Chemical pollution poses a significant threat to aquatic ecosystems by affecting species abundance and destabilising community structure. Pollutants with bioaccumulative properties are transferred through food webs, influencing species across trophic levels and potentially disrupting ecological stability. Whilst the deleterious effects of pollution are well-recognised, the mechanisms by which adaptive foraging may interact with pollutant exposure to mitigate these impacts at the food web scale remain unclear. The present study employs a multi-species mathematical model to evaluate the efficacy of consumers’ adaptive foraging behaviour in maintaining species persistence and biodiversity within food webs of varying structural complexity exposed to a bioaccumulative pollutant. The hypothesis that adaptive foraging can alleviate the negative effects of pollutants was investigated. The results of the study demonstrated that adaptive foraging enables consumers to avoid highly contaminated preys, thereby enhancing community stability. Moreover, adaptive behaviour fosters the persistence of extensive and intricate food webs, even in the presence of pollutants. The findings of the present study offer novel insights into the dynamic interplay between chemical disturbance and behavioural adaptation, thus underscoring the potential of adaptive foraging as a stabilising mechanism in ecotoxicological contexts.
Supplementary Information
The online version contains supplementary material available at 10.1038/s41598-025-23950-8.
Introduction
Anthropogenic pollution is regarded as a significant contributor to global biodiversity loss^1–4^. Compounds of particular concern are those with high bioaccumulation potential and toxicity. This category includes hydrophobic substances, such as certain pesticides, polycyclic aromatic hydrocarbons, polychlorinated biphenyls, numerous metals, and certain pharmaceuticals. These substances have the potential to pose a serious threat to ecosystems^5^. The toxicity of these compounds impacts organisms through two primary direct bioconcentration from the exposure medium and trophic transfer via dietary interactions throughout the different trophic levels of the food web. The relative contribution of direct uptake from the environment versus dietary bioaccumulation to the internal concentration of the contaminant in the organisms is dependent on the hydrophobic properties of the compound, which are commonly characterised by its logKOW (KOW: octanol-water partition ratio). This ratio shifts from dominance of direct uptake from the exposure medium (e.g. water) for compounds with logKOW < 3, to dominance of dietary uptake for compounds with logKOW > 5, particularly in high trophic levels^6^. Consequently, chemicals with a logKOW ≥ 3 or logKOW ≥ 5, depending on the classification, are regarded as posing a high potential for bioaccumulation and biomagnification in food webs^7^.
Bioaccumulative pollutants can disrupt ecological communities by altering species’ biomass and their interactions, which particularly affects secondary consumers and top predators. These impacts extend to foraging behaviour and comprise top-down control mechanisms^8^. In communities exposed to bioaccumulative pollutants, trophic interactions serve as pathways through which both biomass and pollutants are transferred upwards in the food chain. However, if species exhibit adaptive foraging behaviour, i.e. the ability to dynamically adjust prey preferences to maximise energy intake, they may modulate not only energy transfer but also the pathways of pollutant transfer and bioaccumulation through the food web. This has the potential to affect multiple trophic levels and species depending on their foraging decisions. We hypothesise that when species adjust their foraging effort by optimising the per-capita food consumption over time, they may avoid highly polluted preys. This occurs because the abundance of the prey species tends to be suppressed by the deleterious effects of the pollutant, making polluted preys less available. Consequently, less polluted (i.e. healthier) preys should exhibit higher abundance and be preferred by adaptive consumers, thereby fostering species persistence, defined as the proportion of species surviving in the community over time, and biodiversity in food webs, commonly measured by changes in the Shannon-Wiener diversity index (H’), a well-established biodiversity metric reflecting community evenness and richness^9^. Previous studies have demonstrated that foraging adaptation enhances the overall resilience of ecological communities in unpolluted scenarios^10–12^. However, how the relationship between adaptive behaviour and chemical pollution affects the stability of food webs remains to be understood.
Theoretical work demonstrates that the structural properties of food webs, particularly species richness and connectance, influence how ecosystems respond to anthropogenic disturbances^13–17^. Pioneering studies such as those of suggested. Former studies evaluating this issue suggest that greater complexity might lead to destabilisation of ecological networks^16,17^, while subsequent studies produced mixed results^14,18–20^. More recent empirical and modelling research indicates that topological complexity (i.e. high species richness and connectance) may enhance resilience, especially in the context of chemical stress^21–23^. Nevertheless, the ways in which adaptive foraging interacts with food web complexity to modulate community stability to pollutant´s exposure remain unclear. We hypothesise that increased complexity will facilitate more efficient adaptive foraging, thereby reducing pollutant bioaccumulation in top predators and promoting species persistence and biodiversity.
This study examines the effectiveness of adaptive foraging behaviour, captured as optimised changes in prey preference, in promoting the persistence and diversity of species in food webs subjected to the disruptive effects of a bioaccumulative pollutant. Our main research question is whether adaptive foraging behaviour mitigates the impact of chemical pollution and how this effect varies in food webs with different topological features. To address this question, we employed a mathematical model of community dynamics that incorporates a chemical fate and a bioaccumulation component, which allows the simulation of pollutant uptake and transfer across trophic levels. In our model, predators can dynamically change their foraging preference on each of their preys to maximize food intake, thus potentially reducing pollutant uptake via food. Our specific research objectives (i) to evaluate the effect of the exposure to a hydrophobic pollutant on the relationship between adaptive behaviour and the stability of food webs, and (ii) to assess whether species richness and connectance influence the interactive effects of adaptive behaviour and chemical pollution on the stability of complex ecological communities.
Methods
Community model
Species dynamics
Our model for species dynamics is based on the bioenergetic model proposed by Yodzis and Innes^24^ generalised for multi-species systems by Williams and Martinez^25^ and adapted for polluted environments by Garay-Narváez et al.^26^. For all species i, the model expresses change in its biomass density (Bi) over time 1\documentclass[12pt]{minimal}
\usepackage{amsmath}
\usepackage{wasysym}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{amsbsy}
\usepackage{mathrsfs}
\usepackage{upgreek}
\setlength{\oddsidemargin}{-69pt}
\begin{document}$$
:\frac{d{B}{i}}{dt}={\kappa:}{i}{r}{i}{B}{i}\left(1-\frac{{B}{i}}{{K}{i}}\right)-{x}{i}{B}{i}+y{\kappa:}{i}{x}{i}{B}{i}\sum\limits{j\in:{preys}{i}}{F}{i,j}\left(B\right)-\frac{y}{f}\sum\limits_{j\in:{consumers}{i}}\frac{{x}{j}{F}{j,i}{\left(B\right)B}{j}}{{\epsilon:}_{j}},
where the first term represents the logistic growth rate of producers. Specifically, *r*~*i*~ represents the maximum mass-specific production rate, set to 0 for consumer species (i.e., trophic level *TL > 0*), and to 1 for primary producers (*TL = 0*). The carrying capacity of species *i*, *K*~*i*~, is defined as *K*~*TOT*~
*/ n*~*p*~, where *K*~*TOT*~ is the system’s carrying capacity and *n*~*p*~ is the number of producer species in the food web. *K*~*TOT*~ was set to 5^27^. The *κ*~*i*~ function represents the harmful effect of the pollutant on the growth rate of species *i*, calculated 2\documentclass[12pt]{minimal}
\usepackage{amsmath}
\usepackage{wasysym}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{amsbsy}
\usepackage{mathrsfs}
\usepackage{upgreek}
\setlength{\oddsidemargin}{-69pt}
\begin{document}$$
\:{\kappa\:}_{i}=1-\frac{{A}_{i}}{{A}_{i}+{\gamma\:}_{i}{B}_{i}+\eta\:},
$$\end{document}
where *A*~*i*~ is the internal concentration of the pollutant in organisms of species *i*, as outlined in the Pollutant dynamics section. Parameter *γ*~*i*~ is defined as the intrinsic sensitivity of the species to pollutant accumulation, i.e. the growth of species *i* is reduced to a 50% of its value when *A*~*i*~ = *γ*~*i*,~ and was set to 100 for all the species^26^. Parameter *η* guarantees the definition of the function *κ*~*i*~ and was set to 1.
The second term of Eq. (1) represents the exponential decrease of biomass for consumers in the absence of trophic interactions, where *x*~*i*~ is the mass-specific metabolic rate of species *i*, set to 0 for producer species (*TL = 0*), and calculated by Eq. (3) for 3\documentclass[12pt]{minimal}
\usepackage{amsmath}
\usepackage{wasysym}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{amsbsy}
\usepackage{mathrsfs}
\usepackage{upgreek}
\setlength{\oddsidemargin}{-69pt}
\begin{document}$$
\:{x}_{i}=\frac{{a}_{x}}{{a}_{r}}{\left(\frac{{M}_{C}}{{M}_{P}}\right)}^{-0.25},
$$\end{document}
where the allometric constants *a*~*x*~ and *a*~*r*~ were set to 0.314, and 1, respectively (Brose et al. 2006). The body mass of producer species, *M*~*P*~, was set to 1. The body mass of consumer species, *M*~*C*~, was assumed to increase with trophic level considering a constant predator–prey body mass *M*~*C*~
*= 100*^*TL*^^28^.
The third and fourth terms of the equation represent the increase and decrease in biomass, respectively, due to predator-prey interactions. The maximum consumption rate of species *i* when consuming prey species *j*, *y*, is defined as *y = a*~*y*~*/a*~*x*~, where the allometric constant *a*~*y*~ was set to 2.512^28^. The conversion efficiency of consumed resources to biomass of consumer *i*, *ε*~*j*~, was set to 0.45 for herbivores and 0.85 for carnivores^28^. Ingestion efficiency of the consumed resources to biomass of consumer species *i* is represented by *f* and was set to 1 for every species^28^. The set *preys*~*i*~ denotes all species that are preyed upon by species *i*, whereas *consumers*~*i*~ denotes all species that feed on species *i*.
The functional response of predator *i* when consuming prey species *j*, *F*~*i, j*~*(X)*, was calculated as 4\documentclass[12pt]{minimal}
\usepackage{amsmath}
\usepackage{wasysym}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{amsbsy}
\usepackage{mathrsfs}
\usepackage{upgreek}
\setlength{\oddsidemargin}{-69pt}
\begin{document}$$
\:{F}_{i,j}\left(X\right)=\frac{{\alpha\:}_{i,j}{X}_{j}^{h}}{{B}_{0}^{h}+{B}_{0}^{h}{\sum\:}_{k\in\:{consumers}_{j}}{c}_{k,i}{X}_{k}+{\sum\:}_{k\in\:{preys}_{i}}{\alpha\:}_{i,k}{X}_{k}^{h}},
$$\end{document}
where *X* represents the vector of either biomass density or internal pollutant concentration, depending on the context (see Pollutant Dynamics section), *B*~*0*~ is the half-saturation constant, set to 0.5^29^. A type III functional response was assumed (*h = 2*). Predation interference of the species present in the food web on species *i*^28^ was represented as *c*~*k, i*~, where *k* and *i* may or may not be the same species (i.e., denoting intra- and inter-specific interference, respectively). The value of *c*~*k, i*~ was set to 1 for every combination of species, implying equal predation interference across all species. Foraging preferences *α*~*i, j*~ are associated with the preference of predator *i* for prey species *j*, and they are modelled by the following replicator equation^30^.5\documentclass[12pt]{minimal}
\usepackage{amsmath}
\usepackage{wasysym}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{amsbsy}
\usepackage{mathrsfs}
\usepackage{upgreek}
\setlength{\oddsidemargin}{-69pt}
\begin{document}$$
\:\frac{\text{d}{\alpha\:}_{i,j}}{\text{d}t}={v}_{i}\left(Fi{t}_{i,j}-\sum\limits_{k\in\:{preys}_{i}}{\alpha\:}_{i,k}Fi{t}_{i,k}\right),
$$\end{document}6\documentclass[12pt]{minimal}
\usepackage{amsmath}
\usepackage{wasysym}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{amsbsy}
\usepackage{mathrsfs}
\usepackage{upgreek}
\setlength{\oddsidemargin}{-69pt}
\begin{document}$$
\sum\limits_{j\in\:{preys}_{i}}{\alpha\:}_{i,j}=1,
$$\end{document}
where *v*~*i*~ is the adaptation rate of species *i*. *Fit*~*i, j*~ represents the fitness associated with species *i* when consuming species *j*. In our implementation, *Fit*~*i, j*~ is calculated as *Fit*~*i, j =*~
*κ*~*i*~·*x*~*i*~·*y*·*F*~*i, j*~*(B).* This definition reflects the combined effects of species-specific traits and population size on fitness. Since *κ*~*i*~ is constant across all prey *j* on a given predator *i* at each time step, and *x*~*i*~ and *y* are also constants with respect to prey *j*, these multiplicative factors scale all prey fitness values equally for each predator. As a result, the relative differences in fitness used to update foraging preferences depend solely on the values of *F*~*i, j*~*(B)*. The second term (\documentclass[12pt]{minimal}
\usepackage{amsmath}
\usepackage{wasysym}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{amsbsy}
\usepackage{mathrsfs}
\usepackage{upgreek}
\setlength{\oddsidemargin}{-69pt}
\begin{document}$$
\:{\sum\:}_{k\in\:{preys}_{i}}{\alpha\:}_{i,k}Fi{t}_{i,k}
$$\end{document}; Eq. 5) describes the fitness of species *i* when it spreads its foraging effort over every prey species. The initial preference of species *i* for each of its prey species is equal to *1 / n*, where *n* is the number of prey species. Species with a biomass smaller than 10^− 12^ were considered extinct, and their biomasses were instantaneously set to 0.
It is important to note that, since we followed Brose et al.^28^ to parameterise the ecological rates, in which the ecological rates are normalised by the *r*~*i*~ of the smallest producer species and, therefore, the units of parameters and rates are eliminated, in our study parameters and rates are unitless.
#### Pollutant dynamics
The internal concentration of the pollutant in the organisms of species *i*, denoted by *A*~*i*~, was measured as the total concentration of pollutant within the population, without dividing by *B*~*i*~. The dynamics of *A*~*i*~ over time (Eq. 7), *dA*~*i*~*/dt*, were modelled using the equation proposed by Garay-Narváez et al.^26^, based on Kooi et al.^31^. This equation includes the following uptake from the exposure medium (first term), uptake via food (second term), metabolization of pollutant (third term), elimination of pollutant from the population by predation (fourth term), and depuration through excretion and egestion (fifth term).7\documentclass[12pt]{minimal}
\usepackage{amsmath}
\usepackage{wasysym}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{amsbsy}
\usepackage{mathrsfs}
\usepackage{upgreek}
\setlength{\oddsidemargin}{-69pt}
\begin{document}$$
\:\frac{\text{d}{A}_{i}}{\text{d}t}={\omega\:}_{i}C{B}_{i}+y{\kappa\:}_{i}{{x}_{i}B}_{i}\sum\limits_{j\in\:{preys}_{i}}{F}_{i,j}\left(A\right)-{x}_{i}{A}_{i}-\frac{y}{f}\sum\limits_{j\in\:{consumers}_{i}}\frac{{x}_{j}{F}_{j,i}{\left(A\right)B}_{j}}{{\epsilon\:}_{j}}-{\rho\:}_{i}{A}_{i}
$$\end{document}
It is important to note that the dynamics of bioaccumulation of pollutants vary between primary producers and consumers, due to the presence of physiological and trophic differences between these groups. Within the proposed model, the metabolic rate parameter *x*~*i*~ is set to zero for primary producers, reflecting the absence of metabolic loss in the absence of trophic interactions. Thus, fluctuations in the internal pollutant concentration *A*~*i*~ experienced by primary producers are solely attributable to uptake from the exposure medium, loss via predation, and depuration. Conversely, consumers and top predators have non-zero *x*~*i*~, meaning their pollutant dynamics include metabolic processes in addition to uptake, predation loss, and depuration. This distinction explicitly captures different pollutant accumulation pathways across trophic levels.
In Eq. (7), *ω*~*i*~ is the uptake rate of pollutant from the exposure medium of species *i* (excluding dietary pathways) and *C* is the concentration of pollutant in the exposure medium. The functional response *F*~*i, j*~*(A)* describes the feeding rate of predator species *i* on prey species *j* as a function of the internal pollutant concentration *A*~*i*~. The metabolization of pollutants is modelled as proportional to *x*~*i*~, following the approach of Kooi et al.^31^ and Garay-Narváez et al.^26^, based on the assumption, supported by Hendriks et al.^32^, that metabolic activity, as reflected by respiration rates, is a reliable proxy for pollutant turnover in aquatic organisms. Finally, *ρ*~*i*~ is the depuration rate of the pollutant (as the sum of excretion and egestion) for species *i*. Both *ω*~*i*~ and *ρ*~*i*~ were allometrically scaled considering the physicochemical properties of a hydrophobic pollutant (*logK*~*OW*~
*= 3*) and the species *i*’s trophic level, following the equations described by Hendriks et al.^32^. More information about the equations utilised to calculate the pollutant uptake and depuration rates by the modelled species can be found in the Supplementary Information (S1) online.
The concentration dynamics of the pollutant in the exposure medium, *C*, was modelled following Garay-Narváez et al.^26^ assuming a short-term release, chemical excretion from the populations, chemical uptake, and chemical decay by dissipation from the exposure medium (e.g. by microbial degradation, photolysis):8\documentclass[12pt]{minimal}
\usepackage{amsmath}
\usepackage{wasysym}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{amsbsy}
\usepackage{mathrsfs}
\usepackage{upgreek}
\setlength{\oddsidemargin}{-69pt}
\begin{document}$$
\:\frac{\text{d}C}{\text{d}t}=\varPi\:\left(t\right)+\sum\limits_{i}{\rho\:}_{i}{A}_{i}-\sum\limits_{i}{\omega\:}_{i}C{B}_{i}-\psi\:C,
$$\end{document}
The function denoted by *Π(t)* in the first term of Eq. (8) represents the input of the pollutant as a short-term release to the exposure medium (see Eq. 9). The second term of Eq. (8) describes the input of pollutant to the exposure medium from the organisms via depuration processes. The third term represents the loss of pollutant from the exposure medium due to uptake by organisms. The fourth term accounts for the removal of pollutant from degradation, with *ψ* representing the decay rate of the pollutant due to biological and physicochemical processes. Parameter *ψ* has been set to 10^− 4^, corresponding to a relatively persistent compound with a half-life of approximately 6931.5 time steps (*DT*~*50*~ *= ln(2)/10*^*− 4*^).9\documentclass[12pt]{minimal}
\usepackage{amsmath}
\usepackage{wasysym}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{amsbsy}
\usepackage{mathrsfs}
\usepackage{upgreek}
\setlength{\oddsidemargin}{-69pt}
\begin{document}$$
\:\varPi\:\left(t\right)=\frac{{10}^{P}}{\sqrt{\sigma\:\pi\:}}exp\left(-\frac{{\left(t-{t}_{peak}\right)}^{2}}{\sigma\:}\right)
$$\end{document}
The release of pollutant into the exposure medium increases over time until it reaches its peak at time *t*~*peak*~. After reaching the peak, the pollutant continues to be released, but the amount released decreases until the total pollutant concentration reaches *10*^*P*^. In the present study, *t*~*peak*~ is set to occur 4000 time-units after the initial time (*t = 0*). The selection of this value was made to ensure that the systems reached an equilibrium state prior to the introduction of the pollutant, thus isolating the effects of chemical disturbance from the baseline ecological dynamics. Parameter *σ* shapes the release curve. To simulate a rapid pollutant entry into the exposure medium, where releases initiates shortly before *t*~*peak*~ and the majority of the pollutant is released at *t*~*peak*~ (*Π(t*~*peak*~*) ~ 10*^*P*^), we set *σ* to 100, producing a narrow release pulse approximating a flush.
Both biomass density *B*, internal concentration of the pollutant *A* and pollutant concentration in the exposure medium *C*, are nondimensional variables, as it is demonstrated in the dimensional analysis provided by Garay-Narváez et al.^26^.
### Experimental design
The modelling experiment was conducted using a five-level logarithmic pollutant concentration *P* (2, 4, 5, 6, 8) scale, where the actual concentration is given by *10*^*P*^. Furthermore, a series of eleven distinct levels of adaptive behaviour were examined, as indicated by the adaptation rates *v*~*i*~: 10^− 5.6^, 10^− 4.8^, 10^− 4^, 10^− 3.2^, 10^− 2.4^, 10^− 1.6^, 10^− 0.8^, 0.25, 0.5, 0.75 and 1. We considered three levels of species richness (25, 50 and 100) and three levels of connectance (0.15, 0.2 and 0.3). Additionally, we considered the cases *Π(t) = 0* (i.e. no pollution) and *v*~*i*~
*= 0* (i.e. no adaptive behaviour) as controls. The model simulations were conducted with 40 synthetic food webs as replicates for each treatment, generated using the niche model proposed by Williams and Martinez^33^. The parameter values of the dynamic model were computed based on methods described by Brose et al.^28^, Garay-Narváez et al.^26^ and Hendriks et al.^32^.
In this study, topological complexity is defined as the structural characterisation of food webs in terms of their species richness and connectance. When referenced as *interaction between richness and connectance*, this term indicates the combined effect of both metrics included in the linear models to assess non-additive influences on stability variables.
Furthermore, additional simulations were performed on 100 food webs comprising 50 species and a connectance of 0.2, which were subjected to *P =* 5 and *v*~*i*~ ∈ {10^− 5.6^, 10^− 4^, 10^− 2.4^, 10^− 0.8^, 0.5, 1}, to evaluate the impact of adaptive behaviour in a polluted environment. The simulations were conducted in conjunction with the control cases *Π(t) = 0* and *v*~*i*~
*=0*.
The computational code for the model was written in the Julia programming language version 1.10.2^34^. To conduct the model simulations, initial biomass density values for each species in the network were generated randomly from a uniform distribution (*B*~*i*~(0) ∈ [0, 1]). At time *t = 0*, the internal concentration of the pollutant for each species in the network (*A*~*i*~) and pollutant concentration in the exposure medium (*C*) were initialised to zero. Initial prey preference values were calculated as *1 / n*~*i*~, where *n*~*i*~ denotes the number of prey species for the predator species *i*.
We used the *Tsit5* method implemented in the *DifferentialEquations.jl* package^35^ to solve the differential equations, setting an absolute tolerance of 1 × 10^− 6^ and a relative tolerance of 1 × 10^− 3^. The integration time span ranged from 0 to 2 × 10^4^. For each simulation, we recorded time series tracking the trajectory of biomass density, internal concentration of the pollutant, pollutant concentration in the exposure medium, and prey preference. Figure 1 illustrates a typical plot depicting the trajectory of biomass densities for a specific case, demonstrating convergence to a stable state within the chosen integration time.
Fig. 1Biomass dynamics over time for a 25 species food web, with a connectance value of 0.2, a pollutant pulse with a concentration of 5, and a release peak at t = 4,000 (dotted line), without adaptation behaviour (v = 0).
We evaluated the impact of treatments by analysing species persistence and the change in the Shannon-Wiener index as response variables. Species persistence was assessed by computing the ratio of the final number of species to the initial number of species in the food web. H’ is a metric of community diversity that accounts for both species richness and the evenness of biomass distribution among species. Higher values of H’ are indicative of more diverse and balanced communities, which often reflect greater stability and resilience. The change in the Shannon-Wiener index (ΔH’) was determined as the ratio of its value at the end of the simulation to its initial value, thus providing a quantification of shifts in community diversity over time. A ΔH’ value greater than 1 indicates an increase in diversity, suggesting enhanced stability, while a value less than 1 denotes a decline in diversity and possible degradation of community structure. H’ was calculated based on the biomass of each species at the respective time unit^27^.
The impact of the treatments on the different trophic levels was evaluated by grouping the species according to their trophic level, with producers (*TL =* 0), consumers (0 < *TL < max(TL)*) and top predators (*TL = max(TL)*) representing the three categories. The total biomass (∑*B*~*i*~) and the total internal concentration of the pollutant per biomass unit (∑(*A*~*i*~/*B*~*i*~)) were calculated for the three groups of species *i* at the times *t* ∈ {*t*~*peak*~; 8,000; 12,000}. Additionally, species persistence and ΔH’ were computed for each species group, including primary producers, consumers and top predators.
### Statistical analysis
Generalized Linear Models (GLM) with a normal distribution were used to assess interactive effect of the predictor variables on both response variables separately^36,37^. The analyses were performed utilizing the *GLM.jl* package in the Julia programming language^38^. The effects of the interactions between independent variables on the response variables were classified based on the magnitude and direction of the effects, as outlined by Piggott et al.^39^. The interaction categories additive, positive synergistic, negative synergistic, positive antagonistic and negative antagonistic effects. Further details can be found in Piggott et al.^39^. The full model selection results and performance metrics for these analyses are provided in Supplementary Table S1 (species persistence) and Supplementary Table S2 (ΔH’).
Additionally, linear regressions were performed to evaluate the impact of topology (species richness, connectance and the interaction of them) on the response variables of each combination of pollutant and adaptation level. These analyses were also performed utilizing the *GLM.jl* package in Julia.
## Results
The GLM results show that increasing levels of pollutant produce a significant decrease in both species persistence and ΔH’, while adaptive velocity had a significant positive effect on these parameters (see Table 1; Fig. 2). The analysis also demonstrated that the interaction between the concentration level and the adaptation rate had a negative, albeit comparatively weaker, effect on species persistence and ΔH’ than that of the pollutant concentration alone, thus resulting in positive antagonistic effects (see Table 1). The interaction between species richness and pollution stress had a significant positive effect on species persistence and ΔH’ (i.e., positive synergistic effect), compared to the negative effect of species richness alone on both metrics. The same positive effect on species persistence was observed for the interaction of connectance and pollutant concentration (i.e., positive synergistic effect), although both factors exhibited a negative effect on this response variable when considered separately. Conversely, we found a negative but weaker effect for the interaction between connectance and pollutant concentration on ΔH’ (i.e., positive antagonistic effect), compared to the negative effect of connectance alone. Furthermore, we found a positive effect of the interaction between species richness and connectance on species persistence (i.e. positive synergistic), as well as a negative but weaker effect of this interaction on ΔH’, in comparison to the effects of these factors separately (i.e. negative antagonistic). The interaction between species richness and adaptation rate exhibited a positive but weaker effect than adaptation rate alone on both species persistence and ΔH’ (i.e. positive antagonistic). We found a negative but weaker effect than connectance alone of the interaction between connectance and adaptation rate on species persistence (see Table 1; Fig. 2).
Fig. 2Species persistence (a; columns 1–3) and Shannon-Wiener index ratio (b; ΔH’ = final H’ / initial H’; columns 4–6) for food webs consisting of 25 (row 1), 50 (row 2), and 100 (row 3) species, with connectance levels of 0.15 (columns 1 and 4), 0.2 (columns 2 and 5), and 0.3 (columns 3 and 6), under varying pollutant concentrations (y-axis) and adaptation rates (curves). Dot marks indicate the mean value and error bars indicate the standard error. For clarity, only 7 of the 12 adaptation rates modelled in this study are shown.
Table 1Parameter estimates, standard errors (SE), *p*-values and effect’s classification for the GLM fitted to the effects of species richness, connectance, pollutant concentration and adaptive behaviour on species persistence (top) and ΔH’ (bottom) data.Response variableParameterEstimateSE*p*-valueClassificationSpecies persistenceIntercept1.250.0108< 0.001PositiveSpecies Richness (S)−0.006911.33e−4< 0.001NegativeConnectance (C)−1.170.0453< 0.001NegativePollutant (P)−0.1030.00163< 0.001NegativeAdaptive Velocity (A)0.1190.0137< 0.001PositiveS:C0.003930.00053< 0.001Positive synergisticS:P0.0005811.27e−5< 0.001Positive synergisticS:A0.0009911.01e−4< 0.001Positive antagonisticC:P0.02440.00634< 0.001Positive synergisticC:A−0.1810.0503< 0.001Negative antagonisticP:A−0.01570.0012< 0.001Positive antagonisticΔH’Intercept1.120.0108< 0.001PositiveSpecies Richness (S)-0.00231.33e-4< 0.001NegativeConnectance (C)−0.4610.0455< 0.001NegativePollutant (P)−0.04410.00164< 0.001NegativeAdaptive Velocity (A)0.03830.01380.00533PositiveS:C−0.001215.32e−40.0233Negative antagonisticS:P0.0002461.27e−5< 0.001Positive synergisticS:A0.0007721.01e−4< 0.001Positive antagonisticC:P−0.1270.00636< 0.001Negative antagonisticC:A−0.01950.05040.700AdditiveP:A−0.007910.00121< 0.001Positive antagonistic
It was observed that the interaction between adaptation rate - which had a positive effect on species persistence and ΔH’ - and species richness was positive on both response variables. However, these effects were still weaker than that of adaptation alone (i.e., positive antagonistic effect). A negative but weaker effect was observed for the interaction between connectance and adaptation rate, in comparison to the negative effect of connectance alone on species persistence (i.e. negative antagonistic effect; see Table 1; Fig. 2).
Interactions between three and four factors were not considered, as they did not explain more variance than models considering up to two-factor interactions when evaluating the goodness of fit for both response variables in our study (see Supplementary Information S2 online).
The analyses of the impact of adaptive behaviour on the different trophic levels show lower biomass densities of primary producers as adaption rate increases for both unpolluted and polluted scenarios at the three sampled times *t* (Fig. 3a, d, g). Higher biomass densities of consumers were observed as adaption rate increases for both pollution scenarios and sampling times (Fig. 3b, e, h). No differences were observed for the biomass of top predators at different adaptive behaviour levels (Fig. 3c, f, i). Regarding the internal concentration of the pollutant per biomass unit, the results show a decrease on its value as adaption rate increases for both primary producers and consumers, but no differences were found for top predators (Fig. 4). Furthermore, we found that the values of both species persistence and ΔH’ for consumers are higher as adaption rate increases (Fig. 5b, e). Conversely, no differences were observed for either species persistence or ΔH’ of primary producers and top predators as adaption rate increased (Fig. 5a, c, d, f).
Fig. 3Total biomass of producer (left), consumer (middle) and top predator species (right) in time step *t* equal to 4000 (top); 8000 (center) and 12,000 (bottom). Results for pollutant level of 0 (i.e., no pollutant stress) and 5 (i.e., moderate pollutant stress) are shown. Note that *t =* 4000 corresponds to the pollutant release peak, *t*~*peak.*~
Fig. 4Total internal concentration of pollutant per biomass unit (Σ(*A*~*i*~/*B*~*i*~)) of producer (left), consumer (middle) and top predator species (right) in time step *t* equal to 4,000 (top); 8,000 (center) and 12,000 (bottom). Results for pollutant level of 5 are shown.
Fig. 5Species persistence (top) and *ΔH’* (bottom) of the producer (first column), consumer (second column) and top predator species (third column) at *t =* 20,000. Results for pollutant level of 0 and 5 are shown.
When analysing the effect of topological complexity on the stability of food webs with different levels of pollution and adaptation, the linear regressions (see Supplementary Information S3 online) revealed that adaptive behaviour reduces the magnitude of the negative effect of species richness on species persistence when the pollution level is 5 or lower, and the effect was found to be maximal when the adaptation rate was 0.25 (Supplementary Fig. S1.a). In addition, and independent of the adaptation rate, we observed that the magnitude of the negative effect of species richness on species persistence decreased as the pollution level increased. In contrast, adaptation was found to increase the magnitude of the negative effect of connectance on species persistence for pollution levels equal to or less than 5 (Supplementary Fig. S1.b). Furthermore, it was observed that connectance had a greater negative effect on species persistence as pollution levels increased within the range of 0 to 5, but this trend was reversed for pollution levels greater than 5. We found no significant effects of adaptive behaviour on the interaction between species richness and connectance along the different pollution levels for species persistence (Supplementary Fig. S1.c). Conversely, the results show that the positive effect of the species richness-connectance interaction on species persistence is stronger as the pollution level increases, but there is a threshold at *P* = 5, after which it becomes weaker. For ΔH’, we observed a weaker negative effect of species richness with increased levels of pollution, but only for *P ≤* 4 (Supplementary Fig. S1.d). In addition, the results show a negative effect of species richness alone on ΔH’ for each pollution level. We found no effects of the interaction between connectance and adaptive behaviour on ΔH’, but Supplementary Figure S1.e shows a negative effect of connectance alone on ΔH’, which becomes stronger as the pollutant stress increases. Regarding the effect of the interaction between species richness and connectance on ΔH’, we found no relationship between both the topological complexity and adaptive behaviour, and topological complexity alone on ΔH’ (Supplementary Fig. S1.f).
## Discussion
The results of our study indicate that increasing exposure levels to a hydrophobic compound have a significant negative impact on the stability of complex food webs, affecting species persistence and diversity (Supplementary Fig. S2 online). These findings are consistent with those of previous modelling^21,31^ and micro- and mesocosm studies^40–44^. Additionally, our findings corroborate the hypothesis that foraging switching reduces the deleterious impact of a bioaccumulative pollutant on both species persistence and diversity, thereby exerting a stabilising effect on aquatic food webs.
The results of the analyses conducted for each trophic level in this study demonstrate that adaptive behaviour enhances both the persistence of species and the ΔH’ of consumer species (i.e., species with trophic level 0 < *TL* < *max(TL)*) in both unpolluted and polluted scenarios. We expected adaptive behaviour to affect the persistence and ΔH’ of top predators, but our findings show no effect of adaptive behaviour on either of the response variables. However, this discrepancy can be explained by the extinction of all the top predator species at the end of the integration time in the majority of the simulations with *P* = 5.
Furthermore, we found a reduction in the consumers’ internal concentration of the pollutant as the adaptation rate increased. In the model under consideration, an increase in the adaptive foraging behaviour of predators tends to result in a preferential targeting of prey species (including producers) that possess higher biomass and lower bioaccumulated pollutant concentrations. Consequently, when consumer persistence is augmented by adaptive foraging behaviour, predation pressure on primary producers escalates, resulting in a decline in their overall biomass. Given the dependence of bioaccumulation in producers on both exposure and population size, the reduction observed in producer biomass is indicative of a decrease in total pollutant load within the trophic group. Consequently, this effect propagates through the food web via dietary uptake, indirectly reducing the total pollutant burden in consumers. Conversely, the consumption of species with higher biomass and lower bioaccumulated chemical concentrations enables predator species to maintain the metabolic levels required to cope with low and moderate levels of pollutant stress (*P* ≤ 5), which involves the expenditure of energy to metabolise the pollutant. However, despite the stabilising effects observed at the consumer level, the simulations frequently resulted in the extinction of top predator species due to pollution stress, irrespective of the presence of adaptive foraging behaviour. This finding suggests a limitation of their capacity to regulate pollutants and trophic dynamics. This decline in top predators is likely to result in alterations to the overall structure of the food web, which may in turn diminish the potential for pollutant biomagnification at the highest trophic levels. This, in turn, may influence the long-term stability and distribution of pollutants in the ecosystem.
Different outcomes may emerge in the absence of adaptive foraging. Without foraging switching, a larger number of consumer and top predator species tend to persist after the stabilization of species biomasses (i.e. when biomass change rates approach zero). This contrasts with communities exhibiting adaptive foraging, where predator preferences shift towards less polluted and more abundant prey. In non-adaptive systems, predators do not adjust their foraging patterns in response to changes in prey pollution levels or abundance, often resulting in continued consumption of highly polluted or less abundant prey. Consequently, pollutant bioaccumulation may be higher in upper trophic levels, and predation pressure on primary producers is reduced, affecting both biomass distribution and pollutant dynamics within the food web. These alternative dynamics underscore the critical role of adaptive foraging in mitigating pollutant transfer and maintaining stability in complex ecological communities under chemical stress.
It is important to note that our study was conducted using a model hydrophobic contaminant whose bioconcentration and bioaccumulation potential was determined based on its physicochemical properties and the body size and metabolic relationships obtained for the theoretical aquatic food webs tested here. Therefore, it is possible that the use of food webs of different structures or chemicals with different physicochemical properties could result in quantitatively different outcomes. For example, we assumed that the internal sensitivity to the model compound was equivalent (on a per unit biomass basis) for all the species in the food web. This allowed us to simulate a compound with a kind of narcotic mode of action that affects all species in a similar way. However, despite high levels of foraging adaptation, compounds with a specific mode of action, such as herbicides and insecticides, could lead to species extinction and food web collapse. Specifically, herbicides can directly impair or eliminate primary producers, thereby disrupting energy availability for higher trophic levels and consequently causing extinctions throughout the food web. On the other hand, insecticides can target consumer species, reducing the biomass density of herbivore and carnivore populations and destabilising trophic interactions, thereby weakening prey-predator dynamics. These effects can diminish community structure and biodiversity, leading to collapse even in the presence of adaptive behaviour.
Our study shows that increasing either species richness or connectance independently has a negative effect on food web stability. However, when considered together, these topological features enhance stability, providing empirical evidence that challenges the idea that greater complexity leads to reduced stability^16^. Discrepancies between theoretical and empirical studies may stem from factors such as environmental variability, changing interaction strengths and non-random food web architecture. In our model, we accounted specifically for adaptive foraging and structure based on the niche model, both of which contribute to increased stability in complex food webs under pollutant stress.
High species richness raised species persistence and diversity in polluted ecosystems, corroborating previous findings^21^. Similarly, higher connectance lessened the negative impact on species persistence with increasing pollutant stress, although diversity decreased more under increased pollutant stress. The combined effect of species richness and connectance remained positively associated with species persistence even at the highest levels of pollutant concentration, with complexity’s protective effect strengthened at low to moderate levels of pollution. These results emphasize the importance of ecological complexity for biodiversity conservation under pollution stress.
When assessing the interplay between food web complexity and adaptive foraging behaviour, we found that the destabilising effect of species richness in undisturbed communities decreased with increasing adaptive rates. This was observed considering both species persistence and diversity. The detrimental effect of species richness on the species persistence and diversity of adaptive food webs was found to be reduced to a greater extent in the presence of pollution, but only for low and moderate pollution levels (*P ≤ 5*). In contrast, communities with higher connectance demonstrated a more adverse response to pollution in terms of species persistence as adaptation rate (i.e. the speed of foraging switching) increased. This indicates that the effect of connectance on persistence under pollution stress is modulated by the consumers’ adaptive foraging dynamics. Furthermore, our results suggest a strong positive effect of the species richness-connectance interaction on diversity when foraging switching is present, and that this effect is amplified when the adaptation rate is high. However, this phenomenon was observed only when pollution levels reached a certain threshold. Therefore, our study demonstrates that an increase in topological complexity can enhance the stability of food webs in response to chemical perturbations when species’ foraging preferences change over time. These findings are consistent with those of previous studies on the effects of adaptive behaviour on the feasibility of topological complexity, which have investigated both anti-predator strategies of prey^45,46^ and adaptive behaviour of predators, such as foraging switching^10,47–50^. Consequently, the incorporation of adaptive behavioural responses into mathematical models may be necessary to understand the stability patterns of complex food webs.
In summary, our study provides a foundation for predicting the impact of pollutants in complex ecological communities that integrate adaptive foraging behaviour. It shows that smaller and less connected food webs are more susceptible to chemical pollution than larger and more connected ones. This suggests that conservation efforts should be focused on areas with low biodiversity and high risk of contamination. Moreover, we show that the susceptibility of food webs to chemical stress is mitigated by the presence of adaptive changes in prey preference, a phenomenon that is more pronounced in simpler food webs. Therefore, adaptive behaviour serves to mitigate the severity of adverse effects in ecosystems experiencing biodiversity loss and increasing pollution levels. Our study supports the incorporation of the widely occurring phenomenon of adaptive behaviour in multi-species models and highlights the value of adaptive behavioural traits as a means of mitigating the effects of chemical stress and other forms of disturbance in ecological communities.
## Supplementary Information
Below is the link to the electronic supplementary material.
Supplementary Material 1