# Stability Analysis and Optimal Control of a Limnothrissa Miodon Model with Harvesting

Authors: Farikayi K. Mutasa; Brian Jones; Senelani D. Hove-Musekwa; Itai H. Tendaupenyu; Tamuka Nhiwatiwa; Mzime R. Ndebele-Murisa
Journal: Discrete Dynamics in Nature and Society
 (2022)

---

## Abstract

We construct a theoretical, deterministic mathematical model of the dynamics of Limnothrissa miodon with nutrients, phytoplankton, and zooplankton and investigate the effect of harvesting on the population density of Limnothrissa miodon in a lake. For the autonomous model, results from local stability analysis are in agreement with numerical simulations in that the coexistence equilibrium is locally stable, provided certain conditions are satisfied. The coexistence equilibrium is globally unstable if it is feasible. Numerical results show that a stable limit cycle exists for the nonautonomous model. Optimal control results show an optimal harvesting monthly effort of 15394 boat nights which corresponds to 505 fishing units, showing that there is overcapacity in Lake Kariba. A maximum sustainable annual catch of 34669 tonnes is obtained and simulation results show that Limnothrissa miodon abundance is more closely related to nutrient inflow than to harvesting.


---

## Body

## 1. Introduction

Limnothrissa miodon (Boulenger, 1906), also referred to as kapenta, is a major source of protein and income to fishing cooperatives, wholesalers, retailers, and the local community [1]. Management of the Lake Kariba fishery which is shared by both Zimbabwe and Zambia is of paramount importance for the continued survival of the kapenta [2, 3]. Harvesting of kapenta by fishing vessels plays an important role in the dynamics of the sardine Limnothrissa miodon [2, 3] and as a result, it is critical to look into the role of harvesting in its dynamics theoretically and numerically. Harvesting of the natural resource needs to be done in a sustainable manner so that it is not overexploited and becomes extinct. Therefore, mathematical modelling of the Limnothrissa miodon model with harvesting will provide insight into the kapenta fishery dynamics in Lake Kariba. A deterministic model involving nutrients, plankton, and Limnothrissa miodon, as well as harvesting by fishing vessels, has yet to be developed and analysed. We present and analyze a deterministic continuous dynamical system composed of ordinary differential equations that describe the dynamics of Limnothrissa miodon in the presence of nutrients, plankton, and harvesting. The Limnothrissa miodon model will aid us in our understanding of the dynamics of the aquatic ecosystem in the kapenta fishery in Lake Kariba.According to McLachlan [4], the lake productivity is largely due to reloading of nutrients at turnover, the nutrient-rich sediment carried by the inflowing rivers, and the nutrients discharged as runoff into the lake. The lake loses large amounts of nutrients per year since it experiences large outflows of between 50 and 60 km3 as compared to its volume of 160 km3 [5]. Marshall [6] estimated the total mortality of Limnothrissa miodon to be 0.816 month−1 and 1.149 month−1 in 1978 and 1983, respectively. Moreau et al. [7] estimated the total mortality of Limnothrissa miodon to be between 0.37 month−1 and 0.81 month−1 in Lake Tanganyika. For Lake Kariba and Lake Tanganyika, the natural mortality rates for Limnothrissa miodon were 0.44 and 0.30 month−1, respectively, and Pauly’s empirical method was used to estimate the parameters. According to Marshall [6], the fishing effort to mortality relationship of Limnothrissa miodon is 0.731 month−1. From the biomass assessment done in 1994/95, the fishing mortality was estimated to be 2.06 year−1 and 2.18 year−1 for 1993 and 1994, respectively, for kapenta [8]. The Working Group Assessment [9] obtained a fishing mortality of 2.1 year−1, a natural mortality of 3.5 year−1, giving a total mortality of 5.6 year−1 for kapenta in Lake Kariba. Cochrane [10] and Marshall [11] concluded that catches of Limnothrissa miodon vary seasonally with availability of food. The fish are more abundant after turnover in August to September and after the annual flooding of the inflowing rivers in April and May. Marshall [11] also concluded the existence of a possible relationship between the sardine biomass and fishing effort. In their study on the impact of fishing pressure on kapenta production on the Zambian side of Lake Kariba, Chali et al. [2] concluded that both catch per unit effort (CPUE) and kapenta production were declining significantly. Kapenta production and CPUE decreased from 9993 tonnes (t) and 0.145 t/boat/night to 6004 t and 0.085 t/boat/night from 2009 to 2012, respectively. The fishing effort increased from about 68734 to 70706 boat nights from 2009 to 2012, respectively. CPUE is an indirect measure of kapenta abundance in the lake and so a decreasing CPUE shows overexploitation of the sardine. On the Zambian side of the lake, the fishing units increased from 423 in 2009 to 800 in 2013 and this could have led to the decline in the kapenta stocks and overexploitation of the sardine. The low kapenta catches were as a result of a huge rise in the fishing effort which was more than the 500 recommended rigs for the whole lake [12].Machena et al. [13] from their study on the preliminary assessment of the trophic structure of Lake Kariba concluded that the fishing mortality of Limnothrissa miodon could be increased since it is fairly harvested. Magadza [14] in his study estimated a maximum sustainable yield (MSY) of kapenta to be 40000 t and concluded that the catches of kapenta in the fishery increased logistically from 1974 to about 1990, where the highest catch of 37000 t was recorded and then decreased thereafter. Madamombe [15] for the period 1974 to 1999 used the Pella and Tomlinson production model and obtained a maximum sustainable yield of 23336 t with a corresponding effort of 725 rigs. Mudenda [16] estimated MSY with the Schaeffer and Fox models and obtained 23525 t and 24271 t with a corresponding fishing effort of 108109 and 141243 boat nights, respectively. Mudenda [16] concluded that there was overfishing in the fishery ranging from 21%−39%. Tendaupenyu and Pyo [17] for the study period 1988−2009 used a maximum entropy model and estimated a MSY of 25372 t with a corresponding fishing effort of 109731 boat nights and they concluded that there was overcapacity in Lake Kariba and that the kapenta stock was in decline. The analytical model used by Tendaupenyu and Pyo [17] estimated a fishing mortality of 1.21 year−1 compared to the current 0.927 year−1. Tendaupenyu and Pyo [17] used a maximum entropy model and estimated an intrinsic growth rate of 0.42 and a carrying capacity of 240500 t for kapenta in Lake Kariba. Tendaupenyu and Pyo [17] estimated a catchability coefficient of 0.00000153 for kapenta from the Walters and Hilborn model. Dynamical systems have not been used to investigate how kapenta fish populations in Lake Kariba are described and influenced by harvesting. We will be able to quantitatively explain the impact of kapenta fish harvesting by developing and analyzing a mathematical model.The important question is what levels of the fishing effort will lead to a sustained population density of Limnothrissa miodon in the lake and what fishing pressure will result in the extinction of the natural resource. Many fish harvesting models have been formulated and analysed [18–20] and they mainly differ in the harvesting function.The remainder of this paper describes the study area, data collection, model formulation, positivity and existence of solutions, equilibrium states and their stability, and optimal control in Section2. The results which include numerical simulations and analysis are presented in Section 3, a detailed discussion is in Section 4, and conclusions are in Section 5.


## 2. Materials and Methods

### 2.1. Study Area

Lake Kariba is located on the Zambezi River, in a tropical area with seasonal rainfall between latitudes 1628′ to 1804′ S and longitudes 2642′ to 2903′E [15]. The offshore single-species pelagic kapenta fishery is highly mechanized, using light attraction and lift nets from boats, and is licensed [21]. According to the 2011 Frame survey [22], the fishing units on the Zimbabwean and Zambian sides of Lake Kariba have similar technical efficiencies and harvesting systems. Large circular dip nets and mercury light bulbs are used in the fishing gear to attract the kapenta. The lights used below and above the lake’s surface are powered by generators, the number of which varies from rig to rig. To lower or raise the dip net, hydraulic, mechanical, or manual winches are used [8].


### 2.2. Data Collection

The Lake Kariba Fisheries Research Institute collects data on catch, effort in the experimental gillnet, and 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 [8]. The University of Zimbabwe Lake Kariba Research Institute collects data on the lake’s water quality. The data on acoustic surveys of the biomass of kapenta are not often collected because of the lack of equipment and the cost involved and as such a few fishery surveys have been done on Lake Kariba using research vessels. Therefore, data on the biomass of the aquatic ecosystem of Lake Kariba are only limited to the hydroacoustic surveys of 1981, 1982, and 1983 [23]; September 1988 [24]; 1992 [25]; and for the period from 14 to 24 August 2014 [26]. Some data on phytoplankton, zooplankton, and kapenta in Lake Kariba were obtained from Ndebele-Murisa [21]. Some data is missing due to the fact that some of the fish caught in the offshore fishery are sold to buyers by the fishers before they come back to land and therefore are not recorded. As much as 30% of kapenta catch is thought to be sold illegally by boat crews before landing without the knowledge of the owners of the fishing vessels [27].Figure1(a) shows the time series plot of total nitrogen and chlorophyll a in μgl−1 in Lake Kariba from November 2016 to December 2017. The monthly averages of total nitrogen and chlorophyll a in Lake Kariba from April 2014 to December 2017 are illustrated in Figure 1(b). Figures 1(c) and 1(d) show the time series plot from November 2016 to December 2017 and the monthly averages from April 2014 to December 2017 of zooplankton in Lake Kariba in individuals per litre (ind./l), respectively.Figure 1
(a) Time series plot of total nitrogen and chlorophyll a; (b) a plot of monthly averages of total nitrogen and chlorophyll a; (c) time series plot of zooplankton; (d) a plot of monthly averages of zooplankton.
(a)(b)(c)(d)Figure2(a) shows the number of fishing vessels on the Zimbabwean side of Lake Kariba from 2008 to 2018. The time series plot shows a gradual increase in the number of units from 307 in 2008 to 530 in 2018. Figure 2(b) shows the number of fishing units on the Zambian and Zimbabwean sides of the lake from 1994 to 1999.Figure 2
(a) Time series plot of fishing units from 2008 to 2018. (b) Time series plot of fishing units from 1994 to 1999 [27].
(a)(b)Figure3(a) illustrates the time series plot of Limnothrissa miodon catches and the CPUE from November 2016 to December 2017. The catch and CPUE monthly averages for the period 1974 to 2018 are shown in Figure 3(b). Figures 3(a) and 3(b) are for the Zimbabwean side of the lake. Figure 3(b) shows a small peak in the catches of kapenta around March and a big peak around August. The peak in March is attributed to the inflow of nutrients through the inflowing rivers and the flooding of the Zambezi river and the one in August is due to the reloading of nutrients after the mixing of the lake around July.Figure 3
(a) Kapenta catch and CPUE time series plot; (b) a plot of monthly averages of kapenta catch and CPUE; (c) kapenta catch in Lake Kariba time series plot; (d) a plot of effort and CPUE in Lake Kariba.
(a)(b)(c)(d)The catches of kapenta versus effort are shown in Figure3(c). Figure 3(d) shows the effort, catch, and the CPUE of kapenta from 1974 to 2011. Figures 3(c) and 3(d) show data from both the Zambian and Zimbabwean sides of Lake Kariba. For the period 1974 to 2011, the global maximum is at an effort of 104131 boat nights giving a catch of 30943 t in 1990.


