Comparison of three modelling frameworks for aquatic ecosystems: practical aspects and applicability

Freshwater ecosystems are under multiple stressors and it is crucial to find methods to better describe, manage, and sustain aquatic ecosystems. Ecosystem modelling has become an important tool in integrating trophic relationships into food webs, assessing important nodes using network analysis, and making predictions via simulations. Fortunately, several modelling techniques exist, but the question is which approach is relevant and applicable when? In this study, we compare three modelling frameworks (Ecopath, Loop Analysis in R, STELLA software) using a case study of a small aquatic network (8 nodes). The choice of framework depends on the research question and data availability. We approach this topic from a methodological aspect by describing the data requirements and by comparing the applicability and limitations of each modelling approach. Each modelling framework has its specific focus, but some functionalities and outcomes can be compared. The predictions of Loop Analysis as compared to Ecopath’s Mixed Trophic Impact plot are in good agreement at the top and bottom trophic levels, but the middle trophic levels are less similar. This suggests that further comparisons are needed of networks of varying resolution and size. Generally, when data are limiting, Loop Analysis can provide qualitative predictions, while the other two methods provide quantitative results, yet rely on more data.


Introduction
Ecosystem modelling has greatly evolved over the past decades with several software frameworks available to scientists and resource managers (Geary et al., 2020). Ecosystem models have been used as a tool to integrate ecosystem components with processes (e.g., trophic interactions) into a representative model. Here, we focus on network-based trophic models in which the components of a system are described in terms of nodes (e.g., species or functional groups) and their trophic interactions (e.g., diet matrix) (Jordán & Scheuring, 2004;Belgrano et al., 2005). These are simplified versions of an ecosystem, focusing on connecting consumer(s) with the resource(s) (e.g., predator-prey interactions). These interactions can be based on presenceabsence (i.e., connectance web) or weighted (quantitative) connections (Woodward et al., 2005). Each model has limitations and assumptions, under which the question of interest can be examined, which necessarily means that there is no perfect model. However, depending on data availability (e.g., quantitative vs qualitative), there are considerations on which modelling approach is most suitable for the research question in a context-dependent way.
Our objective is to compare three commonly used but rarely compared modelling frameworks (Ecopath, STELLA, Loop Analysis) using a case study of a general lake model. We discuss the data requirements and the potential of incorporating social systems (e.g., socio-ecological models). It is imperative to apply interdisciplinary, system-level thinking to the complex problems facing society today (Richmond, 1993;Saviano et al., 2019). As aquatic ecosystems are under increasing anthropogenic pressure and their biodiversity is threatened globally (Dudgeon et al. 2019;Sala et al., 2000), the need to better describe and understand these systems is urgently needed.
First, we provide a brief overview of the modelling frameworks. In Table 1, a summary table was compiled describing several practical aspects of the three frameworks and their main references. Then, we discuss the preliminary model outcomes using a case study of a freshwater lake by comparing the STELLA and Loop Analysis models to the Ecopath model.

Loop Analysis
Loop Analysis is a qualitative model (Levins, 1974), that provides a method that is useful where species and their natural history are well-known, but not quantified (Dambacher et al., 2003). Using the qualitative modelling framework of Loop Analysis, one can analyze pathways and feedback in the system, making predictions about the response of variables to perturbations. For example, these can be the addition (increased biomass) or deletion (decreased biomass) of other nodes. Based on feedback and pathways, one can qualitatively specify the direction of change. For getting these predictions, Loop Analysis uses differential equations (Bodini, 2000;Bodini & Clerici, 2016;Fábián, 2021).

Data
Lake Balaton is a large (596 km 2 ), shallow (3.25 m average depth), freshwater lake located in Hungary (Istvánovics et al., 2008). A freshwater lake model was created with 8 general functional groups and fishing pressure added, in three modelling frameworks: Ecopath, STELLA, and Loop Analysis. For Lake Balaton, our model ecosystem, an Ecopath model described in Bíró (2002) was taken as baseline. Functional groups were aggregated (Producers = Phytoplankton, Periphyton, Benthic algae; 7 fish species grouped into 'Other Fishes') to be more general. These functional groups are purposefully very general for simplicity and to provide the basis for the methodological comparison. In the next section, we describe the creation of each of the models.

