# Bifurcation Analysis of a 5D Nutrient, Plankton,Limnothrissa miodon Model with Hydrocynus vittatus Predation

Authors: Farikayi K. Mutasa; Brian Jones; Itai H. Tendaupenyu; Tamuka Nhiwatiwa; Mzime R. Ndebele-Murisa
Journal: Journal of Applied Mathematics
 (2022)

---

## Abstract

In this paper, we construct and analyze a theoretical, deterministic5D mathematical model of Limnothrissa miodon with nutrients, phytoplankton, zooplankton, and Hydrocynus vittatus predation. Local stability analysis results agree with the numerical simulations in that the coexistence equilibrium is locally stable provided that certain conditions are satisfied. The coexistence equilibrium is globally stable if certain conditions are met. Existence, stability, and direction of Hopf bifurcations are derived for some parameters. Bifurcation analysis shows that the model undergoes Hopf bifurcation at the coexistence point for the zooplankton growth rate with periodic doubling leading to chaos.


---

## Body

## 1. Introduction

Hydrocynus vittatus (Castelnau, 1861), also referred to as tigerfish, is the major predator of Limnothrissa miodon (Boulenger, 1906), also referred to as kapenta in Lake Kariba [1]. It is therefore important to investigate mathematically the role that tigerfish plays in the dynamics of Limnothrissa miodon. This paper begins by formulating and analyzing a deterministic Limnothrissa miodon model. The model has 5 classes, and these are as follows: concentration of nutrients, population density of phytoplankton, zooplankton population density, density of the Limnothrissa miodon population, and population density of tigerfish. The densities in each class are functions of time and are denoted by Nt, Pt, Zt, Lt, and Rt, respectively. The model is analyzed to determine the effect of predation on the population density of Limnothrissa miodon using qualitative techniques. Numerical simulations are done to illustrate the dynamics of the Limnothrissa miodon model.Mathematical modeling of theLimnothrissa miodon model with tigerfish predation will give us an insight into the dynamics of the kapenta fishery in Lake Kariba. A deterministic model that involves nutrients, phytoplankton, zooplankton, Limnothrissa miodon, and tigerfish has not been formulated and analyzed. In this paper, we formulate and analyze a deterministic, continuous, dynamical system which consists of ordinary differential equations that describe the dynamics of Limnothrissa miodon in the presence of nutrients, phytoplankton, and zooplankton and with tigerfish predation. The Limnothrissa miodon model will help in our understanding of the dynamics of the aquatic ecosystem in the kapenta fishery in Lake Kariba.The major predator in Lake Kariba is tigerfish [2, 3], and after the introduction of Limnothrissa miodon into Lake Kariba, they became a major prey item for tigerfish. Since kapenta inhabit deeper pelagic waters, tigerfish are now also found in that habitat close to the surface. According to Bell-Cross [4], the tigerfish is an efficient and extremely active predator which preys on fish of up to 40% its length. A number of studies have been done to assess the diet of Hydrocynus vittatus [5–8]. Results from the study by Mhlanga [7] showed that Limnothrissa miodon is the dominant food item in the diet of the piscivorous Hydrocynus vittatus. Stable isotope analysis by Marufu et al. [8] showed that Limnothrissa miodon is still the dominant food item consumed by tigerfish. Mhlanga [9] obtained a natural, fishing, total mortality, and exploitation rate of 0.66, 0.335, 0.995, and 0.337, respectively, of tigerfish from the Bumi Basin of Lake Kariba and the Ume River. Balon [10] obtained an instantaneous mortality coefficient for the inshore tigerfish of 0.52, which was similar to the one obtained by Langerman [11] of 0.58. Marshall [12] obtained a correlation of r=0.89 between the abundance of Limnothrissa miodon and Hydrocynus vittatus. Takano and Subramaniam [13] concluded in their study that tigerfish predation and increased fishing pressure are the major factors contributing to the natural mortality of Limnothrissa miodon.Pal and Chatterjee [14] showed the existence of Hopf bifurcations for the phytoplankton growth rate, phytoplankton carrying capacity, time delay, and fish mortality rate, in a plankton-fish model. According to Raw et al. [15], a stable interior point, period-one limit cycles, multiple-period cycles, and chaotic attractors were observed for the zooplankton growth rate bifurcation parameter in a plankton-fish model. They suggested that chaos in plankton-fish dynamics is a result of an excess of predation rate. A hopf bifurcation was observed for the phytoplankton growth rate and harvesting effort [15]. Panja and Jana [16] investigated a plankton-fish model and found that zooplankton consumption rate, fish harvesting rate, and half saturation significantly alter the model stability through a Hopf bifurcation as the parameters are varied.Dynamical systems have not been used to understand how tigerfish predation describes and influences the dynamics of kapenta fish populations in Lake Kariba. By formulating a mathematical model and analyzing it, we will be able to qualitatively explain the impact of predation on the levels of kapenta fish. The qualitative behavior of the solutions of the dynamical system is investigated for a set of parameters through bifurcation analysis.The remainder of this paper presents the materials and methods which describe the study area, data collection in Section2; model formulation, positivity, and existence of solutions in Section 3; equilibrium states and their stability in Section 4; and bifurcation analysis in Section 5. Numerical simulations are presented in Section 6, and concluding remarks in Section 7.


## 2. Materials and Methods

### 2.1. Study Area

Lake Kariba (277 km long; about 5364 km2 in surface area; 160 km3 capacity; 29 m mean depth; and 120 m maximum depth) is located in a tropical area with seasonal rainfall on the Zambezi River between latitudes 1628′ to 1804′S and longitudes 2642′ to 2903′E [17] and was formed by damming the Zambezi River at the Kariba gorge in 1958 and was filled in 1963 [18]. Lake Kariba has an average width of 19.4 km, although the widest portion is 40 km, and is 486 m above sea level, and the shoreline is approximately 2164 km [19]. The lake is almost equally shared by the two riparian countries, Zambia and Zimbabwe, and its catchment area covers 663817 km2 extending over parts of Angola, Zambia, Namibia, Botswana, and Zimbabwe [20]. The offshore single-species pelagic kapenta fishery is highly mechanized and performed by light attraction and lift nets from pontoon rigs and is licence-controlled [21].


### 2.2. Data Collection

The data used in this study was obtained from the Lake Kariba Fisheries Research Institute (LKFRI) and the University of Zimbabwe Lake Kariba Research Station (UZLKRS). The Lake Kariba Fisheries Research Institute collects data on catch, effort in the experimental gillnet, inshore artisanal, and offshore kapenta fisheries. The catch data is measured in metric tonnes (wet weight), and fishing effort is the number of nights fished. The CPUE is the kapenta catch that is landed by a boat after a night of fishing and is measured in tonnes/boat/night. It is an important parameter in fisheries management as it is an indicator of fish abundance and economic performance of the fishery [22]. The University of Zimbabwe Lake Kariba Research Institute collects data on water quality of the lake. Figure 1 shows the monthly average time series of kapenta catch and tigerfish bycatch in tonnes in the Lake Kariba kapenta fishery from 1974 to 2018.Figure 1
Time series plot of kapenta and tigerfish catch.The predator-prey relationship between the tigerfish and kapenta is shown in Figure1. From Figure 1, it is apparent that the tigerfish bycatch and kapenta catch show cyclical behavior and that the tigerfish bycatch generally tracks the peaks in the kapenta catch.


## 2.1. Study Area

Lake Kariba (277 km long; about 5364 km2 in surface area; 160 km3 capacity; 29 m mean depth; and 120 m maximum depth) is located in a tropical area with seasonal rainfall on the Zambezi River between latitudes 1628′ to 1804′S and longitudes 2642′ to 2903′E [17] and was formed by damming the Zambezi River at the Kariba gorge in 1958 and was filled in 1963 [18]. Lake Kariba has an average width of 19.4 km, although the widest portion is 40 km, and is 486 m above sea level, and the shoreline is approximately 2164 km [19]. The lake is almost equally shared by the two riparian countries, Zambia and Zimbabwe, and its catchment area covers 663817 km2 extending over parts of Angola, Zambia, Namibia, Botswana, and Zimbabwe [20]. The offshore single-species pelagic kapenta fishery is highly mechanized and performed by light attraction and lift nets from pontoon rigs and is licence-controlled [21].


## 2.2. Data Collection

