Introduction

The Fussmann (2007) review states, “It is a long-standing debate among ecologists whether chaotic dynamics are likely to occur in food web systems.” Increasing results since the late 1990s, both theoretical and experimental, have supported the conclusion that chaotic and other complex dynamical interactions, are certainly possible and perhaps common (Constantino et al. 1997; Becks et al. 2005; Graham et al. 2007; Beninca et al. 2008; Becks and Arndt 2008, 2013; Vayenas and Pavlou 1999; Massoud et al. 2018). However, bringing together a well-defined experiment and a highly sensitive nonlinear, coupled mathematical model with variables and parameter values appropriate for the experiment is still lacking. Conducting measurements is difficult, especially in field experiments where the environmental conditions are not uniform or constant. Only laboratory-based experiments, which are still quite difficult to perform, have produced microbial population time series of sufficient length and precision to verify the occurrence of deterministic chaotic dynamics (Becks et al. 2005; Graham et al. 2007; Beninca et al. 2008).

Graham et al. (2007) studied a realistic mix of microbes in aerobic bioreactors filled initially with a mixture of wastewater from a treatment plant and simulated wastewater. The main variables recorded as time series were total bacteria, ammonia-oxidizing bacteria, nitrate-oxidizing bacteria and protozoa, along with the concentrations of nitrate, nitrite and total ammonia. Chaotic dynamics were identified by calculating positive Lyapunov exponents using the method of Rosenstein et al. (1993). The Beninca et al. (2008) experiments, which continued for 6.3 years, demonstrated chaotic dynamics in a complex plankton community in a water sample obtained from the Baltic Sea. The experiment was housed in a cylindrical mesocosm containing 90 L of water with a 10 cm sediment layer at the bottom. The behavior of the relevant time series was consistent with Lyapunov exponents averaging about 0.058 per day. Massoud et al. (2018) provided additional quantitative analysis of the Beninca et al. (2008) study dealing with the limits of predictability in data sets displaying deterministic chaos. Both experiments would be difficult to represent mathematically because of the large number of coupled variables involved and their many unspecified couplings.

As noted by Fussmann (2007), a particularly well-defined experiment was that of Becks et al. (2005). That team studied a one predator, two preys, microbial food web in a chemostat, with additional analyses presented during following years (Becks and Arndt 2008,2013). Depending on the feeding (dilution) rate, the experiments produced both classical steady states as well as chaotic dynamics. Fussmann (2007) stated that, “If this food web could be parameterized for a mathematical model that predicts its behavior correctly, ecologists would possess a magnificent system to test dynamical food web theory.” Thus, the main objectives of this paper are to develop a surrogate mathematical model of the Becks et al. (2005) experiments, use this model to analyze the experiments, explain why a surrogate model is the best that can be developed without further experimental measurements and, thereby, motivate additional experiments. Also considered is the question: why might chaotic dynamics occur in microbial food webs? Given the odd mathematical properties of chaotic dynamics, shouldn’t there be some biological advantage for such dynamics to be selected rather than classical dynamics?

The Becks et al. experiments

A conceptual model of the Becks et al. (2005) experiments is shown in Fig. 1. The food web was composed of a nutrient source, two bacteria that consumed nutrient (a rod, r, and a coccus, c), and a ciliate predator, p, that consumed both bacteria. The variable supporting the system was the food supply that was varied by changing the dilution rate (chemostat flow/chemostat volume), having units of inverse time). The four coupled dependent.

Fig. 1
figure 1

Diagram of the Becks et al. (2005) chemostat system. A nutrient solution flowing left to right was consumed by two microbes (rods and cocci), with rods being stronger competitors for nutrient than cocci. A ciliate predator fed on both microbes, but preferred rods over cocci

Full size image

variables were concentrations of nutrient (mg/cc) and each of the three microbe populations (cells/cc). For a fixed set of dilution rates, the three microbe concentrations were measured at a selected set of approximately daily time intervals. Each set of data constituted a time series (concentrations at discrete times), and deterministic chaos was identified by using a computerized version of the analytical procedure developed by Rosenstein et al. (1993) for calculating the largest Lyapunov exponent with the application of the TISEAN package (Hegger et al. 1999). Classical steady states were observed at D = 0.9/day and 0.75/day, chaotic dynamics were observed at D = 0.5/day, and periodic dynamics were observed at a slightly lower D = 0.45/day. (See Becks et al. (2005) for additional details.)

