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

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.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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 omitted; refer to PDF]

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 by $M_S$ per unit designated area and the infected fish population denoted by $M_I$ per unit designated area. Therefore, at time $t$, the total Scomber colias population is $Mt=M_{S}t+M_{I}t$. 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 omitted; refer to PDF]

Let $M_S$ 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 $\alpha$ 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 $\beta$ as follows $$\begin{array}{}\text{(1)}& {\dot{M}}_{S}t={rM}_{S}t1-\frac{M_{S}t+M_{I}t}{L}-\alpha M_{S}t{M}_{I}t-\frac{\beta M_S^{2}tTt}{M_{S}t+\sigma M_{I}t}.\end{array}$$

Let $M_I$ be the density of infected Scomber colias with the rate of transmission of the disease $\alpha$. This population is preyed by the Thunnus thynnus with the response rate $a$. The infected Scomber colias populations died by the death rate $\delta$ as follows $$\begin{array}{}\text{(2)}& {\dot{M}}_{I}t=\alpha M_{S}t{M}_{I}t-\frac{\alpha M_I^{2}tTt}{M_{S}t+\sigma M_{I}t}-\delta M_{I}t.\end{array}$$

The Thunnus thynnus $T$ is a predator of susceptible Scomber colias and infected Scomber colias where ${\eta }_1$ represents the rate of conversion to Thunnus thynnus of susceptible mackerel, ${\eta }_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 $$\begin{array}{}\text{(3)}& \dot{T}t=\frac{{\eta }_1\beta M_S^{2}tTt}{M_{S}t+\sigma M_{I}t}-\frac{{\eta }_{2}a{M}_I^{2}tTt}{M_{S}t+\sigma M_{I}t}-dTt.\end{array}$$

Following the previous assumptions, the biological model is formulated as follows $$\begin{array}{}\text{(4)}& \begin{array}{c}{\dot{M}}_{S}t=r{M}_{S}t1-\frac{M_{S}t+M_{I}t}{L}-\alpha M_{S}t{M}_{I}t-\frac{\beta M_S^{2}tTt}{M_{S}t+\sigma M_{I}t},\hfill \\ {\dot{M}}_{I}t=\alpha M_{S}t{M}_{I}t-\frac{a{M}_I^{2}tTt}{M_{S}t+\sigma M_{I}t}-\delta M_{I}t,\hfill \\ \dot{T}t=\frac{{\eta }_1\beta M_S^{2}tTt}{M_{S}t+\sigma M_{I}t}-\frac{{\eta }_{2}a{M}_I^{2}tTt}{M_{S}t+\sigma M_{I}t}-dTt.\hfill \end{array}\end{array}$$

We assume that all parameters are positive, $M_{S}t\ge 0,M_{I}t\ge 0,Tt\ge 0$. The system of equations (4) has a unique solution in $$\begin{array}{}\text{(5)}& V=M_S,M_I,T/M_S\ge 0,M_I\ge 0\hspace{1em}\text{and}\hspace{1em}T\ge 0,\end{array}$$since the right-hand sides of them are Lipschitz continuous in $V$.

At $M_S=M_I=0,$ the term $\beta M_S^{2}T/M_S+\sigma M_I$ is not defined. As for $M_S$ and $M_I$ are small, $0<M_I/M_S+\sigma M_I\le 1$ and $0<M_S/M_S+\sigma M_I\le 1,$ so it is natural to interpret as zero at $M_S=M_I=0.$ For $a{M}_I^{2}T/M_S+\sigma M_I,$ we have the same interpretation at $M_S=M_I=0.$

3. Biological Model Analysis

3.1. Basic Reproduction Number

We assume that we have the following SI model $$\begin{array}{}\text{(6)}& \begin{array}{c}{\dot{M}}_S=-\alpha M_S{M}_I,\hfill \\ {\dot{M}}_I=\alpha M_S{M}_I-\delta M_I.\hfill \end{array}\end{array}$$

By dividing the second equation by $\delta M_I$, we obtain ${\dot{M}}_I/\delta M_I=\alpha /\delta M_S-1.$ Then, for $\alpha /\delta M_S>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 $M_S$ is such that $\alpha /\delta M_S<1.$ The ratio $\alpha /\delta M_S$ can then be interpreted as the number of contacts that can transmit the disease by infected individuals throughout their period of contagion. If $\alpha /\delta M_S^{\ast }>1$, there will necessarily be an epidemic, whereas in the opposite case, if $\alpha /\delta M_S^{\ast }<1$, only a few individuals will be infected before the spread of the disease does not stop by itself. Let $R_0=\alpha /\delta M_{S_0},$ where $M_{S_0}=L.$ This threshold is the basic reproduction rate.

3.2. Positivity

We will denote by $$\begin{array}{}\text{(7)}& \begin{array}{c}{f}_1=r1-\frac{M_S+M_I}{L}-\alpha M_I-\frac{\beta M_{S}T}{M_S+\sigma M_I},\hfill \\ f_2=\alpha M_S-\frac{a{M}_{I}T}{M_S+\sigma M_I}-\delta ,\hfill \\ f_3=\frac{{\eta }_1\beta M_S^2}{M_S+\sigma M_I}-\frac{{\eta }_{2}a{M}_I^2}{M_S+\sigma M_I}-d.\hfill \end{array}\end{array}$$

Then, the system (4) becomes $$\begin{array}{}\text{(8)}& \begin{array}{c}{\dot{M}}_{S}t=M_{S}t{f}_1{M}_{S}t,M_{I}t,Tt,\hfill \\ {\dot{M}}_{I}t=M_{I}t{f}_2{M}_{S}t,M_{I}t,Tt,\hfill \\ \dot{T}t=Tt{f}_3{M}_{S}t,M_{I}t,Tt,\hfill \end{array}\end{array}$$subject to the initial condition M_S (0) > 0,M_I (0) > 0 and T (0) > 0. The system (8) is defined on the set $$\begin{array}{}\text{(9)}& \Omega =M_S,M_I,T\in {ℝ}^3/M_S\ge 0,M_I\ge 0,T\ge 0.\end{array}$$

Let $Xt=M_{S}t,M_{I}t,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 $$\begin{array}{}\text{(10)}& \begin{array}{c}{M}_{S}t=M_{S}0\exp {\int }_0^t{f}_1{M}_{S}s,M_{I}s,Tsds>0,\hfill \\ M_{I}t=M_{I}0\exp {\int }_0^t{f}_2{M}_{S}s,M_{I}s,Tsds>0,\hfill \\ Tt=T0\exp {\int }_0^t{f}_3{M}_{S}s,M_{I}s,Tsds>0.\hfill \end{array}\end{array}$$

Therefore, all solutions starting from an interior of the first octant remain in it for all future time.

3.3. Boundedness

We define the function $Nt=M_{S}t+M_{I}t+Tt$, and we take the time derivative along the solution of the system $\dot{N}t={\dot{M}}_{S}t+{\dot{M}}_{I}t+\dot{T}t.$ For any positive constant $\rho$, we have $$\begin{array}{}\text{(11)}& \dot{N}t+\rho Nt=r{M}_{S}t-\frac{r{M}_{S}t{M}_{S}t+M_{I}t}{L}-\frac{\beta M_S^{2}tTt+a{M}_I^{2}tTt}{M_{S}t+\sigma M_{I}t}-\delta M_{I}t+\frac{{\eta }_1\beta M_S^{2}tTt-{\eta }_{2}a{M}_I^{2}t\text{T}t}{M_{S}t+\sigma M_{I}t}-dTt+\rho M_{S}t+M_{I}t+Tt=r+\rho M_{S}t+\rho -\delta M_{I}t+\rho -dTt-\frac{1-{\eta }_1\beta M_S^{2}tTt}{M_S+\sigma M_I}-\frac{a+{\eta }_2{M}_I^{2}tTt}{M_S+\sigma M_I}-\frac{r{M}_{S}t{M}_{S}t+M_{I}t}{L}\le r+\rho M_{S}t+\rho -\delta M_{I}t+\rho -dTt.\end{array}$$

We choose $\rho$ as a minimum of $\delta ,d$. On the one hand, we have $$\begin{array}{}\text{(12)}& {\dot{M}}_{S}t+{\dot{M}}_{I}t\le \frac{r}{\alpha }1-\frac{M_{S}t+M_{I}t}{L}.\end{array}$$

Consequently, $\lim {\sup }_{t\to \infty }{M}_{S}t+M_{I}t\le L.$ Then, $\lim {\sup }_{t\to \infty }{M}_{S}t\le L$ and $\lim {\sup }_{t\to \infty }{M}_{I}t\le L.$ That improves that the susceptible and infected preys are always bounded. On the other hand, there exists $t_0$ such that $t\ge t_0,M_S\le 2L$, thus, for $t\ge t_0$$$\begin{array}{}\text{(13)}& \dot{N}t+\rho Nt\le 2L\frac{r}{\alpha }+\rho =\omega .\end{array}$$

Which gives, for $$\begin{array}{}\text{(14)}& Nt\le \frac{\omega }{\rho }+N{t}_0{e}^{-\rho t-t_0}.\end{array}$$

Hence, that $\lim {\sup }_{t\to \infty }Nt\le \omega /\rho$ 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 $$\begin{array}{}\text{(15)}& \begin{array}{c}r1-\frac{M_S+M_I}{L}-\alpha M_I-\frac{\beta M_{S}T}{M_S+\sigma M_I}=0,\hfill \\ \alpha M_S-\frac{a{M}_{I}T}{M_S+\sigma M_I}-\delta =0,\hfill \\ \frac{{\eta }_1\beta M_S^2}{M_S+\sigma M_I}-\frac{{\eta }_{2}a{M}_I^2}{M_S+\sigma M_I}-d=0.\hfill \end{array}\end{array}$$

After calculation, it is obvious that system (15) has five equilibrium points as follows $$\begin{array}{}\text{(16)}& \begin{array}{c}{P}_1=0,0,0,\hfill \\ P_2=\frac{d}{{\eta }_1\beta },0,r\frac{-d+L\beta {\eta }_1}{L{\beta }^2{\eta }_1},\hfill \\ P_3=\frac{\delta }{\alpha },-\frac{Lr}{L\alpha r+l\alpha }\delta -L\alpha ,0,\hfill \\ P_4=L,0,0,\hfill \\ P_5{M}_S,M_I,T.\hfill \end{array}\end{array}$$

