Data
Lake Balaton is a large (596 km2), 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 B1, Appendix D), and fishing yield (Weiperth et al. 2014) were updated to recent values using literature data from 20002020 (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 (Figure 1). In both Ecopath models, total biomass is comparable across trophic levels (measured in t/km2/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. Prebalance (PREBAL) diagnostics (Link 2010) were run to check the model's compliance with basic ecological principles (Appendix B1) and the Mixed Trophic Index (MTI) was obtained (Appendix C, Figure 3). The MTI table is used to compare with the loop analysis predictions.
Stella Model
The STELLA model was created in isee systems software (Figure 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 nonliving 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), RungeKutta 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 analyses (Venables and Ripley 2002; Pinheiro et al. 2013). For simulating sign predictions, we used the R code in Bodini and Clerici (2016). The figure was created with GVEdit Graph File Editor For Graphviz version: 1.02 and Graphviz version: 2.38 (Ellson et al. 2004). 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 preypredator (predatorprey) relationship. 0 means there is no trophic connection between two groups. The diagonal terms of the community matrix are selfeffects of system variables, represented in signed digraphs as links connecting variables with themselves. These links are selfdampening (circleheaded) with selflimiting growth rate, except detritus, because selflimitation was considered only for living groups (Table 2).
Table 2
Community matrix of the Loop Analysis model. Direct trophic interactions are represented in such a way that the elements in the rows impact the elements in the columns. The values can be − 1, 0, or 1. Selfdampening (diagonal elements) is only for living groups.

Pike

OFish

Bream

Mollusca

OBI

Zooplankton

Producers

Detritus

Fishing

Pike

1

1

1

0

0

0

0

0

1

OFish

1

1

0

1

1

1

1

1

1

Bream

1

0

1

1

1

1

0

1

1

Mollusca

0

1

1

1

0

0

1

1

0

OBI

0

1

1

0

1

0

1

1

0

Zooplankton

0

1

1

0

0

1

1

1

0

Producers

0

1

0

1

1

1

1

0

0

Detritus

0

1

1

1

1

1

0

0

0

Fishing

1

1

1

0

0

0

0

0

1

We followed the routine described in Bodini and Clerici (2016) to get the predictions for our network. The loop formula is used for calculating the equilibrium value of the variables following a perturbation, so it can be deduced how does the abundance of a certain variable change (Bodini 2000):
On the left side, xj is the variable with the equilibrium value being calculated and c is the changing parameter (e.g., mortality, fecundity, abundance). On the right side, f is the growth rate, ∂fi /∂c designates whether the growth rate of the ith variable is increasing or decreasing (positive or negative input, respectively), pji(k) is the pathway connecting the variable to the changed biomass variable (where the perturbation enters the system), Fn−k(comp) is the complementary feedback, which buffers or reverses the effects of the pathway and Fn designates the overall feedback of the system, which is a measure of the inertia of the whole system to change (Bodini 2000, Bodini and 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 and Levins 1985b; Raymond et al. 2011). The sign of the coefficients of [A]−1 gives the direction of the expected changes for the variables (Bodini and 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 * N2 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 E, Bodini and Clerici, 2016).
Using this routine, we obtained the matrix of predictions (Appendix E), and the version with directions only (Table 3). Table 3 was turned into a graph with the Graphviz software for better visualization (Figure 4). For direct comparability with Ecopath’s MTI, we used the directions of predictions.
Table 3
Table of predictions by Loop Analysis, only with directions, used for comparison of Loop Analysis and Ecopath’s Mixed Trophic Impact plot.

Pike

OFish

Bream

Mollusca

OBI

Zooplankton

Producers

Detritus

Fishing

Pike

+





+

+

+



0

0

OFish

0

+







0



+

+

Bream

0



+

0

0

0

+



+

Mollusca

0

0

0

+





0



0

OBI

0

0

0



+



0





Zooplankton

0

0

0





+

0



0

Producers

0

+



0

0

0

+



0

Detritus

0

+

+

0

0

0



+

+

Fishing



0

0

0

0

0

0

0

+
