# A Stage-Structured Fishery Model for African Catfish and Nile Tilapia Feeding on Two Food Resources with Harvesting

Authors: L. Nankinga; L. S. Luboobi; J. Y. T. Mugisha; B. Nannyonga; L. Carlsson
Journal: Journal of Applied Mathematics
 (2022)

---

## Abstract

In this paper, a fishery model for African catfish and Nile tilapia is formulated. This model is used to compare financial profit and biomass outtakes in a two-species system versus single species systems. We consider a stage-structured fish population model consisting of the aforementioned fish species together with two food resources. The model dynamics include cannibalism, predator-prey, feeding, reproduction, maturation, development, mortality, and harvesting. We prove consistency of the model in the sense that the solutions will stay bounded and nonnegative over time. Conditions for local stability of fish-free equilibrium point are established. The simulation results reveal asymptotically stable solutions with coexistence of African catfish, Nile tilapia, and two food resources. The major conclusion from our findings is that fisheries should culture both species to maximize the biomass outtake and financial profit.


---

## Body

## 1. Introduction

Fish is a source of animal protein, providing many key nutrients and calories that are needed for physical and mental development [1, 2]. Fish contribute 17 percent of animal protein and 7 percent of all proteins and are crucial for over 3 billion people in developing countries [3].In East Africa, the major species of fish caught are Nile perch (Lates niloticus), Nile tilapia (Oreochromis niloticus), African catfish (Clarias gariepinus), and silver fish (Rastrineobola argentea) to mention but a few. The fishery industry contributes to the gross domestic product (GDP) and has continued to be an important source of foreign exchange earnings through fish exports to the East African region and international markets [4]. Due to the importance of fish to so many people, the demand is very high and has resulted in overfishing and hence a decrease in fish population stock [5, 6]. Therefore, there is a need to develop sustainable fisheries management practices in East Africa.Mathematical modeling plays a significant role in the sustainable management of renewable resources. It describes phenomena leading to better prediction, management, and control techniques. Modeling and analysis of multiple species ecological problems were first done by Volterra [7]. Since then, many unstructured ecological models have been formulated and analyzed to study different phenomena [8–12]. Unstructured models, however, treat all individuals as identical and account for reproduction and death as life history processes that are important for population dynamics.Over the last four decades, there has been a trend towards increased realism of ecological theory. This trend has been supported by introducing physiologically structured population models that distinguish individuals based on physiological characteristics, e.g., size, age, weight, and length. These models include feeding, development, reproduction, and mortality life history processes [13–19].Stage-structured population models, which still carry a lot of interesting information, are simpler versions of physiologically structured population models and have increased in popularity [9, 19–24]. Field and laboratory studies have attributed several interesting features of population dynamics of cannibalistic populations, specifically to the occurrence of cannibalism [13, 19, 25–30]. Nevertheless, we have not encountered any mathematical models that consider several species including cannibalism.In this study, we formulate a stage-structured fishery model that resembles the feeding habits of African catfish and Nile tilapia to compare financial profit and biomass outtakes in a two-species system versus single-species systems. The paper is organized as follows: In Section2, we present the food and feeding habits of African catfish and Nile tilapia. In Section 3, we formulate a stage-structured fishery model with harvesting. In Section 4, we prove the basic properties of the model system (18). Section (5) gives the mathematical analysis of the model. In Section 6, we present the numerical simulations and results of the model equations. We end with a discussion of the results in Section 7.


## 2. Food and Feeding Habits of African Catfish and Nile Tilapia

African catfish (Clarias gariepinus) is a group of domesticated African catfishes of the family of Clariidae, air breathing catfishes. It is widely distributed in African and middle East fresh waters [31–33]. In Uganda, it is found in most natural water bodies, including swamps, streams, rivers, and lakes [34].Nile tilapia (Oreochromis niloticus) is a species of tilapia, a Cichlid fish native to Africa and Israel. As a food resource, the Nile tilapia is the most important fish species in the tropical and subtropical freshwater, often forming the basis of commercially fisheries in many African countries [35]. In Uganda, Nile tilapia is the second most important farmed aquaculture after African catfish [36]. It is native of lakes Albert, Edward, and George and introduced in the 1950s into Lakes Kyoga and Victoria [37].Food and feeding habits of fresh water fish species make up the foundation for the development of a management program of fish capture and culture. According to Otieno et al. [38], fish feed mainly on food items that can fit their mouth and what their stomach can digest. As fish grows, the stomach becomes larger, and its digestive system becomes more developed. African catfish is an opportunistic omnivorous feeder ingesting a wide variety of food items such as algae, macrophytes, zooplankton, insects, fish, detritus, amphibians, and sand grains [39, 40].Tesfahun [41] showed that African catfish preferred to feed on animals while plant items were ingested during pursuing other prey organisms. On top of some common resources, such as zoo plankton, smaller specimens tend to ingest more insects while larger specimens mainly prefer fish prey.Nile tilapia is omnivorous in its feeding habit with plant origin foods as the most dominant food item [41].


## 3. Model Formulation

In this section, we will develop our model system (18). The model we are considering is devoted to the interactions and dynamics of African catfish, Nile tilapia, and two food resources with harvesting. Two food resources are used to capture the difference in feeding habits. African catfish are classified into four stages: small juvenile X1, large juvenile X2, small adult X3, and large adult X4. Nile tilapia are composed of two stages: juvenile Y1 and adult Y2. The stages divide up the the population in body mass intervals. The two food resources are an animal-based food resource RA (e.g., zooplankton) and a plant-based food resource RP (e.g., phytoplankton). African catfish feed on both food resources. In addition, small adult African catfish feed on juvenile Nile tilapia and small juvenile African catfish; large adult African catfish feed on adult Nile tilapia and large juvenile African catfish. Nile tilapia only feed on the plant-based food resources.The variables used in the model are summarized in Table1.Table 1
Variables used in the stage-structured model for African catfish, Nile tilapia, and the two food resources.
VariablesDescriptionX1Biomass of small juvenile African catfish at timetX2Biomass of large juvenile African catfish at timetX3Biomass of small adult African catfish at timetX4Biomass of large adult African catfish at timetY1Biomass of juvenile Nile tilapia at timetY2Biomass of adult Nile tilapia at timetRABiomass of zooplankton (animal items) at timetRPBiomass of phytoplankton (plant items) at timetBasing on the food and feeding pattern of African catfish and Nile tilapia in Section2, the assumed feeding habit pattern is summarized in Table 2.Table 2
Food preferences of African catfish and Nile tilapia: the elements in the first column eat the elements in the first row in a preferred order, where it prefers the food with the lowest number. If no number exists, then it does not eat the corresponding food. The variable descriptions are given in Table1.
ConsumerPreyX1X2Y1Y2RARPX112X212X32134X42134Y11Y21The biomass of small juvenile African catfishX1 increases due to somatic growth at a rate ω1, reproduction at a rate ω4 from large adult African catfish, and reproduction at a rate qν3 from small adult African catfish, where q is the reproduction fraction and ν3 is the maturation rate from small adult African catfish to large adult African catfish. Its biomass decreases due to maturation into the large juvenile stage at a rate ν1, harvesting at a rate h1, cannibalism at a rate c1, and death due to natural causes at a rate μC. The reproduction fraction, q, is a parameter determining the fraction of how much of the energy production should be dedicated to reproduction and how much should be dedicated to maturation.The biomass of large juvenile African catfishX2 increases due to the maturation of small juvenile African catfish into the large juvenile stage at a rate ν1 and the somatic growth at a rate ω2. Its biomass decreases due to maturation of the large juvenile at a rate ν2, cannibalism at a rate c2, harvesting at a rate h2, and death due to natural causes at a rate μC.The biomass of small adult African catfish biomassX3 increases due to maturation of large juvenile African catfish into the small adult stage at a rate ν2 and somatic growth at a rate ω3. Its biomass decreases due to maturation of small adult African catfish at a rate 1−qν3 to large adult African catfish, reproduction at a rate qν3, harvesting at a rate h3, and natural death at a rate μC.The biomass of large adult African catfishX4 increases due to the maturation of small adult African catfish into the large adult stage at a rate 1−qν3. Its biomass decreases due to harvesting at a rate h4 and natural death at a rate μC. Biomass increase due to somatic growth ω4 is lost in reproduction.The biomass of juvenile Nile tilapiaY1 increases due to reproduction of adult Nile tilapia at a rate λ2 and somatic growth at a rate λ1. Its biomass decreases due to maturation into the adult stage at a rate γ, predation at a rate p1, harvesting at a rate f1, and death due to natural causes at a rate μT.The biomass of adult Nile tilapiaY2 increases due to maturation into the adult stage at a rate γ, and it decreases due to predation at a rate p2 harvesting at a rate f2 and natural death at a rate μT. Biomass increase due to somatic growth λ2 is lost in reproduction.According to De Roos et al. [21], in the absence of African catfish and Nile tilapia, the food resource RP is assumed to have semichemostat growth rate GPRP,
(1)GPRP=rPRPmax−RP,where rP is the intrinsic growth rate and RPmax is the carrying capacity of food resource, RP.While the food resourceRA follows the logistic growth rate GARA,
(2)GARA=rARA1−RARAmax,where rA is the intrinsic growth rate and RAmax is the carrying capacity of food resource RA.Due to external influx and spatial conditions such as coverage, not all food resources are available to the fish population, see, e.g., [18] for more details and discussions. That is, the fish cannot reach all food RA; thus, we introduce a threshold constant, Λ; the available animalistic food is defined by
(3)R^A=RA−Λ+≔max0,RA−Λ.This threshold forces the food resourceRA away from zero, thus eliminating the zero growth rate at the origin. A zero growth rate at the origin creates instability in the solution, as well as the possibility of extinction of that food resource. Since the plant food resource, RP, follows a semichemostat growth, it has a nonhorizontal tangent line at the origin; therefore, no analytical nor numerical instability occurs, and hence, we are using a zero threshold constant on the plant food resource. We note here that if we instead had used logistic growth rate models without threshold for the food resources, then cyclic solutions would most probably appear, for a deeper discussion of this and related concepts, see, e.g., [42].