### 2.3. Model Equations

The model has four classes:N denotes nutrient concentration, P denotes phytoplankton population density, Z denotes zooplankton population density, and L denotes Limnothrissa miodon population density. Densities in each class are time-dependent and are denoted by Nt, Pt, Zt, and Lt, respectively. The Limnothrissa miodon model [28] is extended to include harvesting by the fishing vessels. Nutrients are assumed to enter the water body sinusoidally at the rate(1)a1+bsinωt+φ1+csin2ωt+φ2,where a,b,c>0 are constants, φ1 and φ2 are phase angles, and ω is the angular frequency, a constant. φ2 captures the delay in the lake mixing, ω is the frequency of the periodic input of nutrients, and ab and ac are the overall amplitude of the first and second nutrient periodic inputs, respectively. The first term of (1) represents the inflow of nutrients from river inflow and runoff and the second term of (1) captures the nutrient reloading into the lake which occurs around mid-year. The sinusoidal function in (1) is chosen for its simplicity and its parameters are to be determined by fitting. Natural depletion of nutrients occurs at a constant rate μ0. Phytoplankton depletes nutrients at a rate of σ1NP. The rate of growth of phytoplankton is ϕ1σ1NP. It is assumed that the depletion rate of phytoplankton caused by mortality is proportional to P. The change in density of the phytoplankton per unit time per zooplankton as the phytoplankton population density changes is of the modified Holling’s type-I response [29] and is denoted by σ2PZ [29]. The conversion coefficient from phytoplankton to zooplankton is ϕ2. It is assumed that the rate of zooplankton depletion 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 zooplankton density per unit time per Limnothrissa miodon as the zooplankton population density changes. The conversion coefficient from zooplankton to Limnothrissa miodon is phi ϕ3. It is assumed that the rate of Limnothrissa miodon depletion caused by mortality is proportional to L, and the rate of depletion caused by crowding is proportional to L2. Kapenta is harvested at a rate Gt. The model structure improves on the Limnothrissa miodon model in [28] and is shown in Figure 4.Figure 4
Flow diagram of the Limnothrissa miodon model with harvesting.We assume the following:(1)
Catch per unit effort of kapenta is proportional to the fish abundance,Gt/E∝Lt, and can be written as Gt=qELt [18, 19], where q≥0 is a fixed constant of proportionality called the catchability coefficient and E=Et≥0 is the fishing effort per unit time.(2)
Pelagic kapenta have schooling behaviour [18].(3)
Fishing using the rigs involves randomly searching for fish.(4)
Every kapenta in the body of water has an equal chance of being captured.The nonlinear nonautonomous dynamical system with Schaefer harvesting is(2)dNdt=a1+bsinωt+φ1+csin2ωt+φ2−μ0N−σ1NP,dPdt=ϕ1σ1NP−μ1P−σ2PZ,dZdt=ϕ2σ2PZ−μ2Z−σ3ZL,dLdt=ϕ3σ3ZL−μ3L−σ30L2−qEL..Model (2) initial condition is,(3)N0=ψ10,P0=ψ20,Z0=ψ30,L0=ψ40,ψi0>0,i=1,2,3,4,and the mathematically feasible region is defined as(4)Ω=N,P,Z,L∈ℝ4|N≥0,P≥0,Z≥0,L≥0.The crowding of the Limnothrissa miodon population coefficient isσ30>0. The constants of proportionality σii=1,2,3 are positive and depletion rate coefficients are represented by μii=0,1,2,3.


### 2.4. Positivity of Solutions

We prove that, for positive initial data for the model (2), the solutions will remain positive ∀t≥0.Theorem 1.
LetNt≥0,Pt≥0,Zt≥0,Lt≥0 be the initial data. Then, solutions of Nt, Pt, Zt, Lt of model (2) are positive ∀t≥0.Proof.
Based on the first equation of system (2) and Nt∈0,T, it follows that(5)N˙t≥−μ0Nt−σ1NtPt,∀t∈0,T.
Consequently, we obtain(6)Nt≥N0exp∫0t−μ0−σ1Psds≥0,∀t∈0,T.
According to the second equation of system (2),(7)P˙t≥−μ1Pt−σ2PtZt,∀t∈0,T.
The direct integration of (7) yields(8)Pt≥P0exp∫0t−μ1−σ2Zsds≥0,∀t∈0,T.
According to the third equation of system (2),(9)Z˙t≥−μ2Zt−σ3ZtLt,∀t∈0,T.
The direct integration of (9) yields(10)Zt≥Z0exp∫0t−μ2−σ3Lsds≥0,∀t∈0,T.
Considering the variableLt∈0,T, it follows from the fourth equation of system (2) that(11)L˙t≥−Ltμ3+qE+σ30Lt,∀t∈0,T.
The direct integration of (11) yields(12)Lt≥μ3+qEL0e−μ3+qEtμ3+qE+σ30L01−e−μ3+qEt≥0,∀t∈0,T.
As a result, system solutions of (2) with initial conditions (3) remain positive ∀t≥0.


### 2.5. Existence of Solutions

Theorem 2.
A solution of system (2) is feasible.Proof.
We show that all feasible solutions of system (2) are uniformly bounded in Ω⊂ℝ4. Let any solution of system (2) be Nt,Pt,Zt,Lt∈ℝ4, with nonnegative initial conditions.
LetCt=Nt+Pt+Zt+Lt, and then(13)dCdt=a1+bsinωt+φ1+csin2ωt+φ2−μ0N−σ1NP+ϕ1σ1NP−μ1P−σ2PZ+ϕ2σ2PZ−μ2Z−σ3ZL+ϕ3σ3ZL−μ3L−σ30L2−qEL,=a1+bsinωt+φ1+csin2ωt+φ2−μ0N−μ1P−μ2Z−μ3+qEL−σ30L2+σ1ϕ1−1NP+σ2ϕ2−1PZ+σ3ϕ3−1ZL,≤a1+bsinωt+φ1+csin2ωt+φ2−μ0N−μ1P−μ2Z−μ3+qEL,≤a1+bsinwt+φ1+csin2wt+φ2−mCt,where m=minμ0,μ1,μ2,μ3+qE. Thus,(14)dCtdt+mCt≤a1+bsinωt+φ1+csin2ωt+φ2.
Equation (14) is a first-order differential inequality, with the solution given by(15)0<CN,P,Z,L≤abm4sintω+φ1+4bm2ω2sintω+φ1+cm4sin2tω+φ2+cm2ω2sin2tω+φ2m5+5m3ω2+4mω4−abmωm2+4ω2costω+φ1+2cmωm2+ω2cos2tω+φ2+m4+5m2ω2+4ω4m5+5m3ω2+4mω4+C0e−mt,as t⟶∞, (15) becomes(16)0<CN,P,Z,L≤abm4sintω+φ1+4bm2ω2sintω+φ1+cm4sin2tω+φ2+cm2ω2sin2tω+φ2m5+5m3ω2+4mω4−abmωm2+4ω2costω+φ1+2cmωm2+ω2cos2tω+φ2+m4+5m2ω2+4ω4m5+5m3ω2+4mω4.
Therefore, all solutions of system (2) enter the feasible region,(17)Ω=Nt,Pt,Zt,Lt∈ℝ+4:C≤abm4sintω+φ1+4bm2ω2sintω+φ1+cm4sin2tω+φ2+cm2ω2sin2tω+φ2m5+5m3ω2+4mω4−abmωm2+4ω2costω+φ1+2cmωm2+ω2cos2tω+φ2+m4+5m2ω2+4ω4m5+5m3ω2+4mω4+ς,∀ς>0.
This concludes the theorem’s proof.


### 2.6. Equilibrium States for the Autonomous Model

Of biological interest is the coexistence equilibrium. Ifb=c=0 in system (2), then the model is autonomous and has a unique coexistence equilibrium, E∗=N∗,P∗,Z∗,L∗, which is obtained through the solution of the equations:(18)a−μ0N−σ1NP=0,(19)ϕ1σ1N−μ1−σ2Z=0,(20)ϕ2σ2P−μ2−σ3L=0,(21)ϕ3σ3Z−μ3−qE−σ30L=0.Solving forN, P, Z, and L in (18)–(21) gives(22)σ1σ2σ3σ30L∗2+σ1σ2σ3μ3+qE+σ1σ2μ2σ30+μ0σ22ϕ2σ30+μ1σ1σ32ϕ3L∗+μ0μ1ϕ2σ2ϕ3σ3+μ1σ1μ2ϕ3σ3+μ0σ22ϕ2μ3+qE+σ1σ2μ2μ3+qE−aϕ1σ1ϕ2σ2ϕ3σ3=0.If expression (24) is positive, (23) will have a unique positive root:(23)L∗=A12−4σ1σ2σ3σ30A2−A12σ1σ2σ3σ30>0,where(24)A1=σ1σ2σ3μ3+qE+σ1σ2μ2σ30+μ0σ22ϕ2σ30+μ1σ1σ32ϕ3,A2=μ0μ1ϕ2σ2ϕ3σ3+μ1σ1μ2ϕ3σ3+μ0σ22ϕ2μ3+qE+σ1σ2μ2μ3+qE−aϕ1σ1ϕ2σ2ϕ3σ3.And (23) can be written as(25)A12−4σ1σ2σ3σ30A2>A12,σ1σ2σ3σ30A2<0,σ1σ2σ3σ30μ0μ1ϕ2σ2ϕ3σ3+μ1σ1μ2ϕ3σ3+μ0σ22ϕ2μ3+qE+σ1σ2μ2μ3+qE−aϕ1σ1ϕ2σ2ϕ3σ3<0,σ2μ3+qEμ0σ2ϕ2+σ1μ2<ϕ3σ3aϕ1σ1ϕ2σ2−μ0μ1ϕ2σ2+μ1μ2σ1,μ3+qEϕ3σ3<aϕ1σ1ϕ2σ2−μ0μ1ϕ2σ2+μ1μ2σ1σ2μ0σ2ϕ2+σ1μ2.Therefore,L∗ exists whenever(26)aϕ1σ1ϕ2σ2−μ0μ1ϕ2σ2+μ1μ2σ1σ2μ0σ2ϕ2+σ1μ2>μ3+qEϕ3σ3.The equilibriumE∗ is(27)N∗=σ1σ2σ3μ3+qE−μ2σ1σ2σ30−μ0σ22σ30ϕ2+μ1σ1σ32ϕ3+4A3+A422σ12σ32ϕ1ϕ3,P∗=−σ1σ2σ3μ3+qE+μ2σ1σ2σ30−μ0σ22σ30ϕ2−μ1σ1σ32ϕ3+4A3+A422σ1σ22σ30ϕ2,Z∗=σ1σ2σ3μ3+qE−μ2σ1σ2σ30−μ0σ22σ30ϕ2−μ1σ1σ32ϕ3+4A3+A422σ1σ2σ32ϕ3,L∗=−σ1σ2σ3μ3+qE−μ2σ1σ2σ30−μ0σ22σ30ϕ2−μ1σ1σ32ϕ3+4A3+A422σ1σ2σ3σ30,where(28)A3=aσ12σ22σ32σ30ϕ1ϕ2ϕ3,A4=σ1σ2σ3μ3+qE−σ2σ30μ2σ1+μ0σ2ϕ2+μ1σ1σ32ϕ3.From equation array (27), it follows that(29)N∗=σ1σ2σ3μ3+qE−μ2σ1σ2σ30−μ0σ22σ30ϕ2+μ1σ1σ32ϕ3+4A3+A422σ12σ32ϕ1ϕ3,P∗=σ1σ32ϕ1ϕ3N∗−σ2σ3μ3+qE+μ2σ2σ30−μ1σ32ϕ3σ22σ30ϕ2,Z∗=σ2σ30ϕ2P∗+σ3μ3+qE−μ2σ30σ32ϕ3,L∗=σ3ϕ3Z∗−μ3+qEσ30.The coexistence equilibrium,E∗=0.134444,503.68,0.0694424,128.578, for the autonomous model (2), with damped oscillations for the set of default parameters a=618.8, μ0=0.2924, σ1=9.13751, ϕ1=1, μ1=0.974667, σ2=3.655, ϕ2=1, μ2=0.30458, σ3=14.31542, σ30=0.000152, ϕ3=1, μ3=0.834558, q=0.00000918, and E=15250, are shown in Figure 5.Figure 5
Time series plot of (a) nutrients; (b) phytoplankton; (c) zooplankton; (d) Limnothrissa miodon for system (2) for model (2) with b=c=0 and assumed initial condition: N0=10,P0=7,Z0=4,L0=2 using the default parameter values.
(a)(b)(c)(d)


