Abstract
A theoretical model developed by Stone describing a three-level trophic system in the Ocean is analysed. The system consists of two distinct predator-prey networks, linked by competition for nutrients at the lowest level. There is also an interaction at the level of the two preys, in the sense that the presence of one is advantageous to the other when nutrients are low. It is shown that spontaneous oscillations in population numbers are possible, and that they result from a Hopf bifurcation. The limit cycles are analysed using Floquet theory and are found to change from stable to unstable as a solution branch is traversed.
1. Introduction
In a recent paper, Stone [1] presented a mathematical model of a three-level trophic food web involving ocean-dwelling microorganisms. It was his aim to explain the paradoxical nature of Phytoplankton, which, when stressed by nutrient limitation, stimulates the growth of Bacteria, a competitor for the nutrients, by releasing carbon which is then taken up by the Bacteria. Stone [1] argued that this behaviour could be explained in terms of indirect effects. This means that the relationship between two compartments of the model should not be viewed in isolation but rather as part of the whole. His conclusion, drawn from a simple mathematical examination, was that due to the success rate of the grazing Protozoa, which prey on the Bacteria, the Phytoplankton are actually simulating the competitor of a competitor in order to survive. Recently, Hadley and Forbes [2] examined Stones’ [1] model for the case when nutrients are unlimited. They found, using a dynamical systems approach, that there are no limit cycles arising from a Hopf bifurcation in this case. There are, however, oscillatory solutions arising from centre type behaviour. Nevertheless, these were anticipated to be structurally unstable, in the sense that refinements to the model, such as the inclusion of nutrient-dependent growth rates, would replace these oscillations with some other behaviours. This may be an indicator as to why the Phytoplankton do not employ this survival strategy under every circumstance.
Much research has been performed on the stability of interacting populations in biological models (e.g., see May [3], Murdoch et al. [4]), focusing on how populations are drawn to attractors such as limit cycles or steady-states, or if and when extinction occurs. Shertzer et al. [5] found that a simple mechanistic model of a four-species predator-prey interaction agreed with data collected in a chemostat experiment. The model predicted limit cycle, steady-states, and extinction behaviour present in the chemostat although the results differed on the prediction of period and phase. The model was improved to include the ability of prey to evolve defence mechanisms and was found to be accurate in all respects. Hutchison [6] speculated that Plankton communities cannot come to equilibrium but continue to develop to oscillatory solutions or chaos owing to all the external effects they are subject to as their environment changes. Scheffer et al. [7] found that, even in homogeneous and constant environments, various competition and predation models suggest that plankton will never settle to equilibrium. This was also found in experimentation and that chaos quite often resulted, even with low dimensional systems. This has the outcome of making long-term predictions about such systems impossible. However, Verschoor et al. [8] subsequently showed experimentally that bi-and tritrophic food chains with induced defences approached a stable equilibrium without any oscillatory tendency, while those without defences in the algae showed high-amplitude population. For a tritrophic model, van der Stap et al. [9] showed that population stability of phytoplankton occurred when the phytoplankton had a defence mechanism that affected the uptake interaction of the Zooplankton. From these studies, it appears that qualitative behaviour, such as population stability oscillation and chaos, is highly dependent on the number of species present and the conditions under which the results were obtained.
There has been much work done on the effect of interaction functions in dynamical systems. It has been shown by Gross et al. [10] that Holling II, or Michaelis-Menten, type interaction functions can either destabilize or stabilize steady-states, dependent on the form of the interaction function used. A change in stability is achieved by enriching one of the populations. This was found to be an underlying idea in the “paradox of enrichment” (Edwards and Brindley [11]). The “paradox of enrichment’’ was originally attributed to Rosenzweig [12] and further refined by May [13] and Gilpin and Rosenzweig [14]. In this paper we use a Michaelis-Menten term for the growth rates, which are consequently nutrient-dependant. This will be discussed in more detail in Section 2, and a more comprehensive discussion of this idea is given by Murray [15].
Other dynamical structures such as quasiperiodicity have also been found in higher dimension food webs similar to the one being examined here. Ruan [16] examined a three-level model involving Zooplankton, Phytoplankton, and nutrients in limited supply. Modelling nutrients with both instantaneous and delayed recycling, they found that the equilibrium point loses stability when a critical value is reached in nutrient levels and passes via a Hopf-bifurcation into a limit cycle. On the other hand, Wang et al. [17] found in a three-level food web with nonlinear nutrient dependence that there was limit cycle behaviour, quasiperiodicity, and chaos.
In the present paper, we use methods from Dynamical Systems theory to analyse the three-level trophic food web described by Stone [1] (for further discussion on Dynamical Systems theory, see Murray [15] and Edelstein-Keshet [18]). This model displays the interaction between Phytoplankton, Bacteria, Protozoa, Zooplankton, and Nutrients. In particular we are looking for stability of the steady-state populations and possible self-sustained oscillations. The model has been examined previously by Hadley and Forbes [2] in the case when nutrients did not vary. In the present paper, we investigate the possibility that the introduction of variable nutrient concentration in the system will generate instability in the steady-state, causing self-sustained oscillations to occur in the solutions. We find a Hopf bifurcation present in this nonlinear model, leading to oscillatory limit-cycle behaviour in the system. We examine the stability of the limit cycles using Floquet theory.
The model is presented in Section 2 and for convenience, scaled (nondimensional) populations and rates are introduced. In Section 3, we provide an analysis of the model including a determination of the conditions under which a limit-cycle oscillation can arise through Hopf bifurcation. A branch of limit cycles is located numerically in Section 4, and Floquet theory is used to study the change in its stability as the breeding rate of bacteria is altered. The paper concludes in Section 5 with a discussion of the results.
2. The Mathematical Model
The trophic web considered now is that illustrated in Figure 1, from Stones’ [1] original paper, where the direction of the interaction is given by the arrows. There are five interacting components, namely the Bacteria , Phytoplankton , Protozoa , Zooplankton , and Nutrients as shown in Figure 1. The interactions in the diagram are described by the system of five differential equations
It should be noted here that Stone’s model (Figure 1) assumes there is no cross feeding of Protozoa on Phytoplankton or Zooplankton on Bacteria (Stone [1]). This is reflected in (2.1). Here the growth rate of the Bacteria and the reproductive rate of the Phytoplankton, and , respectively, depend on the nutrient concentration . We are using a Michaelis-Menten uptake term for nutrients by Bacteria and Phytoplankton, so that the rates are governed by
These terms reflect the saturation effects seen in the nutrient uptake process, whereby an infinite supply of nutrients will not result in unlimited growth. The terms , , are constants that determine the rate of uptake of the nutrients. In their paper Gross et al. [10] showed a situation where, for Holling II type interaction functions like (2.2), enrichment (increase in a population) destabilized the steady-state solutions of a general food chain similar to the one being considered.
The quantity is the interaction rate between Phytoplankton and Bacteria. The constants and are the mortality rates of the Protozoa and Zooplankton, respectively. The final quantity represents the mass of nutrients contained in each organism. All these quantities are positive.
In their paper, Hadley and Forbes [2] assumed that N did not change and therefore had some fixed value say. For consistency with their results and using (2.2), we obtain the initial value for the reproductive rates
We incorporate (2.3) into the nutrient-dependant reproductive rates by forming the ratio
The nutrient-dependant reproductive rates and are now expressed in the equivalent form
In this way, if , these rates would be the same as in the nutrient-independent case examined by Hadley and Forbes [2]. Here and are constants.
The system (2.1) is now recast in terms of dimensionless variables. In this case, the four populations (B, P, Z, R) are scaled with respect to the quantity , which is a naturally occurring measure of population size. The nutrient concentration is scaled with respect to the quantity which is arbitrary, and so this choice has no implications for the dynamics of the system. Time t is made dimensionless using the quantity , which is a time scale linked roughly to the life-cycle of the phytoplankton. In these nondimensional variables, (2.1) become
The constant is such that .
There are five nondimensional parameter groupings in the system (2.6). These are
The first parameter represents the interaction rate between the Phytoplankton and the Zooplankton. The second quantity represents the growth rate of the Bacteria. The third and fourth parameters and are the mortality rates of the Protozoa and the Zooplankton, respectively. The last parameter represents the interaction rate between Phytoplankton and Bacteria. All these quantities are positive.
3. Analysis of the Model
3.1. Steady-State Populations
We now look for the steady-state solutions which satisfy
We choose , where now represents some arbitrary dimensionless initial concentration of nutrient. This yields the four equilibrium points
where In (3.3) the constants and are the nutrient uptake rates for the bacteria and phytoplankton, respectively. The first steady-state in (3.2) represents the case where only the nutrients remain, and all the species die out. The second and third steady-states are where two species, a predator and prey coupling, survive. The fourth steady-state, which is of most interest to our analysis, is where all the species survive. We will now look to ascertain the stability of these states. It should be mentioned that two more steady-states were found; however they contained negative values for some of the populations and as such were unrealisable in an actual situation.
3.2. Stability of the Steady-States
When the time-dependent populations are close to any of the four steady-states in (3.2), the small-amplitude behaviour may be determined by linearization, in the form
The constant ε represents how close the system is to one of its steady-states. We determine the linearized system near an equilibrium point by substituting (3.4) into the governing equations (2.6) and retaining terms at the first order in ε. This results in the linear matrix system
Here,
are defined for convenience.
We are only interested here in the last steady-state in (3.2), where all the populations survive. These values are substituted into (3.5) to produce the Jacobian matrix
In this expression (3.7), it has proved convenient to define the intermediate quantities
The eigenvalues for the linearized system (3.5), with coefficient matrix (3.7), are found from the characteristic equation
in which it is convenient to define
In order for there to be a Hopf bifurcation present, we need there to be at least one pair of complex conjugate eigenvalues whose real part vanishes. For this to happen, (3.9) must be in the form
We now expand (3.11) into the form
We perform a comparison between the corresponding coefficients of powers of given in (3.9) and (3.12) so that
and then solve these sets of criteria to find and the coefficients and in terms of the quantities in (3.10). This enables us to find the following conditions necessary for a Hopf bifurcation to be present in the system (2.6);
The inequalities (i), (iii), and (iv) are a consequence of the need to ensure that the quantity in (3.13) remains positive, to satisfy the necessary conditions for the Hopf bifurcation. The condition (ii) results from solving (3.13) for the quantities . We solved the above conditions numerically to find a set of parameter values consistent with a Hopf bifurcation. The values for the parameters , , were derived from Stones’ [1] paper. We found that by varying the values of the parameters and from those derived in Stone’s original paper, we were able to find a set of parameter values for which a Hopf bifurcation is possible in the system. The parameters and were kept within the realms of physical possibility; they were in fact less than the values determined from the rate constants used in Stone’s [1] paper. The Hopf curve is shown below in Figure 2, and it indicates the location in the parameter-space of the points for which conditions (3.14) are satisfied, and therefore gives the parameter values at which nonlinear oscillatory limit cycles might be expected to be born.
Both the supercritical and subcritical limit cycles emerge to the right in Figure 2, as stable and unstable structures, respectively. Self sustained oscillations born from a Hopf bifurcation occur for values of on or inside the curve for a given . The Hopf bifurcation first appears at . Each successive value of has two Hopf points, a supercritical (continuous line) and subcritical (dashed line) Hopf point. A stable limit cycle emerges at a supercritical Hopf bifurcation, but the line of subcritical bifurcations at the top of the diagram gives rise to unstable oscillatory behaviour. The supercritical and subcritical branches converge again at , after which the Hopf criteria (3.14) are no longer satisfied for the current parameter values.
Shertzer et al. [5] showed that predator-prey models of systems, where the prey has developed defences to attack, have the best fit with experimentation and that limit cycles are born from the terms modeling the defence mechanism. Our analysis shows that self-sustained oscillations in the unforced system are possible. The oscillations take the form of limit cycles arising out of a Hopf bifurcation (periodic oscillations of species over time in a closed system) which will be investigated further in Section 4.
4. Numerical Results
We used MATLAB to solve the nonlinear system at the parameter values discussed in the previous section. Once we were satisfied that apparently oscillatory solutions were present then, using a shooting algorithm based on Newton’s method, we were able to find if the resultant solutions were in fact periodic limit cycles. The method used will now be described briefly here.
We rescaled (2.6) by introducing a new time variable
where is the period of oscillation and is as yet unknown. We made a guess for the initial conditions and . Using MATLAB, we integrated the set of rescaled equations to find the values after one complete period. If the solution is truly periodic, then the initial conditions should match these values. We used the equilibrium point as the initial value for B(0) as we know a periodic solution set will include this point. We estimated the values of the other variables when using the numerical MATLAB code, and used these as initial values, including a guess for the period which was also estimated from the MATLAB results.
Newton’s method was used to adjust the estimates of , , , , and so that the residual quantities , , , , and are made arbitrarily small. We used a damped version of Newton’s method whereby if the norm of the vector of residuals is not reduced then the length of the correction step used in the method was halved. As limit cycles become more unstable, the method finds it harder to converge to a solution. Integrating the system of equations backwards in time was found to be necessary for highly unstable solutions, since this procedure converts an unstable orbit into a stable one (in negative time).
After we have calculated a limit cycle, we then perform a linear perturbation to the system (2.6) once they have been rescaled with respect to . The resultant coefficient matrix is -periodic in and so we can use Floquet theory to ascertain the stability of the limit cycle. This method is described in Forbes [19]. In brief we found the eigenvalues of a monodromy matrix, which give a measure of how close the perturbation is to the limit cycle. If for all , then the limit cycle is stable; however, if any one of the , then the limit cycle is unstable. The result of our numerical analysis is now considered. We started with the initial set of parameter values
The locations of the Hopf points were obtained from Figure 2.
Figure 2 indicates that as increases from 1.092, the supercritical and subcritical branches of Hopf values diverge along the axis until they reach a minimum and maximum, respectively, at about . These two branches then begin to converge again until they meet at = 10.82, at which point the Hopf criteria (3.14) fail. As previously mentioned, the oscillatory solutions born from these two Hopf points give rise to nonlinear self-sustained oscillations for values of inside the curve in Figure 2, for each value of . The behaviour of these limit cycles at a fixed value of is now studied in more detail, by varying beyond the supercritical point and up to the subcritical value, as determined from Figure 2. The results of these numerical calculations are presented in Figure 3.
Figure 3 illustrates the amplitude of the limit cycle oscillations for the Bacteria, as the reproduction rate is varied. At , there is a supercritical Hopf bifurcation formed (this stable structure travels to the right as increases). This is precisely the value predicted in Figure 2 for , and a stable limit cycle of very small amplitude emerges at this point. We also checked this value by finding the eigenvalues of the Jacobian of the linear perturbation model described in (3.5), where at this value of there exists a purely complex conjugate pair of eigenvalues. The solid line in Figure 3 represents the stable portion of the nonlinear oscillatory solution branch. To follow this branch accurately, we reduced the error tolerance (MATLAB) used for the integration solver to . This was done so as to improve numerical accuracy.
The dashed lines in Figure 3 are the unstable portion of the nonlinear oscillatory solution branch. These limit cycles were found by integrating backward in time in the shooting method algorithm. As was mentioned previously, this was necessary as the limit cycles in this region are highly unstable.
As can be seen from Figure 3, there is a second subcritical Hopf bifurcation at at which an unstable limit cycle oscillation emerges from the steady-state and travels to the right as increases. These unstable solutions are indicated with a dashed line. Again, this value is in agreement with the predictions of our linear perturbation model (3.5) for , as illustrated in Figure 2. Another feature that can be seen from Figure 3 is the fold bifurcation occurring at . This is where the stable and the unstable limit cycle branches, born from the two Hopf bifurcation points, join to create a single unified branch of limit cycles. The Floquet multipliers have been computed for this nonlinear problem, and their behaviour has been monitored carefully along the branch shown in Figure 3. At the fold point, one of these multipliers takes the value one, and then its absolute value increases beyond one as the solutions indicated with a dashed line are traversed, confirming that these limit cycles are genuinely unstable.
The period of the limit cycles in Figure 3 has also been computed, and changes as varies. We found that for the supercritical Hopf point, , the nondimensional period is approximately 7.2236. This equates to about 14.5 days in dimensional terms. Similarly, for the subcritical Hopf point at , the nondimensional period was determined to be 6.7416 or 13.5 days in dimensional terms. At the fold bifurcation, the nondimensional period is 7.1079, which corresponds to 14.2 days.
We now look at two distinct solutions to the system occurring at the same value of . This serves to highlight the nonlinear dynamics of the complex system studied here.
Figure 4 shows that at the value , there are two independent solutions to the system forming two distinct limit cycles, at the same values of the physical parameters. For Bacteria at this value, the stable solution is shown with a solid line. As can be seen from the graph, the amplitude of the stable orbit is much greater than the unstable counterpart, consistently with Figure 3. That two distinct limit cycles can be formed for the same parameter value shows the complexity of the system. We examined the stability of the two limit cycle solutions produced at using Floquet theory. The results of this analysis are illustrated in Figures 5 and 6.
Figure 5 shows the five eigenvalues (red dots) of the monodromy matrix formed when we solve the linear perturbation to the limit cycle. These are shown against the unit circle in the complex eigenvalue plane since this represents the border between stability and instability in Floquet theory. When one of the eigenvalues crosses the unit circle, the solution becomes unstable. Now in Figure 5 we see that one of the eigenvalues is equal to 1. It is known that, for limit cycles, one Floquet multiplier must be equal to one. Physically, this corresponds to the fact that a perturbation tangent to the limit cycle remains on it, representing a neutrally stable event (Seydel [20]). This requirement, in fact, represents a very sensitive test of the numerical accuracy of our method. All the other eigenvalues lie inside the unit circle and therefore the solution is stable. If we then consider Figure 6, we see that the eigenvalues (indicated by red dots) have one pair of complex conjugates with real components that are significantly greater than one. This therefore corresponds to an unstable limit cycle.
The instability of this solution means that as a result of a small perturbation , the resultant solution will diverge from the limit cycle. In Figures 7 and 8, we can see this pronounced change in stability.
Figures 7 and 8 show a perturbation to both the stable and unstable limit cycles formed at . We achieved this perturbation by changing the initial conditions. In Figure 7, we firstly subtracted 0.2 from the initial value of the first species (Bacteria) and subtracted 0.3 from the initial value of the second species (Phytoplankton) and resolved the equations. In the second instance, we added 0.2 to the initial value of the Bacteria and subtracted 0.3 from the initial value of the Phytoplankton. This is equivalent to making a perturbation outside then inside the limit cycle. In Figure 8 we applied the same process by firstly, adding 0.02 to the initial values of each of the four species and the nutrients and resolving the equations, and in the second instance by subtracting 0.02. Each figure shows the limit cycle run for a suitable time as well as the solutions of the perturbations to the limit cycles. From the stable case seen in Figure 7, we see that the perturbation solutions converge to the limit cycle (shown in red). This is as expected for a stable orbit. In Figure 8 the solutions do not converge but instead oscillate and move away from the limit cycle.
5. Discussion
This paper presented an investigation of the solutions to the system proposed by Stone [1]. The system predicts four steady-states (and two more with unphysically realizable negative values for the populations). Of these equilibria, only one predicts long-term survival for all four species. It is around this equilibrium point that we centred our study.
Hadley and Forbes [2] analyzed this same system, although in the special case in which nutrient concentration was not allowed to vary. They showed that the system is degenerate in that case, meaning the unforced equations gave rise to equilibrium points that are centres (excluding the point where all populations vanish). It was suggested that when nutrients are allowed to vary, the degeneracy might be removed and Hopf bifurcations may occur instead of centre behaviour. This has been confirmed in the present investigation. Here, for the fully populated equilibrium point, Hopf bifurcations were found for a range of parameter values. Furthermore, the limit cycles found changed stability as we moved along the solution branch. The Floquet multipliers calculated for these solutions along this branch indicated that there was no period doubling bifurcation present, typically a route to chaos, for the parameter values used.
The periods of the self-sustained oscillations found in this investigation are typically of the order of half a month in dimensional values. This is significantly different to the case in which nutrient supply is constant; for that (degenerate) case, Hadley and Forbes [2] found that centre-type behaviour had periods corresponding to about one day. Thus, while diurnal forcing can result in nonlinear resonances in the case of constant nutrient concentration, it would not be expected to have much effect when nutrients are scarce, as in the present investigation. Here, it is to be expected that seasonal, rather than diurnal, forcing would have a more significant effect on the population dynamics. This is beyond the scope of the present study, however.
It should also be noted that Figure 2 shows the Hopf curve appearing for a small value for and disappearing as increases beyond a certain value. This may be an explanation as to why the Phytoplankton stimulate the Bacteria when nutrients are low. The limit cycles only appear for certain parameter values and particularly when N is relatively small. In effect, the strategy does not work when is abundant. This is consistent with the findings of Hadley and Forbes [2] where no limit cycles were found in that case.
In Section 3 we used the values for the interaction constants stated in Stones’ [1] paper; the constants simulate the system in a stressed state so that phytoplankton releases extracellular organic carbon (EOC). Although our interpretation of the relationship between phytoplankton and bacteria is different to that of Stone [1], we have shown that in our case the system finds a set of stable solutions of increasing amplitude for all populations, which remains stable as the reproductive rate of bacteria increases up until the point a fold bifurcation is reached and stability is lost. We must consider the fact that we have changed the values for the interaction function , between the bacteria and phytoplankton, making it considerably less than that derived from Stones’ [1] dimensional rate constants. In addition, the range of values chosen for the reproductive rate of the bacteria is less than that chosen by Stone [1]. These interaction functions are for the nondimensional system (2.6). A change to is achieved by increasing and a change in can be achieved by either increasing and at the same rate or by decreasing . Although we have changed these values slightly to find a Hopf bifurcation and limit cycles, we would consider it reasonable to assume that they still fall within acceptable ranges for the biological processes they represent.
We have chosen a Michaelis-Menten term for the nutrient uptake. There are of course other possibilities, although any law that limits the rate for large nutrient concentration might reasonably be expected to behave at least qualitatively similarly to the results presented here. We have not considered migration or alternative models for the interaction between the bacteria and phytoplankton. Nevertheless, this study has shed some light on the dynamics of this system. In addition, we have assumed that the mass of nutrient per organism is roughly constant. The dynamics of this system may possibly become more elaborate if this assumption were to be relaxed. However, these are considerations for further study.
Acknowledgment
Comments by an anonymous Referee are gratefully acknowledged.