# Dynamical Systems Analysis of a Five-Dimensional Trophic Food Web Model in the Southern Oceans

Authors: Scott A. Hadley; Lawrence K. Forbes
Journal: Journal of Applied Mathematics
 (2009)

---

## 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.


---

## Body

## 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 Section2 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 Figure1, from Stones’ [1] original paper, where the direction of the interaction is given by the arrows. There are five interacting components, namely the Bacteria (B), Phytoplankton (P), Protozoa (R), Zooplankton (Z), and Nutrients (N) as shown in Figure 1. The interactions in the diagram are described by the system of five differential equationsFigure 1
Stone’s Compartmental model showing interaction between components. The direction of the arrows indicates a direct positive influence by one component on another. HereR, B, N, P, Z represent protozoa, bacteria, nutrients, phytoplankton, and zooplankton, respectively.(2.1)dBdt=rb-erRB+eiBP,dPdt=rpP-ezPZ-eiBP,dZdt=ezPZ-dzZ,dRdt=erRB-drR,1mdNdt=-rbB-rpP+drR+dzZ.
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, rb and rp,  respectively, depend on the nutrient concentration N.  We are using a Michaelis-Menten uptake term for nutrients by Bacteria and Phytoplankton, so that the rates are governed by(2.2)rb(N)=kbNab+bbN,rp(N)=kpNap+bpN.
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 ai, bi, ki  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 quantityei  is the interaction rate between Phytoplankton and Bacteria. The constants dr and dz are the mortality rates of the Protozoa and Zooplankton, respectively. The final quantity  m  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 N0 say. For consistency with their results and using (2.2), we obtain the initial value for the reproductive rates(2.3)rb0=kbN0ab+bbN0,rp0=kpN0ap+bpN0.
We incorporate (2.3) into the nutrient-dependant reproductive rates by forming the ratio(2.4)rb(N)rb0=kbN/(ab+bbN)kbN0/(ab+bbN0).
The nutrient-dependant reproductive rates  rb(N)  and  rp(N)  are now expressed in the equivalent form(2.5)rb(N)=rb0(N/N0)1+cb0(N/N0-1),rp(N)=rp0(N/N0)1+cp0(N/N0-1).
In this way, if =N0, these rates would be the same as in the nutrient-independent case examined by Hadley and Forbes [2].  Here cb0=bbN0/(ab+bbN0) and cp0=bpN0/(ap+bpN0) 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 rp0/er, which is a naturally occurring measure of population size.  The nutrient concentration N  is scaled with respect to the quantity N0 which is arbitrary, and so this choice has no implications for the dynamics of the system.  Time t is made dimensionless using the quantity  1/rp0, which is a time scale linked roughly to the life-cycle of the phytoplankton.  In these nondimensional variables,   (2.1) become(2.6)dBdt=β[N1+cb0(N-1)]B-RB+ηBP,dPdt=[N1+cp0(N-1)]P-αPZ-ηBP,dZdt=αPZ-δZ,dRdt=RB-γR,dNdt=μ[-βNB1+cb0(N-1)-NP1+cp0(N-1)+γR+δZ].
The constant μ is such that 0<μ<1.There are five nondimensional parameter groupings in the system (2.6).  These are(2.7)α=ezer,β=rb0rp0,γ=drrp0,δ=dzrp0,η=eier.
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(B,P,Z,R,N) which satisfy(3.1)dBdt=dPdt=dZdt=dRdt=dNdt=0.
We choose Neq=N0, where now N0 represents  some arbitrary dimensionless initial concentration of nutrient. This yields the four equilibrium points(3.2)(Beq,Peq,Zeq,Req,Neq)=(0,0,0,0,N0),(Beq,Peq,Zeq,Req,Neq)=(γ,0,0,βK,N0),(Beq,Peq,Zeq,Req,Neq)=(0,δα,K1α,0,N0),(Beq,Peq,Zeq,Req,Neq)=(γ,δα,K1-ηγα,βK+ηδα,N0),
where
(3.3)K=N01+cb0(N0-1),K1=N01+cp0(N0-1).
In (3.3) the constants cb0  and  cp0 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(3.4)B(t)=Beq+εB1+O(ε2),P(t)=Peq+εP1+O(ε2),Z(t)=Zeq+εZ1+O(ε2),R(t)=Req+εR1+O(ε2),N(t)=Neq+εN1+O(ε2).
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(3.5)ddt[B1P1Z1R1N1]=[J11ηBeq0-BeqJ15-ηPeqJ22-αPeq0J250αZeqαPeq-δ00Req00Beq-γ0-μβK-μK1μδμγJ55][B1P1Z1R1N1].
Here,(3.6)J11=βK-Req+ηPeq,J15=β(-cb0NeqBeq(1+cb0(Neq-1))2+Beq1+cb0(Neq-1)),J22=K1-αZeq-ηBeq,J25=-cp0NeqPeq(1+cp0(Neq-1))2+Peq1+cp0(Neq-1),J55=-μ(J15+J25)
are defined for convenience.We are only interested here in the last steady-state(γ,δ/α,K1-ηγ/α,βK+ηδ/α,N0) in   (3.2), where all the populations survive.  These values are substituted into (3.5) to produce the Jacobian matrix(3.7)|0γη0-γE15-ηδα0-δ0E250E32000E410000-μβK-μK1μδμγE55|.
In this expression (3.7), it has proved convenient to define the intermediate quantities(3.8)E15=βγ(1-cb0)(1+cb0(N0-1))2,E25=δ(1-cp0)α(1+cp0(N0-1))2,E32=K1-ηγ,E41=βK+ηδα,E55=-μ(E15+E25).
The eigenvalues λ for the linearized system (3.5), with coefficient matrix (3.7), are found from the characteristic equation(3.9)-λ5+Q4λ4+Q3λ3+Q2λ2+Q1λ=0,
in which it is convenient to define(3.10)Q1=-E32E41γδ+E25E41μγ2η-μγK1E25E41-E15E32μδ2ηα-E15E32μβδK,Q2=E41E55γ+E15E41μγ+E32E55δ+E25E32μδ+E55η2δγα-E25μβηγK+E15μδηαK1,Q3=-E41γ-E32δ-E25μK1-η2δγα-E15μβK,Q4=E55.
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(3.11)(λ+ip)(λ-ip)(-λ3+b2λ2+b1λ)=0.
We now expand (3.11) into the form(3.12)-λ5+b2λ4+(b1-p2)λ3+p2b2λ2+p2b1λ=0.
We perform a comparison between the corresponding coefficients of powers of λ  given in (3.9) and (3.12) so that(3.13)b2=Q4,b1-p2=Q3,p2b2=Q2,p2b1=Q1,
and then solve these sets of criteria to find p and the coefficients b1 and b2  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);(3.14)(i)Q2Q4>0,(ii)Q1Q42-Q3Q2Q4-Q22=0,(iii)Q32+4Q1>0,(iv)Q1<0.
The inequalities (i), (iii), and (iv) are a consequence of the need to ensure that the quantity p2 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 Q1,…,Q4. 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 (N0-β) 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.Figure 2
The Hopf curve: The curve shows the location in(N0,β) space at which oscillatory solutions are born. The smooth line shows the supercritical Hopf points and the dashed lines shows the subcritical Hopf points.  Here α=0.4, γ=1.2, δ=0.4, η=0.01, cb0=0.1, cp0=0.1, μ=0.1.Both the supercritical and subcritical limit cycles emerge to the right in Figure2, 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  N0.  The Hopf bifurcation first appears at N0=1.092. Each successive value of N0(N0>1.092) 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 N0=10.82, 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.