For the last equilibrium point $P_5{M}_S,M_I,T$ is given by $$\begin{array}{}\text{(17)}& \begin{array}{c}{M}_S=-\frac{a{M}_{I}r{M}_I-Lr+L\alpha M_I}{ar{M}_I+L\alpha \beta -L\beta \delta },\hfill \\ T=\frac{aL\alpha -r\sigma -1M_I-Lar+\sigma \beta \alpha -\delta a\alpha M_I^{2}r+L\alpha +L\beta \delta \alpha -\delta +ar{M}_I\delta -L\alpha }{a{ar{M}_I+L\beta \alpha -\delta }^2},\hfill \end{array}\end{array}$$where $M_I$ is given in the following calculation.

Since $M_S=a{M}_I-r{M}_I+Lr-L\alpha M_I/ar{M}_I+L\alpha \beta -L\beta \delta \ge 0,$ so we obtain $M_I\le Lr/r+L\alpha$. Hence, there is no economic equilibrium when $M_I>Lr/r+L\alpha .$

Consider the function $g$ defined on an interval $0,Lr/r+L\alpha$ by $$\begin{array}{}\text{(18)}& g{M}_I=M_{I1}{M}_I^3+M_{I2}{M}_I^2+M_{I3}{M}_I+M_{I4}=0,\end{array}$$where $$\begin{array}{}\text{(19)}& \begin{array}{c}{M}_{I1}=a^2{r}^{2}a{\eta }_2-\beta {\eta }_1-L\alpha \beta {\eta }_{1}2r+L\alpha ,\hfill \\ M_{I2}=a^{2}rdr\sigma -1-Ld\alpha +2L\beta {\eta }_{1}r+L\alpha +2L\beta {\eta }_2\alpha -\delta ,\hfill \\ M_{I3}=Laa{r}^{2}d-L\beta {\eta }_1+\alpha -\delta L{\beta }^2{\eta }_2\alpha -\delta -d\beta L\alpha -r2\sigma -1,\hfill \\ M_{I4}=L^{2}d\beta \alpha -\delta ar+\sigma \beta \alpha -\delta .\hfill \end{array}\end{array}$$

If $M_{I4}=0,$ i.e., $\alpha =\delta$, then, equation (18) admits two different solutions: the first one is trivial and the other solution exists if $M_{I3}<0$ and $M_{I1}>0.$ When $\alpha <\delta$, we have $g0=L^{2}d\beta \alpha -\delta ar+\sigma \beta \alpha -\delta <0$ with $ar+\sigma \beta \alpha -\delta >0$ and $$\begin{array}{}\text{(20)}& g\frac{Lr}{r+L\alpha }=M_{I1}{\frac{Lr}{r+L\alpha }}^3+M_{I2}{\frac{Lr}{\text{r}+L\alpha }}^2+M_{I3}\frac{Lr}{r+L\alpha }+M_{I4}=L^2{a{r}^2+\beta \alpha -\delta r+L\alpha }^{2}d\sigma r+L\alpha +Lar{\eta }_2>0.\end{array}$$

Equation (18) admits at least one solution $M_I^{\ast }\in 0,Lr/r+L\alpha$ which implies that system (1) has at least one capital-labor equilibrium $M_S^{\ast },M_I^{\ast },T^{\ast }$ with $$\begin{array}{}\text{(21)}& \begin{array}{c}{M}_S^{\ast }=a{M}_I^{\ast }\frac{Lr-r+L\alpha M_I^{\ast }}{ar{M}_I^{\ast }+L\beta \alpha -\delta },\hfill \\ T^{\ast }=\frac{ar1-\sigma +L\alpha M_I^{\ast }-Lar+\sigma \beta \alpha -\delta a\alpha r+L\alpha M_I^{\ast 2}+ar\delta -L\alpha M_I^{\ast }+L\beta \delta \alpha -\delta }{a{ar{M}_I^{\ast }+L\alpha \beta -L\beta \delta }^2}.\hfill \end{array}\end{array}$$

To determine the uniqueness of the solution $M_I^{\ast }$, we discuss the following cases

(i) If $M_{I1}=M_{I2}=0,$ then the equation (18) has a unique positive solution given by

$$\begin{array}{c}\text{(22)}& M_I^{\ast }=\frac{L^{2}d\beta \delta -\alpha ar+\sigma \beta \alpha -\delta }{Laa{r}^{2}d-L\beta {\eta }_1+\alpha -\delta \beta dr2\sigma -1-Ld\alpha -\beta {\eta }_2\alpha -\delta },d-L\beta {\eta }_1<0,\alpha -\delta >0\hspace{1em}\text{and}\hspace{1em}dr2\sigma -1<Ld\alpha -\beta {\eta }_2\alpha -\delta \hspace{1em}\text{then}\hspace{1em}{M}_I^{\ast }>0.\end{array}$$

(ii) If $M_{I1}=0$ and $M_{I2}\ne 0,$ then equation (18) has a unique positive solution given by $M_I^{\ast }=1/2M_{I2}\sqrt{M_{I3}^2-4M_{I2}{M}_{I4}}-M_{I3}$

(iii) If $M_{I1}\ne 0$ by applying the rule signs of Descartes, equation (18) has a unique positive root $M_I^{\ast },$ if any of the following three conditions is satisfied

$$\begin{array}{}\text{(23)}& \begin{array}{c}{M}_{I2}>0\hspace{1em}\text{and}\hspace{1em}{M}_{I3}>0,\\ M_{I2}>0\hspace{1em}\text{and}\hspace{1em}{M}_{I3}<0,\\ M_{I2}<0\hspace{1em}\text{and}\hspace{1em}{M}_{I3}<0.\end{array}\end{array}$$

To depress the cubic equation, we substitute $M_I=x-M_{I2}/3M_{I1}$ and we solve $g{M}_I=0$. Then, we get $$\begin{array}{}\text{(24)}& \begin{array}{c}{M}_{I1}{x-\frac{M_{I2}}{3M_{I1}}}^3+M_{I2}{x-\frac{M_{I2}}{3M_{I1}}}^2+M_{I3}x-\frac{M_{I2}}{3M_{I1}}+M_{I4}=0,\hfill \\ {x-\frac{M_{I2}}{3M_{I1}}}^3+\frac{M_{I2}}{M_{I1}}{x-\frac{M_{I2}}{3M_{I1}}}^2+\frac{M_{I3}}{M_{I1}}x-\frac{M_{I2}}{3M_{I1}}+\frac{M_{I4}}{M_{I1}}=0,\hfill \\ x^3-\frac{M_{I2}^3}{27M_{I1}^3}-\frac{M_{I2}}{M_{I1}}{x}^2+\frac{M_{I2}^2}{3M_{I1}^2}x+\frac{M_{I2}}{M_{I1}}\frac{M_{I2}^2}{9M_{I1}^2}+x^2-\frac{2M_{I2}}{3M_{I1}}x\hfill \\ +\frac{M_{I3}}{M_{I1}}x-\frac{M_{I2}}{3M_{I1}}+\frac{M_{I4}}{M_{I1}}=0,\hfill \\ x^3+\frac{3M_{I1}{M}_{I3}-M_{I2}^2}{3M_{I1}^2}x+\frac{27M_{I1}^2{M}_{I4}+2M_{I2}^3-9M_{I1}{M}_{I2}{M}_{I3}}{27M_{I1}^3}=0.\hfill \end{array}\end{array}$$

Then, we obtain $gx=x^3+px+q$ where $$\begin{array}{}\text{(25)}& \begin{array}{c}p=\frac{3M_{I1}{M}_{I3}-M_{I2}^2}{3M_{I1}^2},\hfill \\ q=\frac{27M_{I1}^2{M}_{I4}+2M_{I2}^3-9M_{I1}{M}_{I2}{M}_{I3}}{27M_{I1}^3}.\hfill \end{array}\end{array}$$

By using Cardano’s formula, we obtain ${\Delta }_1=q^2+4p^3/27$ then $$\begin{array}{}\text{(26)}& M_I^{\ast }=\sqrt[3]{\frac{-p-\sqrt{q^2+4p^3/27}}{2}}+\sqrt[3]{\frac{-p+\sqrt{q^2+4p^3/27}}{2}}-\frac{M_{I2}}{3M_{I1}}.\end{array}$$

3.5. Local Stability of Equilibrium Points

To ensure the positivity of the population’s biomass, we assure that $$\begin{array}{}\text{(27)}& \begin{array}{c}L\alpha >\delta ,\hfill \\ L\beta {\eta }_1>d,\hfill \\ L\alpha >r,\hfill \end{array}\end{array}$$

We consider the matrix $JS$ that represents the variational matrix of the system at the equilibrium point $P=M_S,M_I,T$$$\begin{array}{}\text{(28)}& JP=\begin{array}{ccc}{J}_{11}& J_{12}& J_{13}\\ J_{21}& J_{22}& J_{23}\\ J_{31}& J_{32}& J_{33}\end{array},\end{array}$$where $$\begin{array}{}\text{(29)}& \begin{array}{c}{J}_{11}=r1-\frac{M_S+M_I}{L}-\alpha M_I-\frac{2\beta M_{S}T{M}_S+\sigma M_I-\beta M_S^{2}T}{{M_S+\sigma M_I}^2},\hfill \\ J_{12}=\frac{1}{L}r{M}_S-\alpha M_S+\frac{\sigma \beta M_S^{2}T}{{M_S+\sigma M_I}^2},\hfill \\ J_{13}=-\frac{\beta M_S^2}{M_S+\sigma M_I},\hfill \\ J_{21}=\alpha M_I-\frac{-a{M}_I^{2}T}{{M_S+\sigma M_I}^2},\hfill \\ J_{22}=\alpha M_S-\frac{2a{M}_{I}T{M}_S+\sigma M_I-a\sigma M_I^{2}T}{{M_S+\sigma M_I}^2}-\delta ,\hfill \\ J_{23}=-\frac{a{M}_I^2}{M_S+\sigma M_I},\hfill \\ J_{31}=\frac{2{\eta }_1\beta M_{S}T{M}_S+\sigma M_I-{\eta }_1\beta M_S^{2}T+{\eta }_{2}a{M}_I^{2}T}{{M_S+\sigma M_I}^2},\hfill \\ J_{32}=\frac{-{\eta }_1\sigma \beta M_S^{2}T-2{\eta }_{2}a{M}_{I}T{M}_S+\sigma M_I-{\eta }_2\sigma a{M}_I^{2}T}{{M_S+\sigma M_I}^2},\hfill \\ J_{33}=\frac{{\eta }_1\beta M_S^2}{M_S+\sigma M_I}-\frac{{\eta }_{2}a{M}_I^2}{M_S+\sigma M_I}-d.\hfill \end{array}\end{array}$$