### 2.7. Stability Analysis for the Autonomous Model

Ifb=c=0 in system (2), then the model has a unique coexistence equilibrium E∗ which is locally stable.Theorem 3.
If the equilibriumE∗ exists, then it is locally asymptotically stable if the conditions in (39) are met.Proof.
The stability of the equilibrium states of model (2) are determined by the Jacobian matrix [30], J, of system (2):(30)J=−μ0−σ1P−σ1N00ϕ1σ1Pϕ1σ1N−μ1−σ2Z−σ2P00ϕ2σ2Zϕ2σ2P−μ2−σ3L−σ3Z00ϕ3σ3Lϕ3σ3Z−μ3−qE−2σ30L.
EvaluatingJ at E∗ results in(31)JE∗=−μ0−σ1c5−σ1c400ϕ1σ1c5ϕ1σ1c4−μ1−σ2c6−σ2c500ϕ2σ2c6ϕ2σ2c5−μ2−σ3c7−σ3c600ϕ3σ3c7ϕ3σ3c6−μ3−qE−2σ30c7,where c4=N∗, c5=P∗, c6=Z∗, and c7=L∗. Then, (31) simplifies to(32)JE∗=−μ0−σ1c5−σ1c400ϕ1σ1c50−σ2c500ϕ2σ2c60−σ3c600ϕ3σ3c7ϕ3σ3c6−μ3−qE−2σ30c7.
The eigenvalues of (32) are the roots of the auxiliary equation(33)Δλ=c6c7σ32λ2+λμ0+c5σ1λ+c4σ1ϕ1ϕ3+λλ2+λμ0+c5σ1λ+c4σ1ϕ1+c5c6λ+μ0+c5σ1σ2ϕ2λ+μ3+qE+2c7σ30−c6σ3ϕ3=0.
The auxiliary (33) can be expressed as follows:(34)Δλ=λ4+a1λ3+a2λ2+a3λ+a4=0,where(35)a1=μ0+μ3+qE+c5σ1+2c7σ30−c6σ3ϕ3,a2=c5qEσ1−c6μ0σ3ϕ3+c5μ3σ1+2c7μ0σ30+c4c5σ12ϕ1−c5c6σ3σ1ϕ3+c5c6σ22ϕ2+c6c7σ32ϕ3+2c5c7σ30σ1+qEμ0+μ0μ3,a3=c4c5σ12ϕ1−c6σ3ϕ3+2c7σ30+qE+μ3+c6c5σ22ϕ2−c6σ3ϕ3+2c7σ30+qE+μ0+μ3+c7σ1σ32ϕ3+c7μ0σ32ϕ3+c52σ1σ22ϕ2,a4=c5c6qEμ0σ22ϕ2+c52c6qEσ1σ22ϕ2+c5c6μ0μ3σ22ϕ2+c52c6μ3σ1σ22ϕ2+2c5c6c7μ0σ30σ22ϕ2−c5c62μ0σ3σ22ϕ2ϕ3+2c52c6c7σ1σ30σ22ϕ2−c52c62σ1σ3σ22ϕ2ϕ3+c4c5c6c7σ12σ32ϕ1ϕ3.
According to the Routh–Hurwitz criterion, all eigenvalues of the auxiliary (34) have negative real parts if(36)a1>0,a1a2−a3>0,a3a1a2−a3−a12a4>0,a4>0.Remark 1.
If the conditionsb=c=0 for E∗ are not satisfied for a given set of parameter values, then the respective steady state has a possibility of oscillatory behaviour for model (2).Theorem 4.
If the conditions in (42) for the Lyapunov function in (40) are met, the equilibrium E∗ is globally asymptotically stable.Proof.
The proof is based on Lyapunov’s second method. LetN−N∗>0, P−P∗>0, Z−Z∗>0, L−L∗>0. Let VN,P,Z,L be a positive definite Lyapunov function [31] such that VN∗,P∗,Z∗,L∗=0 by(37)VN,P,Z,L=ξ1N−N∗−N∗lnNN∗+ξ2P−P∗−P∗lnPP∗+ξ3Z−Z∗−Z∗lnZZ∗+ξ4L−L∗−L∗lnLL∗,where ξi′s, i=1,2,3,4, are positive constants. V is a positive definite function in the set Ψ, except at E∗ where it is zero. The derivative of V with respect to solution of system (2) is(38)V˙=ξ1N−N∗N˙N+ξ2P−P∗P˙P+ξ3Z−Z∗Z˙Z+ξ4L−L∗L˙L,=−ξ1N−N∗μ0+σ1P−aN−ξ2P−P∗μ1+σ2Z−ϕ1σ1N−ξ3Z−Z∗μ2+σ3L−ϕ2σ2P−ξ4L−L∗μ3+qE+σ30L−ϕ3σ3Z.
ThenV˙<0 if(39)μ0+σ1P>aN,μ1+σ2Z>ϕ1σ1N,μ2+σ3L>ϕ2σ2P,μ3+qE+σ30L>ϕ3σ3Z.
Thus,E∗ is globally asymptotically stable in the region bounded by all points (N>N∗, P>P∗, Z>Z∗, L>L∗) in (39).


#### 2.7.1. Fishing Effort

From model (2),(40)dLdt=ϕ3σ3ZL−μ3L−σ30L2−qEL.For the autonomous model (2), (43) at equilibrium is(41)qEL∗=ϕ3σ3Z∗L∗−μ3L∗−σ30L∗2.We let(42)h=ϕ3σ3Z∗L∗−μ3L∗−σ30L∗2,and ∂h/∂L∗=0 yields(43)ϕ3σ3Z∗−μ3−2σ30L∗=0,and ∴(44)h=qEϕ3σ3Z∗−μ3−qE2σ30,and ∂2h/∂L∗2<0 meaning h is the maximum sustainable catch. For sustainability of the kapenta resource, the sustainable fishing effort is(45)qE<ϕ3σ3Z∗−μ3.For the autonomous system, at equilibrium, the kapenta catch is(46)h=qeL∗=0.00000918×15250×128.578=18.00027711μgl−1,and it is shown in Figure 6. The surface (colour green) is for the equation h=ϕ3σ3Z∗L∗−μ3L∗−σ30L∗2 and the surface (colour magenta) is for the equation h=qEL∗.Figure 6
Plot of catch,Zt and Lt at steady state.Assuming the lake is at its capacity of 160 km3 and a carrying capacity of 240500 tonnes obtained by Tendaupenyu and Pyo [17], 18.00027711μgl−1 of Limnothrissa miodon is an average monthly catch of 2889.05 tonnes and an average yearly total of 34668.59 tonnes. From (43),(47)L∗=ϕ3σ3Z∗−μ32σ30.The maximum sustainable catch,h, can be written as(48)h=14σ30ϕ3σ3Z∗−μ32.h is a function of the variables in the food chain and parameters and this result is a true reflection of what happens in the lake as a result of the environmental factors which are not constant. The maximum sustainable catch of 34668.59 t is within the range of the MSY of 25372 t obtained by Tendaupenyu and Pyo [17] and MSY of 40000 t obtained by Magadza [14].


### 2.8. Optimal Harvesting for the Autonomous Model

