# Bioeconomic-Epidemiological Model of Scomber colias Population in the Moroccan Coasts

Authors: Nossaiba Baba; Imane Agmour; Yousef El Foutayeni; Naceur Achtaich
Journal: Journal of Applied Mathematics
 (2021)

---

## Abstract

In this paper, we develop and study a mathematical model for the dynamics of Scomber colias and Thunnus thynnus prey-predator with parasitic helminths. We search to analyze a bioeconomic model in which both susceptible and infected prey populations Scomber colias are exposed to the predator Thunnus thynnus, with varying degrees of exposure. However, the predator feeds preferentially on the most numerous prey types. This implies a kind of switching from the susceptible class to the infected class, and vice versa, as these two types of prey change in numerical superiority. So, the positivity, boundedness, equilibria, stability, and bioeconomic equilibrium are studied. Some numerical simulation of stability is cited. For giving a high yield and keeping the Scomber colias and Thunnus thynnus populations away from extension, we use the Maximum Principle of Pontryagin.


---

## Body

## 1. Introduction

Morocco is one of the largest fish producers in the world, according to the 2018 report from the Food and Agriculture Organization of the United Nations (FAO). In 2018, national fishery production totaled 1371683 tons for a turnover of 11579544 thousand MAD (Moroccan Dirhams). The export volume reached more than 722921 tons for a turnover of 22,531 million MAD (DPM, 2018) [1]. A study published in December 2019 reveals that the main reason for the repression of Moroccan fish is the presence of parasites. The study was carried out by the Hassan II Agronomic and Veterinary Institute (HAVI) and the center specializing in the pathology of aquatic animals at the National Fisheries Research Institute (NFRI). The parasites are the source of major problems with natural fish stocks. Indeed, the parasites determine pathologies slowing the growth and increasing the mortality of their hosts and constitute, consequently, a limiting factor of the success of the productivity in aquaculture. Following the analysis concerning 1678 pieces of fish of different species, 537 of which come from the Atlantic (port of Essaouira and wholesale market of Casablanca) and 1141 pieces from the Mediterranean, the extent of parasitism is almost the same on the Mediterranean coast (31.1%) and the Atlantic (32%). Nematodes (anisakis + acanthocephalans) occupy most of the parasites in the Atlantic and the Mediterranean, with respective prevalence of 21.4% and 24.9%. The plerocercoid larva of the cestode was found in silver saber, with a prevalence of 8.3% at the Atlantic level, as is shown in Figure 1. It should be noted that to be contaminated, the fish must feed on plants or animals carrying parasites. In fact, several animal species can carry parasites and contaminate the fish’s environment through their excrement, which sometimes contains parasite eggs. To develop, these eggs will infect plankton, small crustaceans, snails, etc. If fish eat these organisms, they in turn can carry parasites. These reasons encouraged us to study this phenomenon and to see all the interactions that occur between populations of fish. In this context, many mathematical models have been developed to describe the dynamics of fisheries, and we can refer, for example, [2–7]. Moreover, we can cite [8], in this work, the authors have defined a bioeconomic equilibrium model for Parapeaneus longirostris and small pelagic fish populations in two different areas; the first one is protected against fishing and the second is a free access zone.Figure 1
Infected fish among fish examined in the Atlantic and Mediterranean Sea.They studied the influence of the predator mortality rate variation on the evolution of prey biomass and the profit of coastal trawlers. In [9], the authors have sought to highlight that the increase of the carrying capacity of marine species does not always lead to an increase on the catch levels and on the incomes. They considered a bioeconomic model of Sardina pilchardus, Engraulis encrasicolus, and Xiphias gladius marine species that are exploited by several seiners in the Moroccan Atlantic zone based on the parameters given by NFRI. In [10], the authors studied a model consisting of susceptible and infected prey populations and a predator population. More precisely, they assume that the likelihood of a predator catching a susceptible prey or an infected prey is proportional to the numbers of these two different types of prey species. They assume that the predator will eventually die as a result of eating infected prey. In previous works, in particular, in the works [11–23], the authors did not take into account parasitic helminths which have a negative impact on the evolution and behavior of populations. Therefore, the purpose of this paper is to identify parasitic Helminths in Scomber colias and Thunnus thynnus of great economic importance, which occupy a prominent place in the trophodynamic chain. We search to study a bioeconomic model in which both susceptible and infected prey populations (Scomber colias) are exposed to the predator (Thunnus thynnus), with varying degrees of exposure. However, the predator feeds preferentially on the most numerous prey types. This implies a kind of switching from the susceptible class to the infected class, and vice versa, as these two types of prey change in numerical superiority. Switching may simply come about due to individual behavior changing with varying abundance of that class of prey. The paper is organized as follows. In the next section, we present the biological model of susceptible and infected Scomber colias with the presence of the predators Thunnus thynnus; in other words, we resolve a system of three ordinary differential equations, the first equation describes the natural growth of the susceptible Scomber colias fish population a prey of the Thunnus thynnus fish population, the second equation describes the natural growth of infected Scomber colias fish population a prey of the Thunnus thynnus fish population, and the third equation describes the natural growth of the Thunnus thynnus fish population as a predator of the susceptible and infected Scomber colias. The basic reproduction number, positivity, and boundedness are studied in the first part. The existence of the steady states of this system and their stability are studied using eigenvalue analysis, and we define a bioeconomic equilibrium model for these fish populations exploited by a fishing fleet. In Section 3, we compute some numerical simulations to determine the optimal conditions under which the biological steady state can be attained. In Section 4, we give a numerical simulation of the mathematical model and discussion of the results. Finally, we give a conclusion and some potential perspectives in Section 5.


## 2. Biological Model

Our study is based on an epidemiological model which describes the interaction between the susceptible and infected fish population Scomber colias and their predators Thunnus thynnus. We assume that the disease is transmitted by direct contact with the prey. Due to the presence of the disease, the prey population is divided into two disjoint classes, the susceptible fish population denoted byMS per unit designated area and the infected fish population denoted by MI per unit designated area. Therefore, at time t, the total Scomber colias population is Mt=MSt+MIt. The predator population per unit designated area is denoted by T. If there is no predation and no infection, then the fish population follows a logistic growth model. Both susceptible and infected Scomber colias are preyed upon by Thunnus thynnus, but the Thunnus thynnus preferentially eats the infected fish, because the infection causes the fish to be more vulnerable to predation (as shown in Figure 2).Figure 2
Interaction between susceptible Scomber colias, infected Scomber colias, and Thunnus thynnus.LetMS be the density of susceptible Scomber colias at time t. This population grows according to a logistic equation with growth rate r and carrying capacity L in the absence of the disease. The parameter α represents the rate of transmission multiplied by the area under consideration, measured in units of area. The susceptible Scomber colias is a prey of the Thunnus thynnus with the response rate β as follows
(1)ṀSt=rMSt1−MSt+MItL−αMStMIt−βMS2tTtMSt+σMIt.LetMI be the density of infected Scomber colias with the rate of transmission of the disease α. This population is preyed by the Thunnus thynnus with the response rate a. The infected Scomber colias populations died by the death rate δ as follows
(2)ṀIt=αMStMIt−αMI2tTtMSt+σMIt−δMIt.The Thunnus thynnusT is a predator of susceptible Scomber colias and infected Scomber colias where η1 represents the rate of conversion to Thunnus thynnus of susceptible mackerel, η2 is the rate of conversion to Thunnus thynnus of infected Scomber colias, and d is the death rate of the Thunnus thynnus as follows
(3)Ṫt=η1βMS2tTtMSt+σMIt−η2aMI2tTtMSt+σMIt−dTt.Following the previous assumptions, the biological model is formulated as follows(4)ṀSt=rMSt1−MSt+MItL−αMStMIt−βMS2tTtMSt+σMIt,ṀIt=αMStMIt−aMI2tTtMSt+σMIt−δMIt,Ṫt=η1βMS2tTtMSt+σMIt−η2aMI2tTtMSt+σMIt−dTt.We assume that all parameters are positive,MSt≥0,MIt≥0,Tt≥0. The system of equations (4) has a unique solution in
(5)V=MS,MI,T/MS≥0,MI≥0andT≥0,since the right-hand sides of them are Lipschitz continuous in V.AtMS=MI=0, the term βMS2T/MS+σMI is not defined. As for MS and MI are small, 0<MI/MS+σMI≤1 and 0<MS/MS+σMI≤1, so it is natural to interpret as zero at MS=MI=0. For aMI2T/MS+σMI, we have the same interpretation at MS=MI=0.


## 3. Biological Model Analysis

### 3.1. Basic Reproduction Number