The data used in this study was obtained from the Lake Kariba Fisheries Research Institute (LKFRI) and the University of Zimbabwe Lake Kariba Research Station (UZLKRS). The Lake Kariba Fisheries Research Institute collects data on catch, effort in the experimental gillnet, inshore artisanal, and offshore kapenta fisheries. The catch data is measured in metric tonnes (wet weight), and fishing effort is the number of nights fished. The CPUE is the kapenta catch that is landed by a boat after a night of fishing and is measured in tonnes/boat/night. It is an important parameter in fisheries management as it is an indicator of fish abundance and economic performance of the fishery [22]. The University of Zimbabwe Lake Kariba Research Institute collects data on water quality of the lake. Figure 1 shows the monthly average time series of kapenta catch and tigerfish bycatch in tonnes in the Lake Kariba kapenta fishery from 1974 to 2018.Figure 1
Time series plot of kapenta and tigerfish catch.The predator-prey relationship between the tigerfish and kapenta is shown in Figure1. From Figure 1, it is apparent that the tigerfish bycatch and kapenta catch show cyclical behavior and that the tigerfish bycatch generally tracks the peaks in the kapenta catch.


## 3. Model Formulation

The model has 5 classes:N denoting the concentration of nutrients, P is the population density of phytoplankton, Z is the zooplankton population density, L is the density of the Limnothrissa miodon population, and R is the density of the Hydrocynus vittatus population. The densities in each class are functions of time and are denoted by Nt, Pt, ZtLt, and Rt, respectively. The Limnothrissa miodon model [23] is developed to include predation by Hydrocynus vittatus. It is assumed that nutrients enter the water body at the rate a, where a>0 is a constant and the nutrients are depleted naturally at a constant rate μ0. The nutrients are depleted by phytoplankton at a rate of σ1NP. The growth rate of phytoplankton is γ1NP. It is assumed that the depletion rate of phytoplankton caused by mortality is proportional to P. Phytoplankton is depleted by zooplankton at a rate σ2PZ. The depletion of phytoplankton per unit time by zooplankton is given by σ2PZ and is the modified Holling’s type I response [24], which refers to the change in density of the phytoplankton per unit time per zooplankton as the phytoplankton population density changes. The growth rate of zooplankton is γ2PZ. It is assumed that the depletion rate of zooplankton caused by mortality is proportional to Z. The functional response of zooplankton to the Limnothrissa miodon given by σ3ZL is of the modified Holling’s type I response, which refers to the change in density of the zooplankton per unit time per Limnothrissa miodon as the zooplankton population density changes. The growth rate of Limnothrissa miodon is γ3ZL. It is assumed that the depletion rate of Limnothrissa miodon caused by mortality is proportional to L and its rate of depletion caused by crowding is proportional to L2. Kapenta are harvested at a rate qEL, where q is the catchability coefficient and E is the effort measured as boat nights. Tigerfish search and feed on kapenta; therefore, we use Holling’s type II functional response. The feeding rate saturates at the maximum feeding rate σ4. The feeding rate is half maximal at L=d. The response fL=σ4L/d+L [25] models the fact that the consumption of kapenta is limited by satiation of tigerfish, handling time (killing and eating) and time spent hunting kapenta. The growth rate of Limnothrissa miodon is γ4LR/d+L. The tigerfish are harvested at a rate κηR, where κ is the catchability coefficient and η is the effort. The nonlinear dynamical system is
(1)dNdt=a−μ0N−σ1NP,dPdt=γ1NP−μ1P−σ2PZ,dZdt=γ2PZ−μ2Z−σ3ZL,dLdt=γ3ZL−μ3L−σ30L2−qEL−σ4LRd+L,dRdt=γ4LRd+L−κηR−μ4R,with initial condition(2)N0=ψ10,P0=ψ20,Z0=ψ30,L0=ψ40,R0=ψ50,ψi0>0,i=1,2,3,4,5.and define(3)Ω=N,P,Z,L,R∈ℝ5N≥0,P≥0,Z≥0,L≥0,R≥0,to be the mathematically feasible region. The coefficientσ30 is a positive constant for the crowding of the Limnothrissa miodon population. σ1, σ2, σ3,σ4 are positive constants of proportionality. The μi’s for i=0,1,2,3,4 are the depletion rate coefficients.


### 3.1. Positivity of Solutions

Model system (1) describes the dynamics of an ecosystem and it is necessary to prove that the concentrations of nutrients, and the densities of phytoplankton, zooplankton, and kapenta are positive for all time. For positive initial data for the ecosystem model (1), we prove that the solutions will remain positive ∀t≥0.Theorem 1.
Let the initial data beNt≥0,Pt≥0,Zt≥0,Lt≥0,Rt≥0. Then, solutions of Nt, Pt, Zt, Lt, Rt of system (1) are positive ∀t≥0.Proof.
Considering the variableNt in 0,T, from the first equation of model (1), it follows that
(4)Ṅt≥−μ0Nt−σ1NtPt,∀t∈0,T.
Hence, we obtain(5)Nt≥N0exp∫0t−μ0−σ1Psds≥0,∀t∈0,T.
From the second equation of model (1), it follows that
(6)Ṗt≥−μ1Pt−σ2PtZt,∀t∈0,T.
Direct integration of (6) results in
(7)Pt≥P0exp∫0t−μ1−σ2Zsds≥0,∀t∈0,T.
From the third equation of model (1), it follows that
(8)Żt≥−μ2Zt−σ3ZtLt,∀t∈0,T.
Direct integration of (8) results in
(9)Zt≥Z0exp∫0t−μ2−σ3Lsds≥0,∀t∈0,T.
Considering the variableLt in 0,T, from the fourth equation of model (1), it follows that
(10)L̇t≥−μ3+qE+σ30Lt+σ4Rtd+LtLt,∀t∈0,T.
Direct integration of (10) results in
(11)Lt≥L0exp−∫0tμ3+qE+σ30Ls+σ4Rsd+Lsds≥0,∀t∈0,T.
From the fifth equation of model (1), it follows that
(12)Ṙt≥−μ4+κηRt,∀t∈0,T.
Direct integration of (12) results in
(13)Rt≥R0exp−μ4+κηt≥0,∀t∈0,T.
Therefore, the solutions of system (1) with initial condition (2) remains positive ∀t≥0.


### 3.2. Existence of Solutions

Theorem 2.
A solution of model system (1) is feasible.Proof.
It is necessary to show that system (1) is dissipative; that is, all feasible solutions are uniformly bounded in Ω⊂ℝ5. Let Nt,Pt,Zt,Lt,Rt∈ℝ5 be any solution of system (1) with nonnegative initial conditions.
LetAt=Nt+Pt+Zt+Lt+Rt, then
(14)dAdt=a−μ0N−σ1NP+γ1NP−μ1P−σ2PZ+γ2PZ−μ2Z−σ3ZL+γ3ZL−μ3L−σ30L2−qEL−σ4LRd+L+γ4LRd+L−κηR−μ4R,=a−μ0N−μ1P−μ2Z−μ3+qEL−σ30L2−κηR−μ4R+γ1−σ1NP+γ2−σ2PZ+γ3−σ3ZL+γ4−σ4LRd+L,≤a−μ0N−μ1P−μ2Z−μ3+qEL−κη+μ4R,≤a−uAt,where u=minμ0,μ1,μ2,μ3+qE,κη+μ4. Thus,
(15)dAtdt+uAt≤a.
The solution to Equation (15) is
(16)0<AN,P,Z,L,R≤au1−e−ut+AN0,P0,Z0,L0,R0e−ut,
ast⟶∞, (16) becomes
(17)0<AN,P,Z,L,R≤au.
Therefore, all solutions of the system (1) enter the feasible region,
(18)Ω=Nt,Pt,Zt,Lt,Rt∈ℝ+5:A≤au+ς,∀ς>0.
This completes the proof of the theorem.


## 3.1. Positivity of Solutions

Model system (1) describes the dynamics of an ecosystem and it is necessary to prove that the concentrations of nutrients, and the densities of phytoplankton, zooplankton, and kapenta are positive for all time. For positive initial data for the ecosystem model (1), we prove that the solutions will remain positive ∀t≥0.Theorem 1.
Let the initial data beNt≥0,Pt≥0,Zt≥0,Lt≥0,Rt≥0. Then, solutions of Nt, Pt, Zt, Lt, Rt of system (1) are positive ∀t≥0.Proof.
Considering the variableNt in 0,T, from the first equation of model (1), it follows that
(4)Ṅt≥−μ0Nt−σ1NtPt,∀t∈0,T.
Hence, we obtain(5)Nt≥N0exp∫0t−μ0−σ1Psds≥0,∀t∈0,T.
From the second equation of model (1), it follows that
(6)Ṗt≥−μ1Pt−σ2PtZt,∀t∈0,T.
Direct integration of (6) results in
(7)Pt≥P0exp∫0t−μ1−σ2Zsds≥0,∀t∈0,T.
From the third equation of model (1), it follows that
(8)Żt≥−μ2Zt−σ3ZtLt,∀t∈0,T.
Direct integration of (8) results in
(9)Zt≥Z0exp∫0t−μ2−σ3Lsds≥0,∀t∈0,T.
Considering the variableLt in 0,T, from the fourth equation of model (1), it follows that
(10)L̇t≥−μ3+qE+σ30Lt+σ4Rtd+LtLt,∀t∈0,T.
Direct integration of (10) results in
(11)Lt≥L0exp−∫0tμ3+qE+σ30Ls+σ4Rsd+Lsds≥0,∀t∈0,T.
From the fifth equation of model (1), it follows that
(12)Ṙt≥−μ4+κηRt,∀t∈0,T.
Direct integration of (12) results in
(13)Rt≥R0exp−μ4+κηt≥0,∀t∈0,T.
Therefore, the solutions of system (1) with initial condition (2) remains positive ∀t≥0.