In addition to their main study, Becks et al. (2005) made separate measurements that provided some basis for a mathematical analysis of the experiments. These measurements included: (1) the average volume of the individual microbes, (2) the fact that the rods would out-compete the cocci for nutrient when no predators were present, and (3) that in one measurement, the predators preferred to consume the rods over the cocci at a ratio of 4:1. These supplementary measurements turned out to be vital for even surrogate model development.

Mathematical model development

To formulate mathematically the processes of the Becks et al. experiment, we first generalized the original three-equation chemostat model (nutrient plus two microbes) developed by Kot et al. (1992) to a four-equation model of four dependent variables: nutrient n(t), rods r(t), cocci c(t), and predators p(t) – with t being time. The units of the nutrient concentration are mg/cc, while those of the microbes are cell numbers/cc (cells/cc). In order to conserve mass in the resulting model, an average microbial mass had to be selected for each microbe type. This was done based on measured microbe volume averages published by Becks et al. (2005). We note that modeling based on the mean microbial volume is already an approximation, because behavior could vary during the maturing process. Using the measured microbial volume averages along with an assumed density of 1 g/cc, results in: r mass (mr) = 1.6 × 10–9 mg; c mass (mc) = 8.2 × 10–9 mg, and p mass (mp) = 3.2 × 10–6 mg. With these minor differences, the Kot et al. (1992) equations for a nutrient (mg/cc), a rod (cells/cc) and a predator (cells/cc) would be written as:

$$\begin{gathered} \frac{{{\text{d}}n}}{{{\text{d}}t}} = Dn_{0} - \frac{{\mu_{rn} }}{{Y_{rn} }}\left[ {\frac{{n(rm_{r} )}}{{K_{rn} + n}}} \right] - Dn \hfill \\ \frac{{{\text{d}}(rm_{r} )}}{{{\text{d}}t}} = \mu_{rn} \left[ {\frac{{n(rm_{r} )}}{{K_{rn} + n}}} \right] - \frac{{\mu_{pr} }}{{Y_{pr} }}\left[ {\frac{{(rm_{r} )(pm_{p} )}}{{K_{pr} + (rm_{r} )}}} \right] - D(rm_{r} ) \hfill \\ \frac{{{\text{d}}(pm_{p} )}}{{{\text{d}}t}} = \mu_{pr} \left[ {\frac{{(rm_{r} )(pm_{p} )}}{{K_{pr} + (rm_{r} )}}} \right] - D(pm_{p} ) \hfill \\ \end{gathered}$$
(1)

These three equations are identical to those used by Kot et al. (1992), with their S = n, H = rmr, and P = pmp. The “mx” terms are mean mass of the respective microbe, so r and p are dimensionless numbers of rods and predators per cc, the microbe “concentration” units recorded by Becks et al. (2005). All the other constant terms are various maximum specific growth rates “μxx”, half saturation constants “Kxx” and yield coefficients “Yxx” used to define the Monod kinetics. Model (Monod, 1949).

To extend Eqs. (1) to include one more nutrient-consuming microbe, a coccus, an equation similar to the middle one of Eqs. (1) is added, along with the analogous coupling terms, resulting in the following system of four ordinary differential equations.