We assume that we have the following SI model(6)ṀS=−αMSMI,ṀI=αMSMI−δMI.By dividing the second equation byδMI, we obtain ṀI/δMI=α/δMS−1. Then, for α/δMS>1, each infected individual will contaminate more than one susceptible individual, and the disease will spread to an ever-increasing number of individuals. This will remain the case until the number of susceptible MS is such that α/δMS<1. The ratio α/δMS can then be interpreted as the number of contacts that can transmit the disease by infected individuals throughout their period of contagion. If α/δMS∗>1, there will necessarily be an epidemic, whereas in the opposite case, if α/δMS∗<1, only a few individuals will be infected before the spread of the disease does not stop by itself. Let R0=α/δMS0, where MS0=L. This threshold is the basic reproduction rate.


### 3.2. Positivity

We will denote by(7)f1=r1−MS+MIL−αMI−βMSTMS+σMI,f2=αMS−aMITMS+σMI−δ,f3=η1βMS2MS+σMI−η2aMI2MS+σMI−d.Then, the system (4) becomes
(8)ṀSt=MStf1MSt,MIt,Tt,ṀIt=MItf2MSt,MIt,Tt,Ṫt=Ttf3MSt,MIt,Tt,subject to the initial condition MS (0) > 0,MI (0) > 0 and T (0) > 0. The system (8) is defined on the set
(9)Ω=MS,MI,T∈ℝ3/MS≥0,MI≥0,T≥0.LetXt=MSt,MIt,Tt be the solution of the system (4) at the biological equilibrium. Then, all the solutions of the system (4) are nonnegative. By system of equations (8) with initial condition, we have
(10)MSt=MS0exp∫0tf1MSs,MIs,Tsds>0,MIt=MI0exp∫0tf2MSs,MIs,Tsds>0,Tt=T0exp∫0tf3MSs,MIs,Tsds>0.Therefore, all solutions starting from an interior of the first octant remain in it for all future time.


### 3.3. Boundedness

We define the functionNt=MSt+MIt+Tt, and we take the time derivative along the solution of the system Ṅt=ṀSt+ṀIt+Ṫt. For any positive constant ρ, we have
(11)Ṅt+ρNt=rMSt−rMStMSt+MItL−βMS2tTt+aMI2tTtMSt+σMIt−δMIt+η1βMS2tTt−η2aMI2tTtMSt+σMIt−dTt+ρMSt+MIt+Tt=r+ρMSt+ρ−δMIt+ρ−dTt−1−η1βMS2tTtMS+σMI−a+η2MI2tTtMS+σMI−rMStMSt+MItL≤r+ρMSt+ρ−δMIt+ρ−dTt.We chooseρ as a minimum of δ,d. On the one hand, we have
(12)ṀSt+ṀIt≤rα1−MSt+MItL.Consequently,limsupt→∞MSt+MIt≤L. Then, limsupt→∞MSt≤L and limsupt→∞MIt≤L. That improves that the susceptible and infected preys are always bounded. On the other hand, there exists t0 such that t≥t0,MS≤2L, thus, for t≥t0(13)Ṅt+ρNt≤2Lrα+ρ=ω.Which gives, for(14)Nt≤ωρ+Nt0e−ρt−t0.Hence, thatlimsupt→∞Nt≤ω/ρ independently of the initial conditions. Finally, we conclude that the trajectories of the system are bounded.


### 3.4. Equilibrium Points

For the system (4), we can see that the equilibrium points satisfy the following equations
(15)r1−MS+MIL−αMI−βMSTMS+σMI=0,αMS−aMITMS+σMI−δ=0,η1βMS2MS+σMI−η2aMI2MS+σMI−d=0.After calculation, it is obvious that system (15) has five equilibrium points as follows
(16)P1=0,0,0,P2=dη1β,0,r−d+Lβη1Lβ2η1,P3=δα,−LrLαr+lαδ−Lα,0,P4=L,0,0,P5MS,MI,T.For the last equilibrium pointP5MS,MI,T is given by
(17)MS=−aMIrMI−Lr+LαMIarMI+Lαβ−Lβδ,T=aLα−rσ−1MI−Lar+σβα−δaαMI2r+Lα+Lβδα−δ+arMIδ−LαaarMI+Lβα−δ2,where MI is given in the following calculation.SinceMS=aMI−rMI+Lr−LαMI/arMI+Lαβ−Lβδ≥0, so we obtain MI≤Lr/r+Lα. Hence, there is no economic equilibrium when MI>Lr/r+Lα.Consider the functiong defined on an interval 0,Lr/r+Lα by
(18)gMI=MI1MI3+MI2MI2+MI3MI+MI4=0,where
(19)MI1=a2r2aη2−βη1−Lαβη12r+Lα,MI2=a2rdrσ−1−Ldα+2Lβη1r+Lα+2Lβη2α−δ,MI3=Laar2d−Lβη1+α−δLβ2η2α−δ−dβLα−r2σ−1,MI4=L2dβα−δar+σβα−δ.IfMI4=0, i.e., α=δ, then, equation (18) admits two different solutions: the first one is trivial and the other solution exists if MI3<0 and MI1>0. When α<δ, we have g0=L2dβα−δar+σβα−δ<0 with ar+σβα−δ>0 and
(20)gLrr+Lα=MI1Lrr+Lα3+MI2Lrr+Lα2+MI3Lrr+Lα+MI4=L2ar2+βα−δr+Lα2dσr+Lα+Larη2>0.Equation (18) admits at least one solution MI∗∈0,Lr/r+Lα which implies that system (1) has at least one capital-labor equilibrium MS∗,MI∗,T∗ with
(21)MS∗=aMI∗Lr−r+LαMI∗arMI∗+Lβα−δ,T∗=ar1−σ+LαMI∗−Lar+σβα−δaαr+LαMI∗2+arδ−LαMI∗+Lβδα−δaarMI∗+Lαβ−Lβδ2.To determine the uniqueness of the solutionMI∗, we discuss the following cases
(i)
IfMI1=MI2=0, then the equation (18) has a unique positive solution given by(22)MI∗=L2dβδ−αar+σβα−δLaar2d−Lβη1+α−δβdr2σ−1−Ldα−βη2α−δ,d−Lβη1<0,α−δ>0anddr2σ−1<Ldα−βη2α−δthenMI∗>0.(ii)
IfMI1=0 and MI2≠0, then equation (18) has a unique positive solution given by MI∗=1/2MI2MI32−4MI2MI4−MI3(iii)
IfMI1≠0 by applying the rule signs of Descartes, equation (18) has a unique positive root MI∗, if any of the following three conditions is satisfied(23)MI2>0andMI3>0,MI2>0andMI3<0,MI2<0andMI3<0.To depress the cubic equation, we substituteMI=x−MI2/3MI1 and we solve gMI=0. Then, we get
(24)MI1x−MI23MI13+MI2x−MI23MI12+MI3x−MI23MI1+MI4=0,x−MI23MI13+MI2MI1x−MI23MI12+MI3MI1x−MI23MI1+MI4MI1=0,x3−MI2327MI13−MI2MI1x2+MI223MI12x+MI2MI1MI229MI12+x2−2MI23MI1x+MI3MI1x−MI23MI1+MI4MI1=0,x3+3MI1MI3−MI223MI12x+27MI12MI4+2MI23−9MI1MI2MI327MI13=0.Then, we obtaingx=x3+px+q where
(25)p=3MI1MI3−MI223MI12,q=27MI12MI4+2MI23−9MI1MI2MI327MI13.By using Cardano’s formula, we obtainΔ1=q2+4p3/27 then
(26)MI∗=−p−q2+4p3/2723+−p+q2+4p3/2723−MI23MI1.


### 3.5. Local Stability of Equilibrium Points