## 3.2. Existence of Solutions

Theorem 2.
A solution of model system (1) is feasible.Proof.
It is necessary to show that system (1) is dissipative; that is, all feasible solutions are uniformly bounded in Ω⊂ℝ5. Let Nt,Pt,Zt,Lt,Rt∈ℝ5 be any solution of system (1) with nonnegative initial conditions.
LetAt=Nt+Pt+Zt+Lt+Rt, then
(14)dAdt=a−μ0N−σ1NP+γ1NP−μ1P−σ2PZ+γ2PZ−μ2Z−σ3ZL+γ3ZL−μ3L−σ30L2−qEL−σ4LRd+L+γ4LRd+L−κηR−μ4R,=a−μ0N−μ1P−μ2Z−μ3+qEL−σ30L2−κηR−μ4R+γ1−σ1NP+γ2−σ2PZ+γ3−σ3ZL+γ4−σ4LRd+L,≤a−μ0N−μ1P−μ2Z−μ3+qEL−κη+μ4R,≤a−uAt,where u=minμ0,μ1,μ2,μ3+qE,κη+μ4. Thus,
(15)dAtdt+uAt≤a.
The solution to Equation (15) is
(16)0<AN,P,Z,L,R≤au1−e−ut+AN0,P0,Z0,L0,R0e−ut,
ast⟶∞, (16) becomes
(17)0<AN,P,Z,L,R≤au.
Therefore, all solutions of the system (1) enter the feasible region,
(18)Ω=Nt,Pt,Zt,Lt,Rt∈ℝ+5:A≤au+ς,∀ς>0.
This completes the proof of the theorem.


## 4. Equilibria and Stability Analysis