For model (2) with b=c=0, we formulate a Hamiltonian function and apply Pontryagin’s maximum principle [32].(49)The optima sustainable yieldOSY=total revenueTR−total costTC,where TR and TC are the total revenue from selling harvested kapenta and the total cost of harvesting the kapenta, respectively. For the harvesting function, h=qEL,(50)OSY=αqL−βE.Then,(51)J=∫0∞e−δtαqL−βEtdtis a continuous time inflow of revenue, α is the price per unit of harvested kapenta, β is the utilised cost per unit effort, and δ is annual discount rate in the price of kapenta. Our goal is to maximize J subject to the equations in (18)–(21) and given the control constraints 0<E<Emax. To determine the optimal level of equilibrium, we employ Pontryagin’s maximum principle [32]. For the autonomous model (2) and (54), the Hamiltonian function is(52)H=e−δtαqL−βE+λ1ta−μ0N−σ1NP+λ2tγ1NP−μ1P−σ2PZ+λ3tγ2PZ−μ2Z−σ3ZL+λ4tγ3ZL−μ3L−σ30L2−qEL,where γ1=ϕ1σ1, γ2=ϕ2σ2, γ3=ϕ3σ3, λ1,λ2,λ3 and λ4 are adjoint variables and(53)σt=e−δtαqL2−β−λ4qL2is a switching function [32]. The condition in (54) must be met by the optimal control E that maximises H,(54)E=Emax,σt>0;λ4eδt<α−βqL2,0,σt<0;λ4eδt>α−βqL2,where α−β/qL2 is the net economic revenue from unit kapenta harvest and λ4eδt is the shadow price [32]. The optimal harvesting policy is(55)E=Emax,σt>0,0,σt<0,E∗,σt=0,where E=Emax if the shadow price is less than the net economic revenue on a unit kapenta harvest, E=0 if the shadow price is greater than the net economic revenue on a unit kapenta harvest, and E=E∗ when the shadow price equals the net economic revenue on a unit kapenta harvest [32]. If σt=0, then ∂H/∂E=0, indicating that H does not depend on E, the control variable, which is required for the singular control E∗ to be optimal on the interval 0<E∗<Emax [32]. Substituting σt into (53),(56)λ4=e−δtα−βqL2.The adjoint equations are(57)dλ1dt=−∂H∂N,dλ2dt=−∂H∂P,dλ3dt=−∂H∂Z,dλ4dt=−∂H∂L.From the equations in (60),(58)dλ1dt=−λ1μ0−σ1P−λ2γ1P,dλ2dt=−λ1−σ1N−λ2γ1N−μ1−σ2Z−λ3γ2Z,dλ3dt=−λ2−σ2P−λ3γ2P−μ2−σ3L−λ4γ3L,dλ4dt=−λ3−σ3Z−λ4γ3Z−μ3−2σ30L−qE−αqEe−δt.From (56),(59)dλ4dt=−δλ4.Substitutingσt=0 in (53), we obtain(60)λ4=α−βqL2e−δt.From (59) and the last equation of (58),(61)−δλ4=−λ3−σ3Z−λ4γ3Z−μ3−2σ30L−qE−αqEe−δt.Substituting (60) into (61) results in(62)λ3=e−δtα−β/qL2γ3Z−μ3−2σ30L−qE−δ+αqEσ3Z.Substituting (21) into (62) results in(63)λ3=−e−δtα−β/qL2δ+σ30L−αqEσ3Z.And (63) is written as(64)λ3=A1e−δt,where(65)A1=−α−β/qL2δ+σ30Ls−αqEσ3Z.From (67) and the third equation of (58),(66)−δA1e−δt=−λ2−σ2P−λ3γ2P−μ2−2σ3L−λ4γ3L.Substituting (20) into (66), we obtain(67)λ2=−e−δtδA1−α−β/qL2σ2P.And (67) is written as(68)λ2=A2e−δt,where(69)A2=−δA1−α−β/qL2σ2P.From equation (71) and the second equation of (58),(70)−δA2e−δt=λ1σ1N−λ3γ2Z.Substituting (20) into (66), we obtain(71)λ1=−e−δtδA2−A1γ2Zσ1N.And (67) is written as(72)λ1=A3e−δt,where(73)A3=−δA2−A1γ2Zσ1N.From (72), the first equation of (18), and (58),(74)−δA3e−δt=−λ1−aN−λ2γ1P.(75)A3=γ1PA2a/N+δ.Substituting the values ofA2, A1 into (75), we obtain(76)E=αqL2−βa+δNγ2Pσ2Zδ+σ30L+δδ2+σ3Z+δσ30L+γ1Pσ1N2σ3Z+δδ+σ30Lαq2L2a+δNδ2+γ2Pσ2Z+γ1δPσ1N2.The optimal solutionNδ,Pδ,Zδ,Lδ and the optimal harvesting effort E=Eδ are obtained by solving (18)–(21) with (76). Solving (18)–(21) and (76) using the same default parameter values and with α=0.86, β=0.001, and δ=0.12, we obtain the optimal equilibrium Nδ,Pδ,Zδ,Lδ=0.13448,503.543,0.0695343,128.543 and the optimal effort, E∗=15393.9 boat nights. The optimal number of boats is 15393.9/30.5≈505, in the Lake Kariba fishery. Since the Zimbabwean and Zambian sides of the Lake are in the ratio 0.55:0.45, we therefore recommend 278 boats and 227 on the Zimbabwean and Zambian sides of the kapenta fishery, respectively. The optimal number of fishing units of 505 is in agreement with the recommended 500 in the Lake Kariba fishery [12]. The plot of kapenta for different values of E is shown in Figure 7.Figure 7
Time series of Limnothrissa miodon for varying fishing effort for model system (2) with E=7625 (colour purple), E=15250 (colour magenta), and E∗=15394 (colour blue).The rate of change ofE with α in (77) is positive implying that if the price per unit of harvested kapenta increases the effort increases, for fixed β. This means that as the price rises of kapenta, the number of fishing vessels is likely to increase in the fishery.(77)∂E∂α=βa+δxγ2pσ2wδ+σ30y+δδ2+σ3w+δσ30y+γ1pσ1x2σ3w+δδ+σ30yα2q2y2a+δxδ2+γ2pσ2w+γ1δpσ1x2.The rate of change ofE with β in (78) is negative implying that if the utilised cost per unit effort increases the effort decreases, for fixed α. If the cost of fishing rises, then the number of fishing vessels in the fishery can decline.(78)∂E∂β=−a+δxγ2pσ2wδ+σ30y+δδ2+σ3w+δσ30y+γ1pσ1x2σ3w+δδ+σ30yαq2y2a+δxδ2+γ2pσ2w+γ1δpσ1x2.


## 2.1. Study Area

Lake Kariba is located on the Zambezi River, in a tropical area with seasonal rainfall between latitudes 1628′ to 1804′ S and longitudes 2642′ to 2903′E [15]. The offshore single-species pelagic kapenta fishery is highly mechanized, using light attraction and lift nets from boats, and is licensed [21]. According to the 2011 Frame survey [22], the fishing units on the Zimbabwean and Zambian sides of Lake Kariba have similar technical efficiencies and harvesting systems. Large circular dip nets and mercury light bulbs are used in the fishing gear to attract the kapenta. The lights used below and above the lake’s surface are powered by generators, the number of which varies from rig to rig. To lower or raise the dip net, hydraulic, mechanical, or manual winches are used [8].


## 2.2. Data Collection

The Lake Kariba Fisheries Research Institute collects data on catch, effort in the experimental gillnet, and 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 [8]. The University of Zimbabwe Lake Kariba Research Institute collects data on the lake’s water quality. The data on acoustic surveys of the biomass of kapenta are not often collected because of the lack of equipment and the cost involved and as such a few fishery surveys have been done on Lake Kariba using research vessels. Therefore, data on the biomass of the aquatic ecosystem of Lake Kariba are only limited to the hydroacoustic surveys of 1981, 1982, and 1983 [23]; September 1988 [24]; 1992 [25]; and for the period from 14 to 24 August 2014 [26]. Some data on phytoplankton, zooplankton, and kapenta in Lake Kariba were obtained from Ndebele-Murisa [21]. Some data is missing due to the fact that some of the fish caught in the offshore fishery are sold to buyers by the fishers before they come back to land and therefore are not recorded. As much as 30% of kapenta catch is thought to be sold illegally by boat crews before landing without the knowledge of the owners of the fishing vessels [27].Figure1(a) shows the time series plot of total nitrogen and chlorophyll a in μgl−1 in Lake Kariba from November 2016 to December 2017. The monthly averages of total nitrogen and chlorophyll a in Lake Kariba from April 2014 to December 2017 are illustrated in Figure 1(b). Figures 1(c) and 1(d) show the time series plot from November 2016 to December 2017 and the monthly averages from April 2014 to December 2017 of zooplankton in Lake Kariba in individuals per litre (ind./l), respectively.Figure 1
(a) Time series plot of total nitrogen and chlorophyll a; (b) a plot of monthly averages of total nitrogen and chlorophyll a; (c) time series plot of zooplankton; (d) a plot of monthly averages of zooplankton.
(a)(b)(c)(d)Figure2(a) shows the number of fishing vessels on the Zimbabwean side of Lake Kariba from 2008 to 2018. The time series plot shows a gradual increase in the number of units from 307 in 2008 to 530 in 2018. Figure 2(b) shows the number of fishing units on the Zambian and Zimbabwean sides of the lake from 1994 to 1999.Figure 2
(a) Time series plot of fishing units from 2008 to 2018. (b) Time series plot of fishing units from 1994 to 1999 [27].
(a)(b)Figure3(a) illustrates the time series plot of Limnothrissa miodon catches and the CPUE from November 2016 to December 2017. The catch and CPUE monthly averages for the period 1974 to 2018 are shown in Figure 3(b). Figures 3(a) and 3(b) are for the Zimbabwean side of the lake. Figure 3(b) shows a small peak in the catches of kapenta around March and a big peak around August. The peak in March is attributed to the inflow of nutrients through the inflowing rivers and the flooding of the Zambezi river and the one in August is due to the reloading of nutrients after the mixing of the lake around July.Figure 3
(a) Kapenta catch and CPUE time series plot; (b) a plot of monthly averages of kapenta catch and CPUE; (c) kapenta catch in Lake Kariba time series plot; (d) a plot of effort and CPUE in Lake Kariba.
(a)(b)(c)(d)The catches of kapenta versus effort are shown in Figure3(c). Figure 3(d) shows the effort, catch, and the CPUE of kapenta from 1974 to 2011. Figures 3(c) and 3(d) show data from both the Zambian and Zimbabwean sides of Lake Kariba. For the period 1974 to 2011, the global maximum is at an effort of 104131 boat nights giving a catch of 30943 t in 1990.


## 2.3. Model Equations

The model has four classes:N denotes nutrient concentration, P denotes phytoplankton population density, Z denotes zooplankton population density, and L denotes Limnothrissa miodon population density. Densities in each class are time-dependent and are denoted by Nt, Pt, Zt, and Lt, respectively. The Limnothrissa miodon model [28] is extended to include harvesting by the fishing vessels. Nutrients are assumed to enter the water body sinusoidally at the rate(1)a1+bsinωt+φ1+csin2ωt+φ2,where a,b,c>0 are constants, φ1 and φ2 are phase angles, and ω is the angular frequency, a constant. φ2 captures the delay in the lake mixing, ω is the frequency of the periodic input of nutrients, and ab and ac are the overall amplitude of the first and second nutrient periodic inputs, respectively. The first term of (1) represents the inflow of nutrients from river inflow and runoff and the second term of (1) captures the nutrient reloading into the lake which occurs around mid-year. The sinusoidal function in (1) is chosen for its simplicity and its parameters are to be determined by fitting. Natural depletion of nutrients occurs at a constant rate μ0. Phytoplankton depletes nutrients at a rate of σ1NP. The rate of growth of phytoplankton is ϕ1σ1NP. It is assumed that the depletion rate of phytoplankton caused by mortality is proportional to P. The change in density of the phytoplankton per unit time per zooplankton as the phytoplankton population density changes is of the modified Holling’s type-I response [29] and is denoted by σ2PZ [29]. The conversion coefficient from phytoplankton to zooplankton is ϕ2. It is assumed that the rate of zooplankton depletion 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 zooplankton density per unit time per Limnothrissa miodon as the zooplankton population density changes. The conversion coefficient from zooplankton to Limnothrissa miodon is phi ϕ3. It is assumed that the rate of Limnothrissa miodon depletion caused by mortality is proportional to L, and the rate of depletion caused by crowding is proportional to L2. Kapenta is harvested at a rate Gt. The model structure improves on the Limnothrissa miodon model in [28] and is shown in Figure 4.Figure 4
Flow diagram of the Limnothrissa miodon model with harvesting.We assume the following:(1)
Catch per unit effort of kapenta is proportional to the fish abundance,Gt/E∝Lt, and can be written as Gt=qELt [18, 19], where q≥0 is a fixed constant of proportionality called the catchability coefficient and E=Et≥0 is the fishing effort per unit time.(2)
Pelagic kapenta have schooling behaviour [18].(3)
Fishing using the rigs involves randomly searching for fish.(4)
Every kapenta in the body of water has an equal chance of being captured.The nonlinear nonautonomous dynamical system with Schaefer harvesting is(2)dNdt=a1+bsinωt+φ1+csin2ωt+φ2−μ0N−σ1NP,dPdt=ϕ1σ1NP−μ1P−σ2PZ,dZdt=ϕ2σ2PZ−μ2Z−σ3ZL,dLdt=ϕ3σ3ZL−μ3L−σ30L2−qEL..Model (2) initial condition is,(3)N0=ψ10,P0=ψ20,Z0=ψ30,L0=ψ40,ψi0>0,i=1,2,3,4,and the mathematically feasible region is defined as(4)Ω=N,P,Z,L∈ℝ4|N≥0,P≥0,Z≥0,L≥0.The crowding of the Limnothrissa miodon population coefficient isσ30>0. The constants of proportionality σii=1,2,3 are positive and depletion rate coefficients are represented by μii=0,1,2,3.