Ecopath model
Diet matrix (Appendix A), biomass data (Appendix B-1, Appendix D), and fishing yield (Weiperth et al., 2014) were updated to recent values using literature data from 2000-2020 (manuscript in preparation). 8 nodes describe this lake system, of which six are aggregate functional groups and two nodes are for the main fish species in the lake ( Fig. 1). In both Ecopath models, total biomass is comparable across trophic levels (measured in t/km 2 /year). Total primary production decreased significantly over the two time periods (changing from eutrophic to oligotrophic state, Bernát et al., 2020), being the main ecological difference between the earlier (Bíró, 2002) and this newer Ecopath model. Pre-balance (PREBAL) diagnostics (Link, 2010) were run to check the model's compliance with basic ecological principles (Appendix B-1) and the Mixed Trophic Index (MTI) was obtained (Appendix C, Fig. 3). The MTI table is used to compare with the Loop Analysis predictions.

STELLA model
The STELLA model was created in isee systems software (STELLA, 2021; Fig. 2). We used the same input data (as in Ecopath) for the stocks (biomass of functional groups), flows (diet matrix, annual fishing yield), and converters (mortality) (Appendix D). The model was parametrized starting from the bottom up (i.e., producer stock and output flows, then invertebrate stocks and flows, fishes, and finally detritus), ensuring that the relative contribution of diet groups is accurately represented. Some additional assumptions had to be made in order to be able to run the model: 1) production input is constant (oligotrophic state), 2) natural mortality at 18 °C (average annual water temperature) was retrieved for Sander lucioperca and Abramis brama from FishBase (and their average taken for OFish), 3) for invertebrates, natural losses at the stock outflows were estimated to be highest for Zooplankton (75%), and lower for Mollusca (40%) and OBI (25%) (to balance inflows and outflows), 4) living to non-living flows are unknown (not parametrized) and are not connected to detritus (otherwise the detritus stock would accumulate these), except for producers. Further settings are the following: time step (DT = 1 year), Runge-Kutta 2 integration method.

Loop analysis model
The Loop Analysis model was created in R software (R Development Core Team, 2020). We used MASS 7.3-51.5 and nlme 3.1-148 R packages for the analyzes (Pinheiro et al., 2013;Venables & Ripley, 2002). For simulating sign predictions, we used the R code in Bodini and Clerici (2016). Optionally, a signed digraph figure can be created 1 3 using GVEdit Graph File Editor For Graphviz version: 1.02 and 2.38 (Ellson et al., 2004), but was not applied here. For making the community matrix, Ecopath model's diet matrix was used. The community matrix got a 1 (− 1) value where in the diet matrix was a prey-predator (predator-prey) relationship. 0 means there is no trophic connection between two groups. The diagonal terms of the community matrix are self-effects of system variables, represented in signed digraphs as links connecting variables with themselves.
These links are self-dampening with self-limiting growth rate, except detritus, because self-limitation was considered only for living groups (Table 2). We followed the routine described in Bodini and Clerici (2016) to get the predictions for our network (Table 3). The loop formula is used for calculating the equilibrium value of the variables following a perturbation, so it can be deduced how the abundance of a certain variable change (Bodini, 2000): On the left side, x j is the variable with the equilibrium value being calculated and c is the changing parameter (e.g., mortality, fecundity, abundance). On the right side, is the growth rate, ∂f i /∂c designates whether the growth rate of the i th variable is increasing or decreasing (positive or negative input, respectively), p ji (k) is the pathway connecting the variable to the changed biomass variable (where the perturbation enters the system), F n-k (comp) is the complementary feedback, which buffers or reverses the effects of the pathway and F n designates the overall − 0* 0* 0* 0* 0* 0* 0* + feedback of the system, which is a measure of the inertia of the whole system to change (Bodini, 2000;Bodini & Clerici, 2016). See also Puccia and Levins (1985) for the discussion of the correspondence between matrix algebra and Loop Analysis. A perturbation on variable j (in this case the perturbation is the increase in the biomass of j) has a net effect (the sum of the direct and indirect effects) on variable i given by the j -ith element of the inverse community matrix [A] −1 (see Levins, 1974;Puccia & Levins, 1985;Raymond et al., 2011). The sign of the coefficients of [A] −1 gives the direction of the expected changes for the variables (Bodini & Clerici, 2016). To make predictions, we used a routine that randomly assigns numerical values from a uniform distribution to the coefficients of the community matrix (these coefficients belong to the links of the signed digraph). This was performed 100 * N 2 times, where N is the number of variables in the system. Matrices satisfying the asymptotic Lyapunov criteria were accepted and inverted. The routine of Bodini and Clerici (2016) calculated predictions for the probabilities based on the percentage of positive and negative signs and zeroes in the inverted matrices. They defined a set of rules to make a final table of predictions only from signs (Appendix A in Bodini & Clerici, 2016) which is what we applied in this study (Table 3). Using this routine, from stable matrices we obtained the simulated tables of percentages of ± and 0 and the table of predictions generated from the tables (Table 3).