### 3.1. Biological Assumptions

In our model, we make the following assumptions:(i)
African catfish are assumed to have birth size (in fact, the birth size in our model corresponds to the size at which the fish transform from larval into fry, since this is the stage when the fish starts utilizing the food resource)sb. The size of small juveniles is in the interval sb≤s<s1, the size of large juveniles is in the interval s1≤s<s2, the size of the small adults is in the interval s2≤s<sm, and the size of large adult African catfish is sm(ii)
Nile tilapia are assumed to have birth sizekb. The size of the juveniles is in the interval kb≤s<km, and the size of adults is km(iii)
Large adult African catfish of sizesm and adult Nile tilapia of size km do not grow but invest all energy in reproduction(iv)
Small adult African catfish uses the available energy for growth, maturation, and reproduction(v)
Small juvenile African catfish, large juvenile African catfish, and juvenile Nile tilapia use all available energy for growth and maturation(vi)
The growth rate and reproduction rate depend on food availability(vii)
Maturation rates depend on the net biomass production rates(viii)
Consumption of food follows the Holling type II functional response. The net biomass production rates are assumed to equal the balance between ingestion and mass-specific maintenance rates(ix)
At low food densities, the energy intake of juveniles and adults is insufficient to cover their maintenance requirements. As a result, they do not produce biomass. They only produce biomass when their food intake is sufficient to cover their maintenance demands(x)
The net biomass production rates for juveniles and adults are assumed to equal the balance between ingestion and maintenance rates for African catfish and Nile tilapia, respectively


### 3.2. Vital Rates

We adopt De Roos et al.’s [21] derivation of net biomass production rates of individuals in different stages and the maturation rate from one stage to another.The net biomass production rates of juvenile and adult African catfish are given by(4)ω1=ω2=σCIJA+1−AP−T+,(5)ω3=σCIAF1+1−F1A+1−AP−T+,(6)ω4=σCIAF2+1−F2A+1−AP−T+,where
(7)Fi=YiHC+Yi+1−YiHC+YiXiHC+Xi,i=1,2,A=R^AHC+R^A,P=RPHC+RP.Here,σC is the conversion efficiency of food ingested, IJ and IA are the maximum ingestion rates per unit biomass of juvenile and adult African catfish, respectively. T is the mass specific maintenance rate of African catfish, and HC is the half saturation food level of African catfish.The feeding rateFi is defined in such a way that, for large values of biomass, say Y1, we have F1≈Y1/HC+Y1 and for Y1<<HC it follows that F1≈X1/HC+X1. Note also that 0≤F1≤1.The net biomass production rates of juvenile and adult Nile tilapia are given by(8)λ1=σTψJRPHT+RP−E+,λ2=σTψARPHT+RP−E+.where σT is the conversion efficiency of food ingested by Nile tilapia and ψJ and ψA are the maximum ingestion rates per unit biomass of juvenile and adult Nile tilapia, respectively. HT is the half saturation constant of Nile tilapia, and E is the mass-specific maintenance rate of Nile tilapia.Natural mortality rates of African catfish and Nile tilapia are given byμC and μT, respectively.Cannibalistic mortality rates (the cannibalistic and predation mortality rates are in fact proportional to the biomass of the consumed population stage; in some literature, the biomass of the consumed population is included in these expressions) of small adult African catfish on small juvenile African catfish,c1, and large adult African catfish on large juvenile African catfish, c2, are given by
(9)ci=IAXi+2HC+Xi,i=1,2.Predation rates of small adult African catfish on juvenile Nile tilapia,p1, and adult African catfish on adult Nile tilapia, p2, are given by
(10)pi=IAXi+2HC+Yi,i=1,2.We use the maturation function(11)νω,μ,h,c,z=ω−μ+h+c1−z1−μ+h+c/ω,to define the following maturation rates for African catfish
(12)smalljuvenile⟶largejuvenileν1=νω1,μC,h1,c1,sbs1,largejuvenile⟶smalladultν2=νω2,μC,h2,c2,s1s2,smalladult⟶largeadultν3=νω3,μC,h3,0,s2sm,and for Nile tilapia
(13)juvenile⟶adultγ=νλ1,μT,f1,p1,kbkm,


### 3.3. The Schematic Diagram of the Model

Based on the model description and assumptions made in Section3, the dynamics of the system is represented in Figure 1.Figure 1
Schematic diagram of the dynamics of African catfish,Xi, Nile tilapia, Yi, and two food resources, RA and RP; solid arrows indicate growth, reproduction, mortality, and maturation from one stage to another whereas dashed connections show the flow of energy in the food-web. The remaining parameters in the graph are defined through Equations (4)–(13) (Tables 1 and 3–5).


### 3.4. Model Equations

From the assumptions and interrelation between the variables in Table1 and the parameters in Tables 3–5, as shown in Figure 1, the rates at which the biomass of small juvenile African catfish, large juvenile African catfish, small adult African catfish, large adult African catfish, juvenile Nile tilapia, adult Nile tilapia, and food resources changes are given by the following ordinary differential equations, using the expressions given in Equations (4)–(13)
(14)dX1dt=ω4X4+qν3X3+ω1X1,(15)−ν1+μC+h1+c1X1,(16)dX2dt=ν1X1+ω2X2−ν2+μC+h2+c2X2,(17)dX3dt=ν2X2+ω3X3−ν3+μC+h3X3,(18)dX4dt=1−qν3X3−μC+h4X4,(19)dY1dt=λ2Y2+λ1Y1−γ+μT+f1+p1Y1,(20)dY2dt=γY1−μT+f2+p2Y2,(21)dRAdt=rARA1−RARAmax−FAR^A,(22)dRPdt=rPRPmax−RP−FPRP,where
(23)R^A=RA−Λ+,(24)FA=IJHC+R^AX1+IJHC+R^AX2+1−F1IAHC+R^AX3+1−F2IAHC+R^AX4,(25)FP=1−R^AHC+R^AIJHC+RPX1+1−R^AHC+R^AIJHC+RPX2+1−F11−R^AHC+R^AIAHC+RPX3+1−F21−R^AHC+R^AIAHC+RPX4+ψJHT+RPY1+ψAHT+RPY2,subject to the initial conditions: X10≥0,(26)X20≥0,X30≥0,X40≥0,Y10≥0,Y20≥0,RA0≥0 and RP0≥0.Table 3
Parameters for the food resourcesRA and RP.
Subject and symbolDefinitionValueUnitSourcerAIntrinsic growth rate0.15day-1[57]RAmaxCarrying capacity4×10−3g/L[58]ΛFood threshold constant0.01×RAmaxg/LAssumedrPIntrinsic growth rate0.5day-1[56, 59]RPmaxCarrying capacity1.6×10−2g/LEstimated11Personal communication with biologists A. Larsson and J. Eriksson.Table 4
Parameters for Nile tilapia, denotedY in our model.
ParameterDefinitionValueUnitSourcekbSize at birth0.02g[60]kmMaturation size of adults176g[61]EMaintenance rate2.75×10−6day-1[48, 50]μTNatural mortality rate2.75×10−4day-1[52]HTHalf saturation constant1.26×10−2g/L[53, 54]ψJMaximum ingestion rate for juveniles2.75×10−2day-1[49, 51]ψAMaximum ingestion rate for adults2.33×10−2day-1[23]σTConversion efficiency rate for juveniles0.5—[50]Table 5
Parameters for African catfish, denotedX in our model description. MA is an abbreviation of model assumption.
ParameterDefinitionValueUnitSourcesbBirth size0.02g[55]s1Maturation size to large juveniles10gMAsmMaturation size to small adults190g[46, 62]smaxMaturation size to large adults6000gMAqReproduction fraction0.5—MAμCNatural mortality rate2.65×10−4day-1[52]TMaintenance rate2.65×10−6day-1[48, 50]HCHalf saturation constant1.38×10−2g/L[53, 54]IJMaximum ingestion rate for juveniles2.65×10−2day-1[49, 51]IAMaximum ingestion rate for adults2.25×10−2day-1[23]σCConversion efficiency rate0.5—[50]