## 2.4. Positivity of Solutions

We prove that, for positive initial data for the model (2), the solutions will remain positive ∀t≥0.Theorem 1.
LetNt≥0,Pt≥0,Zt≥0,Lt≥0 be the initial data. Then, solutions of Nt, Pt, Zt, Lt of model (2) are positive ∀t≥0.Proof.
Based on the first equation of system (2) and Nt∈0,T, it follows that(5)N˙t≥−μ0Nt−σ1NtPt,∀t∈0,T.
Consequently, we obtain(6)Nt≥N0exp∫0t−μ0−σ1Psds≥0,∀t∈0,T.
According to the second equation of system (2),(7)P˙t≥−μ1Pt−σ2PtZt,∀t∈0,T.
The direct integration of (7) yields(8)Pt≥P0exp∫0t−μ1−σ2Zsds≥0,∀t∈0,T.
According to the third equation of system (2),(9)Z˙t≥−μ2Zt−σ3ZtLt,∀t∈0,T.
The direct integration of (9) yields(10)Zt≥Z0exp∫0t−μ2−σ3Lsds≥0,∀t∈0,T.
Considering the variableLt∈0,T, it follows from the fourth equation of system (2) that(11)L˙t≥−Ltμ3+qE+σ30Lt,∀t∈0,T.
The direct integration of (11) yields(12)Lt≥μ3+qEL0e−μ3+qEtμ3+qE+σ30L01−e−μ3+qEt≥0,∀t∈0,T.
As a result, system solutions of (2) with initial conditions (3) remain positive ∀t≥0.


## 2.5. Existence of Solutions

Theorem 2.
A solution of system (2) is feasible.Proof.
We show that all feasible solutions of system (2) are uniformly bounded in Ω⊂ℝ4. Let any solution of system (2) be Nt,Pt,Zt,Lt∈ℝ4, with nonnegative initial conditions.
LetCt=Nt+Pt+Zt+Lt, and then(13)dCdt=a1+bsinωt+φ1+csin2ωt+φ2−μ0N−σ1NP+ϕ1σ1NP−μ1P−σ2PZ+ϕ2σ2PZ−μ2Z−σ3ZL+ϕ3σ3ZL−μ3L−σ30L2−qEL,=a1+bsinωt+φ1+csin2ωt+φ2−μ0N−μ1P−μ2Z−μ3+qEL−σ30L2+σ1ϕ1−1NP+σ2ϕ2−1PZ+σ3ϕ3−1ZL,≤a1+bsinωt+φ1+csin2ωt+φ2−μ0N−μ1P−μ2Z−μ3+qEL,≤a1+bsinwt+φ1+csin2wt+φ2−mCt,where m=minμ0,μ1,μ2,μ3+qE. Thus,(14)dCtdt+mCt≤a1+bsinωt+φ1+csin2ωt+φ2.
Equation (14) is a first-order differential inequality, with the solution given by(15)0<CN,P,Z,L≤abm4sintω+φ1+4bm2ω2sintω+φ1+cm4sin2tω+φ2+cm2ω2sin2tω+φ2m5+5m3ω2+4mω4−abmωm2+4ω2costω+φ1+2cmωm2+ω2cos2tω+φ2+m4+5m2ω2+4ω4m5+5m3ω2+4mω4+C0e−mt,as t⟶∞, (15) becomes(16)0<CN,P,Z,L≤abm4sintω+φ1+4bm2ω2sintω+φ1+cm4sin2tω+φ2+cm2ω2sin2tω+φ2m5+5m3ω2+4mω4−abmωm2+4ω2costω+φ1+2cmωm2+ω2cos2tω+φ2+m4+5m2ω2+4ω4m5+5m3ω2+4mω4.
Therefore, all solutions of system (2) enter the feasible region,(17)Ω=Nt,Pt,Zt,Lt∈ℝ+4:C≤abm4sintω+φ1+4bm2ω2sintω+φ1+cm4sin2tω+φ2+cm2ω2sin2tω+φ2m5+5m3ω2+4mω4−abmωm2+4ω2costω+φ1+2cmωm2+ω2cos2tω+φ2+m4+5m2ω2+4ω4m5+5m3ω2+4mω4+ς,∀ς>0.
This concludes the theorem’s proof.


## 2.6. Equilibrium States for the Autonomous Model

Of biological interest is the coexistence equilibrium. Ifb=c=0 in system (2), then the model is autonomous and has a unique coexistence equilibrium, E∗=N∗,P∗,Z∗,L∗, which is obtained through the solution of the equations:(18)a−μ0N−σ1NP=0,(19)ϕ1σ1N−μ1−σ2Z=0,(20)ϕ2σ2P−μ2−σ3L=0,(21)ϕ3σ3Z−μ3−qE−σ30L=0.Solving forN, P, Z, and L in (18)–(21) gives(22)σ1σ2σ3σ30L∗2+σ1σ2σ3μ3+qE+σ1σ2μ2σ30+μ0σ22ϕ2σ30+μ1σ1σ32ϕ3L∗+μ0μ1ϕ2σ2ϕ3σ3+μ1σ1μ2ϕ3σ3+μ0σ22ϕ2μ3+qE+σ1σ2μ2μ3+qE−aϕ1σ1ϕ2σ2ϕ3σ3=0.If expression (24) is positive, (23) will have a unique positive root:(23)L∗=A12−4σ1σ2σ3σ30A2−A12σ1σ2σ3σ30>0,where(24)A1=σ1σ2σ3μ3+qE+σ1σ2μ2σ30+μ0σ22ϕ2σ30+μ1σ1σ32ϕ3,A2=μ0μ1ϕ2σ2ϕ3σ3+μ1σ1μ2ϕ3σ3+μ0σ22ϕ2μ3+qE+σ1σ2μ2μ3+qE−aϕ1σ1ϕ2σ2ϕ3σ3.And (23) can be written as(25)A12−4σ1σ2σ3σ30A2>A12,σ1σ2σ3σ30A2<0,σ1σ2σ3σ30μ0μ1ϕ2σ2ϕ3σ3+μ1σ1μ2ϕ3σ3+μ0σ22ϕ2μ3+qE+σ1σ2μ2μ3+qE−aϕ1σ1ϕ2σ2ϕ3σ3<0,σ2μ3+qEμ0σ2ϕ2+σ1μ2<ϕ3σ3aϕ1σ1ϕ2σ2−μ0μ1ϕ2σ2+μ1μ2σ1,μ3+qEϕ3σ3<aϕ1σ1ϕ2σ2−μ0μ1ϕ2σ2+μ1μ2σ1σ2μ0σ2ϕ2+σ1μ2.Therefore,L∗ exists whenever(26)aϕ1σ1ϕ2σ2−μ0μ1ϕ2σ2+μ1μ2σ1σ2μ0σ2ϕ2+σ1μ2>μ3+qEϕ3σ3.The equilibriumE∗ is(27)N∗=σ1σ2σ3μ3+qE−μ2σ1σ2σ30−μ0σ22σ30ϕ2+μ1σ1σ32ϕ3+4A3+A422σ12σ32ϕ1ϕ3,P∗=−σ1σ2σ3μ3+qE+μ2σ1σ2σ30−μ0σ22σ30ϕ2−μ1σ1σ32ϕ3+4A3+A422σ1σ22σ30ϕ2,Z∗=σ1σ2σ3μ3+qE−μ2σ1σ2σ30−μ0σ22σ30ϕ2−μ1σ1σ32ϕ3+4A3+A422σ1σ2σ32ϕ3,L∗=−σ1σ2σ3μ3+qE−μ2σ1σ2σ30−μ0σ22σ30ϕ2−μ1σ1σ32ϕ3+4A3+A422σ1σ2σ3σ30,where(28)A3=aσ12σ22σ32σ30ϕ1ϕ2ϕ3,A4=σ1σ2σ3μ3+qE−σ2σ30μ2σ1+μ0σ2ϕ2+μ1σ1σ32ϕ3.From equation array (27), it follows that(29)N∗=σ1σ2σ3μ3+qE−μ2σ1σ2σ30−μ0σ22σ30ϕ2+μ1σ1σ32ϕ3+4A3+A422σ12σ32ϕ1ϕ3,P∗=σ1σ32ϕ1ϕ3N∗−σ2σ3μ3+qE+μ2σ2σ30−μ1σ32ϕ3σ22σ30ϕ2,Z∗=σ2σ30ϕ2P∗+σ3μ3+qE−μ2σ30σ32ϕ3,L∗=σ3ϕ3Z∗−μ3+qEσ30.The coexistence equilibrium,E∗=0.134444,503.68,0.0694424,128.578, for the autonomous model (2), with damped oscillations for the set of default parameters a=618.8, μ0=0.2924, σ1=9.13751, ϕ1=1, μ1=0.974667, σ2=3.655, ϕ2=1, μ2=0.30458, σ3=14.31542, σ30=0.000152, ϕ3=1, μ3=0.834558, q=0.00000918, and E=15250, are shown in Figure 5.Figure 5
Time series plot of (a) nutrients; (b) phytoplankton; (c) zooplankton; (d) Limnothrissa miodon for system (2) for model (2) with b=c=0 and assumed initial condition: N0=10,P0=7,Z0=4,L0=2 using the default parameter values.
(a)(b)(c)(d)


## 2.7. Stability Analysis for the Autonomous Model