Model (1) has 6 equilibria:
(a)
The trivial equilibrium is(19)E0=0,0,0,0,0.(b)
The phytoplankton free equilibrium is(20)E1=N1∗,0,0,0,0=aμ0,0,0,0,0.(c)
The zooplankton free equilibrium is(21)E2=N2∗,P2∗,0,0,0,=μ1γ1,aγ1−μ1μ0σ1μ1,0,0,0.E2 is obtained when phytoplankton is taking part in the ecosystem, and zooplankton and Limnothrissa miodon are not taking part in the ecosystem. The phytoplankton population is not enough to support the zooplankton population. E2 exists provided that
(22)aγ1>μ1μ0.Rearranging the inequality (22), we obtain a/μ0>μ1/γ1, and this means that N1∗>N2∗. The nutrients will reach the value a/μ0 at equilibrium in the absence of phytoplankton, which is reduced to the steady state value of μ1/γ1 in the presence of phytoplankton.
(d)
TheLimnothrissa miodon free equilibrium is(23)E3=N3∗,P3∗,Z3∗,0,0,E3=aμ0+σ1P3∗,μ2γ2,γ1N3∗−μ1σ2,0,0.The zooplankton population is insufficient to support the population ofLimnothrissa miodon. When the Limnothrissa miodon population is not present in the ecosystem and both phytoplankton and zooplankton are present, E3 is attained. E3 exists if N3∗>μ1/μ1, i.e. N3∗>N2∗. In order to support the zooplankton population, more nutrients are needed in the ecosystem. From (23),
(24)E3=aγ2γ2μ0+μ2σ1,μ2γ2,aγ1γ2−μ1γ2μ0+μ2σ1σ2γ2μ0+μ2σ1,0,0.E3 exists on condition that
(25)aγ1−μ0μ1σ1μ1>μ2γ2.Inequality (25) can be rearranged to give P2∗>P3∗, meaning that the phytoplankton equilibrium is reduced in the presence of zooplankton.
(e)
The tigerfish free equilibriumE4=N4∗,P4∗,Z4∗,L4∗,0 is obtained by solving the equations:(26)a−μ0N−σ1NP=0,(27)ϕ1σ1N−μ1−σ2Z=0,(28)ϕ2σ2P−μ2−σ3L=0,(29)ϕ3σ3Z−μ3−eq−σ30L=0.Solving forN, P, Z, and L in (26), (27), (28), and (29) gives
(30)σ1σ2σ3σ30L4∗2+σ1σ2σ3μ3+eq+σ1σ2μ2σ30+μ0σ2γ2σ30+μ1σ1σ3γ3L4∗+μ0μ1γ2γ3+μ1σ1μ2γ3+μ0σ2γ2μ3+eq+σ1σ2μ2μ3+eq−aγ1γ2γ3=0.Equation (31) will have a unique positive root if the expression (32) is positive,
(31)L4∗=A12−4σ1σ2σ3σ30A2−A12σ1σ2σ3σ30>0,where
(32)A1=σ1σ2σ3μ3+eq+σ1σ2μ2σ30+μ0σ2γ2σ30+μ1σ1σ3γ3,A2=μ0μ1γ2γ3+μ1σ1μ2γ3+μ0σ2γ2μ3+eq+σ1σ2μ2μ3+eq−aγ1γ2γ3.(31) can be written as
(33)A12−4σ1σ2σ3σ30A2>A12,σ1σ2σ3σ30A2<0,σ1σ2σ3σ30μ0μ1γ2γ3+μ1σ1μ2γ3+μ0σ2γ2μ3+eq+σ1σ2μ2μ3+eq−aγ1γ2γ3<0,σ2μ3+eqμ0γ2+σ1μ2<γ3aγ1γ2−μ0μ1γ2+μ1μ2σ1,μ3+eqγ3<aγ1γ2−μ0μ1γ2+μ1μ2σ1σ2μ0γ2+σ1μ2,μ3+eqγ3<Z3∗.Therefore,L4∗ exists whenever
(34)Z3∗>μ3+eqγ3.The tigerfish free equilibrium is(35)N4∗=μ3+eqσ1σ2σ3−μ2σ1σ2σ30−μ0σ2σ30γ2+μ1σ1σ3γ3+4A3+A422σ1σ3γ1γ3,P4∗=−μ3+eqσ1σ2σ3+μ2σ1σ2σ30−μ0σ2σ30γ2−μ1σ1σ3γ3+4A3+A422σ1σ2σ30γ2,Z4∗=μ3+eqσ1σ2σ3−μ2σ1σ2σ30−μ0σ2σ30γ2−μ1σ1σ3γ3+4A3+A422σ1σ2σ3γ3,L4∗=−μ3+eqσ1σ2σ3−μ2σ1σ2σ30−μ0σ2σ30γ2−μ1σ1σ3γ3+4A3+A422σ1σ2σ3σ30,R4∗=0,where
(36)A3=aσ1σ2σ3σ30γ1γ2γ3,A4=μ3+eqσ1σ2σ3−σ2σ30μ2σ1+μ0γ2+μ1σ1σ3γ3.From equation array (35), it follows that
(37)N4∗=μ3+eqσ1σ2σ3−μ2σ1σ2σ30−μ0σ2σ30γ2+μ1σ1σ3γ3+4A3+A422σ1σ3γ1γ3,P4∗=γ1σ3γ3N4∗−μ3+eqσ2σ3+μ2σ2σ30−μ1σ3γ3σ2σ30γ2,Z4∗=σ30γ2P4∗+μ3+eqσ3−μ2σ30σ3γ3,L4∗=γ3Z4∗−μ3+eqσ30,R4∗=0.(f)
The positive interior equilibrium of system (1) is the one of biological interest. For some set of parameter values, model (1) has a unique positive interior equilibrium, E∗=N∗,P∗,Z∗,L∗,R∗ which is obtained by solving the equations:(38)a−μ0N−σ1NP=0,(39)γ1N−μ1−σ2Z=0,(40)γ2P−μ2−σ3L=0,(41)γ3Z−μ3−σ30L−qE−σ4Rd+L=0,(42)γ4Ld+L−κη−μ4=0.Solving forN, P, Z, L, and R in (38), (39), (40), (41), and (42) gives the positive interior equilibrium
(43)N∗=aγ2−γ4+κη+μ4γ2μ0+μ2σ1−γ4+κη+μ4−dσ1σ3κη+μ4,P∗=μ2−dσ3κη+μ4/−γ4+κη+μ4γ2,Z∗=−γ4+κη+μ4aγ1γ2−μ1γ2μ0+μ2σ1+dμ1σ1σ3κη+μ4σ2γ2μ0+μ2σ1−γ4+κη+μ4−dσ1σ3κη+μ4,L∗=dκη+μ4γ4−κη+μ4,R∗=γ4dγ2μ0+μ2σ1−γ4+κη+μ4−dσ1σ3κη+μ4−γ4+κη+μ4γ3μ1+σ2eq+μ3−dσ2σ30κη+μ4−aγ1γ2γ3−γ4+κη+μ42σ2σ4−γ4+κη+μ42γ2μ0+μ2σ1−γ4+κη+μ4−dσ1σ3κη+μ4.Remark 3.
For some given set of parameter values, the model (1) has no equilibrium points.Theorem 4.
The trivial equilibrium is always stable.Proof.
The auxilliary equation of the jacobian matrix atE0 is
(44)λ+μ0λ+μ1λ+μ2κ+λ+μ4eq+λ+μ3=0.
The eigenvalues of (44) are λ1=−μ0, λ2=−μ1, λ3=−μ2, λ4=−κ−μ4, and λ5=−eq−μ3. All the eigenvalues are negative; therefore, E0 is stable.Theorem 5.
The phytoplankton free equilibrium is locally asymptotically stable ifaγ1<μ0μ1.Proof.
The auxilliary equation of the jacobian matrix atE1 is
(45)λ+μ0−λ+aγ1−μ0μ1μ0λ+μ2λ+eq+μ3λ+κ+μ4=0.
The eigenvalues areλ1=−μ0,λ2=aγ1−μ0μ1/μ0,λ3=−μ2,λ4=−eq−μ3,λ5=−κ−μ4.E1 is asymptotically stable if a/μ0<μ1/γ1, meaning that N1∗<N1∗ for stability of E1.Theorem 6.
The zooplankton free equilibrium is locally asymptotically stable ifa2γ12−4aγ1μ12+4μ0μ13>0, γ2aγ1−μ0μ1/μ1σ1<μ2/γ2, and aγ1>μ0μ1.Proof.
The auxilliary equation of the Jacobian matrix atE2 is
(46)−κ−λ−μ4−eq−λ−μ3γ2aγ1−μ0μ1/μ1σ1−λ−μ2aγ12λμ1σ1+aγ12μ12σ1+γ1λ2μ12σ1−γ1μ0μ13σ1=0.
The eigenvalues areλ1=−eq−μ3, λ2=−κ−μ4, λ3=γ2aγ1−μ0μ1/μ1σ1−μ2, λ4=−a2γ12−4aγ1μ12+4μ0μ13−aγ1/2μ1 and λ5=a2γ12−4aγ1μ12+4μ0μ13−aγ1/2μ1. Hence, E2 is locally asymptotically stable if a2γ12−4aγ1μ12+4μ0μ13>0, γ2aγ1−μ0μ1/μ1σ1<μ2/γ2, and aγ1>μ0μ1.Theorem 7.
The Limnothrissa miodon free equilibrium is locally asymptotically stable if the conditions in (50) are satisfied.Proof.
The auxilliary equation of the jacobian matrix atE3 is
(47)−κ−λ−μ4−b1b2γ1λσ1−b2b3γ2λσ2−b2b3γ2μ0σ2−b22b3γ2σ1σ2−b2λ2σ1−λ3−λ2μ0b3γ3−eq−λ−μ3=0,where b1=aγ2/γ2μ0+μ2σ1, b2=μ2/γ2, and b3=aγ1γ2−μ1γ2μ0+μ2σ1/σ2γ2μ0+μ2σ1. The characteristic Equation (47) can be written as follows:
(48)Δλ=a0λ5+a1λ4+a2λ3+a3λ2+a4λ+a5=0,where
(49)a0=1,a1=b2σ1+eq+κ+μ0+μ3+μ4−b3γ3,a2=b2σ1b1γ1+eq+κ+μ3+μ4+b3b2γ2σ2−γ3b2σ1+κ+μ0+μ4+κ+μ4eq+μ3+μ0eq+κ+μ3+μ4,a3=μ0b2b3γ2σ2+κ+μ4−b3γ3+eq+μ3+b2σ1κ+μ4−b3γ3+eq+μ3+b1γ1−b3γ3+eq+κ+μ3+μ4+b3γ2σ2−b3γ3+b2σ1+eq+κ+μ3+μ4,a4=−b2b1γ1σ1κ+μ4−b3γ3+eq+μ3+b3γ2σ2b2σ1eq+κ+μ3+μ4−b3γ3b2σ1+κ+μ0+μ4+κ+μ4eq+μ3+μ0eq+κ+μ3+μ4,a5=b2b3γ2σ2κ+μ4b2σ1+μ0−b3γ3+eq+μ3.
By the Routh-Hurwitz criterion, it follows that all eigenvalues of the characteristic Equation (48) have negative real parts if
(50)a0>0,a1>0,a1a2−a0a3>0,a3a1a2−a0a3−a12a4>0,a1a3a4a2+a0a3a5a2−a12a42−a02a52−a0a32a4+2a0a1a4a5−a1a5a22>0,a5>0.Theorem 8.
The tigerfish free equilibrium is locally asymptotically stable if the conditions in (54) are satisfied.Proof.
The auxilliary equation of the jacobian matrix atE4 is
(51)λ5+λ4−A1−A2−A4+λ3A1A2+A4A2+A1A4+b4b5γ1σ1+b5b6γ2σ2+b6b7γ3σ3+λ2−A4b4b5γ1σ1−A4b5b6γ2σ2−A4b6b7γ3σ3−A2b4b5γ1σ1−A1b5b6γ2σ2−A2b5b6γ2σ2−A1b6b7γ3σ3−A1A2A4+λA2A4b4b5γ1σ1+A1A2b5b6γ2σ2+A1A4b5b6γ2σ2+A2A4b5b6γ2σ2+A1A4b6b7γ3σ3+b4b5b6b7γ1γ3σ3σ1−A1A2A4b5b6γ2σ2−A4b4b5b6b7γ1γ3σ1σ3=0,where b4=N4∗, b5=P∗, b6=Z∗, b7=L∗, A1=−b5σ−μ0, A2=b6γ3−2b7σ30−eq−μ3, A3=−b7σ4/b7+d, and A4=b7γ4/b7+d−κ−μ4. The characteristic Equation (51) can be written as follows:
(52)Δλ=a0λ5+a1λ4+a2λ3+a3λ2+a4λ+a5=0,where
(53)a0=1,a1=−A1−A2−A4,a2=A1A2+A4A2+A1A4+b4b5γ1σ1+b5b6γ2σ2+b6b7γ3σ3,a3=−A4b4b5γ1σ1−A4b5b6γ2σ2−A4b6b7γ3σ3−A2b4b5γ1σ1−A1b5b6γ2σ2−A2b5b6γ2σ2−A1b6b7γ3σ3−A1A2A4,a4=A2A4b4b5γ1σ1+A1A2b5b6γ2σ2+A1A4b5b6γ2σ2+A2A4b5b6γ2σ2+A1A4b6b7γ3σ3+b4b5b6b7γ1γ3σ3σ1,a5=−A1A2A4b5b6γ2σ2−A4b4b5b6b7γ1γ3σ1σ3.
By the Routh-Hurwitz criterion, it follows that all eigenvalues of the characteristic Equation (52) have negative real parts if
(54)a0>0,a1>0,a1a2−a0a3>0,a3a1a2−a0a3−a12a4>0,a1a3a4a2+a0a3a5a2−a12a42−a02a52−a0a32a4+2a0a1a4a5−a1a5a22>0,a5>0.Theorem 9.
If the equilibriumE∗=N∗,P∗,Z∗,L∗,R∗ exists, then, it is locally-asymptotically stable if the conditions in (60) are satisfied.Proof.
Evaluating the variational matrix atE∗ gives
(55)JE∗=−μ0−c5σ1−c4σ1000c5γ10−c5σ2000c6γ20−c6σ3000c7γ3−eq+c6γ3−μ3−c8σ4d+c7+c7c8σ4d+c72−2c7σ30−c7σ4d+c7000c8γ4d+c7−c7c8γ4d+c720,where c4=N∗, c5=P∗, c6=Z∗, c7=L∗, and c8=R∗. (55) simplifies to
(56)JE∗=B1−λ−c4σ1000c5γ1−λ−c5σ2000c6γ2−λ−c6σ3000c7γ3B2−λB3000B4−λ,where B1=−c5σ1−μ0, B2=c6γ3−c8σ4/c7+d+c7c8σ4/c7+d2−2c7σ30−eq−μ3, B3=−c7σ4/c7+d, and B4=c8γ4/c7+d−c7c8γ4/c7+d2. The eigenvalues of (56) are the roots of the auxiliary equation
(57)detJE∗=0.
The characteristic Equation (57) can be written as follows:
(58)Δλ=a0λ5+a1λ4+a2λ3+a3λ2+a4λ+a5=0,where
(59)a0=1,a1=−B1−B2,a2=B1B2−B3B4+c4c5γ1σ1+c5c6γ2σ2+c6c7γ3σ3,a3=−B2c4c5γ1σ1−B1c5c6γ2σ2−B2c5c6γ2σ2−B1c6c7,a4=−B3B4c4c5γ1σ1+B1B2c5c6γ2σ2−B3B4c5c6γ2σ2+c4c5c6c7γ1γ3σ3σ1,a5=B1B3B4c5c6γ2σ2.
By the Routh-Hurwitz criterion, it follows that all eigenvalues of the characteristic Equation (58) have negative real parts if
(60)a0>0,a1>0,a1a2−a0a3>0,a3a1a2−a0a3−a12a4>0,a1a3a4a2+a0a3a5a2−a12a42−a02a52−a0a32a4+2a0a1a4a5−a1a5a22>0,a5>0.Remark 10.
If the conditions forE∗ in (60) are not satisfied for a given set of parameter values, then the respective steady state will be unstable, and there is a possibility of oscillatory behavior for model (1).Theorem 11.
The equilibriumE∗ is globally-asymptotically stable if the conditions in (63) are satisfied for the Lyapunov function in (61).Proof.
The proof follows Lyapunov’s second method. LetN−N∗>0, P−P∗>0, Z−Z∗>0, L−L∗>0, R−R∗>0. Let VN,P,Z,L,R be a positive Lyapunov function [26, 27], such that VN∗,P∗,Z∗,L∗,R∗=0 by,
(61)VN,P,Z,L,R=β1N−N∗−N∗lnNN∗+β2P−P∗−P∗lnPP∗+β3Z−Z∗−Z∗lnZZ∗+β4L−L∗−L∗lnLL∗,+β5R−R∗−R∗lnRR∗,where βi′s, i=1,2,3,4,5 are positive constants. V is a positive definite function in the set Ψ, except at E∗ where it is zero. The rate of change of V along the solution of system (1) is given by
(62)V̇=β1N−N∗ṄN+β2P−P∗ṖP+β3Z−Z∗ŻZ+β4L−L∗L̇L+β5R−R∗ṘR,=−β1N−N∗μ0+σ1P−aN−β2P−P∗μ1+σ2Z−γ1N−β3Z−Z∗μ2+σ3L−γ2P−β4L−L∗μ3+qE+σ30L+σ4Rd+L−γ3Z−β5R−R∗κη+μ4−γ4Ld+L,=−β1N−N∗μ0+σ1P−aN−μ0−σ1P∗+aN∗−β2P−P∗μ1+σ2Z−γ1N−μ1−σ2Z∗+γ1N∗−β3Z−Z∗μ2+σ3L−γ2P−μ2−σ3L∗+γ2P∗−β4L−L∗μ3+qE+σ30L+σ4Rd+L−γ3Z−μ3−qE−σ30L∗−σ4R∗d+L∗+γ3Z∗−β5R−R∗κη+μ4−γ4Ld+L−κη−μ4+γ4L∗d+L∗,V̇≤−β1aNN∗N−N∗2−N−N∗P−P∗β1σ1−β2γ1−P−P∗Z−Z∗β2σ2−β3γ2−Z−Z∗L−L∗β3σ3−β4γ3−β4σ30L−L∗2−dL−L∗R−R∗β4σ4−β5γ4.
Then,V̇≤0 if
(63)β1σ1≥β2γ1,β2σ2≥β3γ2,β3σ3≥β4γ3,β4σ4≥β5γ4.
Thus, in the region bounded by all points (N>N∗, P>P∗, Z>Z∗, L>L∗, R>R∗) in (63),
E∗ is globally-asymptotically stable.The stability of the periodic solutions is discussed in the bifurcation analysis section.