To ensure the positivity of the population’s biomass, we assure that(27)Lα>δ,Lβη1>d,Lα>r,We consider the matrixJS that represents the variational matrix of the system at the equilibrium point P=MS,MI,T(28)JP=J11J12J13J21J22J23J31J32J33,where
(29)J11=r1−MS+MIL−αMI−2βMSTMS+σMI−βMS2TMS+σMI2,J12=1LrMS−αMS+σβMS2TMS+σMI2,J13=−βMS2MS+σMI,J21=αMI−−aMI2TMS+σMI2,J22=αMS−2aMITMS+σMI−aσMI2TMS+σMI2−δ,J23=−aMI2MS+σMI,J31=2η1βMSTMS+σMI−η1βMS2T+η2aMI2TMS+σMI2,J32=−η1σβMS2T−2η2aMITMS+σMI−η2σaMI2TMS+σMI2,J33=η1βMS2MS+σMI−η2aMI2MS+σMI−d.(1)
At the equilibrium pointP1=0,0,0, the variational matrix is(30)JP1=r000−δ000−d.It has three eigenvalues:λ1=r>0,λ2=−δ<0, and λ3=−d<0. Consequently, the point P1 is unstable.
(2)
AtP2=d/η1β,0,r−d+Lβη1/Lβ2η1 the variational matrix is(31)JP2=0dr−Ldα−drσ+Lrσβη1Lβη1−dη10dα−βδη1βη10rLβη1−dLβrσd−Lβη1Lβ0.It has these eigenvalues:(32)λ1=−β2η1drLβ5η13d−Lβη1<0,λ2=β2η1drLβ5η13d−Lβη1>0,λ3=dα−βδη1βη1>0.Hence, the pointP2 is unstable.
(3)
In the same manner, atP3=δ/α,−Lrδ−Lα/Lαr+lα,0, the variational matrix is(33)JP3=0δr−LαLα−βδ2Lr+LααLrδ+LLαδ−Lrσδ+LLrασr−δ+Lαr+Lα0−aL2r2−δ+Lα2Lαr+LαLrδ+LLαδ−Lrσδ+LLrασ00η1βδ/α2δ/α+σ−Lrδ−Lα/Lαr+Lα−η2a−lrδ−Lα/Lαr+lα2δ/α+σ−Lrδ−Lα/Lαr+Lα−d.It has these eigenvalues:(34)λ1=L7α7r+Lα3rδLα−rδ−LαL15α15r+Lα7,λ2=−L7α7r+Lα3rδLα−rδ−LαL15α15r+Lα7,λ3=λ31+λ32−λ33L2α3δ+r2αδ−r2ασδ+Lr2α2σ+L2rα3σ+2Lrα2δ−Lrα2σδ,where
(35)λ31=L2α2−βδ2η1+dαδ+ar2η2+drασ,λ32=r2δdα−dασ+aδη2−βδη1,λ33=Lrα2βδ2η1−2dαδ+2arδη2−drασ+dασδ,Then, following the conditionsλ1<0,λ2>0 and λ3>0, we can deduce that the point P3 is unstable.
(4)
The variational matrix atP4=L,0,0 is given by(36)JP4=0r−Lα−Lβ0αL−δ000Lβη1−d.Since one of the eigenvaluesλ1=Lα−δ>0 is positive Then, the point P4 is unstable.
(5)
The local stability ofP5MS,MI,T is investigated by considering the variational matrix based on the system (1) with(37)JP5=J11J12−βMS2MS+σMIJ21−TaMSMIMS+σMI2−aMI2MS+σMIJ31J32−d,where
(38)J11=−MSrMS+σMI2+LTσβMILMS+σMI2,J12=−MSr+LαL−σβMS2TMS+σMI2,J21=αMI+aMI2TMS+σMI,J31=Tβη1MS2−aη2MI2σMI−1−MSMIaη2MI−2σβη1MS+σMI2,J32=−Tσβη1MS2+2aη2MSMI+aση2MI2MS+σMI2.(i)
IfMI1=MI2=0, then(39)MI=L2dβδ−αar+σβα−δLaar2d−Lβη1+α−δβdr2σ−1−Ldα−βη2α−δ,T=T11T12+T13aα−δdβLα−r2σ−1−Lβ2η2α−δ+ar2Lβη1−d/aα−δdβr−2rσ+Lα−Lβ2η2α−δ+ar2Lβη1−d2aLβ2α−δα−δdLα−rσ−1−Lβη2α−δ+Lar2η1/α−δdβ−r2σ−1+Lα−Lβ2η2α−δ+ar2Lβη1−d2,MS=Ldar+σβα−δra2r2Lβη1−d−2dσβarα−δ−αβδ−MS1aβα−δLβdα−βη2α−δ−drβ2σ−1MS2,where
(40)T11=Lar+σβα−δα−δβdrσ+Lβη2α−δ+ar2d−Lβη1α−δdβLα−r2σ−1−Lβ2η2α−δ+ar2Lβη1−d,T12=L2d2αβ2ar+σβα−δ2r+Lαα−δ2aα−δ−Lβ2η2α−δ+dβr−2rσ+Lα+ar2Lβη1−d2,T13=Lβα−δ2Lβdαrσ−δ+βδη2+drβδσ−1+Lar2dα−βδη1,MS1=β2arη2+dασLα2+β2δ2α−δ+drσα2+δ2,MS2=α−δLdα−drσ−1−Lβη2α−δ+Lar2η1.In this case, the characteristic polynomialρ3λ3+ρ2λ2+ρ1λ+ρ0 where
(41)ρ3=1,ρ2=TMSMIa+σβMS+σMI2,ρ1=−ρ11+ρ12+ρ13+ρ14+ρ15LMS+σMI3,ρ0=TMSMIρ21+ρ22−ρ23−ρ24−ρ25+ρ26LMS+σMI4,where
(42)ρ11=MSMIαMS+σMI3+TaMSMIr−Lα,ρ12=LTβMS2MI2aση1−aη2+ασ2,ρ13=LT3a2ση2MI4−β2η1MS4,ρ14=TaMSMI3rσ+2Laη2−Lασ,ρ15=LTσβMS3MIα−2βη1,ρ21=aMIMS+σMIβη1MS2+aη2MI2+2σβη1MSMI×r−Lα,ρ22=LaσβMS2MITβη1−5αη2MI,ρ23=3Laσβη2MI3Ta+ασMS,ρ24=LαβMS3MI2aη2+σ2βη1,ρ25=LTa2βη2MSMI2,ρ26=Lβ2η1MS3Ta−ασMS.By considering the previous characteristic polynomial, if one of the coefficients is zero or negative while at least one other coefficient is positive, then there is one or more imaginary roots or roots with a positive real part. If all the coefficients are positive, we calculate the Routh table. The Routh Hurwitz tabular is given byλ3λ2λ11ρ1ρ2ρ0ρ110ρ210 with ρ11=ρ1ρ2−ρ0/ρ23>0,ρ21=ρ0/ρ112>0. So, according to the Routh Hurwitz stability criterion, we can conclude that the equilibrium point in this case is stable.
(ii)
IfMI1=0 and MI2≠0, then(43)MI=12MI2MI32−4MI2MI4−MI3,MS=aMI32−4MI2MI4−MI3Lr−α−rMI32−4MI2MI4−MI32MI2arMI32−4MI2MI42MI2Lβα−δ−arMI3,T=MI3+arσ−1−LaαMI32−4MI2MI4+2LMI2ar+σβα−δaarMI32−4MI2MI4+2LMI2βα−δ−MI32×2MI3+MI32−4MI2MI4Laαr−α−arα+δ−2LMI2βδα−δ.In this case, the characteristic polynomial is as follows:ρ3λ3+ρ2λ2+ρ1λ+ρ0, where
(44)ρ3=1,ρ2=TMSMIa+σβMS+σMI2,ρ1=−ρ11+ρ12+ρ13+ρ14+ρ15LMS+σMI3,ρ0=TMSMIρ21+ρ22−ρ23−ρ24−ρ25+ρ26LMS+σMI4,where
(45)ρ11=MSMIαMS+σMI3+TaMSMIr−Lα,ρ12=LTβMS2MI2aση1−aη2+ασ2,ρ13=LT3a2ση2MI4−β2η1MS4,ρ14=TaMSMI3rσ+2Laη2−Lασ,ρ15=LTσβMS3MIα−2βη1,ρ21=aMIMS+σMIβη1MS2+aη2MI2+2σβη1MSMI×r−Lα,ρ22=LaσβMS2MITβη1−5αη2MI,ρ23=3Laσβη2MI3Ta+ασMS,ρ24=LαβMS3MI2aη2+σ2βη1,ρ25=LTa2βη2MSMI2,ρ26=Lβ2η1MS3Ta−ασMS.The Routh Hurwitz tabular is given byλ3λ2λ11ρ1ρ2ρ0ρ110ρ210 with ρ11=ρ1ρ2−ρ0/ρ23>0,ρ21=ρ0/ρ112>0. Then, according to the Routh Hurwitz stability criterion, we can conclude that the equilibrium point in this case is stable.
(iii)
IfMI1≠0,, then(46)MI=3MI1−p−4/27p3+q2/23+4/27p3+q2−p/23−MI23MI1,MS=amS1−mS2mS3−MI23MI1mS4+mS5,T=t1−t2−t3t4−MI22+3MI1t5−t6−t7Lα−δ3MI1a3LβMI1α−δ−arMI2+t72,where
(47)3LβδMI1α−δ−arMI2δ−LαLα−δ>3arMI1−p−4p3/27+q223+4p3/27+q2−p23,mS1=MI2r+Lα+3LrMI1,mS2=3MI1r+Lα4/27p3+q2−p23+−p−4/27p3+q223,mS3=3MI1−p−4/27p3+q223+4/27p3+q2−p23,mS4=3LβMI1α−δ−arMI2,mS5=3arMI1−p−4p3/27+q223+4/27p3+q2−p23,and
(48)t1=3aMI1−p−4p3/27+q223+4p3/27+q2−p23×Lα−rσ−1,t2=aMI2Lα−rσ−1,t3=3LMI1ar+σβα−δ,t4=aαr+Lα3MI1−p−4p3/27+q223+4p3/27+q2−p23,t5=3LβδMI1α−δ,t6=arMI2δ−Lα,t7=3arMI1−p−4p3/27+q223+4p3/27+q2−p23.In this case, the characteristic polynomial is as follows:ρ3λ3+ρ2λ2+ρ1λ+ρ0 where
(49)ρ3=1,ρ2=TMSMIa+σβMS+σMI2,ρ1=−ρ11+ρ12+ρ13+ρ14+ρ15LMS+σMI3,ρ0=TMSMIρ21+ρ22−ρ23−ρ24−ρ25+ρ26LMS+σMI4,where ρ11+ρ13+ρ15>ρ12+ρ14ρ11<0,ρ15<0,ρ13<0, and
(50)ρ11=MSMIαMS+σMI3+TaMSMIr−Lα,ρ12=LTβMS2MI2aση1−aη2+ασ2,ρ13=LT3a2ση2MI4−β2η1MS4,ρ14=TaMSMI3rσ+2Laη2−Lασ,ρ15=LTσβMS3MIα−2βη1,ρ21=aMIMS+σMIβη1MS2+aη2MI2+2σβη1MSMI×r−Lα,ρ22=LaσβMS2MITβη1−5αη2MI,ρ23=3Laσβη2MI3Ta+ασMS,ρ24=LαβMS3MI2aη2+σ2βη1,ρ25=LTa2βη2MSMI2,ρ26=Lβ2η1MS3Ta−ασMS.The Routh Hurwitz tabular is given byλ3λ2λ11ρ1ρ2ρ0ρ110ρ210 with ρ11=ρ1ρ2−ρ0/ρ23,ρ21=ρ0/ρ112. If p>q with q,p>0 and if t5−t6>t7Lα−δ then ρ2>0. If ρ11+ρ13+ρ15>ρ12+ρ14 then ρ1>0. On the other hand, we have ρ1ρ2>ρ0. We conclude that ρ11>0 and ρ21>0. So, according to the Routh Hurwitz stability criterion, we can conclude that the equilibrium point in this case is stable. Beginning with parameters r=0.9,L=12,α=0.9,β=0.2,σ=0.2,a=0.3,δ=0.2,η1=0.4,η2=0.3, and d=0.2 and initial value 0.1,0.1,0.1, the point P5 tends to the point 0.221;0.784;1.65×10−33. As shown in Figures 3–5. Then, we can deduce that the point P5 is stable. The previous Figures 3–6 show the dynamical behavior of the susceptible populations of Scomber colias and Thunnus thynnus, also the infected population of Scomber colias. It is clear that the biomass of these marine populations converges to the values of the interior equilibrium point. This means biologically that the interior equilibrium point can ensure the existence of the predefined species.Figure 3
Dynamical behavior of the susceptible Scomber colias.Figure 4
Dynamical behavior of the infected Scomber colias.Figure 5
Dynamical behavior of the susceptible Thunnus thynnus.Figure 6
Illustration of the evolution of the susceptible, infected Scomber colias and Thunnus thynnus fish population.Following Figure6, in a finite time the susceptible Scomber colias, infected Scomber colias and Thunnus thynnus converge to the equilibrium point P5.


