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

In this paper, we construct and analyze a theoretical, deterministic 5 D 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 satis ﬁ ed. 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.


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 NðtÞ, PðtÞ, ZðtÞ, LðtÞ, and RðtÞ, 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 the Limnothrissa 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][6][7][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 Section 2; 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.

Materials and Methods
2.1. Study Area. Lake Kariba (277 km long; about 5364 km 2 in surface area; 160 km 3 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 km 2 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].

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.
The predator-prey relationship between the tigerfish and kapenta is shown in Figure 1. 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.

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 NðtÞ, PðtÞ, ZðtÞLðtÞ, and RðtÞ, 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 σ 1 NP: The growth rate of phytoplankton is γ 1 NP: It is assumed that the depletion rate of phytoplankton caused by mortality is proportional to P: Phytoplankton is depleted by zooplankton at a rate σ 2 PZ: The depletion of phytoplankton per unit time by zooplankton is given by σ 2 PZ 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 γ 2 PZ: It is assumed that the depletion rate of zooplankton caused by mortality is proportional to Z: Journal of Applied Mathematics The functional response of zooplankton to the Limnothrissa miodon given by σ 3 ZL 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 γ 3 ZL: 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 L 2 : 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 f ðLÞ = ðσ 4 L/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 γ 4 LR/d + L: The tigerfish are harvested at a rate κηR, where κ is the catchability coefficient and η is the effort. The nonlinear dynamical system is and define 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.
Hence, we obtain From the second equation of model (1), it follows that Direct integration of (6) results in From the third equation of model (1), it follows that Direct integration of (8) results in Considering the variable LðtÞ in ½0, T, from the fourth equation of model (1), it follows that Direct integration of (10) results in From the fifth equation of model (1), it follows that Direct integration of (12) results in Therefore, the solutions of system (1) with initial condition (2) remains positive ∀t ≥ 0:

Equilibria and Stability Analysis
Model (1) has 6 equilibria: (a) The trivial equilibrium is (b) The phytoplankton free equilibrium is (c) The zooplankton free equilibrium is E 2 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. E 2 exists provided that Rearranging the inequality (22), we obtain ða/μ 0 Þ > ðμ 1 / γ 1 Þ, and this means that N * 1 > N * 2 . 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) The Limnothrissa miodon free equilibrium is The zooplankton population is insufficient to support the population of Limnothrissa miodon. When the Limnothrissa miodon population is not present in the ecosystem and both phytoplankton and zooplankton are present, In order to support the zooplankton population, more nutrients are needed in the ecosystem. From (23), ð24Þ E 3 exists on condition that Inequality (25) can be rearranged to give P * 2 > P * 3 , meaning that the phytoplankton equilibrium is reduced in the presence of zooplankton.
Journal of Applied Mathematics Theorem 5. The phytoplankton free equilibrium is locally asymptotically stable if aγ 1 < μ 0 μ 1 .

Theorem 11.
The equilibrium E * 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. Let
Let VðN, P, Z, L, RÞ be a positive Lyapunov function [26,27], such that VðN * , P * , Z * , L * , R * Þ = 0 by, Journal of Applied Mathematics 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 Then 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.

Bifurcation Analysis
The 5D system (1) is written as 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 (2) The auxiliary equation has complex conjugate roots (nonhyperbolicity condition) and if the condition with is satisfied, and there are no other roots with zero real parts and λ i ði = 3, 4, 5:Þ have nonzero real parts if and have negative real parts if is satisfied, and then, there is a Poincare ′ -Andronov-Hopf bifurcation.
The Hopf bifurcation value is 0:63076 for the control parameter σ 4 and is shown in Figure 3  11 Journal of Applied Mathematics −0:113515 + 0:084586i,−0:4748i,+0:4748ig. 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. Figure 3 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 . 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  13 Journal of Applied Mathematics 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 14 Journal of Applied Mathematics where 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. 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). For the 5D system in Equation (64), a local Hopf bifurcation occurs at ðx * ðη c Þ, η c Þ if the Jacobian, J Fðx * ðη c Þ,η c Þ , has a pair of imaginary roots. Using eigenvalues obtained from J Fðx * ðη c Þ,η c Þ , we plot the real part (Re) of the eigenvalues ðλ i ði = 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 J Fðx * ðγ 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 f−20:2857,− 0:0864343 − 0:0701176i,−0:0864343 + 0:0701176i,− 0:620757i,+0:620757ig; therefore, the stable coexistence equilibrium bifurcates into a stable periodic orbit.
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.
The Lyapunov exponents for model (1)   17 Journal of Applied Mathematics 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:

Conclusions
In this paper, we formulated and analyzed a Limnothrissa miodon model with Hydrocynus vittatus predation. Positiv-ity 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.

Data Availability
The data used to support the findings in this study are included in the article.

Conflicts of Interest
There are no conflicts of interest for the authors of this study.