## 5. Bifurcation Analysis

The5D system (1) is written as
(64)ẋ=Fx,η,x∈ℝ5,F∈C∞,where η is a bifurcation parameter. Bifurcation analysis of model (64) is done using the Hopf bifurcation theorem by Guckenheimer and Holmes [28].Theorem 12.
Assume that model (64) has the following characteristics
(1)
The model has a smooth equilibria curve(65)Fx∗η,η=0.(2)
The auxiliary equation(66)Q5λ=λ5+ζ1λ4+ζ2λ3+ζ3λ2+ζ4λ+ζ5=0,
has complex conjugate roots (nonhyperbolicity condition)(67)λ1,2=ληc,λ¯ηc=iω,−iω,ω>0,
and if the condition(68)Ψ=ζ3−ζ1ζ2ζ5ζ2−ζ3ζ4−ζ5−ζ1ζ42=0,
with(69)ω=ζ5−ζ1ζ4ζ3−ζ1ζ2>0,
is satisfied, and there are no other roots with zero real parts andλii=3,4,5. have nonzero real parts if
(70)ζ3−ζ1ω≠0,
and have negative real parts if(71)ζ1>0,ζ3−ζ1ζ2<0,ζ3−ζ1ω>0.(3)
The transversality condition(72)dΨηdηη=ηc≠0,
is satisfied, and then, there is a Poincare′-Andronov-Hopf bifurcation.


## 6. Numerical Simulations

Numerical simulations of the model system (1) are carried out to investigate the dynamics of the Limnothrissa miodon model for the main bifurcation parameter γ2, using the default parameter values in Equation (73)
(73)a=20.3,μ0=0.0096,σ1=0.3,γ1=0.3,μ1=0.032,σ2=0.18,γ2=0.005,μ2=0.08,σ3=0.3,σ30=0.01,γ3=0.36,μ3=0.05,q=0.00027999,E=500,κη=0.16,σ4=0.4,γ4=0.4,d=2,andμ4=0.018.The parameter values in Equation (73) are obtained from published data and others are estimates. A fourth-order Runge-Kutta numerical scheme coded in Wolfram Mathematica is used for the numerical simulations. For model system (1), the units of the variables N, P, Z, L, and R are μgl−1.


### 6.1. Parameterσ4

There is no bifurcation for varyingσ4, but a Hopf bifurcation exists for varying σ4 with γ4=σ4. The predation rate of tigerfish on kapenta is varied from 0.4 to 0.9 with γ4=σ4. Figures 2 and 3(a) show the bifurcation diagrams for the bifurcation parameter σ4.Figure 2
Bifurcation diagram for parameterσ4 versus (a) nutrients; (b) phytoplankton; (c) zooplankton; and (d) Limnothrissa miodon for model system (1) with assumed initial condition: N0=10,P0=7,Z0=4,L0=2,R0=0.5 using the default parameter values.
(a)(b)(c)(d)Figure 3
(a) Bifurcation diagram for tigerfish versusσ4. (b) Plot of real part of eigenvalues of the Jacobian matrix of model (1) for 0.4≤σ4≤0.9.
(a)(b)The Hopf bifurcation value is0.63076 for the control parameter σ4 and is shown in Figure 3(b). From the coefficient criteria in Theorem 12, we obtain σ4=σ4,c=0.63076. Equation (69) gives ω=0.225435,λ1,2=±iω=±0.4748i, and T=2π/ω=27.8714 days. From Theorem 12, Equation (71), we obtain ζ1=18.8633, ζ3−ζ1ζ2=−79.8152 and ζ3−ζ1ω=0.37348. The coefficient criteria in Equation (71) are satisfied and therefore showing existence of a simple Hopf bifurcation at the critical value σ4,c of the bifurcation parameter σ4. The transversality condition in Equation (72) of Theorem 12 is satisfied since dΨσ4/dσ4σ4=σ4,c=−246.512. The eigenvalues at σ4,c are −18.6363,−0.113515−0.08458i,−0.113515+0.084586i,−0.4748i,+0.4748i. The Hopf bifurcation is supercritical since as the bifurcation parameter σ4 is varied with γ4=σ4 the stable positive interior point loses its stability and a stable periodic orbit simultaneously appears.Figure3 shows the effect of varying σ4 for 0.4≤σ4≤0.9 with σ4=γ4. The bifurcation results show that there is a decrease in kapenta population density for 0.4≤σ4<0.63076. Further increase of σ4 for 0.63076<σ4≤0.9 results in a period orbit of period 1 with increasing amplitude. Inefficient tigerfish predation on kapenta results in a decrease in the co-existence equilibrium value for kapenta, and efficient predation increases the chance of periodic behavior in the dynamical system. Therefore, the predation rate σ4, together with the tigerfish growth rate γ4, significantly determines the nature of the predator-prey relationship between tigerfish and kapenta. Simulation results show no chaotic behavior for the bifurcation control parameter σ4.


### 6.2. Parameterγ2