## 3.1. Basic Reproduction Number

We assume that we have the following SI model(6)ṀS=−αMSMI,ṀI=αMSMI−δMI.By dividing the second equation byδMI, we obtain ṀI/δMI=α/δMS−1. Then, for α/δMS>1, each infected individual will contaminate more than one susceptible individual, and the disease will spread to an ever-increasing number of individuals. This will remain the case until the number of susceptible MS is such that α/δMS<1. The ratio α/δMS can then be interpreted as the number of contacts that can transmit the disease by infected individuals throughout their period of contagion. If α/δMS∗>1, there will necessarily be an epidemic, whereas in the opposite case, if α/δMS∗<1, only a few individuals will be infected before the spread of the disease does not stop by itself. Let R0=α/δMS0, where MS0=L. This threshold is the basic reproduction rate.


## 3.2. Positivity

We will denote by(7)f1=r1−MS+MIL−αMI−βMSTMS+σMI,f2=αMS−aMITMS+σMI−δ,f3=η1βMS2MS+σMI−η2aMI2MS+σMI−d.Then, the system (4) becomes
(8)ṀSt=MStf1MSt,MIt,Tt,ṀIt=MItf2MSt,MIt,Tt,Ṫt=Ttf3MSt,MIt,Tt,subject to the initial condition MS (0) > 0,MI (0) > 0 and T (0) > 0. The system (8) is defined on the set
(9)Ω=MS,MI,T∈ℝ3/MS≥0,MI≥0,T≥0.LetXt=MSt,MIt,Tt be the solution of the system (4) at the biological equilibrium. Then, all the solutions of the system (4) are nonnegative. By system of equations (8) with initial condition, we have
(10)MSt=MS0exp∫0tf1MSs,MIs,Tsds>0,MIt=MI0exp∫0tf2MSs,MIs,Tsds>0,Tt=T0exp∫0tf3MSs,MIs,Tsds>0.Therefore, all solutions starting from an interior of the first octant remain in it for all future time.


## 3.3. Boundedness

We define the functionNt=MSt+MIt+Tt, and we take the time derivative along the solution of the system Ṅt=ṀSt+ṀIt+Ṫt. For any positive constant ρ, we have
(11)Ṅt+ρNt=rMSt−rMStMSt+MItL−βMS2tTt+aMI2tTtMSt+σMIt−δMIt+η1βMS2tTt−η2aMI2tTtMSt+σMIt−dTt+ρMSt+MIt+Tt=r+ρMSt+ρ−δMIt+ρ−dTt−1−η1βMS2tTtMS+σMI−a+η2MI2tTtMS+σMI−rMStMSt+MItL≤r+ρMSt+ρ−δMIt+ρ−dTt.We chooseρ as a minimum of δ,d. On the one hand, we have
(12)ṀSt+ṀIt≤rα1−MSt+MItL.Consequently,limsupt→∞MSt+MIt≤L. Then, limsupt→∞MSt≤L and limsupt→∞MIt≤L. That improves that the susceptible and infected preys are always bounded. On the other hand, there exists t0 such that t≥t0,MS≤2L, thus, for t≥t0(13)Ṅt+ρNt≤2Lrα+ρ=ω.Which gives, for(14)Nt≤ωρ+Nt0e−ρt−t0.Hence, thatlimsupt→∞Nt≤ω/ρ independently of the initial conditions. Finally, we conclude that the trajectories of the system are bounded.


## 3.4. Equilibrium Points

