Biodiversity decline causes a loss of functional diversity, which threatens ecosystems through a dangerous feedback loop: this loss may hamper ecosystems' ability to buffer environmental changes, leading to further biodiversity losses. In this context, the increasing frequency of human-induced excessive loading of nutrients causes major problems in aquatic systems. Previous studies investigating how functional diversity influences the response of food webs to disturbances have mainly considered systems with at most two functionally diverse trophic levels. We investigated the effects of functional diversity on the robustness, i.e. resistance, resilience and elasticity, using a tritrophic ---and thus more realistic---plankton food web model. We compared a non-adaptive food chain with no diversity within the individual trophic levels to a more diverse food web with three adaptive trophic levels. The species fitness differences were balanced through trade-offs between defense/growth rate for prey and selectivity/half-saturation constant for predators. We showed that the resistance, resilience and elasticity of tritrophic food webs decreased with larger perturbation sizes and depended on the state of the system when the perturbation occurred. Importantly, we found that a more diverse food web was generally more resistant and resilient but its elasticity was context-dependent. Particularly, functional diversity reduced the probability of a regime shift towards a non-desirable alternative state. The basal-intermediate interaction consistently determined the robustness against a nutrient pulse despite the complex influence of the shape and type of the dynamical attractors. This relationship was strongly influenced by the diversity present and the third trophic level. Overall, using a food web model of realistic complexity, this study confirms the destructive potential of the positive feedback loop between biodiversity loss and robustness, by uncovering mechanisms leading to a decrease in resistance, resilience and potentially elasticity as functional diversity declines.

To study the response of tritrophic food webs with different levels of diversity to a pulse perturbation we used several scripts, libraries and programming languages. Ordinary differential equations were solved numerically in C using the SUNDIALS CVODE solver with relative and absolute tolerances set to 10^-10. Output results were analyzed in Python using several packages among which NumPy, SciPy and Matplotlib.

We run simulations where we perturbed the system with an addition of nutrients and quantified different metrics (minimal and maximal biomasses, return time for each species and trophic level) to investigate the robustness (resistance, resilience and elasticity) of the food web. We used an extinction threshold of 10^-9 microg C/L to prevent numerical problems when state variables are extremely close to 0.

When the attractor was a limit cycle or a chaotic attractor, we investigated how the effect of the perturbation on the dynamics depended on the moment of perturbation i.e. where on the attractor the perturbation was applied. This means that for each point on the attractor, we perturbed the system and calculated the biomasses' minima and maxima as well as the return time following this perturbation.

For this, we sampled the different attractors of the different food webs in a high spatial resolution (points_attractors files). First, starting from an initial condition known to be in the basin of attraction of the relevant attractor, the system was allowed to settle for 10⁵ time units using a large timestep (10^-1 ) such that a point sufficiently close to the attractor could be obtained. If the attractor was a limit cycle, the system was further integrated for approximately one period using a high temporal resolution of 10^-3 This created a set of points S_t on the attractor. Finally, to sample the attractor in a way such that the distances between the sampled points on the attractor did not become very large when the dynamics were moving very fast, the attractor was interpolated and resampled such that the arc length between consecutive sampled points was equal to 1 (in units of biomass). For this set of spatially sampled S_x, the distance between a point on the attractor and the closest point to it in S_x was guaranteed to be smaller than 1 microg C/L. I

We defined the return zone as the volume in phase space inside which the distance to any point on the attractor was less than 5 microg C/L. Note that in the case of chaotic behavior the return time could not be calculated.

For details about the topology of the food webs or the detailed ordinary differential equations see the paper.

In each file the aim is precised at the top of the script. There are two types of file: the ones for producing figures (files starting by figures) and the ones for producing simulations. For the first ones data are available in the folder datasets and if the variable savefig is set to 1 the figure is saved in the folder figures. For the files producing simulations some important steps should be checked before running the code. The corresponding c codes should be compiled and the executable file should be named as precised in the top of the .py and .c files (same name as the c file). We used the following command in the terminal:

gcc -Wall -I/usr/include/gsl -L/usr/local/lib script_name.c -lm -lgsl -lgslcblas -lsundials_cvode -lsundials_nvecserial -g -o script_name

Create also a folder c_ouput_files where all the files written by C will be saved