(1) At the equilibrium point $P_1=0,0,0$, the variational matrix is

$$\begin{array}{}\text{(30)}& J{P}_1=\begin{array}{ccc}r& 0& 0\\ 0& -\delta & 0\\ 0& 0& -d\end{array}.\end{array}$$

It has three eigenvalues: ${\lambda }_1=r>0,{\lambda }_2=-\delta <0$, and ${\lambda }_3=-d<0.$ Consequently, the point $P_1$ is unstable.

(2) At $P_2=d/{\eta }_1\beta ,0,r-d+L\beta {\eta }_1/L{\beta }^2{\eta }_1$ the variational matrix is

$$\begin{array}{}\text{(31)}& J{P}_2=\begin{array}{ccc}0& \frac{dr-Ld\alpha -dr\sigma +Lr\sigma \beta {\eta }_1}{L\beta {\eta }_1}& -\frac{d}{{\eta }_1}\\ 0& \frac{d\alpha -\beta \delta {\eta }_1}{\beta {\eta }_1}& 0\\ \frac{rL\beta {\eta }_1-d}{L\beta }& \frac{r\sigma d-L\beta {\eta }_1}{L\beta }& 0\end{array}.\end{array}$$

It has these eigenvalues: $$\begin{array}{}\text{(32)}& \begin{array}{c}{\lambda }_1=-{\beta }^2{\eta }_1\sqrt{\frac{dr}{L{\beta }^5{\eta }_1^3}d-L\beta {\eta }_1}<0,\hfill \\ {\lambda }_2={\beta }^2{\eta }_1\sqrt{\frac{dr}{L{\beta }^5{\eta }_1^3}d-L\beta {\eta }_1}>0,\hfill \\ {\lambda }_3=\frac{d\alpha -\beta \delta {\eta }_1}{\beta {\eta }_1}>0.\hfill \end{array}\end{array}$$

Hence, the point $P_2$ is unstable.

(3) In the same manner, at $P_3=\delta /\alpha ,-Lr\delta -L\alpha /L\alpha r+l\alpha ,0$, the variational matrix is

$$\begin{array}{}\text{(33)}& J{P}_3=\begin{array}{ccc}0& \delta \frac{r-L\alpha }{L\alpha }& \frac{-\beta {\delta }^{2}Lr+L\alpha }{\alpha Lr\delta +LL\alpha \delta -Lr\sigma \delta +LLr\alpha \sigma }\\ r\frac{-\delta +L\alpha }{r+L\alpha }& 0& -a{L}^2{r}^2\frac{{-\delta +L\alpha }^2}{L\alpha r+L\alpha Lr\delta +LL\alpha \delta -Lr\sigma \delta +LLr\alpha \sigma }\\ 0& 0& \frac{{\eta }_1\beta {\delta /\alpha }^2}{\delta /\alpha +\sigma -Lr\delta -L\alpha /L\alpha r+L\alpha }-\frac{{\eta }_{2}a{-lr\delta -L\alpha /L\alpha r+l\alpha }^2}{\delta /\alpha +\sigma -Lr\delta -L\alpha /L\alpha r+L\alpha }-d\end{array}.\end{array}$$

It has these eigenvalues: $$\begin{array}{}\text{(34)}& \begin{array}{c}{\lambda }_1=L^7{\alpha }^7{r+L\alpha }^3\sqrt{\frac{r\delta L\alpha -r\delta -L\alpha }{L^{15}{\alpha }^{15}{r+L\alpha }^7}},\hfill \\ {\lambda }_2=-L^7{\alpha }^7{r+L\alpha }^3\sqrt{\frac{r\delta L\alpha -r\delta -L\alpha }{L^{15}{\alpha }^{15}{r+L\alpha }^7}},\hfill \\ {\lambda }_3=\frac{{\lambda }_{31}+{\lambda }_{32}-{\lambda }_{33}}{L^2{\alpha }^3\delta +r^2\alpha \delta -r^2\alpha \sigma \delta +L{r}^2{\alpha }^2\sigma +L^{2}r{\alpha }^3\sigma +2Lr{\alpha }^2\delta -Lr{\alpha }^2\sigma \delta },\hfill \end{array}\end{array}$$where $$\begin{array}{}\text{(35)}& \begin{array}{c}{\lambda }_{31}=L^2{\alpha }^2-\beta {\delta }^2{\eta }_1+d\alpha \delta +a{r}^2{\eta }_2+dr\alpha \sigma ,\hfill \\ {\lambda }_{32}=r^2\delta d\alpha -d\alpha \sigma +a\delta {\eta }_2-\beta \delta {\eta }_1,\hfill \\ {\lambda }_{33}=Lr\alpha 2\beta {\delta }^2{\eta }_1-2d\alpha \delta +2ar\delta {\eta }_2-dr\alpha \sigma +d\alpha \sigma \delta ,\hfill \end{array}\end{array}$$

Then, following the conditions ${\lambda }_1<0,{\lambda }_2>0$ and ${\lambda }_3>0$, we can deduce that the point $P_3$ is unstable.

(4) The variational matrix at $P_4=L,0,0$ is given by

$$\begin{array}{}\text{(36)}& J{P}_4=\begin{array}{ccc}0& r-L\alpha & -L\beta \\ 0& \alpha L-\delta & 0\\ 0& 0& L\beta {\eta }_1-d\end{array}.\end{array}$$

Since one of the eigenvalues ${\lambda }_1=L\alpha -\delta >0$ is positive Then, the point $P_4$ is unstable.

(5) The local stability of $P_5{M}_S,M_I,T$ is investigated by considering the variational matrix based on the system (1) with

$$\begin{array}{}\text{(37)}& J{P}_5=\begin{array}{ccc}{J}_{11}& J_{12}& -\frac{\beta M_S^2}{M_S+\sigma M_I}\\ J_{21}& -\frac{Ta{M}_S{M}_I}{{M_S+\sigma M_I}^2}& -\frac{a{M}_I^2}{M_S+\sigma M_I}\\ J_{31}& J_{32}& -d\end{array},\end{array}$$where $$\begin{array}{}\text{(38)}& \begin{array}{c}{J}_{11}=-M_S\frac{r{M_S+\sigma M_I}^2+LT\sigma \beta M_I}{L{M_S+\sigma M_I}^2},\hfill \\ J_{12}=-M_S\frac{r+L\alpha }{L}-\frac{\sigma \beta M_S^{2}T}{{M_S+\sigma M_I}^2},\hfill \\ J_{21}=\alpha M_I+\frac{a{M}_I^{2}T}{M_S+\sigma M_I},\hfill \\ J_{31}=T\frac{\beta {\eta }_1{M}_S^2-a{\eta }_2{M}_I^2\sigma M_I-1-M_S{M}_{I}a{\eta }_2{M}_I-2\sigma \beta {\eta }_1}{{M_S+\sigma M_I}^2},\hfill \\ J_{32}=-\frac{T\sigma \beta {\eta }_1{M}_S^2+2a{\eta }_2{M}_S{M}_I+a\sigma {\eta }_2{M}_I^2}{{M_S+\sigma M_I}^2}.\hfill \end{array}\end{array}$$

(i) If $M_{I1}=M_{I2}=0,$ then

$$\begin{array}{}\text{(39)}& \begin{array}{cc}{M}_I\hfill & =\frac{L^{2}d\beta \delta -\alpha ar+\sigma \beta \alpha -\delta }{Laa{r}^{2}d-L\beta {\eta }_1+\alpha -\delta \beta dr2\sigma -1-Ld\alpha -\beta {\eta }_2\alpha -\delta },\hfill \\ T\hfill & =\frac{T_{11}{T}_{12}+T_{13}a\alpha -\delta d\beta L\alpha -r2\sigma -1-L{\beta }^2{\eta }_2\alpha -\delta +a{r}^{2}L\beta {\eta }_1-d/a{\alpha -\delta d\beta r-2r\sigma +L\alpha -L{\beta }^2{\eta }_2\alpha -\delta +a{r}^{2}L\beta {\eta }_1-d}^2}{a{L{\beta }^2\alpha -\delta \alpha -\delta dL\alpha -r\sigma -1-L\beta {\eta }_2\alpha -\delta +La{r}^2{\eta }_1/\alpha -\delta d\beta -r2\sigma -1+L\alpha -L{\beta }^2{\eta }_2\alpha -\delta +a{r}^{2}L\beta {\eta }_1-d}^2},\hfill \\ M_S\hfill & =\frac{Ldar+\sigma \beta \alpha -\delta r{a}^2{r}^{2}L\beta {\eta }_1-d-2d\sigma \beta ar\alpha -\delta -\alpha \beta \delta -M_{S1}}{a\beta \alpha -\delta L\beta d\alpha -\beta {\eta }_2\alpha -\delta -dr\beta 2\sigma -1M_{S2}},\hfill \end{array}\end{array}$$where $$\begin{array}{}\text{(40)}& \begin{array}{c}{T}_{11}=\frac{Lar+\sigma \beta \alpha -\delta \alpha -\delta \beta dr\sigma +L\beta {\eta }_2\alpha -\delta +a{r}^{2}d-L\beta {\eta }_1}{\alpha -\delta d\beta L\alpha -r2\sigma -1-L{\beta }^2{\eta }_2\alpha -\delta +a{r}^{2}L\beta {\eta }_1-d},\hfill \\ T_{12}=\frac{L^2{d}^2\alpha {\beta }^2{ar+\sigma \beta \alpha -\delta }^{2}r+L\alpha {\alpha -\delta }^2}{a{\alpha -\delta -L{\beta }^2{\eta }_2\alpha -\delta +d\beta r-2r\sigma +L\alpha +a{r}^{2}L\beta {\eta }_1-d}^2},\hfill \\ T_{13}=L\beta {\alpha -\delta }^{2}L\beta d\alpha r\sigma -\delta +\beta \delta {\eta }_2+dr\beta \delta \sigma -1\hfill \\ +La{r}^{2}d\alpha -\beta \delta {\eta }_1,\hfill \\ M_{S1}={\beta }^{2}ar{\eta }_2+d\alpha \sigma L{\alpha }^2+{\beta }^2\delta 2\alpha -\delta +dr\sigma {\alpha }^2+{\delta }^2,\hfill \\ M_{S2}=\alpha -\delta Ld\alpha -dr\sigma -1-L\beta {\eta }_2\alpha -\delta +La{r}^2{\eta }_1.\hfill \end{array}\end{array}$$