For the system (4), we can see that the equilibrium points satisfy the following equations
(15)r1−MS+MIL−αMI−βMSTMS+σMI=0,αMS−aMITMS+σMI−δ=0,η1βMS2MS+σMI−η2aMI2MS+σMI−d=0.After calculation, it is obvious that system (15) has five equilibrium points as follows
(16)P1=0,0,0,P2=dη1β,0,r−d+Lβη1Lβ2η1,P3=δα,−LrLαr+lαδ−Lα,0,P4=L,0,0,P5MS,MI,T.For the last equilibrium pointP5MS,MI,T is given by
(17)MS=−aMIrMI−Lr+LαMIarMI+Lαβ−Lβδ,T=aLα−rσ−1MI−Lar+σβα−δaαMI2r+Lα+Lβδα−δ+arMIδ−LαaarMI+Lβα−δ2,where MI is given in the following calculation.SinceMS=aMI−rMI+Lr−LαMI/arMI+Lαβ−Lβδ≥0, so we obtain MI≤Lr/r+Lα. Hence, there is no economic equilibrium when MI>Lr/r+Lα.Consider the functiong defined on an interval 0,Lr/r+Lα by
(18)gMI=MI1MI3+MI2MI2+MI3MI+MI4=0,where
(19)MI1=a2r2aη2−βη1−Lαβη12r+Lα,MI2=a2rdrσ−1−Ldα+2Lβη1r+Lα+2Lβη2α−δ,MI3=Laar2d−Lβη1+α−δLβ2η2α−δ−dβLα−r2σ−1,MI4=L2dβα−δar+σβα−δ.IfMI4=0, i.e., α=δ, then, equation (18) admits two different solutions: the first one is trivial and the other solution exists if MI3<0 and MI1>0. When α<δ, we have g0=L2dβα−δar+σβα−δ<0 with ar+σβα−δ>0 and
(20)gLrr+Lα=MI1Lrr+Lα3+MI2Lrr+Lα2+MI3Lrr+Lα+MI4=L2ar2+βα−δr+Lα2dσr+Lα+Larη2>0.Equation (18) admits at least one solution MI∗∈0,Lr/r+Lα which implies that system (1) has at least one capital-labor equilibrium MS∗,MI∗,T∗ with
(21)MS∗=aMI∗Lr−r+LαMI∗arMI∗+Lβα−δ,T∗=ar1−σ+LαMI∗−Lar+σβα−δaαr+LαMI∗2+arδ−LαMI∗+Lβδα−δaarMI∗+Lαβ−Lβδ2.To determine the uniqueness of the solutionMI∗, we discuss the following cases
(i)
IfMI1=MI2=0, then the equation (18) has a unique positive solution given by(22)MI∗=L2dβδ−αar+σβα−δLaar2d−Lβη1+α−δβdr2σ−1−Ldα−βη2α−δ,d−Lβη1<0,α−δ>0anddr2σ−1<Ldα−βη2α−δthenMI∗>0.(ii)
IfMI1=0 and MI2≠0, then equation (18) has a unique positive solution given by MI∗=1/2MI2MI32−4MI2MI4−MI3(iii)
IfMI1≠0 by applying the rule signs of Descartes, equation (18) has a unique positive root MI∗, if any of the following three conditions is satisfied(23)MI2>0andMI3>0,MI2>0andMI3<0,MI2<0andMI3<0.To depress the cubic equation, we substituteMI=x−MI2/3MI1 and we solve gMI=0. Then, we get
(24)MI1x−MI23MI13+MI2x−MI23MI12+MI3x−MI23MI1+MI4=0,x−MI23MI13+MI2MI1x−MI23MI12+MI3MI1x−MI23MI1+MI4MI1=0,x3−MI2327MI13−MI2MI1x2+MI223MI12x+MI2MI1MI229MI12+x2−2MI23MI1x+MI3MI1x−MI23MI1+MI4MI1=0,x3+3MI1MI3−MI223MI12x+27MI12MI4+2MI23−9MI1MI2MI327MI13=0.Then, we obtaingx=x3+px+q where
(25)p=3MI1MI3−MI223MI12,q=27MI12MI4+2MI23−9MI1MI2MI327MI13.By using Cardano’s formula, we obtainΔ1=q2+4p3/27 then
(26)MI∗=−p−q2+4p3/2723+−p+q2+4p3/2723−MI23MI1.


## 3.5. Local Stability of Equilibrium Points