Ifb=c=0 in system (2), then the model has a unique coexistence equilibrium E∗ which is locally stable.Theorem 3.
If the equilibriumE∗ exists, then it is locally asymptotically stable if the conditions in (39) are met.Proof.
The stability of the equilibrium states of model (2) are determined by the Jacobian matrix [30], J, of system (2):(30)J=−μ0−σ1P−σ1N00ϕ1σ1Pϕ1σ1N−μ1−σ2Z−σ2P00ϕ2σ2Zϕ2σ2P−μ2−σ3L−σ3Z00ϕ3σ3Lϕ3σ3Z−μ3−qE−2σ30L.
EvaluatingJ at E∗ results in(31)JE∗=−μ0−σ1c5−σ1c400ϕ1σ1c5ϕ1σ1c4−μ1−σ2c6−σ2c500ϕ2σ2c6ϕ2σ2c5−μ2−σ3c7−σ3c600ϕ3σ3c7ϕ3σ3c6−μ3−qE−2σ30c7,where c4=N∗, c5=P∗, c6=Z∗, and c7=L∗. Then, (31) simplifies to(32)JE∗=−μ0−σ1c5−σ1c400ϕ1σ1c50−σ2c500ϕ2σ2c60−σ3c600ϕ3σ3c7ϕ3σ3c6−μ3−qE−2σ30c7.
The eigenvalues of (32) are the roots of the auxiliary equation(33)Δλ=c6c7σ32λ2+λμ0+c5σ1λ+c4σ1ϕ1ϕ3+λλ2+λμ0+c5σ1λ+c4σ1ϕ1+c5c6λ+μ0+c5σ1σ2ϕ2λ+μ3+qE+2c7σ30−c6σ3ϕ3=0.
The auxiliary (33) can be expressed as follows:(34)Δλ=λ4+a1λ3+a2λ2+a3λ+a4=0,where(35)a1=μ0+μ3+qE+c5σ1+2c7σ30−c6σ3ϕ3,a2=c5qEσ1−c6μ0σ3ϕ3+c5μ3σ1+2c7μ0σ30+c4c5σ12ϕ1−c5c6σ3σ1ϕ3+c5c6σ22ϕ2+c6c7σ32ϕ3+2c5c7σ30σ1+qEμ0+μ0μ3,a3=c4c5σ12ϕ1−c6σ3ϕ3+2c7σ30+qE+μ3+c6c5σ22ϕ2−c6σ3ϕ3+2c7σ30+qE+μ0+μ3+c7σ1σ32ϕ3+c7μ0σ32ϕ3+c52σ1σ22ϕ2,a4=c5c6qEμ0σ22ϕ2+c52c6qEσ1σ22ϕ2+c5c6μ0μ3σ22ϕ2+c52c6μ3σ1σ22ϕ2+2c5c6c7μ0σ30σ22ϕ2−c5c62μ0σ3σ22ϕ2ϕ3+2c52c6c7σ1σ30σ22ϕ2−c52c62σ1σ3σ22ϕ2ϕ3+c4c5c6c7σ12σ32ϕ1ϕ3.
According to the Routh–Hurwitz criterion, all eigenvalues of the auxiliary (34) have negative real parts if(36)a1>0,a1a2−a3>0,a3a1a2−a3−a12a4>0,a4>0.Remark 1.
If the conditionsb=c=0 for E∗ are not satisfied for a given set of parameter values, then the respective steady state has a possibility of oscillatory behaviour for model (2).Theorem 4.
If the conditions in (42) for the Lyapunov function in (40) are met, the equilibrium E∗ is globally asymptotically stable.Proof.
The proof is based on Lyapunov’s second method. LetN−N∗>0, P−P∗>0, Z−Z∗>0, L−L∗>0. Let VN,P,Z,L be a positive definite Lyapunov function [31] such that VN∗,P∗,Z∗,L∗=0 by(37)VN,P,Z,L=ξ1N−N∗−N∗lnNN∗+ξ2P−P∗−P∗lnPP∗+ξ3Z−Z∗−Z∗lnZZ∗+ξ4L−L∗−L∗lnLL∗,where ξi′s, i=1,2,3,4, are positive constants. V is a positive definite function in the set Ψ, except at E∗ where it is zero. The derivative of V with respect to solution of system (2) is(38)V˙=ξ1N−N∗N˙N+ξ2P−P∗P˙P+ξ3Z−Z∗Z˙Z+ξ4L−L∗L˙L,=−ξ1N−N∗μ0+σ1P−aN−ξ2P−P∗μ1+σ2Z−ϕ1σ1N−ξ3Z−Z∗μ2+σ3L−ϕ2σ2P−ξ4L−L∗μ3+qE+σ30L−ϕ3σ3Z.
ThenV˙<0 if(39)μ0+σ1P>aN,μ1+σ2Z>ϕ1σ1N,μ2+σ3L>ϕ2σ2P,μ3+qE+σ30L>ϕ3σ3Z.
Thus,E∗ is globally asymptotically stable in the region bounded by all points (N>N∗, P>P∗, Z>Z∗, L>L∗) in (39).


### 2.7.1. Fishing Effort

From model (2),(40)dLdt=ϕ3σ3ZL−μ3L−σ30L2−qEL.For the autonomous model (2), (43) at equilibrium is(41)qEL∗=ϕ3σ3Z∗L∗−μ3L∗−σ30L∗2.We let(42)h=ϕ3σ3Z∗L∗−μ3L∗−σ30L∗2,and ∂h/∂L∗=0 yields(43)ϕ3σ3Z∗−μ3−2σ30L∗=0,and ∴(44)h=qEϕ3σ3Z∗−μ3−qE2σ30,and ∂2h/∂L∗2<0 meaning h is the maximum sustainable catch. For sustainability of the kapenta resource, the sustainable fishing effort is(45)qE<ϕ3σ3Z∗−μ3.For the autonomous system, at equilibrium, the kapenta catch is(46)h=qeL∗=0.00000918×15250×128.578=18.00027711μgl−1,and it is shown in Figure 6. The surface (colour green) is for the equation h=ϕ3σ3Z∗L∗−μ3L∗−σ30L∗2 and the surface (colour magenta) is for the equation h=qEL∗.Figure 6
Plot of catch,Zt and Lt at steady state.Assuming the lake is at its capacity of 160 km3 and a carrying capacity of 240500 tonnes obtained by Tendaupenyu and Pyo [17], 18.00027711μgl−1 of Limnothrissa miodon is an average monthly catch of 2889.05 tonnes and an average yearly total of 34668.59 tonnes. From (43),(47)L∗=ϕ3σ3Z∗−μ32σ30.The maximum sustainable catch,h, can be written as(48)h=14σ30ϕ3σ3Z∗−μ32.h is a function of the variables in the food chain and parameters and this result is a true reflection of what happens in the lake as a result of the environmental factors which are not constant. The maximum sustainable catch of 34668.59 t is within the range of the MSY of 25372 t obtained by Tendaupenyu and Pyo [17] and MSY of 40000 t obtained by Magadza [14].


## 2.7.1. Fishing Effort

From model (2),(40)dLdt=ϕ3σ3ZL−μ3L−σ30L2−qEL.For the autonomous model (2), (43) at equilibrium is(41)qEL∗=ϕ3σ3Z∗L∗−μ3L∗−σ30L∗2.We let(42)h=ϕ3σ3Z∗L∗−μ3L∗−σ30L∗2,and ∂h/∂L∗=0 yields(43)ϕ3σ3Z∗−μ3−2σ30L∗=0,and ∴(44)h=qEϕ3σ3Z∗−μ3−qE2σ30,and ∂2h/∂L∗2<0 meaning h is the maximum sustainable catch. For sustainability of the kapenta resource, the sustainable fishing effort is(45)qE<ϕ3σ3Z∗−μ3.For the autonomous system, at equilibrium, the kapenta catch is(46)h=qeL∗=0.00000918×15250×128.578=18.00027711μgl−1,and it is shown in Figure 6. The surface (colour green) is for the equation h=ϕ3σ3Z∗L∗−μ3L∗−σ30L∗2 and the surface (colour magenta) is for the equation h=qEL∗.Figure 6
Plot of catch,Zt and Lt at steady state.Assuming the lake is at its capacity of 160 km3 and a carrying capacity of 240500 tonnes obtained by Tendaupenyu and Pyo [17], 18.00027711μgl−1 of Limnothrissa miodon is an average monthly catch of 2889.05 tonnes and an average yearly total of 34668.59 tonnes. From (43),(47)L∗=ϕ3σ3Z∗−μ32σ30.The maximum sustainable catch,h, can be written as(48)h=14σ30ϕ3σ3Z∗−μ32.h is a function of the variables in the food chain and parameters and this result is a true reflection of what happens in the lake as a result of the environmental factors which are not constant. The maximum sustainable catch of 34668.59 t is within the range of the MSY of 25372 t obtained by Tendaupenyu and Pyo [17] and MSY of 40000 t obtained by Magadza [14].


## 2.8. Optimal Harvesting for the Autonomous Model

For model (2) with b=c=0, we formulate a Hamiltonian function and apply Pontryagin’s maximum principle [32].(49)The optima sustainable yieldOSY=total revenueTR−total costTC,where TR and TC are the total revenue from selling harvested kapenta and the total cost of harvesting the kapenta, respectively. For the harvesting function, h=qEL,(50)OSY=αqL−βE.Then,(51)J=∫0∞e−δtαqL−βEtdtis a continuous time inflow of revenue, α is the price per unit of harvested kapenta, β is the utilised cost per unit effort, and δ is annual discount rate in the price of kapenta. Our goal is to maximize J subject to the equations in (18)–(21) and given the control constraints 0<E<Emax. To determine the optimal level of equilibrium, we employ Pontryagin’s maximum principle [32]. For the autonomous model (2) and (54), the Hamiltonian function is(52)H=e−δtαqL−βE+λ1ta−μ0N−σ1NP+λ2tγ1NP−μ1P−σ2PZ+λ3tγ2PZ−μ2Z−σ3ZL+λ4tγ3ZL−μ3L−σ30L2−qEL,where γ1=ϕ1σ1, γ2=ϕ2σ2, γ3=ϕ3σ3, λ1,λ2,λ3 and λ4 are adjoint variables and(53)σt=e−δtαqL2−β−λ4qL2is a switching function [32]. The condition in (54) must be met by the optimal control E that maximises H,(54)E=Emax,σt>0;λ4eδt<α−βqL2,0,σt<0;λ4eδt>α−βqL2,where α−β/qL2 is the net economic revenue from unit kapenta harvest and λ4eδt is the shadow price [32]. The optimal harvesting policy is(55)E=Emax,σt>0,0,σt<0,E∗,σt=0,where E=Emax if the shadow price is less than the net economic revenue on a unit kapenta harvest, E=0 if the shadow price is greater than the net economic revenue on a unit kapenta harvest, and E=E∗ when the shadow price equals the net economic revenue on a unit kapenta harvest [32]. If σt=0, then ∂H/∂E=0, indicating that H does not depend on E, the control variable, which is required for the singular control E∗ to be optimal on the interval 0<E∗<Emax [32]. Substituting σt into (53),(56)λ4=e−δtα−βqL2.The adjoint equations are(57)dλ1dt=−∂H∂N,dλ2dt=−∂H∂P,dλ3dt=−∂H∂Z,dλ4dt=−∂H∂L.From the equations in (60),(58)dλ1dt=−λ1μ0−σ1P−λ2γ1P,dλ2dt=−λ1−σ1N−λ2γ1N−μ1−σ2Z−λ3γ2Z,dλ3dt=−λ2−σ2P−λ3γ2P−μ2−σ3L−λ4γ3L,dλ4dt=−λ3−σ3Z−λ4γ3Z−μ3−2σ30L−qE−αqEe−δt.From (56),(59)dλ4dt=−δλ4.Substitutingσt=0 in (53), we obtain(60)λ4=α−βqL2e−δt.From (59) and the last equation of (58),(61)−δλ4=−λ3−σ3Z−λ4γ3Z−μ3−2σ30L−qE−αqEe−δt.Substituting (60) into (61) results in(62)λ3=e−δtα−β/qL2γ3Z−μ3−2σ30L−qE−δ+αqEσ3Z.Substituting (21) into (62) results in(63)λ3=−e−δtα−β/qL2δ+σ30L−αqEσ3Z.And (63) is written as(64)λ3=A1e−δt,where(65)A1=−α−β/qL2δ+σ30Ls−αqEσ3Z.From (67) and the third equation of (58),(66)−δA1e−δt=−λ2−σ2P−λ3γ2P−μ2−2σ3L−λ4γ3L.Substituting (20) into (66), we obtain(67)λ2=−e−δtδA1−α−β/qL2σ2P.And (67) is written as(68)λ2=A2e−δt,where(69)A2=−δA1−α−β/qL2σ2P.From equation (71) and the second equation of (58),(70)−δA2e−δt=λ1σ1N−λ3γ2Z.Substituting (20) into (66), we obtain(71)λ1=−e−δtδA2−A1γ2Zσ1N.And (67) is written as(72)λ1=A3e−δt,where(73)A3=−δA2−A1γ2Zσ1N.From (72), the first equation of (18), and (58),(74)−δA3e−δt=−λ1−aN−λ2γ1P.(75)A3=γ1PA2a/N+δ.Substituting the values ofA2, A1 into (75), we obtain(76)E=αqL2−βa+δNγ2Pσ2Zδ+σ30L+δδ2+σ3Z+δσ30L+γ1Pσ1N2σ3Z+δδ+σ30Lαq2L2a+δNδ2+γ2Pσ2Z+γ1δPσ1N2.The optimal solutionNδ,Pδ,Zδ,Lδ and the optimal harvesting effort E=Eδ are obtained by solving (18)–(21) with (76). Solving (18)–(21) and (76) using the same default parameter values and with α=0.86, β=0.001, and δ=0.12, we obtain the optimal equilibrium Nδ,Pδ,Zδ,Lδ=0.13448,503.543,0.0695343,128.543 and the optimal effort, E∗=15393.9 boat nights. The optimal number of boats is 15393.9/30.5≈505, in the Lake Kariba fishery. Since the Zimbabwean and Zambian sides of the Lake are in the ratio 0.55:0.45, we therefore recommend 278 boats and 227 on the Zimbabwean and Zambian sides of the kapenta fishery, respectively. The optimal number of fishing units of 505 is in agreement with the recommended 500 in the Lake Kariba fishery [12]. The plot of kapenta for different values of E is shown in Figure 7.Figure 7
Time series of Limnothrissa miodon for varying fishing effort for model system (2) with E=7625 (colour purple), E=15250 (colour magenta), and E∗=15394 (colour blue).The rate of change ofE with α in (77) is positive implying that if the price per unit of harvested kapenta increases the effort increases, for fixed β. This means that as the price rises of kapenta, the number of fishing vessels is likely to increase in the fishery.(77)∂E∂α=βa+δxγ2pσ2wδ+σ30y+δδ2+σ3w+δσ30y+γ1pσ1x2σ3w+δδ+σ30yα2q2y2a+δxδ2+γ2pσ2w+γ1δpσ1x2.The rate of change ofE with β in (78) is negative implying that if the utilised cost per unit effort increases the effort decreases, for fixed α. If the cost of fishing rises, then the number of fishing vessels in the fishery can decline.(78)∂E∂β=−a+δxγ2pσ2wδ+σ30y+δδ2+σ3w+δσ30y+γ1pσ1x2σ3w+δδ+σ30yαq2y2a+δxδ2+γ2pσ2w+γ1δpσ1x2.