## 3.1. Biological Assumptions

In our model, we make the following assumptions:(i)
African catfish are assumed to have birth size (in fact, the birth size in our model corresponds to the size at which the fish transform from larval into fry, since this is the stage when the fish starts utilizing the food resource)sb. The size of small juveniles is in the interval sb≤s<s1, the size of large juveniles is in the interval s1≤s<s2, the size of the small adults is in the interval s2≤s<sm, and the size of large adult African catfish is sm(ii)
Nile tilapia are assumed to have birth sizekb. The size of the juveniles is in the interval kb≤s<km, and the size of adults is km(iii)
Large adult African catfish of sizesm and adult Nile tilapia of size km do not grow but invest all energy in reproduction(iv)
Small adult African catfish uses the available energy for growth, maturation, and reproduction(v)
Small juvenile African catfish, large juvenile African catfish, and juvenile Nile tilapia use all available energy for growth and maturation(vi)
The growth rate and reproduction rate depend on food availability(vii)
Maturation rates depend on the net biomass production rates(viii)
Consumption of food follows the Holling type II functional response. The net biomass production rates are assumed to equal the balance between ingestion and mass-specific maintenance rates(ix)
At low food densities, the energy intake of juveniles and adults is insufficient to cover their maintenance requirements. As a result, they do not produce biomass. They only produce biomass when their food intake is sufficient to cover their maintenance demands(x)
The net biomass production rates for juveniles and adults are assumed to equal the balance between ingestion and maintenance rates for African catfish and Nile tilapia, respectively


## 3.2. Vital Rates

We adopt De Roos et al.’s [21] derivation of net biomass production rates of individuals in different stages and the maturation rate from one stage to another.The net biomass production rates of juvenile and adult African catfish are given by(4)ω1=ω2=σCIJA+1−AP−T+,(5)ω3=σCIAF1+1−F1A+1−AP−T+,(6)ω4=σCIAF2+1−F2A+1−AP−T+,where
(7)Fi=YiHC+Yi+1−YiHC+YiXiHC+Xi,i=1,2,A=R^AHC+R^A,P=RPHC+RP.Here,σC is the conversion efficiency of food ingested, IJ and IA are the maximum ingestion rates per unit biomass of juvenile and adult African catfish, respectively. T is the mass specific maintenance rate of African catfish, and HC is the half saturation food level of African catfish.The feeding rateFi is defined in such a way that, for large values of biomass, say Y1, we have F1≈Y1/HC+Y1 and for Y1<<HC it follows that F1≈X1/HC+X1. Note also that 0≤F1≤1.The net biomass production rates of juvenile and adult Nile tilapia are given by(8)λ1=σTψJRPHT+RP−E+,λ2=σTψARPHT+RP−E+.where σT is the conversion efficiency of food ingested by Nile tilapia and ψJ and ψA are the maximum ingestion rates per unit biomass of juvenile and adult Nile tilapia, respectively. HT is the half saturation constant of Nile tilapia, and E is the mass-specific maintenance rate of Nile tilapia.Natural mortality rates of African catfish and Nile tilapia are given byμC and μT, respectively.Cannibalistic mortality rates (the cannibalistic and predation mortality rates are in fact proportional to the biomass of the consumed population stage; in some literature, the biomass of the consumed population is included in these expressions) of small adult African catfish on small juvenile African catfish,c1, and large adult African catfish on large juvenile African catfish, c2, are given by
(9)ci=IAXi+2HC+Xi,i=1,2.Predation rates of small adult African catfish on juvenile Nile tilapia,p1, and adult African catfish on adult Nile tilapia, p2, are given by
(10)pi=IAXi+2HC+Yi,i=1,2.We use the maturation function(11)νω,μ,h,c,z=ω−μ+h+c1−z1−μ+h+c/ω,to define the following maturation rates for African catfish
(12)smalljuvenile⟶largejuvenileν1=νω1,μC,h1,c1,sbs1,largejuvenile⟶smalladultν2=νω2,μC,h2,c2,s1s2,smalladult⟶largeadultν3=νω3,μC,h3,0,s2sm,and for Nile tilapia
(13)juvenile⟶adultγ=νλ1,μT,f1,p1,kbkm,


## 3.3. The Schematic Diagram of the Model

Based on the model description and assumptions made in Section3, the dynamics of the system is represented in Figure 1.Figure 1
Schematic diagram of the dynamics of African catfish,Xi, Nile tilapia, Yi, and two food resources, RA and RP; solid arrows indicate growth, reproduction, mortality, and maturation from one stage to another whereas dashed connections show the flow of energy in the food-web. The remaining parameters in the graph are defined through Equations (4)–(13) (Tables 1 and 3–5).


## 3.4. Model Equations

From the assumptions and interrelation between the variables in Table1 and the parameters in Tables 3–5, as shown in Figure 1, the rates at which the biomass of small juvenile African catfish, large juvenile African catfish, small adult African catfish, large adult African catfish, juvenile Nile tilapia, adult Nile tilapia, and food resources changes are given by the following ordinary differential equations, using the expressions given in Equations (4)–(13)
(14)dX1dt=ω4X4+qν3X3+ω1X1,(15)−ν1+μC+h1+c1X1,(16)dX2dt=ν1X1+ω2X2−ν2+μC+h2+c2X2,(17)dX3dt=ν2X2+ω3X3−ν3+μC+h3X3,(18)dX4dt=1−qν3X3−μC+h4X4,(19)dY1dt=λ2Y2+λ1Y1−γ+μT+f1+p1Y1,(20)dY2dt=γY1−μT+f2+p2Y2,(21)dRAdt=rARA1−RARAmax−FAR^A,(22)dRPdt=rPRPmax−RP−FPRP,where
(23)R^A=RA−Λ+,(24)FA=IJHC+R^AX1+IJHC+R^AX2+1−F1IAHC+R^AX3+1−F2IAHC+R^AX4,(25)FP=1−R^AHC+R^AIJHC+RPX1+1−R^AHC+R^AIJHC+RPX2+1−F11−R^AHC+R^AIAHC+RPX3+1−F21−R^AHC+R^AIAHC+RPX4+ψJHT+RPY1+ψAHT+RPY2,subject to the initial conditions: X10≥0,(26)X20≥0,X30≥0,X40≥0,Y10≥0,Y20≥0,RA0≥0 and RP0≥0.Table 3
Parameters for the food resourcesRA and RP.
Subject and symbolDefinitionValueUnitSourcerAIntrinsic growth rate0.15day-1[57]RAmaxCarrying capacity4×10−3g/L[58]ΛFood threshold constant0.01×RAmaxg/LAssumedrPIntrinsic growth rate0.5day-1[56, 59]RPmaxCarrying capacity1.6×10−2g/LEstimated11Personal communication with biologists A. Larsson and J. Eriksson.Table 4
Parameters for Nile tilapia, denotedY in our model.
ParameterDefinitionValueUnitSourcekbSize at birth0.02g[60]kmMaturation size of adults176g[61]EMaintenance rate2.75×10−6day-1[48, 50]μTNatural mortality rate2.75×10−4day-1[52]HTHalf saturation constant1.26×10−2g/L[53, 54]ψJMaximum ingestion rate for juveniles2.75×10−2day-1[49, 51]ψAMaximum ingestion rate for adults2.33×10−2day-1[23]σTConversion efficiency rate for juveniles0.5—[50]Table 5
Parameters for African catfish, denotedX in our model description. MA is an abbreviation of model assumption.
ParameterDefinitionValueUnitSourcesbBirth size0.02g[55]s1Maturation size to large juveniles10gMAsmMaturation size to small adults190g[46, 62]smaxMaturation size to large adults6000gMAqReproduction fraction0.5—MAμCNatural mortality rate2.65×10−4day-1[52]TMaintenance rate2.65×10−6day-1[48, 50]HCHalf saturation constant1.38×10−2g/L[53, 54]IJMaximum ingestion rate for juveniles2.65×10−2day-1[49, 51]IAMaximum ingestion rate for adults2.25×10−2day-1[23]σCConversion efficiency rate0.5—[50]