To ensure the positivity of the population’s biomass, we assure that(27)Lα>δ,Lβη1>d,Lα>r,We consider the matrixJS that represents the variational matrix of the system at the equilibrium point P=MS,MI,T(28)JP=J11J12J13J21J22J23J31J32J33,where
(29)J11=r1−MS+MIL−αMI−2βMSTMS+σMI−βMS2TMS+σMI2,J12=1LrMS−αMS+σβMS2TMS+σMI2,J13=−βMS2MS+σMI,J21=αMI−−aMI2TMS+σMI2,J22=αMS−2aMITMS+σMI−aσMI2TMS+σMI2−δ,J23=−aMI2MS+σMI,J31=2η1βMSTMS+σMI−η1βMS2T+η2aMI2TMS+σMI2,J32=−η1σβMS2T−2η2aMITMS+σMI−η2σaMI2TMS+σMI2,J33=η1βMS2MS+σMI−η2aMI2MS+σMI−d.(1)
At the equilibrium pointP1=0,0,0, the variational matrix is(30)JP1=r000−δ000−d.It has three eigenvalues:λ1=r>0,λ2=−δ<0, and λ3=−d<0. Consequently, the point P1 is unstable.
(2)
AtP2=d/η1β,0,r−d+Lβη1/Lβ2η1 the variational matrix is(31)JP2=0dr−Ldα−drσ+Lrσβη1Lβη1−dη10dα−βδη1βη10rLβη1−dLβrσd−Lβη1Lβ0.It has these eigenvalues:(32)λ1=−β2η1drLβ5η13d−Lβη1<0,λ2=β2η1drLβ5η13d−Lβη1>0,λ3=dα−βδη1βη1>0.Hence, the pointP2 is unstable.
(3)
In the same manner, atP3=δ/α,−Lrδ−Lα/Lαr+lα,0, the variational matrix is(33)JP3=0δr−LαLα−βδ2Lr+LααLrδ+LLαδ−Lrσδ+LLrασr−δ+Lαr+Lα0−aL2r2−δ+Lα2Lαr+LαLrδ+LLαδ−Lrσδ+LLrασ00η1βδ/α2δ/α+σ−Lrδ−Lα/Lαr+Lα−η2a−lrδ−Lα/Lαr+lα2δ/α+σ−Lrδ−Lα/Lαr+Lα−d.It has these eigenvalues:(34)λ1=L7α7r+Lα3rδLα−rδ−LαL15α15r+Lα7,λ2=−L7α7r+Lα3rδLα−rδ−LαL15α15r+Lα7,λ3=λ31+λ32−λ33L2α3δ+r2αδ−r2ασδ+Lr2α2σ+L2rα3σ+2Lrα2δ−Lrα2σδ,where
(35)λ31=L2α2−βδ2η1+dαδ+ar2η2+drασ,λ32=r2δdα−dασ+aδη2−βδη1,λ33=Lrα2βδ2η1−2dαδ+2arδη2−drασ+dασδ,Then, following the conditionsλ1<0,λ2>0 and λ3>0, we can deduce that the point P3 is unstable.
(4)
The variational matrix atP4=L,0,0 is given by(36)JP4=0r−Lα−Lβ0αL−δ000Lβη1−d.Since one of the eigenvaluesλ1=Lα−δ>0 is positive Then, the point P4 is unstable.
(5)
The local stability ofP5MS,MI,T is investigated by considering the variational matrix based on the system (1) with(37)JP5=J11J12−βMS2MS+σMIJ21−TaMSMIMS+σMI2−aMI2MS+σMIJ31J32−d,where
(38)J11=−MSrMS+σMI2+LTσβMILMS+σMI2,J12=−MSr+LαL−σβMS2TMS+σMI2,J21=αMI+aMI2TMS+σMI,J31=Tβη1MS2−aη2MI2σMI−1−MSMIaη2MI−2σβη1MS+σMI2,J32=−Tσβη1MS2+2aη2MSMI+aση2MI2MS+σMI2.(i)
IfMI1=MI2=0, then(39)MI=L2dβδ−αar+σβα−δLaar2d−Lβη1+α−δβdr2σ−1−Ldα−βη2α−δ,T=T11T12+T13aα−δdβLα−r2σ−1−Lβ2η2α−δ+ar2Lβη1−d/aα−δdβr−2rσ+Lα−Lβ2η2α−δ+ar2Lβη1−d2aLβ2α−δα−δdLα−rσ−1−Lβη2α−δ+Lar2η1/α−δdβ−r2σ−1+Lα−Lβ2η2α−δ+ar2Lβη1−d2,MS=Ldar+σβα−δra2r2Lβη1−d−2dσβarα−δ−αβδ−MS1aβα−δLβdα−βη2α−δ−drβ2σ−1MS2,where
(40)T11=Lar+σβα−δα−δβdrσ+Lβη2α−δ+ar2d−Lβη1α−δdβLα−r2σ−1−Lβ2η2α−δ+ar2Lβη1−d,T12=L2d2αβ2ar+σβα−δ2r+Lαα−δ2aα−δ−Lβ2η2α−δ+dβr−2rσ+Lα+ar2Lβη1−d2,T13=Lβα−δ2Lβdαrσ−δ+βδη2+drβδσ−1+Lar2dα−βδη1,MS1=β2arη2+dασLα2+β2δ2α−δ+drσα2+δ2,MS2=α−δLdα−drσ−1−Lβη2α−δ+Lar2η1.In this case, the characteristic polynomialρ3λ3+ρ2λ2+ρ1λ+ρ0 where
(41)ρ3=1,ρ2=TMSMIa+σβMS+σMI2,ρ1=−ρ11+ρ12+ρ13+ρ14+ρ15LMS+σMI3,ρ0=TMSMIρ21+ρ22−ρ23−ρ24−ρ25+ρ26LMS+σMI4,where
(42)ρ11=MSMIαMS+σMI3+TaMSMIr−Lα,ρ12=LTβMS2MI2aση1−aη2+ασ2,ρ13=LT3a2ση2MI4−β2η1MS4,ρ14=TaMSMI3rσ+2Laη2−Lασ,ρ15=LTσβMS3MIα−2βη1,ρ21=aMIMS+σMIβη1MS2+aη2MI2+2σβη1MSMI×r−Lα,ρ22=LaσβMS2MITβη1−5αη2MI,ρ23=3Laσβη2MI3Ta+ασMS,ρ24=LαβMS3MI2aη2+σ2βη1,ρ25=LTa2βη2MSMI2,ρ26=Lβ2η1MS3Ta−ασMS.By considering the previous characteristic polynomial, if one of the coefficients is zero or negative while at least one other coefficient is positive, then there is one or more imaginary roots or roots with a positive real part. If all the coefficients are positive, we calculate the Routh table. The Routh Hurwitz tabular is given byλ3λ2λ11ρ1ρ2ρ0ρ110ρ210 with ρ11=ρ1ρ2−ρ0/ρ23>0,ρ21=ρ0/ρ112>0. So, according to the Routh Hurwitz stability criterion, we can conclude that the equilibrium point in this case is stable.
(ii)
IfMI1=0 and MI2≠0, then(43)MI=12MI2MI32−4MI2MI4−MI3,MS=aMI32−4MI2MI4−MI3Lr−α−rMI32−4MI2MI4−MI32MI2arMI32−4MI2MI42MI2Lβα−δ−arMI3,T=MI3+arσ−1−LaαMI32−4MI2MI4+2LMI2ar+σβα−δaarMI32−4MI2MI4+2LMI2βα−δ−MI32×2MI3+MI32−4MI2MI4Laαr−α−arα+δ−2LMI2βδα−δ.In this case, the characteristic polynomial is as follows:ρ3λ3+ρ2λ2+ρ1λ+ρ0, where
(44)ρ3=1,ρ2=TMSMIa+σβMS+σMI2,ρ1=−ρ11+ρ12+ρ13+ρ14+ρ15LMS+σMI3,ρ0=TMSMIρ21+ρ22−ρ23−ρ24−ρ25+ρ26LMS+σMI4,where
(45)ρ11=MSMIαMS+σMI3+TaMSMIr−Lα,ρ12=LTβMS2MI2aση1−aη2+ασ2,ρ13=LT3a2ση2MI4−β2η1MS4,ρ14=TaMSMI3rσ+2Laη2−Lασ,ρ15=LTσβMS3MIα−2βη1,ρ21=aMIMS+σMIβη1MS2+aη2MI2+2σβη1MSMI×r−Lα,ρ22=LaσβMS2MITβη1−5αη2MI,ρ23=3Laσβη2MI3Ta+ασMS,ρ24=LαβMS3MI2aη2+σ2βη1,ρ25=LTa2βη2MSMI2,ρ26=Lβ2η1MS3Ta−ασMS.The Routh Hurwitz tabular is given byλ3λ2λ11ρ1ρ2ρ0ρ110ρ210 with ρ11=ρ1ρ2−ρ0/ρ23>0,ρ21=ρ0/ρ112>0. Then, according to the Routh Hurwitz stability criterion, we can conclude that the equilibrium point in this case is stable.
(iii)
IfMI1≠0,, then(46)MI=3MI1−p−4/27p3+q2/23+4/27p3+q2−p/23−MI23MI1,MS=amS1−mS2mS3−MI23MI1mS4+mS5,T=t1−t2−t3t4−MI22+3MI1t5−t6−t7Lα−δ3MI1a3LβMI1α−δ−arMI2+t72,where
(47)3LβδMI1α−δ−arMI2δ−LαLα−δ>3arMI1−p−4p3/27+q223+4p3/27+q2−p23,mS1=MI2r+Lα+3LrMI1,mS2=3MI1r+Lα4/27p3+q2−p23+−p−4/27p3+q223,mS3=3MI1−p−4/27p3+q223+4/27p3+q2−p23,mS4=3LβMI1α−δ−arMI2,mS5=3arMI1−p−4p3/27+q223+4/27p3+q2−p23,and
(48)t1=3aMI1−p−4p3/27+q223+4p3/27+q2−p23×Lα−rσ−1,t2=aMI2Lα−rσ−1,t3=3LMI1ar+σβα−δ,t4=aαr+Lα3MI1−p−4p3/27+q223+4p3/27+q2−p23,t5=3LβδMI1α−δ,t6=arMI2δ−Lα,t7=3arMI1−p−4p3/27+q223+4p3/27+q2−p23.In this case, the characteristic polynomial is as follows:ρ3λ3+ρ2λ2+ρ1λ+ρ0 where
(49)ρ3=1,ρ2=TMSMIa+σβMS+σMI2,ρ1=−ρ11+ρ12+ρ13+ρ14+ρ15LMS+σMI3,ρ0=TMSMIρ21+ρ22−ρ23−ρ24−ρ25+ρ26LMS+σMI4,where ρ11+ρ13+ρ15>ρ12+ρ14ρ11<0,ρ15<0,ρ13<0, and
(50)ρ11=MSMIαMS+σMI3+TaMSMIr−Lα,ρ12=LTβMS2MI2aση1−aη2+ασ2,ρ13=LT3a2ση2MI4−β2η1MS4,ρ14=TaMSMI3rσ+2Laη2−Lασ,ρ15=LTσβMS3MIα−2βη1,ρ21=aMIMS+σMIβη1MS2+aη2MI2+2σβη1MSMI×r−Lα,ρ22=LaσβMS2MITβη1−5αη2MI,ρ23=3Laσβη2MI3Ta+ασMS,ρ24=LαβMS3MI2aη2+σ2βη1,ρ25=LTa2βη2MSMI2,ρ26=Lβ2η1MS3Ta−ασMS.The Routh Hurwitz tabular is given byλ3λ2λ11ρ1ρ2ρ0ρ110ρ210 with ρ11=ρ1ρ2−ρ0/ρ23,ρ21=ρ0/ρ112. If p>q with q,p>0 and if t5−t6>t7Lα−δ then ρ2>0. If ρ11+ρ13+ρ15>ρ12+ρ14 then ρ1>0. On the other hand, we have ρ1ρ2>ρ0. We conclude that ρ11>0 and ρ21>0. So, according to the Routh Hurwitz stability criterion, we can conclude that the equilibrium point in this case is stable. Beginning with parameters r=0.9,L=12,α=0.9,β=0.2,σ=0.2,a=0.3,δ=0.2,η1=0.4,η2=0.3, and d=0.2 and initial value 0.1,0.1,0.1, the point P5 tends to the point 0.221;0.784;1.65×10−33. As shown in Figures 3–5. Then, we can deduce that the point P5 is stable. The previous Figures 3–6 show the dynamical behavior of the susceptible populations of Scomber colias and Thunnus thynnus, also the infected population of Scomber colias. It is clear that the biomass of these marine populations converges to the values of the interior equilibrium point. This means biologically that the interior equilibrium point can ensure the existence of the predefined species.Figure 3
Dynamical behavior of the susceptible Scomber colias.Figure 4
Dynamical behavior of the infected Scomber colias.Figure 5
Dynamical behavior of the susceptible Thunnus thynnus.Figure 6
Illustration of the evolution of the susceptible, infected Scomber colias and Thunnus thynnus fish population.Following Figure6, in a finite time the susceptible Scomber colias, infected Scomber colias and Thunnus thynnus converge to the equilibrium point P5.


## 4. Bioeconomic Model