## 3. Results

### 3.1. Data Fitting

The Fourier series equations (79)–(82) were used in the data fitting.(79)ft=a02+∑n=1∞ancosnt+∑n=1∞bnsinnt,where(80)a0=1π∫−ππftdt,(81)an=1π∫−ππftcosntdt,(82)bn=1π∫−ππftsinntdt.and n=1,2,3,4,… [33].A Fourier series with an order of fit ofn=2 in (83) was used in MATLAB R2016a to estimate the parameters of the total nitrogen data for the period from April 2014 to December 2017. The monthly averages from January to December were used in the data fitting.(83)ft=a0+a1cosωt+b1sinωt+a2cos2ωt+b2sin2ωt.The parameter estimates of (83) are shown in Table 1.Table 1
Fourier coefficients from data fitting.
CoefficientEstimate95% confidence boundsa0618.8533.1,704.6a111.03−166.4,188.5b1120.44.689,236.1a2−44.29−195.6,107.1b236.75−93.15,166.6ω0.61060.4416,0.7795The actual and fitted plots of monthly average total nitrogen are shown in Figure8.Figure 8
Monthly average of total nitrogen (μg / l) and fitted curve.The fitted models in (83) goodness of fit statistics are sum of squared estimate of errors (SSE) =7.615×104, R2=0.6095, and root mean square error (RMSE) =112.7.Factoring trigonometric functions in (83) using Wolfram Mathematica 11.0 results in(84)ft=11.0356.1015+10.9614sin0.0913563+0.6106t−5.21772sin0.878172−1.2212t.and we obtain(85)ft=618.81+0.195385sin0.6106t+0.0913563+0.0930049sin1.2212t+0.878172.Then (85) can be generalised as(86)ft=a1+bsinωt+φ1+csin2ωt+φ2.


### 3.2. Numerical Simulations

Model parameters of system (2) and their interpretations are shown in Table 2. Besides the catchability coefficient q of L, the description, symbol, value, and source of the parameters and coefficients in Table 2 are the same as those used in the study by Mutasa et al. [28], except that the parameters are monthly parameters. Numerical simulations for system (2) are done in Wolfram Mathematica 11.0 using a fourth-order Runge–Kutta numerical scheme. The units of the variables N, P, Z, and L are μgl−1.Table 2
Model parameters and their interpretations.
SymbolDefinitionValueSourceμ0N natural depletion rate coefficient0.2924month−1[34]σ1P uptake rate of N9.13751lμg−1month−1ϕ1P conversion coefficient of N1μ1P natural depletion rate coefficient3.965,0.976−2.44,13.4505month−1[13, 35, 36]σ2Z grazing rate of P18.3−42.7,17.08,3.05−21.045lμg−1month−1[13,37, 38]ϕ2Z conversion coefficient of P grazing0.2−0.75,1.5[39, 40]μ2Z natural depletion rate coefficient0.305,1.6104month−1[13, 41]σ3L grazing rate of Z14.335lμg−1month−1[13]σ30L crowding effect coefficient0.000152lμg−1month−1ϕ3L conversion coefficient of Z grazing1μ3L natural mortality0.2928,0.2745,0.7625−1.11935month−1[7, 9, 13]qL catchability coefficient0.00000153[17]


#### 3.2.1. Nonautonomous Model

For the set of the default parameters, the phase portraits showing the time series of nutrients, phytoplankton, zooplankton, and kapenta with harvesting are shown in Figures9(a)–9(d), respectively. Figure 9 shows stable periodic solutions of period 1 corresponding to the default parameter values. For varying initial conditions, all initial conditions lead to periodic solutions. Different initial conditions lead to the same periodic solution curve and therefore we have a stable limit cycle in the phase space. The time series for zooplankton initially has irregular fluctuations showing that the density of zooplankton is highly variable and eventually stabilizes into a periodic orbit after 40 months. This could be attributed to the initial phase in the lake when the fishery started to operate and then eventually stabilized after some time.Figure 9
Phase portrait showing the time series of (a) nutrients; (b) phytoplankton; (c) zooplankton; (d) Limnothrissa miodon for model (2) and assumed initial condition: N0=9.5,P0=6.5,Z0=3.5,L0=1.5 (time series with colour green); N0=10,P0=7,Z0=4,L0=2 (time series with colour purple); N0=10.5,P0=7.5,Z0=4.5,L0=2.5 (time series with colour magenta); using the default parameter.
(a)(b)(c)(d)Figures10(a)–10(b) show oscillatory dynamics in the phase portraits of the nutrient and phytoplankton; phytoplankton and zooplankton; zooplankton and Limnothrissa miodon; phytoplankton, zooplankton, and Limnothrissa miodon dynamics for system (2), respectively. The phase portraits show that the equilibrium point at E∗=0.134,503.68,0.0694,128.578 (point with colour yellow) for the autonomous model is enclosed by the periodic solutions of the nonautonomous model, but the system does not go to the equilibrium point.Figure 10
Phase portrait showing the dynamics of (a) phytoplankton and nutrients; (b) zooplankton and phytoplankton; (c) Limnothrissa miodon and zooplankton; (d) Limnothrissa miodon, zooplankton, and phytoplankton for model system (2) with assumed initial condition: N0=10,P0=7,Z0=4,L0=2 using the default parameter values.
(a)(b)(c)(d)


#### 3.2.2. Effects of Nutrients and Fishing Effort for the Nonautonomous Model

For the set of the default parameters, Figures11(a) and 11(b) illustrate the effect of varying the inflow of nutrients parameter a in model (2) on the dynamics of Limnothrissa miodon with E=15250. Figures 11(c) and 11(d) illustrate the effect of varying fishing effort in model (2) on the dynamics of Limnothrissa miodon.Figure 11
(a) Time series of Limnothrissa miodon for varyinga for model system (2); (b) 12-month average plot of Limnothrissa miodon for varying a; (c) time series of Limnothrissa miodon for varying fishing effort for model system (2) with E=0, E=22875, and E=38125; (d) 12-month average plot of Limnothrissa miodon to fishing effort and with assumed initial condition: N0=10,P0=7,Z0=4,L0=2 using the default parameter values.
(a)(b)(c)(d)Simulation results from Figures11(a)–11(d) show that the abundance of Limnothrissa miodon is more closely related to the availability of nutrients than to harvesting. Marshall [42] and Paulsen [43] in their studies showed that the abundance of Limnothrissa miodon is mainly due to river inflow and the availability of food compared to harvesting.


## 3.1. Data Fitting

The Fourier series equations (79)–(82) were used in the data fitting.(79)ft=a02+∑n=1∞ancosnt+∑n=1∞bnsinnt,where(80)a0=1π∫−ππftdt,(81)an=1π∫−ππftcosntdt,(82)bn=1π∫−ππftsinntdt.and n=1,2,3,4,… [33].A Fourier series with an order of fit ofn=2 in (83) was used in MATLAB R2016a to estimate the parameters of the total nitrogen data for the period from April 2014 to December 2017. The monthly averages from January to December were used in the data fitting.(83)ft=a0+a1cosωt+b1sinωt+a2cos2ωt+b2sin2ωt.The parameter estimates of (83) are shown in Table 1.Table 1
Fourier coefficients from data fitting.
CoefficientEstimate95% confidence boundsa0618.8533.1,704.6a111.03−166.4,188.5b1120.44.689,236.1a2−44.29−195.6,107.1b236.75−93.15,166.6ω0.61060.4416,0.7795The actual and fitted plots of monthly average total nitrogen are shown in Figure8.Figure 8
Monthly average of total nitrogen (μg / l) and fitted curve.The fitted models in (83) goodness of fit statistics are sum of squared estimate of errors (SSE) =7.615×104, R2=0.6095, and root mean square error (RMSE) =112.7.Factoring trigonometric functions in (83) using Wolfram Mathematica 11.0 results in(84)ft=11.0356.1015+10.9614sin0.0913563+0.6106t−5.21772sin0.878172−1.2212t.and we obtain(85)ft=618.81+0.195385sin0.6106t+0.0913563+0.0930049sin1.2212t+0.878172.Then (85) can be generalised as(86)ft=a1+bsinωt+φ1+csin2ωt+φ2.


## 3.2. Numerical Simulations