## 4. Basic Properties of the Model

We need to show that the model is well-posed in a biologically feasible domain. Since the variablesX1t,X2t,X3t,X4t,Y1t,Y2t,RAt,RPt represent living species, nonnegativity ensures that they never become negative. Boundedness explains that there is a natural restriction to the growth of the population as a consequence of limited resources.


### 4.1. Positivity

We state a theorem which proves that all solutions of the model system (18) are nonnegative for t≥0. This proof could be made by considering cases, in the spirit of [43], but since our model would require 15 such cases, we instead prove the nonnegativity using decoupled ODEs.Let(27)ut=u1t,⋯,u8t,≔X1t,⋯,X4t,Y1t,Y2t,RAt,RPt,we say ut≥0 if all elements in ut are nonnegative.To prove the above statement, we need some preliminary theory. Letℝ+=0,∞.Proposition 1.
Letf:ℝ×ℝ+⟶ℝ+ and g:ℝ⟶ℝ be integrable. All solutions v:ℝ+⟶ℝ to the ordinary differential equations
(28)v′t=fv,t−gtvt,t>0,v0=v0,satisfies either case (1) vt>0 for all t>0 if v0>0 or case (2) vt≥0 for all t>0 if v0=0.Proof.
Rewriting the equation we get after multiplication with the integrating factorht=e∫0tgsds>0(29)ddthtvt=htft,v≥0,integrating both sides gives
(30)∫0tddshsvsds=htvt−h0v0≥0,or equivalently
(31)vt≥v0e−∫0tgsds.
From this inequality, the two statements above follows instantly.Theorem 2.
The solutions,u, to system (18) satisfy ut≥0 for all t≥0 if the initial data is nonnegative, i.e., u0≥0. If for some j∈1,⋯,6, the initial data uj0>0 then ujt>0 for all t>0.The ecological interpretation of this theorem is that, on the one hand, if one species is present in the model at the start of simulations, then its biomass is always positive. On the other hand, if we have high harvest pressure or high cannibalism, the species biomass might decline exponentially towards extinction.Proof (of Theorem2).
That all solutions are nonnegative to system (18) follows as a direct consequence of case (2) in Proposition 1 since all coupled ODEs in system (18), studied as single de-coupled ODEs can be written in the form given in Proposition 1 including the fulfillment of the hypotheses. In the same way, applying case (1) in Proposition 1 proves the second part of the statement.


### 4.2. Boundedness

Theorem 3.
All solutionsut of the model system (18) with initial conditions u0≥0 are uniformly bounded within the region Ω=ut≥0⊂ℝ+8.Proof.
LetNt be the total biomass of system (18) at time t such that
(32)Nt=X1t+X2t+X3t+X4t+Y1t+Y2t+RAt+RPt.
Taking the derivatives ofNt gives
(33)N′t=X1′t+X2′t+X3′t+X4′t+Y1′t+Y2′t+RA′t+RP′t,which along the solution path of the model system (18) yields
(34)N′t=rARA1−RARAmax−FAR^A+rPRPmax−RP−FPRP+ω1−μC−h1−c1X1+ω2−μC−h2−c2X2+ω3−μC−h3X3+ω4−μC−h4X4+λ1−μT−f1−p1Y1+λ2−μT−f2−p2Y2.
Substituting Equations (24) and (25) into Equation (34), we get (for convenience, we now use X5=Y1 and X6=Y2)
(35)N′t=rARA1−rARAmax+rPRAmax−RP−∑i=16πiXi,where
(36)π1=1−σCIJA+1−AP+T+h1+μC+c1,π2=1−σCIJA+1−AP+T+h2+μC+c2,π3=1−σCIAF1+1−F1A+1−AP+T+h3+μC,π4=1−σCIAF2+1−F2A+1−AP+T+h4+μC,π5=1−σTψJP.
For anyβ>0, we add βNt to both sides of Equation (35) to get
(37)dNtdt+βNt=RArA1−RARAmax+β+rPRPmax−rP−βRP−∑i=16πi−βXi.
Sinceπii=1,2,3,4,5,6 are nonnegative continuous functions, each function bounded by a positive constant and rP>0, for example, πi≥T1+μC+c1>0. Thus, we define β=minπ1,π2,π3,π4,π5,π6,rP. This now implies that
(38)dNtdt+βNt≤rPRPmax,(39)+RArA1−RARAmax+β.
We study the right hand side(40)fRA=rPRPmax+RArA1−RARAmax+β,which is a second degree polynomial in RA, and hence, attains its maximum value at
(41)RA=RAmaxrA+β2rA.
SubstitutingRA from Equation (41) into inequality (38) yields
(42)dNtdt+βNt≤rPRPmax+RAmaxrA+β24rA.
Using the product rule of derivatives on the left hand side and integrating from0 to t on both sides give the bound
(43)Nt≤Φβ1−e−βt+N0e−βt,where N0 is the total biomass at time t=0 and
(44)Φ=rPRPmax+RAmaxrA+β24rA>0.
By Theorem2, we get the bound
(45)0≤Nt≤maxN0,Φβ,for all t≥0.
Hence, the solutions of the model inℝ+8 are confined in the region
(46)Ω=X1,X2,X3,X4,Y1,Y2,RA,RP∈ℝ+8,
:0≤Nt≤maxN0,Φ/β. Therefore, we conclude that the model is biologically and mathematically well-posed in the region Ω.


## 4.1. Positivity

We state a theorem which proves that all solutions of the model system (18) are nonnegative for t≥0. This proof could be made by considering cases, in the spirit of [43], but since our model would require 15 such cases, we instead prove the nonnegativity using decoupled ODEs.Let(27)ut=u1t,⋯,u8t,≔X1t,⋯,X4t,Y1t,Y2t,RAt,RPt,we say ut≥0 if all elements in ut are nonnegative.To prove the above statement, we need some preliminary theory. Letℝ+=0,∞.Proposition 1.
Letf:ℝ×ℝ+⟶ℝ+ and g:ℝ⟶ℝ be integrable. All solutions v:ℝ+⟶ℝ to the ordinary differential equations
(28)v′t=fv,t−gtvt,t>0,v0=v0,satisfies either case (1) vt>0 for all t>0 if v0>0 or case (2) vt≥0 for all t>0 if v0=0.Proof.
Rewriting the equation we get after multiplication with the integrating factorht=e∫0tgsds>0(29)ddthtvt=htft,v≥0,integrating both sides gives
(30)∫0tddshsvsds=htvt−h0v0≥0,or equivalently
(31)vt≥v0e−∫0tgsds.
From this inequality, the two statements above follows instantly.Theorem 2.
The solutions,u, to system (18) satisfy ut≥0 for all t≥0 if the initial data is nonnegative, i.e., u0≥0. If for some j∈1,⋯,6, the initial data uj0>0 then ujt>0 for all t>0.The ecological interpretation of this theorem is that, on the one hand, if one species is present in the model at the start of simulations, then its biomass is always positive. On the other hand, if we have high harvest pressure or high cannibalism, the species biomass might decline exponentially towards extinction.Proof (of Theorem2).
That all solutions are nonnegative to system (18) follows as a direct consequence of case (2) in Proposition 1 since all coupled ODEs in system (18), studied as single de-coupled ODEs can be written in the form given in Proposition 1 including the fulfillment of the hypotheses. In the same way, applying case (1) in Proposition 1 proves the second part of the statement.


## 4.2. Boundedness