The basic idea of this section is to define a bioeconomic model of the susceptible Scomber colias, the infected Scomber colias, and the Thunnus thynnus population exploited by three fishing fleets, and we seek to maximize the profit of each fishing fleets. The mathematical formula of catches is given byCij=qjEijxj, for all i,j=1,2,3; the catches Cij of species j by fishing fleets i are proportional to the fishing effort Eij exerted by fishing fleets i to exploit species j and the instantaneous biomass of the stock, where qj is the coefficient of catchability of species j, and xj is the biomass of population j, more precisely, x1≔MS,x2≔MI,x3≔T. The total catches of the species j by all fishing fleets are denoted by Cj=C1j+C2j+C3j. We note that Ej=∑i=13Eij represents the total fishing effort devoted to the species j by all fishing fleets, and Ei=Ei1,Ei2,Ei3T is the vector of fishing effort which must be provided by the fishing fleets i to catch the susceptible Scomber colias, infected Scomber colias and Thunnus thynnus. Then, the bioeconomic model is presented as follows
(51)ṀS=rMS1−MS+MIL−αMSMI−βMS2TMS+σMI−q1E1MS,ṀI=αMSMI−aMI2TMS+σMI−δMI−q2E2MI,Ṫ=η1βMS2TMS+σMI−η2aMI2TMS+σMI−dT−q3E3T.The biomass at biological equilibrium is the solution of the following system(52)r1−MS+MIL=αMI+βMS2TMS+σMI+q1E1,αMS=aMITMS+σMI+δ+q2E2,η1βMS2MS+σMI=η2aMI2MS+σMI+d+q3E3,MS=d+E3q3a11−a12E1+a13E3+a14E22aβη1r+Lα+drσ+rσE3q3,MI=δη1−a12/σE1+a14E2+a17E32arβη1+2Laαβη1+2adrσβη1+2arσβη1E3q3,T=−−m36α−δ−m35α−δE3+a16E2+a15E2−m31+E1m32−E2m33−E3m342Lα+drσ+rσE3q3+ra2η12β−δ+a12/σE1−a14/σE2−Lσβ2δq3E3,where
(53)a11=2ardσ+1+Lα+σδη1,a12=Laσβη1q1,a13=σq32ar+Lσβ2δη1,a14=Lσ2β2η1q2d+E3q3,a15=2arσβη1q2d+E3q3,a16=2aβη1q2r+Lα,a17=Lσβ2δη1q3,m31=da11+σδη1,m32=a12d+1,m33=a14d+η1,m34=q32ar2dσ+1+Lα+Lσ2β2δη1d+1+σδη1−q3a12E1+q3a14E2+q3a13E3,m35=q3a11−q3a12E1+q3a14E2−q3a13E3,m36=da11−da12E1+da14E2−da13E3.The profitPi of each fishing fleet i is equal to the total revenue TRi minus the total cost TCi. Mathematically, the profit for each fishing fleet is represented by Pi=TRi−TCi. Now, we give the expressions of total revenue and total costs of each fishing fleet. We use the fact that the total revenue TR depends linearly on the catch, as follows: TR=price×Catches. With the previous notations, we have
(54)TRi=p1Ci1+p2Ci2+p3Ci3=p1q1Ei1MS+p2q2Ei2MI+p3q3Ei3T=p1q1Ei1d+E3q3a11−a12E1+a13E3+a14E22aβη1r+Lα+drσ+rσE3q3+p2q2Ei2δη1−a12/σE1+a14E2+a17E32aβη1rdσ+1+Lα+rσE3q3+p3q3Ei3σm36α−δ+m35α−δE3−a16E2−a15E21×E1m32−m31−E2m33−E3m34−2a2βη12σδLσβ2E3q3+1−E1a12+E2a14rσd+E3q3+1+Lα,where pjj=1,2,3 is the price per unit biomass of the susceptible Scomber colias, infected Scomber colias and Thunnus thynnus species, with p1>p2. Now, we calculate the expression of the total effort cost, for that, we shall assume that TCi=ci,Ei, where TCi is the total cost of the fishing fleets i. Therefore, the expression of the profit is represented by the following function PiE=TRi−TCi.(55)P1E=p1C11+p2C12+p3C13−c11E11−c12E12−c13E13=p1q1E11d+E3q3a11−a12E1+a13E3+a14E22aβη1r+Lα+drσ+rσE3q3−c11E11+p2q2E12δη1−a12/σE1+a14E2+a17E32arβη1+2Laαβη1+2adrσβη1+2arσβη1E3q3−c12E12+p3q3E13σm31−E1m32+E2m33+E3m34m36α−δ−E2a15+a16+E3m35α−δ2a2βη12σδLσβ2E3q3+1−E1a12+E2a14rdσ+σd+E3q3+Lα−c13E13,P2E=p1C21+p2C22+p3C23−c21E21−c22E22−c23E23=p1q1E21d+E3q3a11−a12E1+a13E3+a14E22aβη1r+Lα+drσ+rσE3q3−c21E21+p2q2E22δη1−a12/σE1+a14E2+a17E32arβη1+2Laαβη1+2adrσβη1+2arσβη1E3q3−c22E22+p3q3E23σm31−E1m32+E2m33+E3m34m36α−δ−E2a15+a16+E3m35α−δ2a2βη12σδLσβ2E3q3+1−E1a12+E2a14rdσ+σd+E3q3+Lα−c23E23,and
(56)P3E=p1C31+p2C32+p3C33−c31E31−c32E32−c33E33=p1q1E31d+E3q3a11−a12E1+a13E3+a14E22aβη1r+Lα+drσ+rσE3q3−c31E31+p2q2E32δη1−a12/σE1+a14E2+a17E32arβη1+2Laαβη1+2adrσβη1+2arσβη1E3q3−c32E32+p3q3E33σm31−E1m32+E2m33+E3m34m36α−δ−E2a15+a16+E3m35α−δ2a2βη12σδLσβ2E3q3+1−E1a12+E2a14rdσ+σd+E3q3+Lα−c33E33.


## 5. Optimal Control