Model parameters of system (2) and their interpretations are shown in Table 2. Besides the catchability coefficient q of L, the description, symbol, value, and source of the parameters and coefficients in Table 2 are the same as those used in the study by Mutasa et al. [28], except that the parameters are monthly parameters. Numerical simulations for system (2) are done in Wolfram Mathematica 11.0 using a fourth-order Runge–Kutta numerical scheme. The units of the variables N, P, Z, and L are μgl−1.Table 2
Model parameters and their interpretations.
SymbolDefinitionValueSourceμ0N natural depletion rate coefficient0.2924month−1[34]σ1P uptake rate of N9.13751lμg−1month−1ϕ1P conversion coefficient of N1μ1P natural depletion rate coefficient3.965,0.976−2.44,13.4505month−1[13, 35, 36]σ2Z grazing rate of P18.3−42.7,17.08,3.05−21.045lμg−1month−1[13,37, 38]ϕ2Z conversion coefficient of P grazing0.2−0.75,1.5[39, 40]μ2Z natural depletion rate coefficient0.305,1.6104month−1[13, 41]σ3L grazing rate of Z14.335lμg−1month−1[13]σ30L crowding effect coefficient0.000152lμg−1month−1ϕ3L conversion coefficient of Z grazing1μ3L natural mortality0.2928,0.2745,0.7625−1.11935month−1[7, 9, 13]qL catchability coefficient0.00000153[17]


### 3.2.1. Nonautonomous Model

For the set of the default parameters, the phase portraits showing the time series of nutrients, phytoplankton, zooplankton, and kapenta with harvesting are shown in Figures9(a)–9(d), respectively. Figure 9 shows stable periodic solutions of period 1 corresponding to the default parameter values. For varying initial conditions, all initial conditions lead to periodic solutions. Different initial conditions lead to the same periodic solution curve and therefore we have a stable limit cycle in the phase space. The time series for zooplankton initially has irregular fluctuations showing that the density of zooplankton is highly variable and eventually stabilizes into a periodic orbit after 40 months. This could be attributed to the initial phase in the lake when the fishery started to operate and then eventually stabilized after some time.Figure 9
Phase portrait showing the time series of (a) nutrients; (b) phytoplankton; (c) zooplankton; (d) Limnothrissa miodon for model (2) and assumed initial condition: N0=9.5,P0=6.5,Z0=3.5,L0=1.5 (time series with colour green); N0=10,P0=7,Z0=4,L0=2 (time series with colour purple); N0=10.5,P0=7.5,Z0=4.5,L0=2.5 (time series with colour magenta); using the default parameter.
(a)(b)(c)(d)Figures10(a)–10(b) show oscillatory dynamics in the phase portraits of the nutrient and phytoplankton; phytoplankton and zooplankton; zooplankton and Limnothrissa miodon; phytoplankton, zooplankton, and Limnothrissa miodon dynamics for system (2), respectively. The phase portraits show that the equilibrium point at E∗=0.134,503.68,0.0694,128.578 (point with colour yellow) for the autonomous model is enclosed by the periodic solutions of the nonautonomous model, but the system does not go to the equilibrium point.Figure 10
Phase portrait showing the dynamics of (a) phytoplankton and nutrients; (b) zooplankton and phytoplankton; (c) Limnothrissa miodon and zooplankton; (d) Limnothrissa miodon, zooplankton, and phytoplankton for model system (2) with assumed initial condition: N0=10,P0=7,Z0=4,L0=2 using the default parameter values.
(a)(b)(c)(d)


### 3.2.2. Effects of Nutrients and Fishing Effort for the Nonautonomous Model

For the set of the default parameters, Figures11(a) and 11(b) illustrate the effect of varying the inflow of nutrients parameter a in model (2) on the dynamics of Limnothrissa miodon with E=15250. Figures 11(c) and 11(d) illustrate the effect of varying fishing effort in model (2) on the dynamics of Limnothrissa miodon.Figure 11
(a) Time series of Limnothrissa miodon for varyinga for model system (2); (b) 12-month average plot of Limnothrissa miodon for varying a; (c) time series of Limnothrissa miodon for varying fishing effort for model system (2) with E=0, E=22875, and E=38125; (d) 12-month average plot of Limnothrissa miodon to fishing effort and with assumed initial condition: N0=10,P0=7,Z0=4,L0=2 using the default parameter values.
(a)(b)(c)(d)Simulation results from Figures11(a)–11(d) show that the abundance of Limnothrissa miodon is more closely related to the availability of nutrients than to harvesting. Marshall [42] and Paulsen [43] in their studies showed that the abundance of Limnothrissa miodon is mainly due to river inflow and the availability of food compared to harvesting.


## 3.2.1. Nonautonomous Model

For the set of the default parameters, the phase portraits showing the time series of nutrients, phytoplankton, zooplankton, and kapenta with harvesting are shown in Figures9(a)–9(d), respectively. Figure 9 shows stable periodic solutions of period 1 corresponding to the default parameter values. For varying initial conditions, all initial conditions lead to periodic solutions. Different initial conditions lead to the same periodic solution curve and therefore we have a stable limit cycle in the phase space. The time series for zooplankton initially has irregular fluctuations showing that the density of zooplankton is highly variable and eventually stabilizes into a periodic orbit after 40 months. This could be attributed to the initial phase in the lake when the fishery started to operate and then eventually stabilized after some time.Figure 9
Phase portrait showing the time series of (a) nutrients; (b) phytoplankton; (c) zooplankton; (d) Limnothrissa miodon for model (2) and assumed initial condition: N0=9.5,P0=6.5,Z0=3.5,L0=1.5 (time series with colour green); N0=10,P0=7,Z0=4,L0=2 (time series with colour purple); N0=10.5,P0=7.5,Z0=4.5,L0=2.5 (time series with colour magenta); using the default parameter.
(a)(b)(c)(d)Figures10(a)–10(b) show oscillatory dynamics in the phase portraits of the nutrient and phytoplankton; phytoplankton and zooplankton; zooplankton and Limnothrissa miodon; phytoplankton, zooplankton, and Limnothrissa miodon dynamics for system (2), respectively. The phase portraits show that the equilibrium point at E∗=0.134,503.68,0.0694,128.578 (point with colour yellow) for the autonomous model is enclosed by the periodic solutions of the nonautonomous model, but the system does not go to the equilibrium point.Figure 10
Phase portrait showing the dynamics of (a) phytoplankton and nutrients; (b) zooplankton and phytoplankton; (c) Limnothrissa miodon and zooplankton; (d) Limnothrissa miodon, zooplankton, and phytoplankton for model system (2) with assumed initial condition: N0=10,P0=7,Z0=4,L0=2 using the default parameter values.
(a)(b)(c)(d)


## 3.2.2. Effects of Nutrients and Fishing Effort for the Nonautonomous Model

For the set of the default parameters, Figures11(a) and 11(b) illustrate the effect of varying the inflow of nutrients parameter a in model (2) on the dynamics of Limnothrissa miodon with E=15250. Figures 11(c) and 11(d) illustrate the effect of varying fishing effort in model (2) on the dynamics of Limnothrissa miodon.Figure 11
(a) Time series of Limnothrissa miodon for varyinga for model system (2); (b) 12-month average plot of Limnothrissa miodon for varying a; (c) time series of Limnothrissa miodon for varying fishing effort for model system (2) with E=0, E=22875, and E=38125; (d) 12-month average plot of Limnothrissa miodon to fishing effort and with assumed initial condition: N0=10,P0=7,Z0=4,L0=2 using the default parameter values.
(a)(b)(c)(d)Simulation results from Figures11(a)–11(d) show that the abundance of Limnothrissa miodon is more closely related to the availability of nutrients than to harvesting. Marshall [42] and Paulsen [43] in their studies showed that the abundance of Limnothrissa miodon is mainly due to river inflow and the availability of food compared to harvesting.


## 4. Discussion

We developed and examined the dynamics of a mathematical model that includes nutrients, phytoplankton, zooplankton, and Limnothrissa miodon in this paper. The phytoplankton growth rate, phytoplankton mortality, phytoplankton grazing, zooplankton growth rate, zooplankton mortality, grazing on zooplankton, and Limnothrissa miodon mortality are assumed to be Holling type-I forms. Theoretical analysis, including the positivity and existence of model solutions (2), is investigated. We obtained the autonomous model’s critical points and examined their stabilities. The equilibrium points’ local and global stability conditions are established. The autonomous and nonautonomous models were numerically simulated. The study’s highlights are as follows:(i)
The application of a dynamical system to the food chain, nutrients⟶ phytoplankton ⟶ zooplankton ⟶ Limnothrissa miodon, is used to demonstrate the qualitative behavior of Limnothrissa miodon in the presence of harvesting. Given a set of parameter values and initial conditions, the phase portraits are damped oscillations for the autonomous model and show periodic behavior for the nonautonomous model.(ii)
For the autonomous model’s local stability analysis, the theoretical results agree with the numerical simulations in that the coexistence equilibrium is stable if certain conditions are satisfied. If the equilibriumE∗ is feasible, it is globally unstable. Numerical results show that the nonautonomous model in Lake Kariba has a stable limit cycle.(iii)
Data fitting results for the period 2015 to 2017 show that there are two notable peaks in April and August in Lake Kariba which correspond to the inflow of nutrients from river inflow and runoff and nutrient reloading, respectively.(iv)
The rate of nutrient inflow has a positive effect on the Limnothrissa miodon abundance. Numerical simulations of system (2) show that Limnothrissa miodon abundance is negatively affected by harvesting.(v)
Abundance of Limnothrissa miodon has been shown to vary seasonally in relation to food abundance in Lake Kariba.(vi)
Optimal control results from the autonomous model show 505 optimal fishing units, showing that currently there is overcapacity in Lake Kariba.(vii)
A maximum sustainable annual catch of 34669 tonnes per year in Lake Kariba is obtained for the autonomous model.(viii)
The Limnothrissa miodon will flourish in the lake for model system (2) as long as the harvesting rate is less than 34669 tonnes per year, with parameter values as shown in Table 2. As a result, the harvesting rate is critical to Lake Kariba’s productivity.(ix)
Simulation results show that kapenta abundance is more closely related to inflow of nutrients than to harvesting.


## 5. Conclusions

The basic kapenta model [28] is extended to include harvesting by fishing vessels. The kapenta fishery in Lake Kariba provides a significant source of income for many people in Zimbabwe and Zambia. As a result, fishery management is critical to the fishery’s long-term survival, as it is a major source of income for fishing cooperatives, wholesalers, retailers, and the Lake Kariba community. Many people use kapenta as a relish and a source of protein. Lake Kariba kapenta fishing is in jeopardy due to dwindling fish populations. Kapenta harvesting must be done in a sustainable manner so that the resource does not become extinct. As a result, mathematical modelling of the kapenta model with harvesting provided us with insight into the dynamics of the Lake Kariba kapenta fishery. We developed a mathematical model for kapenta harvesting with the goal of ensuring that sustainable levels of kapenta are maintained. Dynamical systems have not been used to investigate how kapenta fish populations in a lake are described and influenced by harvesting. We were able to qualitatively explain the impact of harvesting by developing a mathematical model and analyzing it. Once we fully comprehend the dynamics, we will be able to intervene in order to maintain sustainable levels of kapenta fish. Currently there is overcapacity in Lake Kariba, and the kapenta is being overexploited. Based on our optimal control results, we recommend 505 fishing vessels, 278 for the Zimbabwean side and 227 for the Zambian side of Lake Kariba.



---

*Source: 8504993-2022-05-05.xml*