In this case, the characteristic polynomial ${\rho }_3{\lambda }^3+{\rho }_2{\lambda }^2+{\rho }_1\lambda +{\rho }_0$ where $$\begin{array}{}\text{(41)}& \begin{array}{c}{\rho }_3=1,\hfill \\ {\rho }_2=\frac{T{M}_S{M}_{I}a+\sigma \beta }{{M_S+\sigma M_I}^2},\hfill \\ {\rho }_1=-\frac{{\rho }_{11}+{\rho }_{12}+{\rho }_{13}+{\rho }_{14}+{\rho }_{15}}{L{M_S+\sigma M_I}^3},\hfill \\ {\rho }_0=T{M}_S{M}_I\frac{{\rho }_{21}+{\rho }_{22}-{\rho }_{23}-{\rho }_{24}-{\rho }_{25}+{\rho }_{26}}{L{M_S+\sigma M_I}^4},\hfill \end{array}\end{array}$$where $$\begin{array}{}\text{(42)}& \begin{array}{c}{\rho }_{11}=M_S{M}_I\alpha {M_S+\sigma M_I}^3+Ta{M}_S{M}_{I}r-L\alpha ,\hfill \\ {\rho }_{12}=LT\beta M_S^2{M}_I^{2}a\sigma {\eta }_1-a{\eta }_2+\alpha {\sigma }^2,\hfill \\ {\rho }_{13}=LT3a^2\sigma {\eta }_2{M}_I^4-{\beta }^2{\eta }_1{M}_S^4,\hfill \\ {\rho }_{14}=Ta{M}_S{M}_I^{3}r\sigma +2La{\eta }_2-L\alpha \sigma ,\hfill \\ {\rho }_{15}=LT\sigma \beta M_S^3{M}_I\alpha -2\beta {\eta }_1,\hfill \\ {\rho }_{21}=a{M}_I{M}_S+\sigma M_I\beta {\eta }_1{M}_S^2+a{\eta }_2{M}_I^2+2\sigma \beta {\eta }_1{M}_S{M}_I\times r-L\alpha ,\hfill \\ {\rho }_{22}=La\sigma \beta M_S^2{M}_{I}T\beta {\eta }_1-5\alpha {\eta }_2{M}_I,\hfill \\ {\rho }_{23}=3La\sigma \beta {\eta }_2{M}_I^{3}Ta+\alpha \sigma M_S,\hfill \\ {\rho }_{24}=L\alpha \beta M_S^3{M}_{I}2a{\eta }_2+{\sigma }^2\beta {\eta }_1,\hfill \\ {\rho }_{25}=LT{a}^2\beta {\eta }_2{M}_S{M}_I^2,\hfill \\ {\rho }_{26}=L{\beta }^2{\eta }_1{M}_S^{3}Ta-\alpha \sigma M_S.\hfill \end{array}\end{array}$$

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 $\begin{array}{c}{\lambda }^3\\ {\lambda }^2\\ \lambda \\ 1\end{array}\begin{array}{cc}1\hfill & {\rho }_1\hfill \\ {\rho }_2\hfill & {\rho }_0\hfill \\ {\rho }_{11}\hfill & 0\hfill \\ {\rho }_{21}\hfill & 0\hfill \end{array}$ with ${\rho }_{11}={\rho }_1{\rho }_2-{\rho }_0/{\rho }_2^3>0,{\rho }_{21}={\rho }_0/{\rho }_{11}^2>0.$ So, according to the Routh Hurwitz stability criterion, we can conclude that the equilibrium point in this case is stable.

(ii) If $M_{I1}=0$ and $M_{I2}\ne 0,$ then

$$\begin{array}{}\text{(43)}& \begin{array}{c}{M}_I=\frac{1}{2M_{I2}}\sqrt{M_{I3}^2-4M_{I2}{M}_{I4}}-M_{I3},\hfill \\ M_S=\frac{a\sqrt{M_{I3}^2-4M_{I2}{M}_{I4}}-M_{I3}Lr-\alpha -r\sqrt{M_{I3}^2-4M_{I2}{M}_{I4}}-M_{I3}}{2M_{I2}ar\sqrt{M_{I3}^2-4M_{I2}{M}_{I4}}2M_{I2}L\beta \alpha -\delta -ar{M}_{I3}},\hfill \\ T=\frac{M_{I3}+ar\sigma -1-La\alpha \sqrt{M_{I3}^2-4M_{I2}{M}_{I4}}+2L{M}_{I2}ar+\sigma \beta \alpha -\delta }{a{ar\sqrt{M_{I3}^2-4M_{I2}{M}_{I4}}+2L{M}_{I2}\beta \alpha -\delta -M_{I3}}^2}\times \hfill \\ 2M_{I3}+\sqrt{M_{I3}^2-4M_{I2}{M}_{I4}}La\alpha r-\alpha -ar\alpha +\delta -2L{M}_{I2}\beta \delta \alpha -\delta .\hfill \end{array}\end{array}$$

In this case, the characteristic polynomial is as follows: ${\rho }_3{\lambda }^3+{\rho }_2{\lambda }^2+{\rho }_1\lambda +{\rho }_0$, where $$\begin{array}{}\text{(44)}& \begin{array}{c}{\rho }_3=1,\hfill \\ {\rho }_2=\frac{T{M}_S{M}_{I}a+\sigma \beta }{{M_S+\sigma M_I}^2},\hfill \\ {\rho }_1=-\frac{{\rho }_{11}+{\rho }_{12}+{\rho }_{13}+{\rho }_{14}+{\rho }_{15}}{L{M_S+\sigma M_I}^3},\hfill \\ {\rho }_0=T{M}_S{M}_I\frac{{\rho }_{21}+{\rho }_{22}-{\rho }_{23}-{\rho }_{24}-{\rho }_{25}+{\rho }_{26}}{L{M_S+\sigma M_I}^4},\hfill \end{array}\end{array}$$where $$\begin{array}{}\text{(45)}& \begin{array}{c}{\rho }_{11}=M_S{M}_I\alpha {M_S+\sigma M_I}^3+Ta{M}_S{M}_{I}r-L\alpha ,\hfill \\ {\rho }_{12}=LT\beta M_S^2{M}_I^{2}a\sigma {\eta }_1-a{\eta }_2+\alpha {\sigma }^2,\hfill \\ {\rho }_{13}=LT3a^2\sigma {\eta }_2{M}_I^4-{\beta }^2{\eta }_1{M}_S^4,\hfill \\ {\rho }_{14}=Ta{M}_S{M}_I^{3}r\sigma +2La{\eta }_2-L\alpha \sigma ,\hfill \\ {\rho }_{15}=LT\sigma \beta M_S^3{M}_I\alpha -2\beta {\eta }_1,\hfill \\ {\rho }_{21}=a{M}_I{M}_S+\sigma M_I\beta {\eta }_1{M}_S^2+a{\eta }_2{M}_I^2+2\sigma \beta {\eta }_1{M}_S{M}_I\times r-L\alpha ,\hfill \\ {\rho }_{22}=La\sigma \beta M_S^2{M}_{I}T\beta {\eta }_1-5\alpha {\eta }_2{M}_I,\hfill \\ {\rho }_{23}=3La\sigma \beta {\eta }_2{M}_I^{3}Ta+\alpha \sigma M_S,\hfill \\ {\rho }_{24}=L\alpha \beta M_S^3{M}_{I}2a{\eta }_2+{\sigma }^2\beta {\eta }_1,\hfill \\ {\rho }_{25}=LT{a}^2\beta {\eta }_2{M}_S{M}_I^2,\hfill \\ {\rho }_{26}=L{\beta }^2{\eta }_1{M}_S^{3}Ta-\alpha \sigma M_S.\hfill \end{array}\end{array}$$

The Routh Hurwitz tabular is given by $\begin{array}{c}{\lambda }^3\\ {\lambda }^2\\ \lambda \\ 1\end{array}\begin{array}{cc}1\hfill & {\rho }_1\hfill \\ {\rho }_2\hfill & {\rho }_0\hfill \\ {\rho }_{11}\hfill & 0\hfill \\ {\rho }_{21}\hfill & 0\hfill \end{array}$ with ${\rho }_{11}={\rho }_1{\rho }_2-{\rho }_0/{\rho }_2^3>0,{\rho }_{21}={\rho }_0/{\rho }_{11}^2>0.$ Then, according to the Routh Hurwitz stability criterion, we can conclude that the equilibrium point in this case is stable.

(iii) If $M_{I1}\ne 0,$, then

$$\begin{array}{}\text{(46)}& \begin{array}{c}{M}_I=\frac{3M_{I1}\sqrt[3]{-p-\sqrt{4/27p^3+q^2}/2}+\sqrt[3]{\sqrt{4/27p^3+q^2}-p/2}-M_{I2}}{3M_{I1}},\hfill \\ M_S=\frac{a{m}_{S1}-m_{S2}{m}_{S3}-M_{I2}}{3M_{I1}{m}_{S4}+m_{S5}},\hfill \\ T=\frac{t_1-t_2-t_3{t_4-M_{I2}}^2+3M_{I1}{t}_5-t_6-t_{7}L\alpha -\delta }{3M_{I1}a{3L\beta M_{I1}\alpha -\delta -ar{M}_{I2}+t_7}^2},\hfill \end{array}\end{array}$$where $$\begin{array}{c}\text{(47)}& \frac{3L\beta \delta M_{I1}\alpha -\delta -ar{M}_{I2}\delta -L\alpha }{L\alpha -\delta }>3ar{M}_{I1}\sqrt[3]{\frac{-p-\sqrt{4p^3/27+q^2}}{2}}+\sqrt[3]{\frac{\sqrt{4p^3/27+q^2}-p}{2}},\begin{array}{c}{m}_{S1}=M_{I2}r+L\alpha +3Lr{M}_{I1},\hfill \\ m_{S2}=3M_{I1}r+L\alpha \sqrt[3]{\frac{\sqrt{4/27p^3+q^2}-p}{2}}+\sqrt[3]{\frac{-p-\sqrt{4/27p^3+q^2}}{2}},\hfill \\ m_{S3}=3M_{I1}\sqrt[3]{\frac{-p-\sqrt{4/27p^3+q^2}}{2}}+\sqrt[3]{\frac{\sqrt{4/27p^3+q^2}-p}{2}},\hfill \\ m_{S4}=3L\beta M_{I1}\alpha -\delta -ar{M}_{I2},\hfill \\ m_{S5}=3ar{M}_{I1}\sqrt[3]{\frac{-p-\sqrt{4p^3/27+q^2}}{2}}+\sqrt[3]{\frac{\sqrt{4/27p^3+q^2}-p}{2}},\hfill \end{array}\end{array}$$and $$\begin{array}{}\text{(48)}& \begin{array}{c}{t}_1=3a{M}_{I1}\sqrt[3]{\frac{-p-\sqrt{4p^3/27+q^2}}{2}}+\sqrt[3]{\frac{\sqrt{4p^3/27+q^2}-p}{2}}\times L\alpha -r\sigma -1,\hfill \\ t_2=a{M}_{I2}L\alpha -r\sigma -1,\hfill \\ t_3=3L{M}_{I1}ar+\sigma \beta \alpha -\delta ,\hfill \\ t_4=a\alpha r+L\alpha 3M_{I1}\sqrt[3]{\frac{-p-\sqrt{4p^3/27+q^2}}{2}}+\sqrt[3]{\frac{\sqrt{4p^3/27+q^2}-p}{2}},\hfill \\ t_5=3L\beta \delta M_{I1}\alpha -\delta ,\hfill \\ t_6=ar{M}_{I2}\delta -L\alpha ,\hfill \\ t_7=3ar{M}_{I1}\sqrt[3]{\frac{-p-\sqrt{4p^3/27+q^2}}{2}}+\sqrt[3]{\frac{\sqrt{4p^3/27+q^2}-p}{2}}.\hfill \end{array}\end{array}$$