The interior equilibrium,E∗=0.602831,112.216,0.82694,1.6036,0.825832 for model (1), with damped oscillations for the set of default parameters are shown in Figure 4(a) for Limnothrissa miodon and 5(a) for zooplankton, Limnothrissa miodon, and tigerfish.Figure 4
Time series plot ofLimnothrissa miodon for (a) γ2=0.005; (b) γ2=0.0095; (c) γ2=0.075; and (d) γ2=0.16 for model (1) with assumed initial condition: N0=9.8,P0=6.8,Z0=3.8,L0=1.8,R0=0.3; N0=10,P0=7,Z0=4,L0=2,R0=0.5; N0=10.2,P0=7.2,Z0=4.2,L0=2.2,R0=0.7, using the default parameter values.
(a)(b)(c)(d)Figure 5
Phase portrait of zooplankton,Limnothrissa miodon and tigerfish (a) γ2=0.005, (b) γ2=0.0095, (c) γ2=0.075, and (d) γ2=0.16 for model (1) with assumed initial condition: N0=9.8,P0=6.8,Z0=3.8,L0=1.8,R0=0.3 (time series with blue color); N0=10,P0=7,Z0=4,L0=2,R0=0.5 (time series with magenta color); N0=10.2,P0=7.2,Z0=4.2,L0=2.2,R0=0.7 (time series with purple color), using the default parameter values.
(a)(b)(c)(d)The positive interior equilibriumE∗ for model system (1) is permanent for the set of default parameter values. For the main bifurcation parameter, model (1) bifurcates into a limit cycle of period 1 at γ2,c=0.00818058. Time series plots and phase portraits for γ2=0.0095 are shown in Figures 4(b) and 5(b). Time series plots and phase portraits for γ2=0.075 are shown in Figures 4(c) and 5(c). The phase portraits show periodic doubling. Time series plots and phase portraits for γ2=0.16 are shown in Figures 4(d) and 5(d). The phase portraits show a chaotic attractor at γ2=0.16. Varying the zooplankton growth rate γ2 results in some interesting dynamics for model (1). The stability of model (1) changes from a stable coexistence equilibrium into a stable periodic orbit and then periodic doubling enroute to chaos. Therefore, we can conclude that γ2 plays an important role in the dynamics of model (1). The trajectories in Figures 4 and 5 are aperiodic as they are showing erratic behavior. The irregular fluctuations can be attributed to the nonlinearity of the model (1). By visual inspection of the phase portraits in Figures 4(d) and 5(d) and from the magnitude of the maximal Lyapunov exponent, 0.107838, it can be seen that the trajectories are sensitive to initial conditions and, therefore, it is impossible to predict the long-term behavior of model system (1) for γ2=0.16. The trajectories are random and bounded in phase space and they converge to a strange attractor, which has a complex (fractal) geometric structure.We use the Hopf bifurcation theorem by Guckenheimer and Holmes [28] together with conditions by Douskos and Markellos [29] to characterize the bifurcation of the 5D model (1). For the main control bifurcation parameter γ2, we obtain the auxiliary equation
(74)Q5λ=λ5+ζ1λ4+ζ2λ3+ζ3λ2+ζ4λ+ζ5=0,where ζ1=A1/A0; ζ2=A2/A0; ζ3=A3/A0; ζ4=A4/A0; and ζ5=A5/A0, and
(75)A0=−1.γ25−70.1351γ24−1844.6γ23−21561.9γ22−94515.4γ2,A1=564.448γ25+29688.γ24+520444.γ23+3.03997×106γ22−17409.7γ2−15909.2,A2=−1086.15γ25−57143.6γ24−1.00231×106γ23−5.86456×106γ22−56805.4γ2−2166.84,A3=200890.γ25+7.1036×106γ24+6.38216×107γ23+1.79792×107γ22−1.06833×106γ2+1199.88,A4=−42672.5γ25−1.54096×106γ24−1.4681×107γ23−1.36957×107γ22−20141.7γ2+38.9047,A5=−428.171γ25−22520.4γ24−394799.γ23−2.30622×106γ22+10655.1γ2−7.61994.Using the coefficient criteria of Hopf bifurcations without using eigenvalues and solving Equation (68) numerically for the bifurcation parameter γ2 in Wolfram Mathematica 11, we obtain γ2=γ2,c=0.00818058. Equation (69) gives ω=0.385339,λ1,2=±iω=±0.620757i, and T=2π/ω=16.3056 days. From Theorem 5.8, Equation (71), we obtain ζ1=20.4586, ζ3−ζ1ζ2=−71.7456, and ζ3−ζ1ω=0.251287. The plots of γ2 versus ζ1, ζ3−ζ1ζ2, and ζ5 are shown in Figures 6(a)–6(c), respectively. The coefficient criteria in Equation (71) are satisfied and therefore showing existence of a simple Hopf bifurcation at the critical value γ2,c of the bifurcation parameter γ2. The transversality condition is satisfied since dΨγ2/dγ2γ2=γ2,c=−64347.8 and Figure 6(d) show that for 0.005≤γ2≤0.01, the condition in Equation (72) of Theorem 5.9 is satisfied. The Hopf bifurcation is supercritical since as the bifurcation parameter γ2 is varied, the stable positive interior point loses its stability, and a stable limit cycle simultaneously appears.Figure 6
Plots for parameterγ2 versus (a) ζ1; (b) ζ3−ζ1ζ2; (c) ζ3−ζ1ω; and (d) dΨγ2/dγ2 for model system (1) with assumed initial condition: N0=10,P0=7,Z0=4,L0=2,R0=0.5 using the default parameter values in (73).
(a)(b)(c)(d)Bifurcation diagrams for model (1) were plotted using MATLAB R2016a code. ODE45 solver, which is based on an explicit Runge-Kutta ((4), (5)) formula and the Dormand-Prince method, was used for the numerical solution of the ordinary differential equations in (1). The bifurcation diagrams in Figures 7 and 8(a) show the change in stability for model (1) from a positive interior equilibrium into a limit cycle and period-doubling enroute to chaos for the control parameter, 0<γ2≤0.18. The zooplankton growth rate parameter γ2, therefore, has a tremendous effect on the dynamics of model system (1).Figure 7
Bifurcation diagram for parameterγ2 versus (a) nutrients; (b) phytoplankton; (c) zooplankton; and (d) Limnothrissa miodon for model system (1) with assumed initial condition: N0=10,P0=7,Z0=4,L0=2,R0=0.5 using the default parameter values in (73).
(a)(b)(c)(d)Figure 8
Bifurcation diagram for tigerfish versusγ2 of model (1) for 0<γ2≤0.18 and other default parameter values and with assumed initial condition: N0=10,P0=7,Z0=4,L0=2,R0=0.5.For the5D system in Equation (64), a local Hopf bifurcation occurs at x∗ηc,ηc if the Jacobian, JFx∗ηc,ηc, has a pair of imaginary roots. Using eigenvalues obtained from JFx∗ηc,ηc, we plot the real part (Re) of the eigenvalues λii=1,2,3,4,5. against the bifurcation parameter and find the bifurcation point where the curve crosses the axis of the bifurcation parameter as this is where the real part of an eigenvalue of an equilibrium passes through zero. The set of eigenvalues of JFx∗γ2,c,γ2,c are shown in the complex plane in Figure 9(a). The Hopf bifurcation occurs at the point where Reλ=0. Using the FindRoot command in Wolfram Mathematica 11, which searches for a numerical root starting from some initial root, the bifurcation value for the control parameter γ2 is 0.00818058 and is shown in Figure 9(b). Model (1) loses its stability whenever γ2>0.00818058 for the default set of parameter values. The coexistence equilibrium enters into a Hopf bifurcation at γ2=0.00818058. The eigenvalues at γ2,c are −20.2857,−0.0864343−0.0701176i,−0.0864343+0.0701176i,−0.620757i,+0.620757i; therefore, the stable coexistence equilibrium bifurcates into a stable periodic orbit.Figure 9
(a) Complex plane for0.002≤γ2≤0.18. (b) Plot of real part of the eigenvalue λ2 of the Jacobian matrix of model (1) for 0.002≤γ2≤0.02.
(a)(b)The Lyapunov spectrum of model system (1) consists of 5 eigenvalues (λj,j=1,2,3,4,5) called Lyapunov exponents. The rate of separation of infinitesimally close trajectories can be characterized by λj. The Lyapunov spectrum for model (1) is found using code in Wolfram Mathematica 11 and is shown in Figures 10 and 11.Figure 10
Lyapunov spectrum for model (1) for (a) γ2=0.055; (b) γ2=0.095 and other default parameter values and with assumed initial condition: N0=10,P0=7,Z0=4,L0=2,R0=0.5.
(a)(b)Figure 11
Lyapunov spectrum for model (1) for γ2=0.16 and other default parameter values and with assumed initial condition: N0=10,P0=7,Z0=4,L0=2,R0=0.5.The Lyapunov exponents for model (1) with γ2=0.055 and γ2=0.095 are (-0.0256886, -0.0426366, -0.0563857, -0.0809535, and -34.1637) and (0.000382075, -0.0373874, -0.0946365, -0.103307, and -20.1639) which correspond to the positive interior equilibrium and the limit cycle, respectively. The largest Lyapunov exponent for the model with γ2=0.055 is less than zero, indicating that nearby trajectories converge, and therefore, the interior equilibrium is a stable fixed point. For γ2=0.095, the maximal Lyapunov exponent is equal to zero, meaning that nearby trajectories converge to a closed orbit, and therefore, the periodic orbit (limit cycle) is stable.The Lyapunov exponents for model (1) with γ2=0.16 are (0.107838, -0.0039426, -0.0135559, -0.190155, and -5.87422). The maximal Lyapunov exponent is λ1=0.107838, and since λ1>0, this is an indication of divergence of nearby trajectories, and therefore, the dynamical system (1) is unstable and chaotic at γ2=0.16. For the dynamical system (1), a strange attractor exists in phase space at γ2=0.16.