Theorem 3.
All solutionsut of the model system (18) with initial conditions u0≥0 are uniformly bounded within the region Ω=ut≥0⊂ℝ+8.Proof.
LetNt be the total biomass of system (18) at time t such that
(32)Nt=X1t+X2t+X3t+X4t+Y1t+Y2t+RAt+RPt.
Taking the derivatives ofNt gives
(33)N′t=X1′t+X2′t+X3′t+X4′t+Y1′t+Y2′t+RA′t+RP′t,which along the solution path of the model system (18) yields
(34)N′t=rARA1−RARAmax−FAR^A+rPRPmax−RP−FPRP+ω1−μC−h1−c1X1+ω2−μC−h2−c2X2+ω3−μC−h3X3+ω4−μC−h4X4+λ1−μT−f1−p1Y1+λ2−μT−f2−p2Y2.
Substituting Equations (24) and (25) into Equation (34), we get (for convenience, we now use X5=Y1 and X6=Y2)
(35)N′t=rARA1−rARAmax+rPRAmax−RP−∑i=16πiXi,where
(36)π1=1−σCIJA+1−AP+T+h1+μC+c1,π2=1−σCIJA+1−AP+T+h2+μC+c2,π3=1−σCIAF1+1−F1A+1−AP+T+h3+μC,π4=1−σCIAF2+1−F2A+1−AP+T+h4+μC,π5=1−σTψJP.
For anyβ>0, we add βNt to both sides of Equation (35) to get
(37)dNtdt+βNt=RArA1−RARAmax+β+rPRPmax−rP−βRP−∑i=16πi−βXi.
Sinceπii=1,2,3,4,5,6 are nonnegative continuous functions, each function bounded by a positive constant and rP>0, for example, πi≥T1+μC+c1>0. Thus, we define β=minπ1,π2,π3,π4,π5,π6,rP. This now implies that
(38)dNtdt+βNt≤rPRPmax,(39)+RArA1−RARAmax+β.
We study the right hand side(40)fRA=rPRPmax+RArA1−RARAmax+β,which is a second degree polynomial in RA, and hence, attains its maximum value at
(41)RA=RAmaxrA+β2rA.
SubstitutingRA from Equation (41) into inequality (38) yields
(42)dNtdt+βNt≤rPRPmax+RAmaxrA+β24rA.
Using the product rule of derivatives on the left hand side and integrating from0 to t on both sides give the bound
(43)Nt≤Φβ1−e−βt+N0e−βt,where N0 is the total biomass at time t=0 and
(44)Φ=rPRPmax+RAmaxrA+β24rA>0.
By Theorem2, we get the bound
(45)0≤Nt≤maxN0,Φβ,for all t≥0.
Hence, the solutions of the model inℝ+8 are confined in the region
(46)Ω=X1,X2,X3,X4,Y1,Y2,RA,RP∈ℝ+8,
:0≤Nt≤maxN0,Φ/β. Therefore, we conclude that the model is biologically and mathematically well-posed in the region Ω.


## 5. Analysis of the Model

### 5.1. Existence of Equilibrium Points

The points at which the derivatives of the system (18) are zero will be referred to as equilibrium points or steady state solutions. The system has six possible feasible equilibrium points. The vector below indicates the values of X1,X2,X3,X4,Y1,Y2,RA,RP.E1 The fish-free with plant based food resource equilibrium point 0,0,0,0,0,0,0,RPmax.E2 The fish-free with two food resources equilibrium point 0,0,0,0,0,0,RAmax,RPmax.E3 The African catfish free equilibrium point 0,0,0,0,Y1a°,Y2a°,0,RPa°.E4 The African catfish free equilibrium point 0,0,0,0,Y1a°a°,Y2a°a°,RAa°a°,RPa°a°.E5 The Nile tilapia free equilibrium point X1††,X2††,X3††,X4††,0,0,RA††,RP††.E6 The coexistence equilibrium point X1∗∗,X2∗∗,X3∗∗,X4∗∗,Y1∗∗,Y2∗∗,RA∗∗,RP∗∗.


### 5.2. Local Stability Analysis

In this subsection, we investigate the local stability of each of the equilibrium solutions of the model system (18) obtained by calculating the Jacobian matrix J corresponding to each of the equilibrium points.Theorem 4.
The equilibrium points with zero animalistic food resourceE1 and E3 are all unstable.Proof.
The Jacobian matrix evaluated at equilibrium pointE10,0,0,0,0,0,0,RPmax is given by
(47)JE1=b110qν3ω40000ν1b220000000ν2b330000000b43b4400000000b55λ2000000γb6600000000rA00000000−rP,where
b11=ω1−ν1+μC+h1, b22=ω2−ν2+μC+h2, b33=ω3−ν3+μC+h3, b43=1−qν3, b44=−μC+h4, b55=λ1−γ+μT+f1, and b66=−μT+f2. We observe that the seventh and eighth columns contain only diagonal terms which form two eigenvalues rA and −rP out of eight eigenvalues. Since rA is positive, we thus have a nonstable equilibrium point.Similarly, using the method above, the equilibrium pointE3 is unstable.Theorem 5.
The fish-free equilibrium pointE20,0,0,0,0,0,RAmax,RPmax is locally asymptotically stable if RC<1 and RT<1 and unstable otherwise.Proof.
The local stability of the equilibrium pointE2 is investigated using the eigenvalues or trace determinant criteria of the Jacobian matrix. An equilibrium point is locally asymptotically stable if the Jacobian matrix evaluated at that point has negative eigenvalues or has a negative trace and a positive determinant.
The Jacobian matrix evaluated at equilibrium pointE20,0,0,0,0,0,RAmax,RPmax is given by
(48)JE2=a110a13a140000a21a220000000a32a330000000a34a4400000000a55a56000000a65a6600000000a7700000000a88where
a11=ω1−ν1+μC+h1, a13=qν3, a14=ω4, a21=ν1, a22=ω2−ν2+μC+h2, a32=ν2, a33=ω3−ν3+μC+h3, a43=1−qν3, a44=−μC+h4, a55=λ1−γ+μT+f1, a56=λ2, J65=γ, a66=−μT+f2, a77=−rA, and a88=−rP.
Observe that if the harvesting is high enough, thenaii<0 for i=1,2,⋯,6. We need to show that all eigenvalues of JE2 are negative. As the seventh and eighth columns contain only diagonal terms which form the two eigenvalues −rA and −rP, the other six eigenvalues can be obtained from the submatrix J1E2 formed by excluding the seventh and eighth rows and columns of JE2. Hence, we have
(49)J1E2=a110a13a1400a21a2200000a32a3300000a43a44000000a55a540000a65a66.
The eigenvalues of the matrixJ1E2 are the roots of the characteristic equation detJ1E2−λI=0. That is
(50)deta11−λ0a13a1400a21a22−λ00000a32a33−λ00000a43a44−λ000000a55−λa560000a65a66−λ
=0, which is equal to solving (since this is a block diagonal matrix)detAdetB=0 where
(51)A=a11−λ0a13a14a21a22−λ000a32a33−λ000a43a44−λ,B=a55−λa56a65a66−λ.
From(52)deta55−λa56a65a66−λ=0.
We get a second degree characteristic equation given by(53)a55−λa66−λ−a56a65=0,that is
(54)λ2−a55+a66λ+a55a66−a56a65=0.
Using the Routh-Hurwitz criterion, the solutions have negative real parts if(55)−a55+a66>0,a55a66−a56a65>0,
Now(56)−a55+a66>0,is equivalent to
(57)γ+2μT+f1−λ1+f2>0.
That is, this inequality is satisfied if(58)f1+f2>λ1−γ−2μT,a55a66−a56a65>0,is equivalent to
(59)γ+μT+f1−λ1μT+f2−λ2γ>0.
That is, iff1>λ1−γ−μT, we get
(60)1−λ2γγ+μT+f1−λ1μT+f2>0.
This inequality is satisfied if and only if(61)RT=def.λ2γγ+μT+f1−λ1μT+f2<1.
Here,RT is the recovery potential of Nile tilapia based on growth in biomass.
Considering(62)deta11−λ0a13a14a21a22−λ000a32a33−λ000a43a44−λ=0.
We get a fourth degree characteristic equation given by(63)A0λ4+A1λ3+A2λ2+A3λ+A4=0,where
(64)A0=1,A1=−a11+a22+a33+a44,A2=a33a22+a33a11+a44a33+a44a22+a44a11+a22a11,A3=−a13a21a32+a44a33a22+a44a33a11+a33a22a11+a44a22a11,A4=−a14a21a32a43+a44a13a21a32+a44a33a22a11.
Using Routh-Hurwitz criteria [44, 45], the roots of the characteristic Equation (63) are negative or have negative real parts if
(65)A1>0,A1A2−A0A3>0,A1A2−A0A3A3−A12A4>0,A4>0.
Criteria 1. That is,A1>0 is true, provided a11+a22+a33+a44<0, which is true if the harvesting rates are high enough.
Criteria 2. We have to show thatA1A2−A0A3>0; that is, we need to check if the left hand side is greater than zero; the left hand side is
(66)A1A2−A0A3=−a11+a22+a33+a44,(67)a11a22+a11a33+a11a44+a22a33+a22a44+a33a44−1−a11a22a33−a21a13a32−a11a22a44−a11a33a44−a22a33a44=a112a22−a112a33−a112a44−a11a222−4a11a22a33−2a11a22a44−a11a332−2a11a33a44−a11a442−a222a33−a222a44−a22a332−2a22a33a44−a22a442−a332a44−a33a442+a21a13a32.
All terms in the right hand side of Equation (67) are greater than zero (since aii<0,i=1,2,3,4) with the possible exception of +a21a13a32. From the definitions, +a21a13a32 is equivalent to qν1ν2ν3 which is also greater than zero, hence showing that A1A2−A0A3>0.
Criteria 3. The analytical proof of(68)A1A2−A0A3A3−A12A4>0,is very hard. Instead, we observe the graphs in Figure 2 where they show that this criterion is satisfied, since at the harvesting rate zero, the graph of A1A2−A0A3A3−A12A4 against harvesting rate starts at a positive value and is increasing. Figure 2(a) shows the values of this expression for a larger set of harvesting rates. Figure 2(b) shows that this expression is positive even for small values of the harvesting rate.
Finally, Criteria 4:A4>0 is the same as
(69)−a14a21a32a43+a44a13a21a32+a44a33a22a11>0.
This is equivalent to(70)RC=def.MN<1,where
(71)M=qμC+h4+1−qω4μC+h4,(72)N=ν1ν2ν3ν1+μC+h1−ω1ν2+μC+h2−ω2ν3+μC+h3−ω3.
The functionRC can be viewed as the recovery potential of African catfish based on growth in biomass. The recovery potential represents the net generational biomass production (per unit body mass) in a virgin environment [24]. Thus, by the Routh-Hurwitz criterion, the fish-free equilibrium is asymptotically locally stable iff RC<1 and RT<1, defined by Equations (61) and (71) which completes the proof.Figure 2
Graphical representation of Routh-Hurwitz Criteria 3, where they–axis represents the expression A1A2−A0A3A3−A12A4 and the x–axis represents the harvesting rate. Other parameter values used can be found in Tables 3–5.
(a)(b)The Jacobian matrices evaluated at the equilibrium pointsE4, E5, and E6 lead to characteristic polynomials of degree eight whose eigenvalues are too complex to be obtained analytically. Therefore, we revert to numerical simulations, see Section 6.