In this case, the characteristic polynomial is as follows: ${\rho }_3{\lambda }^3+{\rho }_2{\lambda }^2+{\rho }_1\lambda +{\rho }_0$ where $$\begin{array}{}\text{(49)}& \begin{array}{c}{\rho }_3=1,\hfill \\ {\rho }_2=\frac{T{M}_S{M}_{I}a+\sigma \beta }{{M_S+\sigma M_I}^2},\hfill \\ {\rho }_1=-\frac{{\rho }_{11}+{\rho }_{12}+{\rho }_{13}+{\rho }_{14}+{\rho }_{15}}{L{M_S+\sigma M_I}^3},\hfill \\ {\rho }_0=T{M}_S{M}_I\frac{{\rho }_{21}+{\rho }_{22}-{\rho }_{23}-{\rho }_{24}-{\rho }_{25}+{\rho }_{26}}{L{M_S+\sigma M_I}^4},\hfill \end{array}\end{array}$$where ${\rho }_{11}+{\rho }_{13}+{\rho }_{15}>{\rho }_{12}+{\rho }_{14}$${\rho }_{11}<0,{\rho }_{15}<0,{\rho }_{13}<0$, and $$\begin{array}{}\text{(50)}& \begin{array}{c}{\rho }_{11}=M_S{M}_I\alpha {M_S+\sigma M_I}^3+Ta{M}_S{M}_{I}r-L\alpha ,\hfill \\ {\rho }_{12}=LT\beta M_S^2{M}_I^{2}a\sigma {\eta }_1-a{\eta }_2+\alpha {\sigma }^2,\hfill \\ {\rho }_{13}=LT3a^2\sigma {\eta }_2{M}_I^4-{\beta }^2{\eta }_1{M}_S^4,\hfill \\ {\rho }_{14}=Ta{M}_S{M}_I^{3}r\sigma +2La{\eta }_2-L\alpha \sigma ,\hfill \\ {\rho }_{15}=LT\sigma \beta M_S^3{M}_I\alpha -2\beta {\eta }_1,\hfill \\ {\rho }_{21}=a{M}_I{M}_S+\sigma M_I\beta {\eta }_1{M}_S^2+a{\eta }_2{M}_I^2+2\sigma \beta {\eta }_1{M}_S{M}_I\times r-L\alpha ,\hfill \\ {\rho }_{22}=La\sigma \beta M_S^2{M}_{I}T\beta {\eta }_1-5\alpha {\eta }_2{M}_I,\hfill \\ {\rho }_{23}=3La\sigma \beta {\eta }_2{M}_I^{3}Ta+\alpha \sigma M_S,\hfill \\ {\rho }_{24}=L\alpha \beta M_S^3{M}_{I}2a{\eta }_2+{\sigma }^2\beta {\eta }_1,\hfill \\ {\rho }_{25}=LT{a}^2\beta {\eta }_2{M}_S{M}_I^2,\hfill \\ {\rho }_{26}=L{\beta }^2{\eta }_1{M}_S^{3}Ta-\alpha \sigma M_S.\hfill \end{array}\end{array}$$

The Routh Hurwitz tabular is given by $\begin{array}{c}{\lambda }^3\\ {\lambda }^2\\ \lambda \\ 1\end{array}\begin{array}{cc}1\hfill & {\rho }_1\hfill \\ {\rho }_2\hfill & {\rho }_0\hfill \\ {\rho }_{11}\hfill & 0\hfill \\ {\rho }_{21}\hfill & 0\hfill \end{array}$ with ${\rho }_{11}={\rho }_1{\rho }_2-{\rho }_0/{\rho }_2^3,{\rho }_{21}={\rho }_0/{\rho }_{11}^2.$ If $p>q$ with $q,p>0$ and if $t_5-t_6>t_{7}L\alpha -\delta$ then ${\rho }_2>0.$ If ${\rho }_{11}+{\rho }_{13}+{\rho }_{15}>{\rho }_{12}+{\rho }_{14}$ then ${\rho }_1>0$. On the other hand, we have ${\rho }_1{\rho }_2>{\rho }_0.$ We conclude that ${\rho }_{11}>0$ and ${\rho }_{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,\alpha =0.9,\beta =0.2,\sigma =0.2,a=0.3,\delta =0.2,{\eta }_1=0.4,{\eta }_2=0.3$, and $d=0.2$ and initial value $\mathrm{0.1,0.1,0.1}$, the point $P_5$ tends to the point $0.221;0.784;1.65\times 10^{-33}.$ As shown in Figures 3–5. Then, we can deduce that the point $P_5$ 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 omitted; refer to PDF]

Following Figure 6, in a finite time the susceptible Scomber colias, infected Scomber colias and Thunnus thynnus converge to the equilibrium point $P_5$.

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 by $C_{ij}=q_j{E}_{ij}{x}_j,$ for all $i,j=1,2,3$; the catches $C_{ij}$ of species $j$ by fishing fleets $i$ are proportional to the fishing effort $E_{ij}$ exerted by fishing fleets $i$ to exploit species $j$ and the instantaneous biomass of the stock, where $q_j$ is the coefficient of catchability of species $j$, and $x_j$ is the biomass of population $j,$ more precisely, $x_1≔M_S,x_2≔M_I,x_3≔T.$ The total catches of the species $j$ by all fishing fleets are denoted by $C_j=C_{1j}+C_{2j}+C_{3j}$. We note that $E_j={\sum }_{i=1}^3{E}_{ij}$ represents the total fishing effort devoted to the species $j$ by all fishing fleets, and $E^i={E_{i1},E_{i2},E_{i3}}^T$ 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 $$\begin{array}{}\text{(51)}& \begin{array}{c}{\dot{M}}_S=r{M}_{S}1-\frac{M_S+M_I}{L}-\alpha M_S{M}_I-\frac{\beta M_S^{2}T}{M_S+\sigma M_I}-q_1{E}_1{M}_S,\hfill \\ {\dot{M}}_I=\alpha M_S{M}_I-\frac{a{M}_I^{2}T}{M_S+\sigma M_I}-\delta M_I-q_2{E}_2{M}_I,\hfill \\ \dot{T}=\frac{{\eta }_1\beta M_S^{2}T}{M_S+\sigma M_I}-\frac{{\eta }_{2}a{M}_I^{2}T}{M_S+\sigma M_I}-dT-q_3{E}_{3}T.\hfill \end{array}\end{array}$$

The biomass at biological equilibrium is the solution of the following system $$\begin{array}{c}\text{(52)}& \begin{array}{c}r1-\frac{M_S+M_I}{L}=\alpha M_I+\frac{\beta M_S^{2}T}{M_S+\sigma M_I}+q_1{E}_1,\hfill \\ \alpha M_S=\frac{a{M}_{I}T}{M_S+\sigma M_I}+\delta +q_2{E}_2,\hfill \\ \frac{{\eta }_1\beta M_S^2}{M_S+\sigma M_I}=\frac{{\eta }_{2}a{M}_I^2}{M_S+\sigma M_I}+d+q_3{E}_3,\hfill \end{array}\begin{array}{c}{M}_S=\frac{d+E_3{q}_3{a}_{11}-a_{12}{E}_1+a_{13}{E}_3+a_{14}{E}_2}{2a\beta {\eta }_{1}r+L\alpha +dr\sigma +r\sigma E_3{q}_3},\hfill \\ M_I=\frac{\delta {\eta }_1-a_{12}/\sigma E_1+a_{14}{E}_2+a_{17}{E}_3}{2ar\beta {\eta }_1+2La\alpha \beta {\eta }_1+2adr\sigma \beta {\eta }_1+2ar\sigma \beta {\eta }_1{E}_3{q}_3},\hfill \\ T=-\frac{-m_{36}\alpha -\delta -m_{35}\alpha -\delta E_3+a_{16}{E}_2+a_{15}{E}_2-m_{31}+E_1{m}_{32}-E_2{m}_{33}-E_3{m}_{34}}{2L\alpha +dr\sigma +r\sigma E_3{q}_3+r{a}^2{\eta }_1^2\beta -\delta +a_{12}/\sigma E_1-a_{14}/\sigma E_2-L\sigma {\beta }^2\delta q_3{E}_3},\hfill \end{array}\end{array}$$where $$\begin{array}{}\text{(53)}& \begin{array}{c}{a}_{11}=2ard\sigma +1+L\alpha +\sigma \delta {\eta }_1,\hfill \\ a_{12}=La\sigma \beta {\eta }_1{q}_1,\hfill \\ a_{13}=\sigma q_{3}2ar+L\sigma {\beta }^2\delta {\eta }_1,\hfill \\ a_{14}=L{\sigma }^2{\beta }^2{\eta }_1{q}_{2}d+E_3{q}_3,\hfill \\ a_{15}=2ar\sigma \beta {\eta }_1{q}_{2}d+E_3{q}_3,\hfill \\ a_{16}=2a\beta {\eta }_1{q}_{2}r+L\alpha ,\hfill \\ a_{17}=L\sigma {\beta }^2\delta {\eta }_1{q}_3,\hfill \\ m_{31}=d{a}_{11}+\sigma \delta {\eta }_1,\hfill \\ m_{32}=a_{12}d+1,\hfill \\ m_{33}=a_{14}d+{\eta }_1,\hfill \\ m_{34}=q_{3}2ar2d\sigma +1+L\alpha +L{\sigma }^2{\beta }^2\delta {\eta }_{1}d+1\hfill \\ +\sigma \delta {\eta }_1-q_3{a}_{12}{E}_1+q_3{a}_{14}{E}_2+q_3{a}_{13}{E}_3,\hfill \\ m_{35}=q_3{a}_{11}-q_3{a}_{12}{E}_1+q_3{a}_{14}{E}_2-q_3{a}_{13}{E}_3,\hfill \\ m_{36}=d{a}_{11}-d{a}_{12}{E}_1+d{a}_{14}{E}_2-d{a}_{13}{E}_3.\hfill \end{array}\end{array}$$