## 3.1. Steady-State Populations

We now look for the steady-state solutions(B,P,Z,R,N) which satisfy(3.1)dBdt=dPdt=dZdt=dRdt=dNdt=0.
We choose Neq=N0, where now N0 represents  some arbitrary dimensionless initial concentration of nutrient. This yields the four equilibrium points(3.2)(Beq,Peq,Zeq,Req,Neq)=(0,0,0,0,N0),(Beq,Peq,Zeq,Req,Neq)=(γ,0,0,βK,N0),(Beq,Peq,Zeq,Req,Neq)=(0,δα,K1α,0,N0),(Beq,Peq,Zeq,Req,Neq)=(γ,δα,K1-ηγα,βK+ηδα,N0),
where
(3.3)K=N01+cb0(N0-1),K1=N01+cp0(N0-1).
In (3.3) the constants cb0  and  cp0 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(3.4)B(t)=Beq+εB1+O(ε2),P(t)=Peq+εP1+O(ε2),Z(t)=Zeq+εZ1+O(ε2),R(t)=Req+εR1+O(ε2),N(t)=Neq+εN1+O(ε2).
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(3.5)ddt[B1P1Z1R1N1]=[J11ηBeq0-BeqJ15-ηPeqJ22-αPeq0J250αZeqαPeq-δ00Req00Beq-γ0-μβK-μK1μδμγJ55][B1P1Z1R1N1].
Here,(3.6)J11=βK-Req+ηPeq,J15=β(-cb0NeqBeq(1+cb0(Neq-1))2+Beq1+cb0(Neq-1)),J22=K1-αZeq-ηBeq,J25=-cp0NeqPeq(1+cp0(Neq-1))2+Peq1+cp0(Neq-1),J55=-μ(J15+J25)
are defined for convenience.We are only interested here in the last steady-state(γ,δ/α,K1-ηγ/α,βK+ηδ/α,N0) in   (3.2), where all the populations survive.  These values are substituted into (3.5) to produce the Jacobian matrix(3.7)|0γη0-γE15-ηδα0-δ0E250E32000E410000-μβK-μK1μδμγE55|.
In this expression (3.7), it has proved convenient to define the intermediate quantities(3.8)E15=βγ(1-cb0)(1+cb0(N0-1))2,E25=δ(1-cp0)α(1+cp0(N0-1))2,E32=K1-ηγ,E41=βK+ηδα,E55=-μ(E15+E25).
The eigenvalues λ for the linearized system (3.5), with coefficient matrix (3.7), are found from the characteristic equation(3.9)-λ5+Q4λ4+Q3λ3+Q2λ2+Q1λ=0,
in which it is convenient to define(3.10)Q1=-E32E41γδ+E25E41μγ2η-μγK1E25E41-E15E32μδ2ηα-E15E32μβδK,Q2=E41E55γ+E15E41μγ+E32E55δ+E25E32μδ+E55η2δγα-E25μβηγK+E15μδηαK1,Q3=-E41γ-E32δ-E25μK1-η2δγα-E15μβK,Q4=E55.
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(3.11)(λ+ip)(λ-ip)(-λ3+b2λ2+b1λ)=0.
We now expand (3.11) into the form(3.12)-λ5+b2λ4+(b1-p2)λ3+p2b2λ2+p2b1λ=0.
We perform a comparison between the corresponding coefficients of powers of λ  given in (3.9) and (3.12) so that(3.13)b2=Q4,b1-p2=Q3,p2b2=Q2,p2b1=Q1,
and then solve these sets of criteria to find p and the coefficients b1 and b2  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);(3.14)(i)Q2Q4>0,(ii)Q1Q42-Q3Q2Q4-Q22=0,(iii)Q32+4Q1>0,(iv)Q1<0.
The inequalities (i), (iii), and (iv) are a consequence of the need to ensure that the quantity p2 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 Q1,…,Q4. 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 (N0-β) 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.Figure 2
The Hopf curve: The curve shows the location in(N0,β) space at which oscillatory solutions are born. The smooth line shows the supercritical Hopf points and the dashed lines shows the subcritical Hopf points.  Here α=0.4, γ=1.2, δ=0.4, η=0.01, cb0=0.1, cp0=0.1, μ=0.1.Both the supercritical and subcritical limit cycles emerge to the right in Figure2, 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  N0.  The Hopf bifurcation first appears at N0=1.092. Each successive value of N0(N0>1.092) 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 N0=10.82, 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(4.1)τ=t2πPt,
where Pt is the period of oscillation and is as yet unknown. We made a guess for the initial conditions [B(0),P(0),Z(0),R(0),N(0)] and Pt.  Using MATLAB, we integrated the set of rescaled equations to find the values [B(2π),P(2π),Z(2π),R(2π),N(2π)] after one complete period. If the solution is truly periodic, then the initial conditions should match these values.  We used the equilibrium point Beq=γ 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 B(t)=γ using the numerical MATLAB code, and used these as initial values, including a guess for the period Pt which was also estimated from the MATLAB results.Newton’s method was used to adjust the estimates ofP(0),  Z(0), R(0),  N(0), and Pt  so that the residual quantities B(0)-B(2π), P(0)-P(2π), Z(0)-Z(2π), R(0)-R(2π), and N(0)-N(2π)  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 2π-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 |λi|<1 for all i, then the limit cycle is stable; however, if any one of the |λi|>1, then the limit cycle is unstable.  The result of our numerical analysis is now considered.  We started with the initial set of parameter values(4.2)α=0.4,β=0.34,γ=1.2,δ=0.4,η=0.01,cp0=0.1,cb0=0.1,N0=2.
The locations of the Hopf points were obtained from Figure 2.Figure2 indicates that as N0 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 N0=2.  These two branches then begin to converge again until they meet at  N0 = 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 N0.  The behaviour of these limit cycles at a fixed value of N0  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
The amplitude of the limit cycles ofB(t) formed at successive values of β  for fixed N0=2. Here α=0.4, γ=1.2, δ=0.4, η=0.01, cb0=0.1, cp0=0.1, μ=0.1.Figure3 illustrates the amplitude of the limit cycle oscillations for the Bacteria, as the reproduction rate  β is varied.  At β=0.3375, there is a supercritical Hopf bifurcation formed (this stable structure travels to the right as N0  increases). This is precisely the value predicted in Figure 2 for N0=2, 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 10-9. This was done so as to improve numerical accuracy.The dashed lines in Figure3 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 Figure3, there is a second subcritical Hopf bifurcation at β=0.3675 at which an unstable limit cycle oscillation emerges from the steady-state and travels to the right as N0 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  N0=2, as illustrated in Figure 2.     Another feature that can be seen from Figure 3 is the fold bifurcation occurring at β=0.369. 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 Figure3 has also been computed, and changes as β varies. We found that for the supercritical Hopf point, β=0.3375, the nondimensional period is approximately 7.2236.  This equates to about 14.5 days in dimensional terms. Similarly, for the subcritical Hopf point at  β=0.375, 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.Figure4 shows that at the value β=0.368, 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 β=0.368 using Floquet theory.  The  results of this analysis are illustrated in   Figures 5 and 6.Figure 4
The unstable (dashed line) and stable (solid line) limit cycles for the solutionB(t) at the same value of β=0.368.Figure 5
The eigenvalues of the stable limit cycle displayed on the unit circle.Figure 6
The eigenvalues of the unstable limit cycle displayed on the unit circle.Figure5 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  |λ|=1 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.Figure 7
The stable limit cycle (red) and the perturbation solutions (blue) for perturbations above and below the limit cycle.Figure 8
The unstable limit cycle (red) and the perturbation solutions (blue) for perturbations above and below the limit cycle.Figures7 and 8 show a perturbation to both the stable and unstable limit cycles formed at β=0.368. 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 Figure2 shows the Hopf curve appearing for a small value for N and disappearing as N 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 N is abundant.  This is consistent with the findings of Hadley and Forbes [2] where no limit cycles were found in that case.In Section3 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 rb and a change in η can be achieved by either increasing er and ez at the same rate or by decreasing ei.  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.



---

*Source: 575047-2010-02-16.xml*