## 6.1. Parameterσ4

There is no bifurcation for varyingσ4, but a Hopf bifurcation exists for varying σ4 with γ4=σ4. The predation rate of tigerfish on kapenta is varied from 0.4 to 0.9 with γ4=σ4. Figures 2 and 3(a) show the bifurcation diagrams for the bifurcation parameter σ4.Figure 2
Bifurcation diagram for parameterσ4 versus (a) nutrients; (b) phytoplankton; (c) zooplankton; and (d) Limnothrissa miodon for model system (1) with assumed initial condition: N0=10,P0=7,Z0=4,L0=2,R0=0.5 using the default parameter values.
(a)(b)(c)(d)Figure 3
(a) Bifurcation diagram for tigerfish versusσ4. (b) Plot of real part of eigenvalues of the Jacobian matrix of model (1) for 0.4≤σ4≤0.9.
(a)(b)The Hopf bifurcation value is0.63076 for the control parameter σ4 and is shown in Figure 3(b). From the coefficient criteria in Theorem 12, we obtain σ4=σ4,c=0.63076. Equation (69) gives ω=0.225435,λ1,2=±iω=±0.4748i, and T=2π/ω=27.8714 days. From Theorem 12, Equation (71), we obtain ζ1=18.8633, ζ3−ζ1ζ2=−79.8152 and ζ3−ζ1ω=0.37348. The coefficient criteria in Equation (71) are satisfied and therefore showing existence of a simple Hopf bifurcation at the critical value σ4,c of the bifurcation parameter σ4. The transversality condition in Equation (72) of Theorem 12 is satisfied since dΨσ4/dσ4σ4=σ4,c=−246.512. The eigenvalues at σ4,c are −18.6363,−0.113515−0.08458i,−0.113515+0.084586i,−0.4748i,+0.4748i. The Hopf bifurcation is supercritical since as the bifurcation parameter σ4 is varied with γ4=σ4 the stable positive interior point loses its stability and a stable periodic orbit simultaneously appears.Figure3 shows the effect of varying σ4 for 0.4≤σ4≤0.9 with σ4=γ4. The bifurcation results show that there is a decrease in kapenta population density for 0.4≤σ4<0.63076. Further increase of σ4 for 0.63076<σ4≤0.9 results in a period orbit of period 1 with increasing amplitude. Inefficient tigerfish predation on kapenta results in a decrease in the co-existence equilibrium value for kapenta, and efficient predation increases the chance of periodic behavior in the dynamical system. Therefore, the predation rate σ4, together with the tigerfish growth rate γ4, significantly determines the nature of the predator-prey relationship between tigerfish and kapenta. Simulation results show no chaotic behavior for the bifurcation control parameter σ4.


## 6.2. Parameterγ2

The interior equilibrium,E∗=0.602831,112.216,0.82694,1.6036,0.825832 for model (1), with damped oscillations for the set of default parameters are shown in Figure 4(a) for Limnothrissa miodon and 5(a) for zooplankton, Limnothrissa miodon, and tigerfish.Figure 4
Time series plot ofLimnothrissa miodon for (a) γ2=0.005; (b) γ2=0.0095; (c) γ2=0.075; and (d) γ2=0.16 for model (1) with assumed initial condition: N0=9.8,P0=6.8,Z0=3.8,L0=1.8,R0=0.3; N0=10,P0=7,Z0=4,L0=2,R0=0.5; N0=10.2,P0=7.2,Z0=4.2,L0=2.2,R0=0.7, using the default parameter values.
(a)(b)(c)(d)Figure 5
Phase portrait of zooplankton,Limnothrissa miodon and tigerfish (a) γ2=0.005, (b) γ2=0.0095, (c) γ2=0.075, and (d) γ2=0.16 for model (1) with assumed initial condition: N0=9.8,P0=6.8,Z0=3.8,L0=1.8,R0=0.3 (time series with blue color); N0=10,P0=7,Z0=4,L0=2,R0=0.5 (time series with magenta color); N0=10.2,P0=7.2,Z0=4.2,L0=2.2,R0=0.7 (time series with purple color), using the default parameter values.
(a)(b)(c)(d)The positive interior equilibriumE∗ for model system (1) is permanent for the set of default parameter values. For the main bifurcation parameter, model (1) bifurcates into a limit cycle of period 1 at γ2,c=0.00818058. Time series plots and phase portraits for γ2=0.0095 are shown in Figures 4(b) and 5(b). Time series plots and phase portraits for γ2=0.075 are shown in Figures 4(c) and 5(c). The phase portraits show periodic doubling. Time series plots and phase portraits for γ2=0.16 are shown in Figures 4(d) and 5(d). The phase portraits show a chaotic attractor at γ2=0.16. Varying the zooplankton growth rate γ2 results in some interesting dynamics for model (1). The stability of model (1) changes from a stable coexistence equilibrium into a stable periodic orbit and then periodic doubling enroute to chaos. Therefore, we can conclude that γ2 plays an important role in the dynamics of model (1). The trajectories in Figures 4 and 5 are aperiodic as they are showing erratic behavior. The irregular fluctuations can be attributed to the nonlinearity of the model (1). By visual inspection of the phase portraits in Figures 4(d) and 5(d) and from the magnitude of the maximal Lyapunov exponent, 0.107838, it can be seen that the trajectories are sensitive to initial conditions and, therefore, it is impossible to predict the long-term behavior of model system (1) for γ2=0.16. The trajectories are random and bounded in phase space and they converge to a strange attractor, which has a complex (fractal) geometric structure.We use the Hopf bifurcation theorem by Guckenheimer and Holmes [28] together with conditions by Douskos and Markellos [29] to characterize the bifurcation of the 5D model (1). For the main control bifurcation parameter γ2, we obtain the auxiliary equation
(74)Q5λ=λ5+ζ1λ4+ζ2λ3+ζ3λ2+ζ4λ+ζ5=0,where ζ1=A1/A0; ζ2=A2/A0; ζ3=A3/A0; ζ4=A4/A0; and ζ5=A5/A0, and
(75)A0=−1.γ25−70.1351γ24−1844.6γ23−21561.9γ22−94515.4γ2,A1=564.448γ25+29688.γ24+520444.γ23+3.03997×106γ22−17409.7γ2−15909.2,A2=−1086.15γ25−57143.6γ24−1.00231×106γ23−5.86456×106γ22−56805.4γ2−2166.84,A3=200890.γ25+7.1036×106γ24+6.38216×107γ23+1.79792×107γ22−1.06833×106γ2+1199.88,A4=−42672.5γ25−1.54096×106γ24−1.4681×107γ23−1.36957×107γ22−20141.7γ2+38.9047,A5=−428.171γ25−22520.4γ24−394799.γ23−2.30622×106γ22+10655.1γ2−7.61994.Using the coefficient criteria of Hopf bifurcations without using eigenvalues and solving Equation (68) numerically for the bifurcation parameter γ2 in Wolfram Mathematica 11, we obtain γ2=γ2,c=0.00818058. Equation (69) gives ω=0.385339,λ1,2=±iω=±0.620757i, and T=2π/ω=16.3056 days. From Theorem 5.8, Equation (71), we obtain ζ1=20.4586, ζ3−ζ1ζ2=−71.7456, and ζ3−ζ1ω=0.251287. The plots of γ2 versus ζ1, ζ3−ζ1ζ2, and ζ5 are shown in Figures 6(a)–6(c), respectively. The coefficient criteria in Equation (71) are satisfied and therefore showing existence of a simple Hopf bifurcation at the critical value γ2,c of the bifurcation parameter γ2. The transversality condition is satisfied since dΨγ2/dγ2γ2=γ2,c=−64347.8 and Figure 6(d) show that for 0.005≤γ2≤0.01, the condition in Equation (72) of Theorem 5.9 is satisfied. The Hopf bifurcation is supercritical since as the bifurcation parameter γ2 is varied, the stable positive interior point loses its stability, and a stable limit cycle simultaneously appears.Figure 6
Plots for parameterγ2 versus (a) ζ1; (b) ζ3−ζ1ζ2; (c) ζ3−ζ1ω; and (d) dΨγ2/dγ2 for model system (1) with assumed initial condition: N0=10,P0=7,Z0=4,L0=2,R0=0.5 using the default parameter values in (73).
(a)(b)(c)(d)Bifurcation diagrams for model (1) were plotted using MATLAB R2016a code. ODE45 solver, which is based on an explicit Runge-Kutta ((4), (5)) formula and the Dormand-Prince method, was used for the numerical solution of the ordinary differential equations in (1). The bifurcation diagrams in Figures 7 and 8(a) show the change in stability for model (1) from a positive interior equilibrium into a limit cycle and period-doubling enroute to chaos for the control parameter, 0<γ2≤0.18. The zooplankton growth rate parameter γ2, therefore, has a tremendous effect on the dynamics of model system (1).Figure 7
Bifurcation diagram for parameterγ2 versus (a) nutrients; (b) phytoplankton; (c) zooplankton; and (d) Limnothrissa miodon for model system (1) with assumed initial condition: N0=10,P0=7,Z0=4,L0=2,R0=0.5 using the default parameter values in (73).
(a)(b)(c)(d)Figure 8
Bifurcation diagram for tigerfish versusγ2 of model (1) for 0<γ2≤0.18 and other default parameter values and with assumed initial condition: N0=10,P0=7,Z0=4,L0=2,R0=0.5.For the5D system in Equation (64), a local Hopf bifurcation occurs at x∗ηc,ηc if the Jacobian, JFx∗ηc,ηc, has a pair of imaginary roots. Using eigenvalues obtained from JFx∗ηc,ηc, we plot the real part (Re) of the eigenvalues λii=1,2,3,4,5. against the bifurcation parameter and find the bifurcation point where the curve crosses the axis of the bifurcation parameter as this is where the real part of an eigenvalue of an equilibrium passes through zero. The set of eigenvalues of JFx∗γ2,c,γ2,c are shown in the complex plane in Figure 9(a). The Hopf bifurcation occurs at the point where Reλ=0. Using the FindRoot command in Wolfram Mathematica 11, which searches for a numerical root starting from some initial root, the bifurcation value for the control parameter γ2 is 0.00818058 and is shown in Figure 9(b). Model (1) loses its stability whenever γ2>0.00818058 for the default set of parameter values. The coexistence equilibrium enters into a Hopf bifurcation at γ2=0.00818058. The eigenvalues at γ2,c are −20.2857,−0.0864343−0.0701176i,−0.0864343+0.0701176i,−0.620757i,+0.620757i; therefore, the stable coexistence equilibrium bifurcates into a stable periodic orbit.Figure 9
(a) Complex plane for0.002≤γ2≤0.18. (b) Plot of real part of the eigenvalue λ2 of the Jacobian matrix of model (1) for 0.002≤γ2≤0.02.
(a)(b)The Lyapunov spectrum of model system (1) consists of 5 eigenvalues (λj,j=1,2,3,4,5) called Lyapunov exponents. The rate of separation of infinitesimally close trajectories can be characterized by λj. The Lyapunov spectrum for model (1) is found using code in Wolfram Mathematica 11 and is shown in Figures 10 and 11.Figure 10
Lyapunov spectrum for model (1) for (a) γ2=0.055; (b) γ2=0.095 and other default parameter values and with assumed initial condition: N0=10,P0=7,Z0=4,L0=2,R0=0.5.
(a)(b)Figure 11
Lyapunov spectrum for model (1) for γ2=0.16 and other default parameter values and with assumed initial condition: N0=10,P0=7,Z0=4,L0=2,R0=0.5.The Lyapunov exponents for model (1) with γ2=0.055 and γ2=0.095 are (-0.0256886, -0.0426366, -0.0563857, -0.0809535, and -34.1637) and (0.000382075, -0.0373874, -0.0946365, -0.103307, and -20.1639) which correspond to the positive interior equilibrium and the limit cycle, respectively. The largest Lyapunov exponent for the model with γ2=0.055 is less than zero, indicating that nearby trajectories converge, and therefore, the interior equilibrium is a stable fixed point. For γ2=0.095, the maximal Lyapunov exponent is equal to zero, meaning that nearby trajectories converge to a closed orbit, and therefore, the periodic orbit (limit cycle) is stable.The Lyapunov exponents for model (1) with γ2=0.16 are (0.107838, -0.0039426, -0.0135559, -0.190155, and -5.87422). The maximal Lyapunov exponent is λ1=0.107838, and since λ1>0, this is an indication of divergence of nearby trajectories, and therefore, the dynamical system (1) is unstable and chaotic at γ2=0.16. For the dynamical system (1), a strange attractor exists in phase space at γ2=0.16.