The profit $P_i$ of each fishing fleet $i$ is equal to the total revenue ${\text{TR}}_i$ minus the total cost ${\text{TC}}_i$. Mathematically, the profit for each fishing fleet is represented by $P_i={\text{TR}}_i-{\text{TC}}_i.$ Now, we give the expressions of total revenue and total costs of each fishing fleet. We use the fact that the total revenue $\text{TR}$ depends linearly on the catch, as follows: $\text{TR}=\text{price}\times \text{Catches}$. With the previous notations, we have $$\begin{array}{}\text{(54)}& {\text{TR}}_i=p_1{C}_{i1}+p_2{C}_{i2}+p_3{C}_{i3}=p_1{q}_1{E}_{i1}{M}_S+p_2{q}_2{E}_{i2}{M}_I+p_3{q}_3{E}_{i3}T=p_1{q}_1{E}_{i1}\frac{d+E_3{q}_3{a}_{11}-a_{12}{E}_1+a_{13}{E}_3+a_{14}{E}_2}{2a\beta {\eta }_{1}r+L\alpha +dr\sigma +r\sigma E_3{q}_3}+p_2{q}_2{E}_{i2}\frac{\delta {\eta }_1-a_{12}/\sigma E_1+a_{14}{E}_2+a_{17}{E}_3}{2a\beta {\eta }_{1}rd\sigma +1+L\alpha +r\sigma E_3{q}_3}+p_3{q}_3{E}_{i3}\frac{\sigma m_{36}\alpha -\delta +m_{35}\alpha -\delta E_3-a_{16}{E}_2-a_{15}{E}_2}{1}\times \frac{E_1{m}_{32}-m_{31}-E_2{m}_{33}-E_3{m}_{34}}{-2a^2\beta {\eta }_1^2\sigma \delta L\sigma {\beta }^2{E}_3{q}_3+1-E_1{a}_{12}+E_2{a}_{14}r\sigma d+E_3{q}_3+1+L\alpha },\end{array}$$where ${p_j}_{j=1,2,3}$ is the price per unit biomass of the susceptible Scomber colias, infected Scomber colias and Thunnus thynnus species, with $p_1>p_2.$ Now, we calculate the expression of the total effort cost, for that, we shall assume that ${\text{TC}}_i=c^i,E^i,$ where ${\text{TC}}_i$ is the total cost of the fishing fleets $i$. Therefore, the expression of the profit is represented by the following function $P_{i}E={\text{TR}}_i-$${\text{TC}}_i.$$$\begin{array}{c}\text{(55)}& P_{1}E=p_1{C}_{11}+p_2{C}_{12}+p_3{C}_{13}-c_{11}{E}_{11}-c_{12}{E}_{12}-c_{13}{E}_{13}\hfill =p_1{q}_1{E}_{11}\frac{d+E_3{q}_3{a}_{11}-a_{12}{E}_1+a_{13}{E}_3+a_{14}{E}_2}{2a\beta {\eta }_{1}r+L\alpha +dr\sigma +r\sigma E_3{q}_3}-c_{11}{E}_{11}\hfill \\ +p_2{q}_2{E}_{12}\frac{\delta {\eta }_1-a_{12}/\sigma E_1+a_{14}{E}_2+a_{17}{E}_3}{2ar\beta {\eta }_1+2La\alpha \beta {\eta }_1+2adr\sigma \beta {\eta }_1+2ar\sigma \beta {\eta }_1{E}_3{q}_3}-c_{12}{E}_{12}\hfill \\ +p_3{q}_3{E}_{13}\frac{\sigma m_{31}-E_1{m}_{32}+E_2{m}_{33}+E_3{m}_{34}{m}_{36}\alpha -\delta -E_2{a}_{15}+a_{16}+E_3{m}_{35}\alpha -\delta }{2a^2\beta {\eta }_1^2\sigma \delta L\sigma {\beta }^2{E}_3{q}_3+1-E_1{a}_{12}+E_2{a}_{14}rd\sigma +\sigma d+E_3{q}_3+L\alpha }\hfill \\ -c_{13}{E}_{13},\hfill \\ P_{2}E=p_1{C}_{21}+p_2{C}_{22}+p_3{C}_{23}-c_{21}{E}_{21}-c_{22}{E}_{22}-c_{23}{E}_{23}\hfill \\ =p_1{q}_1{E}_{21}\frac{d+E_3{q}_3{a}_{11}-a_{12}{E}_1+a_{13}{E}_3+a_{14}{E}_2}{2a\beta {\eta }_{1}r+L\alpha +dr\sigma +r\sigma E_3{q}_3}-c_{21}{E}_{21}\hfill \\ +p_2{q}_2{E}_{22}\frac{\delta {\eta }_1-a_{12}/\sigma E_1+a_{14}{E}_2+a_{17}{E}_3}{2ar\beta {\eta }_1+2La\alpha \beta {\eta }_1+2adr\sigma \beta {\eta }_1+2ar\sigma \beta {\eta }_1{E}_3{q}_3}-c_{22}{E}_{22}\hfill \\ +p_3{q}_3{E}_{23}\frac{\sigma m_{31}-E_1{m}_{32}+E_2{m}_{33}+E_3{m}_{34}{m}_{36}\alpha -\delta -E_2{a}_{15}+a_{16}+E_3{m}_{35}\alpha -\delta }{2a^2\beta {\eta }_1^2\sigma \delta L\sigma {\beta }^2{E}_3{q}_3+1-E_1{a}_{12}+E_2{a}_{14}rd\sigma +\sigma d+E_3{q}_3+L\alpha }\hfill \\ -c_{23}{E}_{23},\hfill \end{array}$$and $$\begin{array}{c}\text{(56)}& P_{3}E=p_1{C}_{31}+p_2{C}_{32}+p_3{C}_{33}-c_{31}{E}_{31}-c_{32}{E}_{32}-c_{33}{E}_{33}\hfill =p_1{q}_1{E}_{31}\frac{d+E_3{q}_3{a}_{11}-a_{12}{E}_1+a_{13}{E}_3+a_{14}{E}_2}{2a\beta {\eta }_{1}r+L\alpha +dr\sigma +r\sigma E_3{q}_3}-c_{31}{E}_{31}\hfill \\ +p_2{q}_2{E}_{32}\frac{\delta {\eta }_1-a_{12}/\sigma E_1+a_{14}{E}_2+a_{17}{E}_3}{2ar\beta {\eta }_1+2La\alpha \beta {\eta }_1+2adr\sigma \beta {\eta }_1+2ar\sigma \beta {\eta }_1{E}_3{q}_3}-c_{32}{E}_{32}\hfill \\ +p_3{q}_3{E}_{33}\frac{\sigma m_{31}-E_1{m}_{32}+E_2{m}_{33}+E_3{m}_{34}{m}_{36}\alpha -\delta -E_2{a}_{15}+a_{16}+E_3{m}_{35}\alpha -\delta }{2a^2\beta {\eta }_1^2\sigma \delta L\sigma {\beta }^2{E}_3{q}_3+1-E_1{a}_{12}+E_2{a}_{14}rd\sigma +\sigma d+E_3{q}_3+L\alpha }\\ -c_{33}{E}_{33}.\hfill \end{array}$$

5. Optimal Control

$$\begin{array}{}\text{(57)}& \begin{array}{c}{\dot{M}}_S=r{M}_{S}1-\frac{M_S+M_I}{L}-\alpha M_S{M}_I-\frac{\beta M_S^{2}T}{M_S+\sigma M_I}-C_1{M}_S,\hfill \\ {\dot{M}}_I=\alpha M_S{M}_I-\frac{a{M}_I^{2}T}{M_S+\sigma M_I}-\delta M_I-C_2{M}_I,\hfill \\ \dot{T}=\frac{{\eta }_1\beta M_S^{2}T}{M_S+\sigma M_I}-\frac{{\eta }_{2}a{M}_I^{2}T}{M_S+\sigma M_I}-dT-C_{3}T.\hfill \end{array}\end{array}$$

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 control $C_1$, $C_2$, and $C_3$. It is very important to show that all the control variables are nonnegative. The objective function $J{C}_1,C_2,C_3$ is defined as $$\begin{array}{}\text{(58)}& {\int }_{t_0}^{t_f}{e}^{-\gamma t}{p}_1{C}_1+C_2+C_3-v_1{C}_1^2+C_2^2+C_3^2-\frac{c_1{C}_1+C_2+C_3}{q_1{M}_S+q_2{M}_I+q_{3}T}dt.\end{array}$$