$$\begin{gathered} \frac{{{\text{d}}n}}{{{\text{d}}t}} = Dn_{0} - \frac{{\mu_{rn} }}{{Y_{rn} }}\left[ {\frac{{n(rm_{r} )}}{{K_{rn} + n}}} \right] - \frac{{\mu_{cn} }}{{Y_{cn} }}\left[ {\frac{{n(cm_{c} )}}{{K_{cn} + n}}} \right] - Dn \hfill \\ \frac{{{\text{d}}(rm_{r} )}}{{{\text{d}}t}} = \mu_{rn} \left[ {\frac{{n(rm_{r} )}}{{K_{rn} + n}}} \right] - \left\{ {\frac{{\mu_{pr} }}{{Y_{pr} }}\left[ {\frac{{(rm_{r} )(pm_{p} )}}{{K_{pr} + (rm_{r} )}}} \right]} \right\} - D(rm_{r} ) \hfill \\ \frac{{{\text{d}}(cm_{c} )}}{{{\text{d}}t}} = \mu_{cn} \left[ {\frac{{n(cm_{c} )}}{{K_{cn} + n}}} \right] - \left\{ {\frac{{\mu_{pc} }}{{Y_{pc} }}\left[ {\frac{{(cm_{c} )(pm_{p} )}}{{K_{pc} + (cm_{c} )}}} \right]} \right\} - D(cm_{c} ) \hfill \\ \frac{{{\text{d}}(pm_{p} )}}{{{\text{d}}t}} = \mu_{pr} \left[ {\frac{{(rm_{r} )(pm_{p} )}}{{K_{pr} + (rm_{r} )}}} \right] + \mu_{pc} \left[ {\frac{{(cm_{c} )(pm_{p} )}}{{K_{pc} + (cm_{c} )}}} \right] - D(pm_{p} ) \hfill \\ \end{gathered}$$
(2)

The parameters involved are maximum specific growth rates (μxx), half saturation constants (Kxx) and yield coefficients (Yxx), a total of twelve. Equations (2) are direct generalizations of the Kot et al. (1992) model, and their mathematical validity was checked by showing that a mass balance was maintained, and when the Kot et al. (1992) initial conditions and parameter values were used, with one microbe forced to die out, output equivalent to Fig. 3 of Kot et al. (1992) was reproduced. (The braces {} in Eqs. 2 and 3 are used to identify terms used later in the mathematical development.) As written, however, Eqs. (2) do not include information from the supplemental experiments of Becks et al. (2005) or information on how a chemostat operates, such as potential nutrient recycling from dying biomass. Using Eqs. (2) alone, we were not able to produce chaotic dynamics unless we used a sinusoidal dilution rate as Kot et al. (1992) did.

In expanding Eqs. (2) to simulate the Becks et al. supplemental experiments, the specific death rate of predators and their biomass recycling to nutrient are likely to be important, because their average mass is about 1000 times greater than that of each feeding microbe. Moreover, populations of feeding microbes decrease mainly due to consumption by predators, while nothing consumes the dying predators. When predators die, their bodies simply break down with remains consumed or flushed out of the well-mixed chemostat. Because of these considerations, natural death and biomass recycling of the feeding microbes were assumed to be negligible relative to predators, and this was supported also by numerical experiments. The specific death rate for predators and the nutrient recycling terms will be identified in the final system of equations developed below.

As observed in the Becks et al. (2005) supplemental experiments, in the absence of predators r would out-compete c for nutrient, and at an identical population of r and c (4 × 106 cells/cc), p consumed r cells over c cells in the ratio of 4:1. We decided to incorporate the additional information into Eqs. (2) as follows: (1) When r and c are low, p chooses them on an equal basis even though in general r is preferred over c (Starving organisms are not choosy?), (2) At high r and c, as observed in the experiments, p chooses r four times more than c, and (3) in competition for nutrient with no predators, the cocci die out first. In their more recent paper, Becks and Arndt (2013) discussed the possibility that a p preference change with r and c concentrations could be an important but unmeasured process in the Becks et al. (2005) study.