(57)ṀS=rMS1−MS+MIL−αMSMI−βMS2TMS+σMI−C1MS,ṀI=αMSMI−aMI2TMS+σMI−δMI−C2MI,Ṫ=η1βMS2TMS+σMI−η2aMI2TMS+σMI−dT−C3T.The fundamental problem from the economic point of view of the exploitation of renewable resources is to determine the optimal trade-off between present and future catches. The aim of this section is the profit-making aspect of Scomber colias and Thunnus thynnus. It is a thorough study of the optimal catches policy and the profit earned by catches, focusing on quadratic costs and conservation of this fish population by constraining the latter to always stay above a critical threshold. The reason for using quadratic costs is that it allows us to derive an analytical expression for the optimal catches; the resulting solution is different from the solution obtained case of linear cost function (the bang-bang solution). It is assumed that price is a function, which decreases with increasing biomass. Thus, to maximize the total discounted profit from the fishery, we use the controlC1, C2, and C3. It is very important to show that all the control variables are nonnegative. The objective function JC1,C2,C3 is defined as
(58)∫t0tfe−γtp1C1+C2+C3−v1C12+C22+C32−c1C1+C2+C3q1MS+q2MI+q3Tdt.Subject to the system (4), where p1 is the constant price per unit biomass of harvested population, v1 is an economic constant, γ is the instantaneous annual discount rate, and C1,C2,C3 represent the harvesting rates on infected, susceptible Scomber colias and Thunnus thynnus, respectively. The coefficient c1 is the constant fishing cost per unit effort. We are seeking to find an optimal control C1∗t,C2∗t, and C3∗t such that
(59)JC1∗,C2∗,C3∗=maxJC1,C2,C3:C1,C2,C3∈U,where U=Ci:0≤Cit≤1,Lebesguemeasurablet=t0,tffori=1,2,3 is the control set. To show the existence of the optimal control with the initial conditions t=0, we state and prove Theorems 1 and 2 below. This will also help us to analyze the properties of the system (57) with positive initial conditions. For all t>0, since the model describes infected, susceptible Scomber colias and Thunnus thynnus populations. Using the optimal control in the system (57) to see the existence of optimal control with the necessary conditions satisfying the Pontryagin’s Maximum Principle [24]. We applied Pontryagin’s Maximum Principle to convert equations (57)–(59) into a problem of maximizing Lagrange, L, with respect to C1, C2, and C3 to find the maximum value of the Lagrangian. According to Mwantobe et al. [25], this could be effectuated by considering Hamiltonian, H.
(60)H=p1C1+C2+C3−v1C12+C22+C32−c1C1+C2+C3q1MS+q2MI+q3T+θ1MStr1−MSt+MItL−αMIt−βMS2tTtMSt+σMIt−C1+θ2αMStMIt−aMI2TtMSt+σMIt−δMIt−C2MIt+θ3η1βMS2tTtMSt+σMIt−η2aMI2tTtMSt+σMIt−dTt−C3Tt,with θ1, θ2, and θ3 the adjoint variables associated to MS, MI, and T, respectively. To ensure the existence of the previous optimal control, we use the following theorem.Theorem 1 (Sufficient conditions).
The optimal control problem given by (58), along with the state equations of system (57), admits a control Ci∗∈Ufori=1,2,3 such that
(61)JC1∗,C2∗,C3∗=maxJC1,C2,C3:C1,C2,C3∈U.Proof.
See the theorems in Lashari et al. [26], Lashari and Zaman [27], and Lenhart and Workman [28].Theorem 2 (Necessary conditions).
Given optimal controlsC1∗,C2∗,C3∗ and solutions MS∗,MI∗, and T∗, there exist θ1,θ2, and θ3 the adjoint variables satisfying the following system
(62)dθ1dt=−∂H∂MS=−c1q1C1+C2+C3Tq3+q1MS+q2MI2−θ1r2MS−L+MIL+αMI+TβMSMS+2σMIMS+σMI2+C1−θ2MIαMS+σMI2+TaMIMS+σMI2−θ3Tβη1MSMS+2σMIMS+σMI2+Taη2MI2MS+σMI2,dθ2dt=−∂H∂MI=−c1q1C1+C2+C3Tq3+q1MS+q2MI2−θ1TσβMS2MS+σMI2−rMSL−αMS−θ2αMS−TaMI2MS+σMIMS+σMI2−δ+C2+θ3Tσβη1MS2MS+σMI2+Taη2MI2MS+σMIMS+σMI2,dθ3dt=−∂H∂T=−c1q1C1+C2+C3Tq3+q1MS+q2MI2+θ1βMS2MS+σMI+θ2aMI2MS+σMI−θ3βη1MS2MS+σMI−aη2MI2MS+σMI−d+C3,with the transversality conditions θitf=0,i=1,2,3. Where
(63)C1∗=minmax0,c1−θ1q1MS2∗+T∗p1q3+p1q1MS∗+p1q2MI∗−θ1q2MS∗MI∗−Tθ1q3MS∗−1T∗q3+q1MS∗+q2MI,1,C2∗=minmax0,c1−θ2q2MI2∗+T∗p1q3+p1q1MS∗+p1q2MI∗−θ2q1MS∗MI∗−T∗θ2q3MI∗−1T∗q3+q1MS∗+q2MI∗,1,C3∗=minmax0,c1+T∗p1q3+p1q1MS∗+p1q2MI∗−T2∗θ3q3−T∗θ3q1MS−Tθ3q2MI∗−1T∗q3+q1MS∗+q2MI∗,1.
According to Lashari et al. [26], if x,u is an optimal solution of an optimal control problem, then there exists a nontrivial vector function θ=θ1,⋯,θn that satisfies the following equation
(64)dxdt=∂Ht,x,u,θ∂θ;0=∂Ht,x,u,θ∂u;θ̇=∂Ht,x,u,θ∂x.Proof.
Consider the system of differential equations in (65) governing the adjoint variables θ1,θ2,θ3 This is obtained by differentiating the Hamiltonian H in equation (60) with respect to MS,MI, and T. According to Fleming and Rishel [29], these are the state variables, by applying the first and third equations in equation (64) into equation (60)
(65)dθ1dt=−∂H∂MS=−c1q1C1+C2+C3Tq3+q1MS+q2MI2−θ1r2MS−L+MIL+αMI+TβMSMS+2σMIMS+σMI2+C1−θ2MIαMS+σMI2+TaMIMS+σMI2−θ3Tβη1MSMS+2σMIMS+σMI2+Taη2MI2MS+σMI2,dθ2dt=−∂H∂MI=−c1q1C1+C2+C3Tq3+q1MS+q2MI2−θ1TσβMS2MS+σMI2−rMSL−αMS−θ2αMS−TaMI2MS+σMIMS+σMI2−δ+C2+θ3Tσβη1MS2MS+σMI2+Taη2MI2MS+σMIMS+σMI2,dθ3dt=−∂H∂T=−c1q1C1+C2+C3Tq3+q1MS+q2MI2+θ1βMS2MS+σMI+θ2aMI2MS+σMI−θ3βη1MS2MS+σMI−aη2MI2MS+σMI−d+C3.
With the transversality conditions(66)θ1tf=θ2tf=θ3tf=0.
To evaluate the optimal control of the control variable set, whereCi=0,1. Let MS=MS∗,MI=MI∗, and T=T∗, and applying the second equation in equation (64), and differentiating the Hamiltonian, H, in equation (60) with respect to the control variables C1,C2, and C3(67)∂H∂C1=−c1+θ1q1MS2−p1q1MS+Tθ1q3MS−Tp1q3−p1q2MI+θ1q2MSMI+2C1v1Tq3+q1MS+q2MITq3+q1MS+q2MI=0,∂H∂C2=−c1+θ2q2MI2−Tp1q3−p1q1MS−p1q2MI+θ2q1MSMI+2C2v1Tq3+q1MS+q2MI+Tθ2q3MITq3+q1MS+q2MI=0,∂H∂C3=−c1−Tp1q3−p1q1MS−p1q2MI+T2θ3q3+2C3v1Tq3+q1MS+q2MI+Tθ3q1MS+Tθ3q2MITq3+q1MS+q2MI=0.
By applying the optimal control to the control variable set,Ci∗=0,1 for i=1,2,3 into equation (11)
(68)C1∗=minmax0,c1−θ1q1MS2∗+T∗p1q3+p1q1MS∗+p1q2MI∗−θ1q2MS∗MI∗−Tθ1q3MS∗−1T∗q3+q1MS∗+q2MI,1,C2∗=minmax0,c1−θ2q2MI2∗+T∗p1q3+p1q1MS∗+p1q2MI∗−θ2q1MS∗MI∗−T∗θ2q3MI∗−1T∗q3+q1MS∗+q2MI∗,1,C3∗=minmax0,c1+T∗p1q3+p1q1MS∗+p1q2MI∗−T2∗θ3q3−T∗θ3q1MS−Tθ3q2MI∗−1T∗q3+q1MS∗+q2MI∗,1.
This shows that the uniqueness of the optimal control of the model has been achieved for smalltf based on prior boundedness of the state variables as well as adjoint variables.


## 6. Numerical Simulation

To study the optimal control problem numerically, we use the forward-backward Rung-Kutta sweep method. We wrote a code in MATLAB based on this method. The results are given in the following graphs. Note that all the parameters are taken as follows:r=0.9, L=12, α=0.9, β=0.2, σ=0.2, a=0.3, δ=0.2, η1=0.4, η2=0.3, d=0.2, C1=0.2, c1=0.15, C3=0.1, and C2=variable.Following Figures7–9, we can show that when the constant effort harvesting C2 on the susceptible Scomber colias is very low, i.e., C2=0.01, then susceptible and infected Scomber colias and Thunnus thynnus populations with constant harvesting effort exhibit oscillations with very long period and the densities of both populations are very high; however, under optimal harvesting, both populations are kept in a low density, but away from extension, with periodic solutions which have very short periods. In Figures 10–12, we show that when the constant harvesting on the infected Scomber colias is medium, C2=0.17, then the susceptible, infected Scomber Colias and Thunnus Thynnus populations with constant harvesting efforts keep oscillating, but in this case, the density decreases a little bit and also the period decreases in a drastic manner and the susceptible, infected Scomber Colias and Thunnus Thynnus populations have a similar behavior as the previous case.Figure 7
Susceptible Scomber colias population with optimal harvesting and effort harvesting, with very low effort.Figure 8
Infected Scomber colias population with optimal harvesting and effort harvesting, with very low effort.Figure 9
Thunnus thynnus population with optimal harvesting and effort harvesting, with very low effort.Figure 10
Susceptible Scomber colias population with optimal harvesting and effort harvesting, with medium effort.Figure 11
Infected Scomber colias population with optimal harvesting and effort harvesting, with medium effort.Figure 12
Thunnus thynnus population with optimal harvesting and effort harvesting, with medium effort.Now, forC2=0.4, i.e., for a high value of the constant harvesting on the infected Scomber colias then all the fish populations with constant effort harvesting go extinct, however, the susceptible, infected Scomber Colias and Thunnus Thynnus populations exhibit periodic solutions under optimal effort harvesting as we can see in Figures 13–15. Hence, the optimal effort harvesting might reach a maximum as high as 0.8, which could never have been achieved through constant effort harvesting.Figure 13
Susceptible Scomber colias population with optimal harvesting and effort harvesting, with high effort.Figure 14
Infected Scomber colias population with optimal harvesting and effort harvesting, with high effort.Figure 15
Thunnus thynnus population with optimal harvesting and effort harvesting, with high effort.


## 7. Conclusion

In the literature, we found that the authors consider either epidemic models or bioeconomic models but without treating the disease of fish populations. In this work, we have developed and studied a new model concept that combines the bioeconomic model, and the epidemiological model (bioeconomic-epidemiological model) of prey-predator marine populations is developed in which both susceptible and infected prey populations (Scomber colias) are exposed to the predator (Thunnus thynnus), with varying degrees of exposure. We resolved the bioeconomic model at a steady state, and we were able to calculate the profits of each fishing fleet. The optimal harvesting policy is discussed using the Maximum Principle of Pontryagin to give a high yield and keep both population Scomber colias and Thunnus thynnus away from extension.



---

*Source: 8892388-2021-03-20.xml*