Subject to the system (4), where $p_1$ is the constant price per unit biomass of harvested population, $v_1$ is an economic constant, $\gamma$ is the instantaneous annual discount rate, and $C_1,C_2,C_3$ represent the harvesting rates on infected, susceptible Scomber colias and Thunnus thynnus, respectively. The coefficient $c_1$ is the constant fishing cost per unit effort. We are seeking to find an optimal control $C_1^{\ast }t,C_2^{\ast }t$, and $C_3^{\ast }t$ such that $$\begin{array}{}\text{(59)}& J{C}_1^{\ast },C_2^{\ast },C_3^{\ast }=\max J{C}_1,C_2,C_3\text{:}{C}_1,C_2,C_3\in U,\end{array}$$where $U=C_i:0\le C_{i}t\le 1,\text{Lebesguemeasurable}t=t_0,t_f\text{for}i=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 $C_1$, $C_2$, and $C_3$ to find the maximum value of the Lagrangian. According to Mwantobe et al. [25], this could be effectuated by considering Hamiltonian, $H$. $$\begin{array}{c}\text{(60)}& H=p_1{C}_1+C_2+C_3-v_1{C}_1^2+C_2^2+C_3^2-\frac{c_1{C}_1+C_2+C_3}{q_1{M}_S+q_2{M}_I+q_{3}T}\hfill +{\theta }_1{M}_{S}tr1-\frac{M_{S}t+M_{I}t}{L}-\alpha M_{I}t-\frac{\beta M_S^{2}tTt}{M_{S}t+\sigma M_{I}t}-C_1\hfill \\ +{\theta }_2\alpha M_{S}t{M}_{I}t-\frac{a{M}_I^{2}Tt}{M_{S}t+\sigma M_{I}t}-\delta M_{I}t-C_2{M}_{I}t\hfill \\ +{\theta }_3\frac{{\eta }_1\beta M_S^{2}tTt}{M_{S}t+\sigma M_{I}t}-\frac{{\eta }_{2}a{M}_I^{2}tTt}{M_{S}t+\sigma M_{I}t}-dTt-C_{3}Tt,\hfill \end{array}$$with ${\theta }_1$, ${\theta }_2$, and ${\theta }_3$ the adjoint variables associated to $M_S$, $M_I$, 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 $C_i^{\ast }\in U$for$i=1,2,3$ such that $$\begin{array}{}\text{(61)}& J{C}_1^{\ast },C_2^{\ast },C_3^{\ast }=\max J{C}_1,C_2,C_3\text{:}{C}_1,C_2,C_3\in U.\end{array}$$

Proof.

See the theorems in Lashari et al. [26], Lashari and Zaman [27], and Lenhart and Workman [28].

Theorem 2 (Necessary conditions).

Given optimal controls $C_1^{\ast },C_2^{\ast },C_3^{\ast }$ and solutions $M_S^{\ast },M_I^{\ast }$, and $T^{\ast }$, there exist ${\theta }_1,{\theta }_2,$ and ${\theta }_3$ the adjoint variables satisfying the following system $$\begin{array}{}\text{(62)}& \begin{array}{c}\frac{d{\theta }_1}{dt}=-\frac{\partial H}{\partial M_S}=-\frac{c_1{q}_1{C}_1+C_2+C_3}{{T{q}_3+q_1{M}_S+q_2{M}_I}^2}-{\theta }_1\frac{r2M_S-L+M_I}{L}+\alpha M_I+T\beta M_S\frac{M_S+2\sigma M_I}{{M_S+\sigma M_I}^2}+C_1-{\theta }_2\frac{M_I{\alpha {M_S+\sigma M_I}^2+TaM}_I}{{M_S+\sigma M_I}^2}-{\theta }_3\frac{T\beta {\eta }_1{M}_S{M}_S+2\sigma M_I}{{M_S+\sigma M_I}^2}+\frac{Ta{\eta }_2{M}_I^2}{{M_S+\sigma M_I}^2},\hfill \\ \frac{d{\theta }_2}{dt}=-\frac{\partial H}{\partial M_I}=-\frac{c_1{q}_1{C}_1+C_2+C_3}{{T{q}_3+q_1{M}_S+q_2{M}_I}^2}-{\theta }_1\frac{T\sigma \beta M_S^2}{{M_S+\sigma M_I}^2}-\frac{r{M}_S}{L}-\alpha M_S-{\theta }_2\alpha M_S-\frac{Ta{M}_{I}2M_S+\sigma M_I}{{M_S+\sigma M_I}^2}-\delta +C_2+{\theta }_3\frac{T\sigma \beta {\eta }_1{M}_S^2}{{M_S+\sigma M_I}^2}+\frac{Ta{\eta }_2{M}_{I}2M_S+\sigma M_I}{{M_S+\sigma M_I}^2},\hfill \\ \frac{d{\theta }_3}{dt}=-\frac{\partial H}{\partial T}=-\frac{c_1{q}_1{C}_1+C_2+C_3}{{T{q}_3+q_1{M}_S+q_2{M}_I}^2}+{\theta }_1\frac{\beta M_S^2}{M_S+\sigma M_I}+{\theta }_2\frac{a{M}_I^2}{M_S+\sigma M_I}-{\theta }_3\frac{\beta {\eta }_1{M}_S^2}{M_S+\sigma M_I}-\frac{a{\eta }_2{M}_I^2}{M_S+\sigma M_I}-d+C_3,\hfill \end{array}\end{array}$$with the transversality conditions ${\theta }_i{t}_f=0,i=1,2,3$. Where $$\begin{array}{}\text{(63)}& \begin{array}{c}{C}_1^{\ast }=\min \max 0,\frac{c_1-{\theta }_1{q}_1{M}_S^{2\ast }+T^{\ast }{p}_1{q}_3+p_1{q}_1{M}_S^{\ast }+p_1{q}_2{M}_I^{\ast }-{\theta }_1{q}_2{M}_S^{\ast }{M}_I^{\ast }-T{\theta }_1{q}_3{M}_S^{\ast }-1}{T^{\ast }{q}_3+q_1{M}_S^{\ast }+q_2{M}_I},1,\\ C_2^{\ast }=\min \max 0,\frac{c_1-{\theta }_2{q}_2{M}_I^{2\ast }+T^{\ast }{p}_1{q}_3+p_1{q}_1{M}_S^{\ast }+p_1{q}_2{M}_I^{\ast }-{\theta }_2{q}_1{M}_S^{\ast }{M}_I^{\ast }-T^{\ast }{\theta }_2{q}_3{M}_I^{\ast }-1}{T^{\ast }{q}_3+q_1{M}_S^{\ast }+q_2{M}_I^{\ast }},1,\\ C_3^{\ast }=\min \max 0,\frac{c_1+T^{\ast }{p}_1{q}_3+p_1{q}_1{M}_S^{\ast }+p_1{q}_2{M}_I^{\ast }-T^{2\ast }{\theta }_3{q}_3-T^{\ast }{\theta }_3{q}_1{M}_S-T{\theta }_3{q}_2{M}_I^{\ast }-1}{T^{\ast }{q}_3+q_1{M}_S^{\ast }+q_2{M}_I^{\ast }},1.\end{array}\end{array}$$

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 $\theta ={\theta }_1,\cdots ,{\theta }_n$ that satisfies the following equation $$\begin{array}{}\text{(64)}& \frac{dx}{dt}=\frac{\partial Ht,x,u,\theta }{\partial \theta };0=\frac{\partial Ht,x,u,\theta }{\partial u};\dot{\theta }=\frac{\partial Ht,x,u,\theta }{\partial x}.\end{array}$$

Proof.

Consider the system of differential equations in (65) governing the adjoint variables ${\theta }_1,{\theta }_2,{\theta }_3$ This is obtained by differentiating the Hamiltonian $H$ in equation (60) with respect to $M_S,M_I$, 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) $$\begin{array}{}\text{(65)}& \begin{array}{c}\frac{d{\theta }_1}{dt}=-\frac{\partial H}{\partial M_S}=-\frac{c_1{q}_1{C}_1+C_2+C_3}{{T{q}_3+q_1{M}_S+q_2{M}_I}^2}-{\theta }_1\frac{r2M_S-L+M_I}{L}+\alpha M_I+T\beta M_S\frac{M_S+2\sigma M_I}{{M_S+\sigma M_I}^2}+C_1-{\theta }_2\frac{M_I{\alpha {M_S+\sigma M_I}^2+TaM}_I}{{M_S+\sigma M_I}^2}-{\theta }_3\frac{T\beta {\eta }_1{M}_S{M}_S+2\sigma M_I}{{M_S+\sigma M_I}^2}+\frac{Ta{\eta }_2{M}_I^2}{{M_S+\sigma M_I}^2},\hfill \\ \frac{d{\theta }_2}{dt}=-\frac{\partial H}{\partial M_I}=-\frac{c_1{q}_1{C}_1+C_2+C_3}{{T{q}_3+q_1{M}_S+q_2{M}_I}^2}-{\theta }_1\frac{T\sigma \beta M_S^2}{{M_S+\sigma M_I}^2}-\frac{r{M}_S}{L}-\alpha M_S-{\theta }_2\alpha M_S-\frac{Ta{M}_{I}2M_S+\sigma M_I}{{M_S+\sigma M_I}^2}-\delta +C_2+{\theta }_3\frac{T\sigma \beta {\eta }_1{M}_S^2}{{M_S+\sigma M_I}^2}+\frac{Ta{\eta }_2{M}_{I}2M_S+\sigma M_I}{{M_S+\sigma M_I}^2},\hfill \\ \frac{d{\theta }_3}{dt}=-\frac{\partial H}{\partial T}=-\frac{c_1{q}_1{C}_1+C_2+C_3}{{T{q}_3+q_1{M}_S+q_2{M}_I}^2}+{\theta }_1\frac{\beta M_S^2}{M_S+\sigma M_I}+{\theta }_2\frac{a{M}_I^2}{M_S+\sigma M_I}-{\theta }_3\frac{\beta {\eta }_1{M}_S^2}{M_S+\sigma M_I}-\frac{a{\eta }_2{M}_I^2}{M_S+\sigma M_I}-d+C_3.\hfill \end{array}\end{array}$$

With the transversality conditions $$\begin{array}{}\text{(66)}& {\theta }_1{t}_f={\theta }_2{t}_f={\theta }_3{t}_f=0.\end{array}$$

To evaluate the optimal control of the control variable set, where $C_i=0,1$. Let $M_S=M_S^{\ast },M_I=M_I^{\ast }$, and $T=T^{\ast },$ and applying the second equation in equation (64), and differentiating the Hamiltonian, $H$, in equation (60) with respect to the control variables $C_1,C_2$, and $C_3$$$\begin{array}{}\text{(67)}& \begin{array}{c}\frac{\partial H}{\partial C_1}=-\frac{c_1+{\theta }_1{q}_1{M}_S^2-p_1{q}_1{M}_S+T{\theta }_1{q}_3{M}_S-T{p}_1{q}_3-p_1{q}_2{M}_I+{\theta }_1{q}_2{M}_S{M}_I+2C_1{v}_{1}T{q}_3+q_1{M}_S+q_2{M}_I}{T{q}_3+q_1{M}_S+q_2{M}_I}=0,\hfill \\ \frac{\partial H}{\partial C_2}=-\frac{c_1+{\theta }_2{q}_2{M}_I^2-T{p}_1{q}_3-p_1{q}_1{M}_S-p_1{q}_2{M}_I+{\theta }_2{q}_1{M}_S{M}_I+2C_2{v}_{1}T{q}_3+q_1{M}_S+q_2{M}_I+T{\theta }_2{q}_3{M}_I}{T{q}_3+q_1{M}_S+q_2{M}_I}=0,\hfill \\ \frac{\partial H}{\partial C_3}=-\frac{c_1-T{p}_1{q}_3-p_1{q}_1{M}_S-p_1{q}_2{M}_I+T^2{\theta }_3{q}_3+2C_3{v}_{1}T{q}_3+q_1{M}_S+q_2{M}_I+T{\theta }_3{q}_1{M}_S+T{\theta }_3{q}_2{M}_I}{T{q}_3+q_1{M}_S+q_2{M}_I}=0.\hfill \end{array}\end{array}$$