## 7. Conclusions

In this paper, we formulated and analyzed aLimnothrissa miodon model with Hydrocynus vittatus predation. Positivity and existence of solutions for the model were shown. Local stability analysis results agree with the numerical simulations in that the coexistence equilibrium is locally stable provided that certain conditions are satisfied. The coexistence equilibrium is globally stable if certain conditions are met. In bifurcation analysis, the Hopf bifurcation theorem together with certain conditions for a 5D system was used to find the bifurcation point, angular frequency, period, pair of imaginary eigenvalues, stability, and direction for a bifurcation control parameter. The eigenvalue method was also used to find the bifurcation point of a bifurcation control parameter. Results from the Hopf bifurcation theorem and from the eigenvalue method for the 5D system were the same. Therefore, either method can be used to find the bifurcation point for a given control parameter for the 5D system. The model changes in stability from a positive interior equilibrium to a limit cycle and period-doubling enroute to chaos for the zooplankton growth rate control parameter, indicating that zooplankton growth rate has a tremendous effect on the dynamics of the model. Therefore, a variation in the zooplankton growth rate can significantly alter the dynamics of Limnothrissa miodon in Lake Kariba. The Lyapunov spectrum of exponents was used to show the convergence or divergence of nearby trajectories and to check for the existence of chaos in the aquatic environment. A supercritical simple Hopf bifurcation exists for varying the biological Hydrocynus vittatus predation on Limnothrissa miodon parameter only if the parameter is equal to the Hydrocynus vittatus growth rate coefficient. Hydrocynus vittatus predation on Limnothrissa miodon changes a stable coexistence equilibrium into a stable periodic orbit if the predation rate coefficient is the same as the Hydrocynus vittatus growth rate coefficient. We therefore conclude that Hydrocynus vittatus predation on Limnothrissa miodon significantly alters the dynamics of Limnothrissa miodon whenever the Hydrocynus vittatus predation rate coefficient is the same as the Hydrocynus vittatus growth rate coefficient. Hopf bifurcation for varying a predation rate coefficient and the Hopf bifurcation leading to chaos for the zooplankton growth rate control parameter is in agreement with findings from other authors. The periodic orbit obtained in bifurcation analysis reflects what actually happens in Lake Kariba as evidenced from actual data shown in Figure 5.1, which shows the cyclical behavior of the Hydrocynus vittatus bycatch and Limnothrissa miodon catch. Therefore, the Hopf bifurcation results reflect what really happens in the Lake Kariba kapenta fishery. Bifurcation results show that a higher zooplankton growth rate, which implies efficient grazing on phytoplankton, increases the chance of chaos in the dynamical system. Inefficient grazing of phytoplankton by zooplankton implies a lower zooplankton growth rate resulting in simple dynamics in the model system. Therefore, the phytoplankton-zooplankton oscillations and the nature of the zooplankton predation play an important role in model dynamics. An increase in the parameter γ2 leads to destabilization of the dynamical system (1). The chaotic solutions for model system (1) for the bifurcation parameter γ2 could not be validated with actual data for nutrients, phytoplankton, zooplankton, kapenta, and tigerfish from the Lake Kariba kapenta fishery, since the sample size of the data available is small and the data was collected at irregular time intervals. Therefore, the unavailability of data is a limitation in validating the deterministic chaos results for the bifurcation control parameter γ2. In certain ecological systems, chaotic dynamics are expected to contribute to the unpredictability and irregularity of ecological time series [30]. It is debatable whether or not this chaotic behavior observed in Figures 4(d) and 5(d) occurs in the kapenta fishery. The fundamental issue stems from the fact that in most ecological systems, there is a strong stochastic disturbance from environmental factors such as temperature and weather, and this makes determining whether the irregular structure in the data is due to chaotic dynamics or stochastic perturbations difficult [30]. Phytoplankton dynamics usually show erratic and eruptive “busts and blooms” and have similar characteristics of deterministic chaos [31]. According to Stone and Ezrati [32], it is reasonable to have oscillating and chaotic dynamics in nonlinear deterministic ecological systems with growth processes that are strong and without neglecting the possibility of stochastic processes influencing the variability which arises in nature. For future studies, we intend to model the dynamics of Limnothrissa miodon with lake water temperature.



---

*Source: 1095441-2022-07-20.xml*