A simple way to impose condition (3), is to set the value of \(\mu_{rn}\), the maximum specific growth rate of r on n, equal to \(k\mu_{cn}\), with \(k > 1\). Then, for \(k\) being sufficiently large, r will always outcompete c at low cell numbers. Incorporating conditions 1 and 2 leads to a new modification of the Monod kinetics model; it results in a kinetic expression that incorporates a generalization of the supplemental results of the Becks et al. (2005) experiments (p preference ratio for r vs. c varying from 1:1 to 4:1 rather than staying constant at 4:1. It also decreases the number of free parameters in the overall model, which we view as positive.

Based on the Becks et al. (2005, Fig. 1) data, the r and c concentrations are ranging from about 1 × 105 cells/cc to 2 × 106 cells/cc. On a cell numbers basis -see braced terms in the 2nd and 3rd members of Eqs. (2)- the respective consumption rates of p on r and p on c are expressed as follows:

$$m_{r} \frac{{{\text{d}}r_{p} }}{{{\text{d}}t}} = \frac{{\mu_{pr} rm_{r} pm_{p} }}{{Y_{pr} (K_{pr} + m_{r} r)}}\;{\text{and}}\;m_{c} \frac{{{\text{d}}c_{p} }}{{{\text{d}}t}} = \frac{{\mu_{pc} cm_{c} pm_{p} }}{{Y_{pc} (K_{pc} + m_{c} c)}}$$
(3)

Noting that the minimum values of r and c are approximately 1 × 105 cells/cc, we set drp/dt = dcp/dt (preference ratio of 1:1), and obtain from (3):

$$\frac{{\mu_{pr} (1 {\times} 10^{5} )pm_{p} }}{{Y_{pr} (K_{pr} + 1 {\times} 10^{5} m_{r} )}} \, = \frac{{\mu_{pc} (1 {\times} 10^{5} )pm_{p} }}{{Y_{pc} (K_{pc} + 1 {\times} 10^{5} m_{c} )}}$$
(4)

Assuming that for maximum values of r and c (approximately r = c = 2 × 106 cells/cc) drp/dt = 4dcp/dt (preference ratio of 4:1), we obtain from (3):

$$\frac{{\mu_{pr} \left( {2 \times 10^{6} } \right)pm_{p} }}{{Y_{pr} \left( {K_{pr} + 2 \times 10^{6} m_{r} } \right)}} = \frac{{4\mu_{pc} \left( {2 \times 10^{6} } \right)pm_{p} }}{{Y_{pc} \left( {K_{pc} + 2 \times 10^{6} m_{c} } \right)}}$$
(5)

We will modify the numerator and denominator of drp/dt, given in Eqs. (3), in order to achieve the equalities specified in Eqs. (4) and (5), with a linear variation in between. The maximum specific growth rate of p on r will become μpr(m1r + i1), and the denominator will become Kpr + m2r + mrr, with m1, i1 and m2 constants. Substituting this in drp/dt yields

$$\frac{{{\text{d}}r_{p} }}{{{\text{d}}t}} = \frac{{\mu_{pr} \left( {m_{1} r + i_{1} } \right)r\left( {pm_{p} } \right)}}{{Y_{pr} (K_{pr} + m_{2} r + m_{r} r)}}$$
(6)

Thus, we are changing \(\mu_{pr}\) and \(K_{pr}\) from constants to linear functions of r. To keep the overall units consistent, m1 has the units of (cells/cc)−1, i1 is dimensionless, and m2 has the units of mg. Equation (6) is a new semi-empirical relationship expressing a modified form of the Monod kinetics equation, motivated by the Becks et al. experiment, and is given to consider a p preference change for r relative to c. Then to satisfy Eqs. (4) and (5), the following conditions must be met by Eq. (6) (note the terms that cancel):

$$\begin{gathered} \frac{{\mu_{pr} (10^{5} m_{1} + i_{1} )}}{{Y_{pr} (K_{pr} + 10^{5} m_{2} + 10^{5} m_{r} )}} \, = \frac{{\mu_{pc} }}{{Y_{pc} (K_{pc} + 10^{5} m_{c} )}} \, ,\,{\text{and}} \hfill \\ \frac{{\mu_{pr} (2 \times 10^{6} m_{1} + i_{1} )}}{{Y_{pr} (K_{pr} + 2 \times 10^{6} m_{2} + 2 \times 10^{6} m_{r} )}} \, = \frac{{4\mu_{pc} }}{{Y_{pc} (K_{pc} + 2 \times 10^{6} m_{c} )}} \, \hfill \\ \end{gathered}$$
(7)

To obtain the equalities in Eqs. (7), we will first modify the numerators followed by the denominators to maintain overall equality. This will also require introducing the equalities μpr = μpc, Ypr = Ypc and Kpr = Kpc. (Little is known concerning the true values of these parameters and equating them reduces the number of unknown parameters in the overall model.) Accordingly, the introduced parameters m1 and i1 must satisfy the following conditions: 105m1 + i1 = 1, and 2 × 106m1 + i1 = 4, also setting μpr = μpc. These conditions yield m1 = 1.579 × 10–6 (cells/cc)−1 and i1 = 0.8421. The corresponding conditions on m2 are as follows:

$$\begin{gathered} Y_{pr} \left( {K_{pr} + 10^{5} m_{2} + 10^{5} m_{r} } \right) \, = \, Y_{pc} \left( {K_{pc} + 10^{5} m_{c} } \right),\;{\text{and}}\; \hfill \\ Y_{pr} \left( {K_{pr} + 2 \times 10^{6} m_{2} + 2 \times 10^{6} m_{r} } \right) \, = Y_{pc} \left( {K_{pc} + 2 \times 10^{6} m_{c} } \right) \hfill \\ \end{gathered}$$
(8)

If we set Ypr = Ypc and Kpr = Kpc, Eqs. (8) become,

$$\begin{gathered} 10^{5} m_{2} + 10^{5} m_{r} \, = \, 10^{5} m_{c} ,\;{\text{and}} \hfill \\ 2 \times 10^{6} m_{2} + 2 \times 10^{6} m_{r} \, = 2 \times 10^{6} m_{c} \hfill \\ \end{gathered}$$
(9)

It can be seen from both relationships that \(m_{2} = m_{c} - m_{r} = 8.2 \times 10^{ - 9} - 1.6 \times 10^{ - 9} = 6.6 \times 10^{ - 9}\) mg. Thus, after dividing through by the microbial masses, Eqs. (2) being adapted to the Becks et al. (2005) supplemental experiments are written as follows:

$$\begin{gathered} \frac{{{\text{d}}n}}{{{\text{d}}t}} = Dn_{0} - \frac{{\{ k\mu_{cn} \} }}{{Y_{rn} }}\left[ {\frac{{n(rm_{r} )}}{{K_{rn} + n}}} \right] - \frac{{\mu_{cn} }}{{Y_{cn} }}\left[ {\frac{{n(cm_{c} )}}{{K_{cn} + n}}} \right] - Dn \, + \, pm_{p} \delta_{p} ({\text{EF}}) \hfill \\ \frac{{{\text{d}}r}}{{{\text{d}}t}} = \{ k\mu_{cn} \} \left[ {\frac{n(r)}{{K_{rn} + n}}} \right] - \frac{{\{ \mu_{pc} \} ((1.58 \times 10^{ - 6} )r + 0.842)}}{{\{ Y_{pc} \} }}\left[ {\frac{{r(pm_{p} )}}{{\{ K_{pc} \} + (6.6 \times 10^{ - 9} )r + (rm_{r} )}}} \right] - Dr \hfill \\ \frac{{{\text{d}}c}}{{{\text{d}}t}} = \mu_{cn} \left[ {\frac{n(c)}{{K_{cn} + n}}} \right] - \frac{{\mu_{pc} }}{{Y_{pc} }}\left[ {\frac{{c(pm_{p} )}}{{K_{pc} + (cm_{c} )}}} \right] - Dc \hfill \\ \frac{{{\text{d}}p}}{{{\text{d}}t}} = \{ \mu_{pc} \} ((1.579 \times 10^{ - 6} )r + 0.842)\left[ {\frac{{p(rm_{r} )}}{{\{ K_{pc} \} + (6.6 \times 10^{ - 9} )r + (rm_{r} )}}} \right] + \mu_{pc} \left[ {\frac{{p(cm_{c} )}}{{K_{pc} + (cm_{c} )}}} \right] - Dp - p\delta_{p} \, \hfill \\ \end{gathered}$$
(10)

In these equations, δp is the specific death rate for predators, which is recycled to nutrient at the efficiency “EF ≤ 1” in the 1st equation of set (10). For clarity, the modified rod kinetic parameters are placed within braces {}. Further clarification of the mathematical manipulations resulting in Eqs. (10) may be found in the Supplemental Material.

In the final model, the predators have been made to prefer rods over cocci with an increasing rod population, so the rods are disadvantaged and would tend to die out. Based on numerical experiments, this is prevented by setting the k factor in the first of Eqs. (10) to 1.5. As presented in Molz et al. (2019), dimensionless forms of Eqs. (10) were developed in order to aid the numerical simulations of the Becks et al. experiments and the understanding of parameter interactions. These dimensionless equations may also serve as a basis for further study of microbial experiments.

Summary of modeling results

Equations (10) were solved using MATLAB software (ODE45) subject to the approximate initial conditions of n(t = 0) = 0.03 mg/cc, r(t = 0) = 4.2 × 106 cells/cc, c(t = 0) = 1. × 106 cells/cc and p(t = 0) = 3000 cells/cc. A set of parameters that produced chaotic dynamics at the measured dilution rate of 0.5/d is listed in Table 1. Based on parameter values given in Kot et al. (1992) and value ranges given in Kravchenko et al. (2004), Table 1 values appear reasonable. Other than the dilution rate and mean microbe masses, we arrived at the remaining parameter values by “trial-and-error” numerical simulation experiments. Faybishenko et al. (2018) presented additional mathematical details consistent with chaotic behavior. The biological details of the Becks et al. (2005) experiments were not considered.

Table 1 Parameter values used in Eqs. (10) that yield chaotic dynamics with the dilution rate utilized in Becks et al. (2005)
Full size table

The results of simulations, based on the parameter values listed in Table 1, are shown in Figs. 2 and 3. Listed in Table 2 are the ranges of maximum Lyapunov exponents given in Becks et al. (2005), which are reasonably close to those calculated from our simulated time series.

Fig. 2
figure 2

System-space (phase-space) plot of r(t), c(t) and p(t) for the parameters listed in Table 1 for D = 0.5/d (0.0208/h). A strange phase-space attractor is evident

Full size image
Fig. 3
figure 3

Irregular, non-periodic, dynamics of n(t), r(t), c(t) and p(t) associated with the system-space plot in Fig. 2. By counting apparent concentration peaks, the mean frequency of the simulated data appears lower than that observed in Becks et al. (2005), but still within 50%

Full size image
Table 2 Lyapunov exponents used to identify the presence of deterministic chaos based on the time series of concentrations shown in Fig. 3 (Calculations were conducted using the R package “fractal” version: 2.0–1)
Full size table

As in the Becks et al. experiments, after obtaining chaotic dynamics at a dilution rate D = 0.5/day (0.0208/hr), we decreased D to 0.45/day (0.01875/hr) with other parameters as listed in Table 1. It is shown in Fig. 4 that the output starts in an irregular manner, and then becomes steady state with the cocci dying out. In the experiments, the results were periodic with all microbes surviving.

Fig. 4
figure 4

Resulting nutrient and microbe dynamics with D = 0.45/d (0.01875/h.), as was done in the Becks et al. (2005) experiments

Full size image

For a D value of 0.9/d (0.0375/hr), both simulations and experiments produced a classical steady state with one microbe dying out (See Fig. 5 for initial model results). With this more classical result, it should be relatively easy to “fit” the model to the Becks et al. results. To illustrate this, we adjusted three parameters listed in Table 1. The maximum specific growth rate for the cocci was changed from 0.1248/h to.

Fig. 5
figure 5

Resulting nutrient and microbe kinetics with D = 0.9/d (0.0375/h) using the parameter values in Table 1

Full size image

0.1/hr., and the mean masses of the rods and predators were changed respectively from 1.6 × 10–9 mg and 3.2 × 10–6 mg to 2.2 × 10–9 mg and 2.23 × 10–6 mg – not exceedingly large changes. The resulting simulations shown in Fig. 6 are essentially identical to the Becks et al. (2005) experimental results for D = 0.9/d (0.0375/h). However, with the same parameter values the chaotic dynamics could not be reproduced at a D of 0.5/day (0.0208/hr).

Fig. 6
figure 6

Simulated nutrient and microbe dynamics that match the Becks et al. (2005, Fig. 1) results for D = 0.9/d (0.0375/h)

Full size image

We also solved Eqs. (10) with D = 0.75/d (0.03125/h.), which produced a classical steady state with the rods dying out, as shown in Fig. 7; Becks et al. achieved a steady state also, but with all microbes surviving.

Fig. 7
figure 7

Resulting nutrient and microbe dynamics with D = 0.75/d (0.03125/h). Both the model and experiments produced classical steady states, but all microbes survived in the experiments while the rods died out in the model

Full size image

Discussion and conclusions

The present study yields several interesting parallels between the experimental and modeling results. The nature and magnitudes of the model and experimental outputs were similar. For the chaotic case with a dilution rate of 0.5/d (0.0208/h) in both the experiment and the mathematical model, the Lyapunov exponents were in reasonable agreement (Table 2) with chaotic dynamics resulting. For a dilution rate of 0.9/d (0.0375/h) in both the modeling and experiments and a moderate tuning of three unmeasured parameter values (+ 38% for mr, − 43% for mp and -2% for μcn), the model outputs were identical to the experimental results. So, we have a well-defined experiment and a surrogate model that produced both chaotic and classical time series output patterns depending on the dilution rate. This lends further support to the conclusion of Becks et al. (2005) that chaotic dynamics in combination with classical dynamics occurred in their experiments.

From our viewpoint, the key methodology for obtaining internal chaotic dynamics was to cause the preference of p for r vs. c to vary with the r concentration, and to make r more competitive for nutrient than c, as well as to recycle some dying p biomass. This was consistent with the experimental details, and it was vital to the model development that Becks et al. (1995) performed the supplemental measurements related to mean microbial volume, predator food preferences and relative competitiveness of the rods and cocci for nutrients. In this area, our surrogate model analysis yields a prediction: the occurrence of chaotic dynamics requires a predator preference for r versus c to increase significantly with increases in r and c populations, and this could be checked experimentally by extending the Becks et al. (2005) supplemental experiment on predator preference change. Becks et al. (2005) measured this preference at single r and c values and got a ratio of 4:1, but it could be done at a variety of values to see if our linear relationship is valid, or if there is another one affecting the overall dynamics. It would illustrate how potential couplings can be complex beyond any of the classical kinetic expressions such as the Monod expression. The applicable experimental technique was developed by Christofferson et al. (1997.)

What the results are telling us is that future experiments should be designed to have a mix of modeling and experimental testing so that the parameter values or functions for the actual experimental kinetics are better defined from the outset. Such precision is required mathematically in nonlinear, coupled models producing both classical and chaotic dynamics, as in the present analysis. This is a significant experimental challenge, but our results provide clear motivation in that direction. The question is: can one measure and specify precisely the kinetics of interacting components of a microbial system so that a corresponding mathematical model reproduces the chaotic and classical dynamics observed in the experiment? In considering the Becks et al. (2005) experiment, one would presume that the microbes involved respond to their changing environment in a nearly continuous manner directed by their individual genetic codes. In different environments, such responses would be expected to have extensive flexibility and variability derived from millions of years of evolution. How does one specify such kinetics before performing the various experiments? Becks et al. (2005) moved in that direction by performing short term supplemental experiments along with their main experiments, which is probably one reason why Fussmann (2007) recommended those experiments for theoretical study.

The mathematical model represented by Eqs. (10) may be generalized further in several ways. Additional microbes may be added, and the interaction kinetics may be modified further to represent what is occurring in a given experiment. Such modification is not difficult in a formal mathematical sense; the challenge is to base such modifications on experimental evidence and motivation. Another interesting idea would be to perform experiments involving a spatial variable. This leads one to consider biofilms (Davey and O’Toole 2000; Picioreanu et al. 2005). It is observed that nicely structured and controlled biofilms can be grown on electrodes, and electrical measurements can be made very precisely (Yoho et al. 2014).

Returning to the question posed in the introduction: If chaotic dynamics can occur in place of classical dynamics in microbial food webs, what might be a biological advantage? A possible answer is suggested by a throretical analysis of the chaotic double pendulum (Mori et al. 1991). These scientists found that when the pendulum entered its chaotic mode, the system energy dissipation rate due to friction fell into a local minimum, and this was unexpected. An analogous measurement could be made on chemical reactions displaying chaotic dynamics, such as the Belousov–Zhabotinskii reaction system (Hudson and Mankin 1981). If a system energy dissipation minimum (increased system efficiency) is associated with chaotic biological systems, then such systems would be favored by natural selection. As presented by Zotin (2012), the minimum energy dissipation principle developed by Prigogine applies to thermodynamic systems at or near a steady state. It is conceivable that a deterministic chaotic attractor (e.g., Figure 2) may be viewed as a system generalization of a classical steady state as suggested by Molz and Faybishenko (2013). So, the minimum energy dissipation principle applied potentially to a strange attractor should be tested experimentally. Any potential energy dissipation principles are not included in our present surrogate model. How this might be done could become a subject of future research.