## 5.1. Existence of Equilibrium Points

The points at which the derivatives of the system (18) are zero will be referred to as equilibrium points or steady state solutions. The system has six possible feasible equilibrium points. The vector below indicates the values of X1,X2,X3,X4,Y1,Y2,RA,RP.E1 The fish-free with plant based food resource equilibrium point 0,0,0,0,0,0,0,RPmax.E2 The fish-free with two food resources equilibrium point 0,0,0,0,0,0,RAmax,RPmax.E3 The African catfish free equilibrium point 0,0,0,0,Y1a°,Y2a°,0,RPa°.E4 The African catfish free equilibrium point 0,0,0,0,Y1a°a°,Y2a°a°,RAa°a°,RPa°a°.E5 The Nile tilapia free equilibrium point X1††,X2††,X3††,X4††,0,0,RA††,RP††.E6 The coexistence equilibrium point X1∗∗,X2∗∗,X3∗∗,X4∗∗,Y1∗∗,Y2∗∗,RA∗∗,RP∗∗.


## 5.2. Local Stability Analysis

In this subsection, we investigate the local stability of each of the equilibrium solutions of the model system (18) obtained by calculating the Jacobian matrix J corresponding to each of the equilibrium points.Theorem 4.
The equilibrium points with zero animalistic food resourceE1 and E3 are all unstable.Proof.
The Jacobian matrix evaluated at equilibrium pointE10,0,0,0,0,0,0,RPmax is given by
(47)JE1=b110qν3ω40000ν1b220000000ν2b330000000b43b4400000000b55λ2000000γb6600000000rA00000000−rP,where
b11=ω1−ν1+μC+h1, b22=ω2−ν2+μC+h2, b33=ω3−ν3+μC+h3, b43=1−qν3, b44=−μC+h4, b55=λ1−γ+μT+f1, and b66=−μT+f2. We observe that the seventh and eighth columns contain only diagonal terms which form two eigenvalues rA and −rP out of eight eigenvalues. Since rA is positive, we thus have a nonstable equilibrium point.Similarly, using the method above, the equilibrium pointE3 is unstable.Theorem 5.
The fish-free equilibrium pointE20,0,0,0,0,0,RAmax,RPmax is locally asymptotically stable if RC<1 and RT<1 and unstable otherwise.Proof.
The local stability of the equilibrium pointE2 is investigated using the eigenvalues or trace determinant criteria of the Jacobian matrix. An equilibrium point is locally asymptotically stable if the Jacobian matrix evaluated at that point has negative eigenvalues or has a negative trace and a positive determinant.
The Jacobian matrix evaluated at equilibrium pointE20,0,0,0,0,0,RAmax,RPmax is given by
(48)JE2=a110a13a140000a21a220000000a32a330000000a34a4400000000a55a56000000a65a6600000000a7700000000a88where
a11=ω1−ν1+μC+h1, a13=qν3, a14=ω4, a21=ν1, a22=ω2−ν2+μC+h2, a32=ν2, a33=ω3−ν3+μC+h3, a43=1−qν3, a44=−μC+h4, a55=λ1−γ+μT+f1, a56=λ2, J65=γ, a66=−μT+f2, a77=−rA, and a88=−rP.
Observe that if the harvesting is high enough, thenaii<0 for i=1,2,⋯,6. We need to show that all eigenvalues of JE2 are negative. As the seventh and eighth columns contain only diagonal terms which form the two eigenvalues −rA and −rP, the other six eigenvalues can be obtained from the submatrix J1E2 formed by excluding the seventh and eighth rows and columns of JE2. Hence, we have
(49)J1E2=a110a13a1400a21a2200000a32a3300000a43a44000000a55a540000a65a66.
The eigenvalues of the matrixJ1E2 are the roots of the characteristic equation detJ1E2−λI=0. That is
(50)deta11−λ0a13a1400a21a22−λ00000a32a33−λ00000a43a44−λ000000a55−λa560000a65a66−λ
=0, which is equal to solving (since this is a block diagonal matrix)detAdetB=0 where
(51)A=a11−λ0a13a14a21a22−λ000a32a33−λ000a43a44−λ,B=a55−λa56a65a66−λ.
From(52)deta55−λa56a65a66−λ=0.
We get a second degree characteristic equation given by(53)a55−λa66−λ−a56a65=0,that is
(54)λ2−a55+a66λ+a55a66−a56a65=0.
Using the Routh-Hurwitz criterion, the solutions have negative real parts if(55)−a55+a66>0,a55a66−a56a65>0,
Now(56)−a55+a66>0,is equivalent to
(57)γ+2μT+f1−λ1+f2>0.
That is, this inequality is satisfied if(58)f1+f2>λ1−γ−2μT,a55a66−a56a65>0,is equivalent to
(59)γ+μT+f1−λ1μT+f2−λ2γ>0.
That is, iff1>λ1−γ−μT, we get
(60)1−λ2γγ+μT+f1−λ1μT+f2>0.
This inequality is satisfied if and only if(61)RT=def.λ2γγ+μT+f1−λ1μT+f2<1.
Here,RT is the recovery potential of Nile tilapia based on growth in biomass.
Considering(62)deta11−λ0a13a14a21a22−λ000a32a33−λ000a43a44−λ=0.
We get a fourth degree characteristic equation given by(63)A0λ4+A1λ3+A2λ2+A3λ+A4=0,where
(64)A0=1,A1=−a11+a22+a33+a44,A2=a33a22+a33a11+a44a33+a44a22+a44a11+a22a11,A3=−a13a21a32+a44a33a22+a44a33a11+a33a22a11+a44a22a11,A4=−a14a21a32a43+a44a13a21a32+a44a33a22a11.
Using Routh-Hurwitz criteria [44, 45], the roots of the characteristic Equation (63) are negative or have negative real parts if
(65)A1>0,A1A2−A0A3>0,A1A2−A0A3A3−A12A4>0,A4>0.
Criteria 1. That is,A1>0 is true, provided a11+a22+a33+a44<0, which is true if the harvesting rates are high enough.
Criteria 2. We have to show thatA1A2−A0A3>0; that is, we need to check if the left hand side is greater than zero; the left hand side is
(66)A1A2−A0A3=−a11+a22+a33+a44,(67)a11a22+a11a33+a11a44+a22a33+a22a44+a33a44−1−a11a22a33−a21a13a32−a11a22a44−a11a33a44−a22a33a44=a112a22−a112a33−a112a44−a11a222−4a11a22a33−2a11a22a44−a11a332−2a11a33a44−a11a442−a222a33−a222a44−a22a332−2a22a33a44−a22a442−a332a44−a33a442+a21a13a32.
All terms in the right hand side of Equation (67) are greater than zero (since aii<0,i=1,2,3,4) with the possible exception of +a21a13a32. From the definitions, +a21a13a32 is equivalent to qν1ν2ν3 which is also greater than zero, hence showing that A1A2−A0A3>0.
Criteria 3. The analytical proof of(68)A1A2−A0A3A3−A12A4>0,is very hard. Instead, we observe the graphs in Figure 2 where they show that this criterion is satisfied, since at the harvesting rate zero, the graph of A1A2−A0A3A3−A12A4 against harvesting rate starts at a positive value and is increasing. Figure 2(a) shows the values of this expression for a larger set of harvesting rates. Figure 2(b) shows that this expression is positive even for small values of the harvesting rate.
Finally, Criteria 4:A4>0 is the same as
(69)−a14a21a32a43+a44a13a21a32+a44a33a22a11>0.
This is equivalent to(70)RC=def.MN<1,where
(71)M=qμC+h4+1−qω4μC+h4,(72)N=ν1ν2ν3ν1+μC+h1−ω1ν2+μC+h2−ω2ν3+μC+h3−ω3.
The functionRC can be viewed as the recovery potential of African catfish based on growth in biomass. The recovery potential represents the net generational biomass production (per unit body mass) in a virgin environment [24]. Thus, by the Routh-Hurwitz criterion, the fish-free equilibrium is asymptotically locally stable iff RC<1 and RT<1, defined by Equations (61) and (71) which completes the proof.Figure 2
Graphical representation of Routh-Hurwitz Criteria 3, where they–axis represents the expression A1A2−A0A3A3−A12A4 and the x–axis represents the harvesting rate. Other parameter values used can be found in Tables 3–5.
(a)(b)The Jacobian matrices evaluated at the equilibrium pointsE4, E5, and E6 lead to characteristic polynomials of degree eight whose eigenvalues are too complex to be obtained analytically. Therefore, we revert to numerical simulations, see Section 6.


