Model analysis
For all simulations, we varied maximum nutrient availability\(N_{\max}\) from \(0.1\) to \(65\) µg P·L-1 in steps
of \(0.1\), ranging from oligotrophic to eutrophic conditions (Carlson
and Simpson 1996). Equilibrium densities were examined along a grid of\(p_{Z}\) and \(N_{\max}\) levels. We investigated the equilibrium
dynamics of the study system for the case of (1) fixed food preference
of \(Z\) for \(P_{E}\) vs. \(F\), covering the full range of possible
preference levels within the interval of [0,1], and (2) adaptive
food preference of Z.
For the fixed preference case , we numerically solved the system
of Eqs. 1-5 for transects along \(N_{\max}\) for fixed values of\(p_{Z}\). Transects were calculated for \(p_{Z}\) values from 0 to 1 in
steps of 0.01 using a standard Runge-Kutta solver (\(ode45\)) in MATLAB
R2018b (MathWorks). For each run, we integrated the system for 20000
time units, starting from the initial conditions\({N\left(0\right)=8,P}_{E}\left(0\right)=0.472\times 10^{-6}\ \),\(P_{I}\left(0\right)=8.68\times 10^{-6}\),\(F\left(0\right)=2.4\times 10^{-6}\), and\(Z\left(0\right)=0.21\). For the forward calculations, the
transects started at \(N_{\max}=0.1\). The system dynamics were
followed along increasing values of \(N_{\max}\) by using the solution
reached at the end of the previous run to be the initial condition for
the next run with slightly increased \(N_{\max}\) by 0.1. For each run,
the mean biomass of each state variable was calculated using the last
20% of total time steps. In case the system exhibited oscillatory
dynamics, we also calculated the standard deviation of the mean. A
species was considered to be extinct if its mean biomass over the last
20% of time steps was less than \(0.0001\) µg P·L-1.
In case of extinction at the end of the previous run, 0.0001 was added
to the initial biomass of the corresponding population in the following
run, allowing for re-invasion. The continuation was performed until
reaching \(N_{\max}=65\) µg P·L-1. Following the
same procedure, system dynamics were also continued for fixed values of\(N_{\max}\) along increasing values of \(p_{Z}\). All transects were
also followed in the opposite directions, i.e. along decreasing values
of \(N_{\max}\) and \(p_{Z}\), results were identical to the forward
calculations. Bifurcation analysis and specifically the detection of a
Hopf bifurcation was performed using MatCont 7p2 (Dhooge et al. 2006), a
software package usable within MATLAB.
For the adaptive preference case , system dynamics for Eqs. 1-6
were followed along increasing \(N_{\max}\) , using the same procedure
as for the fixed preference case, with starting values\(N_{\max}=0.1\) and \(p_{Z}(0)=0.5\), increasing \(N_{\max}\) by
steps of 0.1. System dynamics of Eqs 1-6 tended to reach equilibrium
dynamics faster than for the non-adaptive case (Eqs. 1-5), therefore,
for the adaptive preference case, system dynamics were only followed for\(10000\) time units.
Furthermore, we calculated the energy flow for different prey
preference and nutrient levels. The energy flow can be assumed to be
equivalent to the flow of the limiting nutrient (Rooney et al. 2006).
Given that the biomasses in our model system are expressed in terms of
the limiting nutrient, energy flow corresponds to the biomass flow along
each trophic link. We evaluated two measures of energy flow along each
trophic link: (1) net-energy gain of the respective consumer biomass,
and (2) gross energy flow. For example, the net energy gain of \(Z\)along the \(P_{E}-Z\) link is given by\(g_{P_{E}Z}=\frac{{e_{p}\cdot p}_{Z}\cdot toc\cdot P_{E}^{2}\cdot Z}{1+h_{P_{E}}\cdot p_{Z}\cdot toc\cdot P_{E}^{2}+h_{F}\cdot\left(1-p_{Z}\right)\cdot toc\cdot F^{2}}\),
while the gross energy flow along the \(P_{E}-Z\) link is given by\(\frac{g_{P_{E}Z}}{e_{P}}\). In case of oscillatory dynamics, we used
the mean biomass of all state variables for the calculation of energy
flow. The relative contribution of fungi to net energy gain of
zooplankton was calculated by\(\frac{g_{\text{FZ}}}{g_{P_{E}Z}+g_{\text{FZ}}}\times 100\). We also
calculated the transfer efficiency of the mycoloop by calculating the
ratio between the net energy gain of zooplankton from fungi and the net
energy gain of the inedible phytoplankton host from nutrient uptake\(\frac{g_{\text{FZ}}}{g_{NP_{I}}}\times 100\).