Results and discussion
The STELLA and Loop Analysis models were compared to the Ecopath model (which is the most comprehensive framework).

Comparison: ecopath and loop analysis
The Mixed Trophic Impact plot from Ecopath (Fig. 3A) could be compared with the predictions of Loop Analysis (Fig. 3B). The MTI quantifies how an infinitesimal increase of any of the impacting groups is predicted to have on the impacted groups (Christensen et al., 2005), while Loop Analysis gives qualitative predictions. The advantage of an MTI plot is also to highlight interactions (and system components) whose importance otherwise might not have been realized (Ulanowicz & Puccia, 1990), including direct and indirect interactions (e.g. competition) in the system (Christensen et al., 2005). For the MTI, a net impact matrix is constructed of all system components in which the element q ij indicates the net impact of i upon j. Specifically, the net impact (q ij ) will equal the difference between g ij (the amount that i serves as a prey item for j) and f ji (the detrimental impact the consumer i has on the resource j). Also, the indirect effects are taken into account (if i has an effect on j, and j influences k, then there is an indirect pathway between i and k), see Ulanowicz & Puccia, 1990. Fishing pressure is regarded as "predator", and for detritus g i,j is set to zero. MTI values fall within the range of -1 to + 1, inclusive. The MTI routine based on the net impact matrix by Ulanowicz and Puccia (1990) is automatically implemented in Ecopath and the MTI plot (and numerical values) can be accessed under the Ecopath output tab.
For the correlation, we excluded the self-effects due to methodological differences between MTI and Loop Analysis. MTI predicts self-limiting growth (negative effect), so if a group's abundance or biomass increases, it will not increase infinitely. In contrast, in the prediction matrix of Loop Analysis, the affecting groups are increasing as the representation of the perturbation.
We found a positive correlation between the percentage of positive predictions in Loop Analysis (LA + (%)) and the quantitative values of MTI (r = 0.62, Fig. 4). According to both methods, the groups with the most positive impact on  (Fig. 3). The strongest negative impacts in Loop Analysis matched those of the MTI plot (i.e., Pike-OFish, Producers-Detritus; Fig. 3). The predictions of Loop Analysis as compared to the outcome of the MTI are in good agreement at the top and bottom trophic levels, but the middle trophic levels (e.g., OBI) are less similar. Most of the directions of predictions are similar, but their strength in the MTI plot is stronger (Fig. 3). The quantitative method results in stronger strength of impacts (Fig. 3A) as compared to the qualitative method (Fig. 3B). The difference between the qualitative and quantitative methods is apparent if we observe the Fishing-Pike effect. Fishing has a singlestep direct negative impact on the group Pike (Figs. 1 and 4B), while in the quantified network the interaction strength makes the effect more subtle (Fig. 3A). Some of the visibly outlier points of the network (e.g.,  are nodes that represent the most distant groups within the network (having the most trophic steps between them), which means that their predictability is most difficult due to pathway redundancy (Ulanowicz, 1980;Ulanowicz & Puccia, 1990).

Comparison: ecopath and STELLA
Simple dynamic changes can be visualized (graphical and numerical outputs) in the STELLA model, for example, one could easily double the predation rate of pike on bream, which would change system dynamics. However, such predictions from the STELLA framework are not as comprehensive as Ecopath's MTI, therefore we did not compare them with each other. When focusing on trophic connections between the living groups, we found it possible to recreate Ecopath's trophic structure with our simple assumptions, meaning that the predicted biomass magnitudes are meaningful and comparable to Ecopath's biomass values shown in Fig. 1. For example, for fishes, biomass (t/km 2 /year) in the steady-state are comparable between the two frameworks (Bream: 16.88,OFish: 17.34 in STELLA at the end of the simulation; and Bream: 14.03, OFish: 16.64 in Ecopath). The STELLA model initially overestimated the biomass of pike, therefore an additional fixed loss of 2.4 t/km 2 /year was subtracted. The primary production is constant, which represents the overall production in the lake during an oligotrophic state (minor interannual deviations are negligible as long as the lake remains in an oligotrophic state). Ultimately, trophic efficiency (called "ecotrophic efficiency" in Ecopath, which is the proportion of the production utilized in the system (e.g. accumulating or moving onto higher trophic levels)) in our model was estimated to be high for producers (~ 90%) and invertebrates (~ 70%), and lower for fish and zooplankton (~ 20-30%) (meaning that primary production and invertebrates are readily consumed, while zooplankton and fish species are mainly lost via natural mortality and diseases), following best practices detailed in Heymans et al. (2016). For the invertebrate groups, both models showed deviation from the expected value (Ecopath gave high values for the Mollusca group (Appendix B-1), and STELLA predicted high values for the Zooplankton group 17.77 vs the input of 9.7 t/km 2 /year in Ecopath. This difference probably comes from the simplicity of processes (flows) of the STELLA model at the second trophic level. Fath et al. (2007) discuss the applicability of STELLA numerical simulation as an input to an Ecopath model. To test this, we plugged in the biomass results from the end of numerical simulation (when steady state was reached) into Ecopath (all other settings remained the same as in the original Ecopath model), resulting in comparable PREBAL diagnostics to the original Ecopath model (Appendix B-2). While the Ecopath framework allows for some uncertainty if data is limited (e.g., the software can estimate one input parameter per group), the STELLA model needs to be fully parameterized in order to run, which in some cases leads to inevitable oversimplification of the system (e.g., omission of living to detritus flows due to lack of quantitative data).

Conclusion
In our freshwater lake case study, we found that the output of Ecopath's MTI plot and the Loop Analysis model are in good agreement. Since Loop Analysis relies only on network topology (assuming no knowledge of interaction strengths in the input; Novak et al., 2011), and Ecopath needs more details of the network (quantitative interaction strength; Ulanowicz & Puccia, 1990), MTI can provide more robust results. Bondavalli et al. (2009) present a hybrid approach, based on a quantitative model of a low-resolution lagoon (7 nodes), in which the net impact matrix was amended with an additional variable (capturing exchanges between the ecosystem components and the surrounding environment), after which the loop analysis procedure was performed. The resolution of the food web can influence many aspects of network analysis (e.g., Martinez et al. 1991;Giacomuzzo & Jordán, 2021). The connectance of the network influences the proportion of correct predictions, thus the reliability of predictions diminishes as network size and connectance increase (Novak et al., 2011). As the number of groups increases within a network, the predictability of Loop Analysis decreases due to pathway redundancy (Ulanowicz, 1980;Ulanowicz & Puccia, 1990), thus it is best recommended for small networks (< 24 nodes, Novak et al., 2011).
Our findings agree with Fath et al. (2007), suggesting transferability between the two frameworks, STELLA and Ecopath. Our recommendation is, that when a system is well documented with quantitative data and the processes are clear, STELLA models can be a great way to better understand the system as a whole (Gertseva et al., 2004;Power et al., 1995;Xuan & Chang, 2014), to highlight important feedback loops (Hayes, 2012;Richmond, 1994), and even to raise awareness about an environmental problem (Jørgensen & De Bernardi, 1997, 1998, and educate the public with its interactive interface (isee Exchange, 2021). Cellier (2008) details further advantages and disadvantages of using the STELLA software for dynamic modelling. Due to the details of parametrization, STELLA is best applicable for small networks (< 10 nodes).
Both frameworks, Loop Analysis and STELLA can be used to complement an existing Ecopath model, for example, if time dynamics predictions are needed, but Ecosim models are not yet available. A limitation of our study is that we only used one small network. It would be important to compare these methods using several networks (of different size, resolution, etc.). Discrepancies between quantitative and qualitative approaches are not surprising and expected to some extent (e.g., less detail and/or information loss is inevitable in qualitative methods compared to quantitative methods). In general, data availability (e.g., lack of information on nodes or on interaction strength) influences food web properties and sampling effects need to be taken into consideration (Berlow et al., 2004;Goldwasser & Roughgarden, 1997). However, not every research question requires numerically parametrized models. For example, it can be qualitatively shown that the effect of an external factor (e.g., aquaculture) negatively impacts certain groups (primary producers, zooplankton, and deposit-feeders), and positively impacts predators and scavengers (Forget et al., 2020). Modelling is thus important to bring out complex, network-level interactions, which might not be evident simply from single parts. The challenge now is to extend these modelling frameworks to social-ecological systems (Martone et al., 2017;Niquil et al., 2021;Rodriguez et al., 2021).    (Link, 2010). Groups are 1 = Pike, 2 = OFish, 3 = Bream, 4 = Mollusca, 5 = OBI, 6 = Zooplantkon, 7 = Producers. Blue bars indicate values estimated by the model (e.g., biomass data was not feasible to estimate from survey data). Biomass for Mollusca is estimated by Ecopath to be high and Zooplankton biomass (based on survey data) is low. Biomass data were obtained from the best available data from years 2000-2020 for Lake Balaton, Hungary (modified from Bíró, 2002) Fig. 3). 1-8 are coded the same as in the MTI plot (Fig. 3), and 9 = Fishing Impacting group