## 6. Numerical Simulations

In this section, we carry out numerical simulations using MATLAB to investigate coexistence and stability of African catfish, Nile tilapia, and the two types of food resources at steady state. We further investigate proportional harvesting in a system that has both African catfish and Nile tilapia, only African catfish, and only Nile tilapia, respectively, in terms of biomass yield and financial profits at steady state.


### 6.1. Model Parameterization

The parameter values are obtained from the literature review when possible. The maturation size of African catfish is 28-30 cm in total body length [46], and the maturation size of Nile tilapia is 20.8 cm [47].We used the standard quarter-power scaling law [48–50] to estimate: the maximum ingestion rate with proportionality constant for the maximum ingestion rate 0.01 [49, 51], maintenance rate with proportionality constant 0.01 [48, 50], and mortality rate with proportionality constant 0.001 [52]. Following Lessmark [53] and Byström and García-Berthou [54], we used uniform half saturation constants for Nile tilapia and African catfish, evaluated by the mean weight of, respectively, species.In the model, the birth size of African catfish is 0.02 grams [55]. The maturation size of large juveniles and large adults is assumed to be 30 cm and 100 cm, respectively. Since our model is based on weight, we use the following relation between length, L (cm), and weight, W (grams)
(73)W=c1Lc2,where c1=0.002 and c2=3.232 for the African catfish [10] and c1=0.019 and c2=3.01 for the Nile tilapia [56]. Conversion efficiency rate is assumed to be independent of food type and fish species and is set to 0.5 [49, 50]. The maximum ingestion rates of adult individuals are assumed to be 0.85 times the ingestion rates of juveniles for both species [23]. Finally, we assume an equal natural mortality rate at all stages.The parameter values used in the model are summarized in Tables3–5.


### 6.2. Numerical Stability Analysis ofE4, E5, and E6

The remaining equilibrium points to investigate the stability areE4, E5, and E6.The stability of equilibrium pointE40,0,0,0,Y1a°a°,Y2a°a°,RAa°a°,RPa°a° occurs when 0.51≤hi≤0.77, i=1,2,3,4. We tested the stability by running the model with different harvesting rates h=0.77,0.77,0.77,0.77,f1,f2, using one hundred different initial values of the form X1,X2,X3,X4,Y1,Y2,RA,RP, taken randomly from the parameter space. In all these simulations, we ended up with the equilibrium solution E4; we therefore conclude that E4 is a stable equilibrium point provided that the harvesting rate of the catfish population is high enough and f1 and f2 are smaller than 0.005.In an analogous manner, we established the stability ofE5 when both f1 and f2 exceed 0.01 and hi≤0.005 for i=1,2,3,4.Finally, when all harvesting rates are low, the simulations always converge to the equilibrium pointE6.


### 6.3. Simulation Results

Our simulations are stable in the sense that when positive initial values are used, the solutions approach the positive steady state (all simulations shown in other figures in this paper have reached steady state) as indicated by Figure3 which shows the dynamics of the solutions to system (18). This signifies coexistence and stability of African catfish, Nile tilapia, and the two different food resources at steady state.Figure 3
Dynamics of African catfish, Nile tilapia, and two different food resources. They–axis shows the biomass in grams per litre (in a logarithmic scale) for each fish compartment together with the food resources. The initial values are chosen positive; all harvest rates are set to 0.005 per day; other parameter values can be found in Tables 3–5.Recently, Lundström et al. [23] adopted the model introduced by De Roos et al. [21], in which Lundström et al. compared the properties between a stage-structured fish population model with an age-structured population model. They investigated how harvesting strategies that qualify for pretty good yield can account for conservation and concluded that equal harvesting rate of juveniles and adults is always a good strategy. We present some of these properties, using equal harvesting in all stages in our model: Figure 4 for African catfish and Nile tilapia and Figure 5 Nile tilapia only.Figure 4
Simulation results, using system (18), including both African catfish and Nile tilapia. The figures show the dependence of impact on size structure (red curve), impact on biomass (blue curve), total biomass (green curve), and yield (black curve) with respect to equal harvesting rates, i.e., h⟶=h,h,h,h,h,h. The initial values are chosen positive, and other parameter values can be found in Tables 3–5.Figure 5
Simulation results, using system (18), including only the Nile tilapia. The figures show the dependence of impact on size structure (red curve), impact on biomass (blue curve), total biomass (green curve), and yield (black curve) with respect to equal harvesting rates, i.e., h⟶=h,h. The initial values are chosen positive, and other parameter values can be found in Tables 3–5.One result by Lundström et al. [23] is that when studying the impact of harvesting on population size structure, they found a positive change in the fraction of juveniles compared to adults in the population. In contrast to this, we find that the impact of population size structure may have a negative dependence on the harvest rate, see Figure 4. This negative impact occurs because African catfish is cannibalistic and predates on Nile tilapia, in that adults feed on juveniles, and this reduces the fraction of juveniles in the system. When we simulate the single fish species Nile tilapia (see Figure 5), our model shows similar results as in [23]; in particular, the impact on size structure is positively dependent on the harvesting rate, which signifies that increasing the uniform harvesting rate there is an increase in the fraction of juveniles.In fisheries, one can use different mesh sizes in the nets; we therefore have modeled system (18) under a specific type of harvesting rate. In Figures 6 and 7, we consider unbalanced harvesting rates on juveniles and adults for both African catfish and Nile tilapia, i.e., harvesting pressure h⟶=1−ph,1−ph,ph,ph,1−ph,ph, where the proportion 0≤p≤1. The figures display the biomass yield (Figure 6) and financial profit (Figure 7) with respect to different harvesting rates, h, and proportions, p. The maximum sustainable biomass yield and maximum sustainable financial profits are attained when p is large; that is, we mainly harvest adults.Figure 6
Simulation results, using system (18), including both African catfish and Nile tilapia with proportional harvesting rates, i.e., h⟶p=1−ph,1−ph,ph,ph,1−ph,ph where 0≤p≤1. The initial values are chosen positive, and other parameter values can be found in Tables 3–5.Figure 7
Simulation results, using system (18), including both African catfish and Nile tilapia with proportional harvesting rates, i.e., h⟶p=1−ph,1−ph,ph,ph,1−ph,ph where 0≤p≤1. The initial values are chosen positive, and other parameter values can be found in Tables 3–5.In fisheries, one can harvest the fish stock by catch and release, for example, with a bottle trap or a fyke net; we therefore model system (18) with variable harvesting rates on all fish compartments. Table 6 shows the results after using pattern search to find optimal harvest rates that maximize biomass yield and financial profit for both species in the system, only Nile tilapia in the system, and only African catfish in the system. Our finding is that harvesting in a system that has both Nile tilapia and African catfish gives the highest biomass yield and financial profit (the optimal values for biomass yield and financial profit in Table 6 are found at different harvesting rates), compared to harvesting in a system having Nile tilapia only or African catfish only.Table 6
Simulation results using optimal harvest rates found in MATLAB, using patternSearch. The columns present the maximum sustainable biomass yield (grams per litre and day) and maximum sustainable financial profit (USD per litre and day). The rows indicate the fish stock scenarios.
Maximum biomass yieldMaximum financial profitOnly Nile tilapia0.00240.0074Only African catfish0.00140.0057Both species0.00270.0081