By applying the optimal control to the control variable set, $C_i^{\ast }=0,1$ for $i=1,2,3$ into equation (11) $$\begin{array}{}\text{(68)}& \begin{array}{c}{C}_1^{\ast }=\min \max 0,\frac{c_1-{\theta }_1{q}_1{M}_S^{2\ast }+T^{\ast }{p}_1{q}_3+p_1{q}_1{M}_S^{\ast }+p_1{q}_2{M}_I^{\ast }-{\theta }_1{q}_2{M}_S^{\ast }{M}_I^{\ast }-T{\theta }_1{q}_3{M}_S^{\ast }-1}{T^{\ast }{q}_3+q_1{M}_S^{\ast }+q_2{M}_I},1,\\ C_2^{\ast }=\min \max 0,\frac{c_1-{\theta }_2{q}_2{M}_I^{2\ast }+T^{\ast }{p}_1{q}_3+p_1{q}_1{M}_S^{\ast }+p_1{q}_2{M}_I^{\ast }-{\theta }_2{q}_1{M}_S^{\ast }{M}_I^{\ast }-T^{\ast }{\theta }_2{q}_3{M}_I^{\ast }-1}{T^{\ast }{q}_3+q_1{M}_S^{\ast }+q_2{M}_I^{\ast }},1,\\ C_3^{\ast }=\min \max 0,\frac{c_1+T^{\ast }{p}_1{q}_3+p_1{q}_1{M}_S^{\ast }+p_1{q}_2{M}_I^{\ast }-T^{2\ast }{\theta }_3{q}_3-T^{\ast }{\theta }_3{q}_1{M}_S-T{\theta }_3{q}_2{M}_I^{\ast }-1}{T^{\ast }{q}_3+q_1{M}_S^{\ast }+q_2{M}_I^{\ast }},1.\end{array}\end{array}$$

This shows that the uniqueness of the optimal control of the model has been achieved for small $t_f$ 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$, $\alpha =0.9$, $\beta =0.2$, $\sigma =0.2$, $a=0.3$, $\delta =0.2$, ${\eta }_1=0.4$, ${\eta }_2=0.3$, $d=0.2$, $C_1=0.2$, $c_1=0.15$, $C_3=0.1$, and $C_2=\text{variable}$.

Following Figures 7–9, we can show that when the constant effort harvesting $C_2$ on the susceptible Scomber colias is very low, i.e., $C_2=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, $C_2=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 omitted; refer to PDF]

Now, for $C_2=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 omitted; refer to PDF][figure omitted; refer to PDF][figure omitted; refer to PDF]

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.

References

[1] S. Dahani, N. Bouchriti, I. Benabbes, A. Boudakkou, A. Chiaar, "Occurrence Des parasites dans les poissons collectés au niveau du littoral Marocain," European Scientific Journal ESJ, vol. 15 no. 36,DOI: 10.19044/esj.2019.v15n36p497, 2019.

[2] M. Bentounsi, I. Agmour, N. Achtaich, Y. El Foutayeni, "The impact of price on the profits of fishermen exploiting tritrophic prey- predator fish populations," International Journal of Differential Equations, vol. 2018,DOI: 10.1155/2018/2381483, 2018.

[3] A. Laurec, J. C. Leguen, "Dynamique des populations marines exploitées," Concepts et modeles, vol. 45, 1981.

[4] M. Bentounsi, I. Agmour, N. Achtaich, Y. El Foutayeni, "Stability analysis of a biological model of a marine resources allowing density dependent migration," International Frontier Science Letters, vol. 12, pp. 22-34, DOI: 10.18052/www.scipress.com/IFSL.12.22, 2017.

[5] H. B. Shi, S. Ruan, "Spatial, temporal and spatiotemporal patterns of diffusive predator-prey models with mutual interference," IMA Journal of Applied Mathematics, vol. 80 no. 5, pp. 1534-1568, DOI: 10.1093/imamat/hxv006, 2015.

[6] Y. Sekerci, S. Petrovskii, "Global warming can lead to depletion of oxygen by disrupting phytoplankton photosynthesis: a mathematical modelling approach," Geosciences, vol. 8 no. 6,DOI: 10.3390/geosciences8060201, 2018.

[7] Y. Sekerci, S. Petrovskii, "Pattern formation in a model oxygen-plankton system," Computation, vol. 6 no. 4,DOI: 10.3390/computation6040059, 2018.

[8] N. Baba, I. Agmour, N. Achtaich, Y. El Foutayeni, "The mathematical study for mortality coefficients of small pelagic species," Communications in Mathematical Biology and Neuroscience, vol. 2019,DOI: 10.28919/cmbn/4117, 2019.

[9] I. Agmour, M. Bentounsi, N. Achtaich, Y. El Foutayeni, "Carrying capacity influence on the incomes of seiners exploiting marine species in the Atlantic coast of Morocco," Mathematical Biosciences, vol. 305, pp. 10-17, DOI: 10.1016/j.mbs.2018.08.012, 2018.

[10] D. Greenhalgh, Q. J. A. Khan, F. A. Al-Kharousi, "Eco-epidemiological model with fatal disease in the prey," Nonlinear Analysis: Real World Applications, vol. 53,DOI: 10.1016/j.nonrwa.2019.103072, 2020.

[11] D. Riad, K. Hattaf, N. Yousfi, "Dynamics of capital-labour model with Hattaf-Yousfi functional response," British Journal of Mathematics & Computer Science, vol. 18 no. 5,DOI: 10.9734/BJMCS/2016/28640, 2016.

[12] L. S. Pontryagin, V. G. Boltyansk, K. V. Gamkrelidre, E. F. Mishchenko, The Mathematical Theory of Optimal Processes, 1964.

[13] I. Agmour, M. Bentounsi, N. Baba, Y. El Foutayeni, N. Achtaich, C. Aouiti, "Impact of wind speed on fishing effort," Modeling Earth Systems and Environment, vol. 6 no. 2, pp. 1007-1015, DOI: 10.1007/s40808-020-00736-7, 2020.

[14] Y. El Foutayeni, M. Bentounsi, I. Agmour, N. Achtaich, "Bioeconomic model of zooplankton-phytoplankton in the central area of Morocco," Modeling Earth Systems and Environment, vol. 6 no. 1, pp. 461-469, DOI: 10.1007/s40808-019-00693-w, 2020.

[15] I. Agmour, M. Bentounsi, N. Achtaich, Y. El Foutayeni, "Catchability coefficient influence on the fishermen’s net economic revenues," Communications in Mathematical Biology and Neuroscience, vol. 2018,DOI: 10.28919/cmbn/3512, 2018.

[16] I. Agmour, N. Achtaich, Y. El Foutayeni, "Stability analysis of a competing fish populations model with the presence of a predator," International Journal of Nonlinear Science, vol. 26, pp. 108-121, 2018.

[17] M. Bentounsi, I. Agmour, N. Achtaich, Y. El Foutayeni, "The Hopf bifurcation and stability of delayed predator-prey system," Computational and Applied Mathematics, vol. 37 no. 5, pp. 5702-5714, DOI: 10.1007/s40314-018-0658-7, 2018.

[18] I. Agmour, M. Bentounsi, N. Achtaich, Y. El Foutayeni, "Optimization of the two fishermen’s profits exploiting three competing species where prices depend on harvest," International Journal of Differential Equations, vol. 2017,DOI: 10.1155/2017/3157294, 2017.

[19] M. Bentounsi, I. Agmour, N. Achtaich, Y. El Foutayeni, "Intrinsic growth rates influence on the net economic rents of fishermen," International Journal of Dynamical Systems and Differential Equations, vol. 9 no. 4, pp. 362-379, DOI: 10.1504/IJDSDE.2019.103735, 2019.

[20] H. El Bouanani, Y. El Foutayeni, M. Khaladi, "A new method for solving non-linear complementarity problems," International Journal of Nonlinear Science, vol. 19, pp. 81-90, 2015.

[21] Y. El Foutayeni, H. El Bouanani, M. Khaladi, "The linear complementarity problem and a method to find all its solutions," Information in Sciences and Computing, vol. 3, 2015.

[22] Y. El Foutayeni, M. Khaladi, "Equilibrium points and their stability properties of a multiple delays model," Differential Equations and Dynamical Systems, vol. 28 no. 2, pp. 255-272, DOI: 10.1007/s12591-016-0321-y, 2020.

[23] Y. El Foutayeni, H. El Bouanani, M. Khaladi, "An (m+1)-step iterative method of convergence order (m+2) for linear complementarity problems," Journal of Applied Mathematics and Computing, vol. 54 no. 1-2, pp. 229-242, DOI: 10.1007/s12190-016-1006-y, 2017.

[24] L. Pontryagin, V. G. Boltyanskii, E. Mishchenko, The Mathematical Theory of Optimal Processes, 1962.

[25] P. M. Mwamtobe, S. Abelman, J. M. Tchuenche, A. Kasambara, "Optimal (control of) intervention strategies for malaria epidemic in Karonga District, Malawi," Abstract and Applied Analysis, vol. 2014,DOI: 10.1155/2014/594256, 2014.

[26] A. A. Lashari, K. Hattaf, G. Zaman, X. Z. Li, "Backward bifurcation and optimal control of a vector-borne disease," Applied Mathematics & Information Sciences, vol. 7 no. 1, pp. 301-309, DOI: 10.12785/amis/070138, 2013.

[27] A. A. Lashari, G. Zaman, "Optimal control of a vector borne disease with horizontal transmission," Nonlinear Analysis: Real World Applications, vol. 13 no. 1, pp. 203-212, DOI: 10.1016/j.nonrwa.2011.07.026, 2012.

[28] S. Lenhart, J. T. Workman, Optimal Control Applied to Biological Models, 2007.

[29] W. H. Fleming, R. W. Rishel, Deterministic and Stochastic Optimal Control, 1975.