## 6.1. Model Parameterization

The parameter values are obtained from the literature review when possible. The maturation size of African catfish is 28-30 cm in total body length [46], and the maturation size of Nile tilapia is 20.8 cm [47].We used the standard quarter-power scaling law [48–50] to estimate: the maximum ingestion rate with proportionality constant for the maximum ingestion rate 0.01 [49, 51], maintenance rate with proportionality constant 0.01 [48, 50], and mortality rate with proportionality constant 0.001 [52]. Following Lessmark [53] and Byström and García-Berthou [54], we used uniform half saturation constants for Nile tilapia and African catfish, evaluated by the mean weight of, respectively, species.In the model, the birth size of African catfish is 0.02 grams [55]. The maturation size of large juveniles and large adults is assumed to be 30 cm and 100 cm, respectively. Since our model is based on weight, we use the following relation between length, L (cm), and weight, W (grams)
(73)W=c1Lc2,where c1=0.002 and c2=3.232 for the African catfish [10] and c1=0.019 and c2=3.01 for the Nile tilapia [56]. Conversion efficiency rate is assumed to be independent of food type and fish species and is set to 0.5 [49, 50]. The maximum ingestion rates of adult individuals are assumed to be 0.85 times the ingestion rates of juveniles for both species [23]. Finally, we assume an equal natural mortality rate at all stages.The parameter values used in the model are summarized in Tables3–5.


## 6.2. Numerical Stability Analysis ofE4, E5, and E6

The remaining equilibrium points to investigate the stability areE4, E5, and E6.The stability of equilibrium pointE40,0,0,0,Y1a°a°,Y2a°a°,RAa°a°,RPa°a° occurs when 0.51≤hi≤0.77, i=1,2,3,4. We tested the stability by running the model with different harvesting rates h=0.77,0.77,0.77,0.77,f1,f2, using one hundred different initial values of the form X1,X2,X3,X4,Y1,Y2,RA,RP, taken randomly from the parameter space. In all these simulations, we ended up with the equilibrium solution E4; we therefore conclude that E4 is a stable equilibrium point provided that the harvesting rate of the catfish population is high enough and f1 and f2 are smaller than 0.005.In an analogous manner, we established the stability ofE5 when both f1 and f2 exceed 0.01 and hi≤0.005 for i=1,2,3,4.Finally, when all harvesting rates are low, the simulations always converge to the equilibrium pointE6.


## 6.3. Simulation Results

Our simulations are stable in the sense that when positive initial values are used, the solutions approach the positive steady state (all simulations shown in other figures in this paper have reached steady state) as indicated by Figure3 which shows the dynamics of the solutions to system (18). This signifies coexistence and stability of African catfish, Nile tilapia, and the two different food resources at steady state.Figure 3
Dynamics of African catfish, Nile tilapia, and two different food resources. They–axis shows the biomass in grams per litre (in a logarithmic scale) for each fish compartment together with the food resources. The initial values are chosen positive; all harvest rates are set to 0.005 per day; other parameter values can be found in Tables 3–5.Recently, Lundström et al. [23] adopted the model introduced by De Roos et al. [21], in which Lundström et al. compared the properties between a stage-structured fish population model with an age-structured population model. They investigated how harvesting strategies that qualify for pretty good yield can account for conservation and concluded that equal harvesting rate of juveniles and adults is always a good strategy. We present some of these properties, using equal harvesting in all stages in our model: Figure 4 for African catfish and Nile tilapia and Figure 5 Nile tilapia only.Figure 4
Simulation results, using system (18), including both African catfish and Nile tilapia. The figures show the dependence of impact on size structure (red curve), impact on biomass (blue curve), total biomass (green curve), and yield (black curve) with respect to equal harvesting rates, i.e., h⟶=h,h,h,h,h,h. The initial values are chosen positive, and other parameter values can be found in Tables 3–5.Figure 5
Simulation results, using system (18), including only the Nile tilapia. The figures show the dependence of impact on size structure (red curve), impact on biomass (blue curve), total biomass (green curve), and yield (black curve) with respect to equal harvesting rates, i.e., h⟶=h,h. The initial values are chosen positive, and other parameter values can be found in Tables 3–5.One result by Lundström et al. [23] is that when studying the impact of harvesting on population size structure, they found a positive change in the fraction of juveniles compared to adults in the population. In contrast to this, we find that the impact of population size structure may have a negative dependence on the harvest rate, see Figure 4. This negative impact occurs because African catfish is cannibalistic and predates on Nile tilapia, in that adults feed on juveniles, and this reduces the fraction of juveniles in the system. When we simulate the single fish species Nile tilapia (see Figure 5), our model shows similar results as in [23]; in particular, the impact on size structure is positively dependent on the harvesting rate, which signifies that increasing the uniform harvesting rate there is an increase in the fraction of juveniles.In fisheries, one can use different mesh sizes in the nets; we therefore have modeled system (18) under a specific type of harvesting rate. In Figures 6 and 7, we consider unbalanced harvesting rates on juveniles and adults for both African catfish and Nile tilapia, i.e., harvesting pressure h⟶=1−ph,1−ph,ph,ph,1−ph,ph, where the proportion 0≤p≤1. The figures display the biomass yield (Figure 6) and financial profit (Figure 7) with respect to different harvesting rates, h, and proportions, p. The maximum sustainable biomass yield and maximum sustainable financial profits are attained when p is large; that is, we mainly harvest adults.Figure 6
Simulation results, using system (18), including both African catfish and Nile tilapia with proportional harvesting rates, i.e., h⟶p=1−ph,1−ph,ph,ph,1−ph,ph where 0≤p≤1. The initial values are chosen positive, and other parameter values can be found in Tables 3–5.Figure 7
Simulation results, using system (18), including both African catfish and Nile tilapia with proportional harvesting rates, i.e., h⟶p=1−ph,1−ph,ph,ph,1−ph,ph where 0≤p≤1. The initial values are chosen positive, and other parameter values can be found in Tables 3–5.In fisheries, one can harvest the fish stock by catch and release, for example, with a bottle trap or a fyke net; we therefore model system (18) with variable harvesting rates on all fish compartments. Table 6 shows the results after using pattern search to find optimal harvest rates that maximize biomass yield and financial profit for both species in the system, only Nile tilapia in the system, and only African catfish in the system. Our finding is that harvesting in a system that has both Nile tilapia and African catfish gives the highest biomass yield and financial profit (the optimal values for biomass yield and financial profit in Table 6 are found at different harvesting rates), compared to harvesting in a system having Nile tilapia only or African catfish only.Table 6
Simulation results using optimal harvest rates found in MATLAB, using patternSearch. The columns present the maximum sustainable biomass yield (grams per litre and day) and maximum sustainable financial profit (USD per litre and day). The rows indicate the fish stock scenarios.
Maximum biomass yieldMaximum financial profitOnly Nile tilapia0.00240.0074Only African catfish0.00140.0057Both species0.00270.0081


## 7. Discussion

The major aim of this study was to compare financial profit and biomass outtakes in a two-species system versus single-species systems. To facilitate this comparison, a nonlinear deterministic stage-structured fishery model for the dynamics of African catfish, Nile tilapia, and two food resources has been developed. African catfish is classified into four stages to capture the cannibalistic behavior, while Nile tilapia is categorized into two stages; in addition, the fish feed on two separate food resources in a preferred order. We established the region where the model system (18) is biologically feasible and mathematically well posed. The analytical results revealed that the fish-free two food resource equilibrium point is asymptotically stable if the recovery potentials RC<1 and RT<1. Numerical simulations indicated one positive equilibrium solution for each set of parameters regardless of positive initial values.We established the two-species system to be better in both financial profit and biomass outtakes, compared to a single-species fish stock. Considering equal harvesting rates on both African catfish and Nile tilapia, results show that the impact of harvesting on population size structure is negative. The negative impact signifies a decrease in the fraction of juveniles in the system. This is because African Catfish is cannibalistic and predates on Nile tilapia.Applying optimal harvest rates for harvesting both Nile tilapia and African catfish, African catfish only, and Nile tilapia only to maximize biomass yield and financial profits, our findings to some extent agree with [63, 64], in that, small fisheries using ponds or small lakes should implement both African catfish and Nile tilapia to harvest more biomass yield and gain more financial profits.We emphasize that our findings are based on the model system (18); to prove these results, one should compare two species versus one species in fisheries, alternatively corroborating our findings, using different models. This study can be extended by analyzing a fishery model where there are more than two fish species and more food resources to make the model more realistic. Another direction for future research would be to investigate under what conditions the hydra effect occurs in cannibalistic predator prey models, cf. [65, 66].



---

*Source: 4112015-2022-01-17.xml*