A Complex Dynamic of an Eco-Epidemiological Mathematical Model with Migration

Abstract

In this paper, we propose an eco-epidemiological mathematical model in order to describe the effect of migration on the dynamics of a prey–predator population. The functional response of the predator is governed by the Holling type II function. First, from the perspective of mathematical results, we develop results concerning the existence, uniqueness, positivity, boundedness, and dissipativity of solutions. Besides, many thresholds have been computed and used to investigate the local and global stability results by using the Routh–Hurwitz criterion and Lyapunov principle, respectively. We have also established the appearance of limit cycles resulting from the Hopf bifurcation. Numerical simulations are performed to explore the effect of migration on the dynamic of prey and predator populations.

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

Generally speaking, spatial or geographic heterogeneity plays an important role in the transmission process of many infectious diseases when species interact. This propagation of infectious diseases has been known to be an important regulating factor for human and animal population sizes. The study of the effect of infectious disease propagation in species in interaction has attracted some attention in ecology and epidemiology due to its seriousness threats all over the world. In most of the mathematical models, environment has been considered as homogeneous. However, in reality, environment is heterogeneous, and it can be considered as a set of different localities connected by migration [1, 2]. In particular, for prey–predator population, infectious diseases coupled with prey–predator model produce a complex dynamic given the multitude of species living in the environment and interacting with each other. This complexity is increased in the presence of migration in each species. We cannot ignore this factor because it is common in the ecological system and plays a major role in the natural regulation of populations.

After the pioneer work of Kermack and McKendricK [3], epidemiological models have drawn consideration of many biologists and ecologists. The infectious diseases in the community of prey–predator models has been studied extensively in literature [1, 4–14]. In [14], Tewa et al. proposed and investigated a prey–predator model with an SIS infectious disease affecting preys or predators or both. They showed that the disease can disappear from the community, persist in one or two populations of the community. Savadogo et al. [8] proposed and analyzed an eco-epidemiological model describing the effect of predation in the dynamic of propagation of disease. Results concerning the local and global stability have been analyzed according to Routh–Hurwitz criterion and Lyapunov principle, respectively. They also established the Hopf-bifurcation to highlight the periodic fluctuation with extinction or persistence of the disease in the preys and predators communities. In [5], Biswas et al. proposed and investigated a cannibalistic eco-epidemiological model with disease in the predator species. They conclude that cannibalism is process of regulation and govern the disease transmission in the predators community. Many works have been investigated in literature on ecological models with disease in prey only [15–18]. For instance, Greenhalgh et al. [17], studied an eco-epidemiological model with fatal disease in the prey population. In [16], Chattopadhyay and Arino studied a prey–predator model with disease in the prey population. They showed persistence and extinction conditions of the populations and also determined conditions for which the system enters a Hopf-bifurcation.

Furthermore, it is well-known that migration can be described as an important demographic process that occurs in all living beings. It can be occurred for many reasons such as education, employment, marriages, war for human beings. For animal species, migration is usually due to the climate change, for habitat, looking for food, predation, cannibalism, reproduction, urbanization, deforestation, etc. Many researches are mainly focused on the effect of infectious diseases in the community of prey–predator models under the influence of migration [4, 11, 13, 19–21]. Indeed, Kant and Kumar [11] proposed and studied a prey–predator model with disease in both species by taking into account the process of migration in only prey population. In their mathematical analysis, existence, positivity and stability of equilibria has been investigated. Moreover, the epidemiological thresholds are computed and used to determine the conditions of the disease propagation. In [4], Arora and Vivek proposed and investigated a prey–predator model when disease spread among prey populations with migration in both species. Holling type II functional response is used for interaction between prey and predator species. Their mathematical analysis has permitted to establish existence of equilibria and the local and global stability. Chowdhury [13] investigated an eco-epidemiological model describing the effect of fast migration on prey–predator model between two different patches. Mathematically, the author proved the asymptotic stability of the unique fast equilibrium point and the aggregated model.

It is in this line of thought that, in this paper, we propose and study an eco-epidemiological prey–predator model to study the effect of migration on the dynamic of prey and predator. Indeed, motivated by the works of Arora and Kumar [4], Kant and Kumar [11], and Tewa et al. [14], our main goal in this work, is to analyze the effect of migration on the dynamic of prey–predator model in the presence of an SIS infectious diseases. We found inspiration in the work of an eco-epidemiological model studied by Tewa et al. [14], by taking into account the migration. Besides, we established the conditions of existence when our model admits at least coexistence equilibrium. Based upon thorough of mathematical analysis of our model under consideration, all the equilibria point of the system are adequately characterized, and their stability analysis are investigated following the Routh–Hurwitz criteria and Lyapunov principle. Also, we established the existence of the limit cycles of the system studied arising from Hopf-bifurcation. Finally, in the numerical simulation, bifurcation diagrams and phase portraits are given, and some complex and rich dynamic behaviors, such as limit cycle, periodic solutions are found. The results show that variation of migration parameters may affect prey and predator density, thereby controlling species density can effectively maintain the ecological balance.

The remaining part of this paper is structured as follows after the statement of the problem. Section 2 is devoted to the formulation of the eco-epidemiological model. In Section 3, a mathematical analysis of the model is established, including well-posedness, stabilities analysis, and Hopf-bifurcation. We perform some numerical simulations to support our main results in Section 4. The paper ends with a conclusion and discussion in Section 5.

2. Mathematical Formulation of the Eco-Epidemiological Model

In this section, our goal is to establish an eco-epidemiological model in order to study the effect of migration on the dynamic of preys and predators population. Let’s denote by $X$ and $Y$ as the susceptible and infectious prey density, respectively, such that $Ht=Xt+Yt$ is the total prey population at any time $t>0.$ It should be pointed out that the density dependence affects the birth and the death of the populations. Therefore, we need to separate the effects of the density dependence. Denoting by $b$ and $\mu$ the natural birth and death rates parameters, respectively, the parameter $\theta$ is such that $b-r\theta H/K$ is the birth rate of susceptible prey and $\mu +1-\theta rH/K$ is the mortality rate [22].

We list the following key assumptions useful in the mathematical formulation of our prey–predator system.

(H1): When there is no predator, the prey population growth logistically $r1-H/KH,$ where, $H$ represents the total population of preys, $K$ denotes the carrying capacity, and $r=b-\mu$ is the intrinsic growth rate;

(H2): The Holling function response of type II is used to represent the process of predation and is defined by $k_{1}H/1+k_{2}H,$ where $k_1$ and $k_2$ denote, respectively, predator search and satiety rates;

(H3): The disease is transmitted by contact between an infected and susceptible prey by standard incidence $\beta XY/H$;

(H4): The disease is not genetically inherited. The infected population can recover or become susceptible to diseases;

(H5): Only susceptible prey is capable of reproducing and contributing to their carrying capacity.

According to the above assumptions and the interaction diagram of Figure 1, the dynamics of the global eco-epidemiological model is given by the following set of differential equations:$$\begin{array}{}\text{(1)}& \begin{array}{cc}\dot{H}=rH1-\frac{H}{K}-k_1\frac{X+qY}{1+k_{2}H}P-m_{1}X+m_{2}Y,\hfill & H0>0,\\ \dot{X}=b-\frac{r\theta H}{K}H-\mu +\frac{1-\theta rH}{K}X-\beta \frac{XY}{H}+\lambda Y-\frac{k_{1}XP}{1+k_{2}H}-m_{1}X,& X0>0,\\ \dot{Y}=\beta \frac{XY}{H}-\lambda Y-\mu +\frac{1-\theta rH}{K}Y-\frac{k_{1}qYP}{1+k_{2}H}-m_{2}Y,\hfill & Y0\ge 0,\\ \dot{P}=\frac{\omega X+qYP}{1+k_{2}H}-\gamma P-m_{3}P,\hfill & P0>0,\end{array}\end{array}$$where

(1) $\gamma$ denotes the natural mortality rate of predators;

(2) $\beta$ is the contact rate between susceptible and infectious preys;

(3) $\lambda$ denotes the recovery rate for infectious preys;

(4) $P$ denotes the predator population;

(5) $e=\omega /k_1$ represents the conversion rate of prey biomass into predatory biomass, with $0<e<1$;

(6) $q>1$ represents the fact that the infected preys are more easy to be caught than the healthy preys;

(7) $m_1$, $m_2$, and $m_3$ denote the migration rates of susceptible prey, infected prey, and predators, respectively, where $m_1>m_2$.

By setting $I=Y/H$ and $S=X/H=1-I$ the proportions of infected and susceptible preys, respectively, in the prey population [1, 14], for any time $t>0,$ system (1) becomes:$$\begin{array}{}\text{(2)}& \begin{array}{cc}\dot{H}=r1-\frac{H}{K}-\frac{k_{1}P1+Iq-1}{1+k_{2}H}-m_1+m_2-m_{1}IH={\pi }_{1}H,I,P,& H0=H_0>0,\\ \dot{I}=\beta 1-I-b+\lambda -\frac{\theta rH}{K}-\frac{k_{1}1-Iq-1P}{1+k_{2}H}-m_{2}I={\pi }_{2}H,I,P,\hfill & I0=I_0>0,\\ \dot{P}=\frac{\omega 1+Iq-1H}{1+k_{2}H}-\gamma -m_{3}P={\pi }_{3}H,I,P,\hfill & P0=P_0>0.\end{array}\end{array}$$

[figure(s) omitted; refer to PDF]

3. Mathematical Investigation of Model

This section deals with mathematical analysis of system (2) [1].

3.1. Existence, Positivity, and Boundedness Properties

For system (2) to be ecologically and epidemiologically meaningful, it is important to prove that all its state variables are nonnegative for all time. Then, we rewrite model (2) in the following form:$$\begin{array}{}\text{(3)}& \dot{Z}t=GZt,\text{where}Zt={Ht,It,Pt}^T,\end{array}$$and $G:{\mathbb{R}}_{+}^3\to {\mathbb{R}}^3$ is defined by:$$\begin{array}{}\text{(4)}& GZ=\begin{array}{c}{G}_{1}H,I,P\\ \\ G_{2}H,I,P\\ \\ G_{3}H,I,P\end{array}=\begin{array}{c}r1-\frac{H}{K}-\frac{k_{1}P1+Iq-1}{1+k_{2}H}-m_1+m_2-m_{1}IH\\ \\ \beta 1-I-b+\lambda -\frac{\theta rH}{K}-\frac{k_{1}1-Iq-1P}{1+k_{2}H}-m_{2}I\\ \\ \frac{\omega 1+Iq-1H}{1+k_{2}H}-\gamma -m_{3}P\end{array}.\end{array}$$

The following results hold for model (2) [1, 8].

Theorem 1.

The nonnegative orthant ${\mathbb{R}}_{+}^3$ is positively invariant by system (2).

Theorem 2.

System (2) admits a unique global solution $Ht,It,Pt$ defined on interval $0,T_{\text{max}}$. Moreover, the set $H\ge 0,I\ge 0,P\ge 0$ is positively invariant for system (2).

Proof.

Indeed,

(1) The theorem of Cauchy–Lipschitz assures the existence and uniqueness of local solution of system (2) on $0,T_{\text{max}}$ given the regularity of the functions involved in the model.

(2) Now let us show that $H\ge 0,I\ge 0,P\ge 0$ is positively invariant under the flow of system (2). Thus, one has$$\begin{array}{}\text{(5)}& \begin{array}{cc}{\pi }_{1}K,0,P=-\frac{k_{1}P}{1+k_{2}K}-m_{1}K\le 0,\hfill & \text{for}P\ge 0,\\ {\pi }_{2}0,1,P=-b+\lambda +m_2\le 0,\hfill & \text{for}P\ge 0,\\ {\pi }_{3}H,I,0=0,\hfill & \text{for}H\ge 0,I\ge 0.\end{array}\end{array}$$Consequently, the set $H\ge 0,I\ge 0,P\ge 0$ is positively invariant for system (2).

The boundedness of system (2) is given by the following theorem.

Theorem 3.

The closed set defined by:$$\begin{array}{}\text{(6)}& D_1=H,I,P\in {\mathbb{R}}_{+}^3/0\le H\le K,0\le I\le 1,0\le eHt+It+Pt\le \frac{erK+4\zeta -m_1}{4\zeta },\end{array}$$is a compact forward invariant set for system (2). Moreover, this set is absorbing for $r>0$.

Proof.

Indeed, adding the three equations of system (2), we obtain:$$\begin{array}{}\text{(7)}& Wt=eHt+It+Pt.\end{array}$$

The time derivatives $Wt$ along the trajectories of system (2) is given by:$$\begin{array}{}\text{(8)}& \begin{array}{ccc}\dot{W}& =& er1-\frac{H}{K}-\frac{k_{1}P1+Iq-1}{1+k_{2}H}-m_1+m_2-m_{1}IH\\ & +& \beta 1-I-b+\lambda -\frac{\theta rH}{K}-\frac{k_{1}1-Iq-1P}{1+k_{2}H}-m_{2}I\\ & +& \frac{\omega 1+Iq-1H}{1+k_{2}H}-\gamma -m_{3}P.\end{array}\end{array}$$

From Equation (8), we get:$$\begin{array}{}\text{(9)}& \dot{W}+\zeta W\le er1-\frac{H}{K}H+e\zeta -e{m}_1,\end{array}$$where $\zeta =\text{min}{m}_3+\gamma ;m_2+b_1+\lambda -\beta >0.$ Let us consider:$$\begin{array}{}\text{(10)}& fH=er1-\frac{H}{K}H+e\zeta -e{m}_1,\end{array}$$then the maximum value of $fH$ at $H=K/2$ is $erK+4\zeta -m_1/4\zeta .$ Thus, from Equation (10) we get the following differential inequality:$$\begin{array}{}\text{(11)}& \dot{W}+\zeta W\le \frac{erK+4\zeta -m_1}{4\zeta }.\end{array}$$

By integrating the differential inequality Equation (7) and by using the theory of Birkoff and Rota [23] yields:$$\begin{array}{}\text{(12)}& Wt\le \frac{erK+4\zeta -m_1}{4\zeta }1-e^{-\zeta t}+W0e^{-\zeta t}.\end{array}$$

Then,$$\begin{array}{}\text{(13)}& \underset{t⟶+\mathrm{\infty }}{\text{lim}}Wt\le \frac{erK+4\zeta -m_1}{4\zeta }.\end{array}$$

Thus, the system (2) are bounded. Particularly, $Wt\le erK+4\zeta -m_1/4\zeta$ if $W0\le erK+4\zeta -m_1/4\zeta .$ As result, $D_1$ is positively invariant. If $Wt>erK+4\zeta -m_1/4\zeta$ either the solution enters $D_1$ in finite time or $Wt$ approaches $erK+4\zeta -m_1/4\zeta$ and $It$ approach zero. Therefore, $D_1$ is an attractive set. Thus, system (2) is ecologically and epidemiologically well-posed.

3.2. Stability Analysis of Trivial Equilibria

In this subsection, we discuss the existence and the local stability of each trivial equilibrium point. Let’s consider the following epidemiological and ecological thresholds, with clear and distinct biological meaning [14, 24]:

(i) ${\mathcal{R}}_0=\beta /b+\lambda +m_2$ is the threshold that determines when disease die out or persist;

(ii) ${\mathcal{R}}_1=\beta /b+\lambda -\theta r+m_2$ is the threshold that determines the local stability of $E_0=K,0,0$ and $E_2=H^0,0,P^0$ when there is no disease;

(iii) ${\mathcal{R}}_{ec}=\omega K/\gamma 1+k_{2}K$ represents the ecological threshold of a prey–predator model without disease and migration.

The trivial equilibria is obtained by setting the right-hand sides of system (2) to zero. The explicit expressions of the trivial equilibria are given by the following proposition:

Proposition 1.

(i) $E_0=0,0,0$ is a trivial equilibrium point of system (2) and is always admissible.

(ii) $E_1=K,0,0$ is a prey free equilibrium point of system (2) and exists unconditionally.

(iii) $E_2=H^0,0,P^0$ is a disease free equilibrium, where$$\begin{array}{}\text{(14)}& H^0=\frac{\gamma +m_3}{\omega -k_2{m}_3+\gamma },P^0=\frac{er1+k_{2}K\gamma {\mathcal{R}}_{ec}-1-m_3-m_{1}K\omega -k_2{m}_3+\gamma }{K{\omega -k_2{m}_3+\gamma }^2}.\end{array}$$

This equilibrium is ecologically admissible if:$$\begin{array}{}\text{(15)}& \omega >k_2\gamma +m_3,{\mathcal{R}}_{ec}>1,m_3<\gamma {\mathcal{R}}_{ec}-1\text{and}{m}_1<\frac{r\gamma 1+k_{2}K\gamma {\mathcal{R}}_{ec}-1-m_3}{K\omega -k_2\gamma +m_3}.\end{array}$$

Now, we are in position to investigate the local stability of trivial equilibria.

Proposition 2.

(i) $E_0$ is always unstable,

(ii) $E_1=K,0,0$ is locally asymptotically stable if and only if:$$\begin{array}{}\text{(16)}& \begin{array}{c}{\mathcal{R}}_1<1,\\ \frac{r\gamma +m_3{1+k_{2}K}^2}{r{1+k_{2}K}^2-k_1}<{\mathcal{R}}_{ec}<1+\frac{r+m_3}{\gamma }.\end{array}\end{array}$$

(iii) $E_2=H^0,0,P^0$ is locally asymptotically stable if and only if:$$\begin{array}{}\text{(17)}& {\mathcal{R}}_{ec}>1+\frac{\omega m_{1}K+m_3}{\gamma }+\frac{\gamma +m_{3}K\theta r-m_1{k}_{2}q-1}{\gamma q-11+k_{2}K},\end{array}$$$$\begin{array}{}\text{(18)}& m_1>\frac{r\gamma +m_3}{K1-k_2\gamma +m_3}.\end{array}$$

Proof.

Let us consider the variational matrix of system (2):$$\begin{array}{}\text{(19)}& JH,I,P=\begin{array}{ccc}{B}_{11}& B_{12}& B_{13}\\ B_{21}& B_{22}& B_{23}\\ B_{31}& B_{32}& B_{33}\end{array},\end{array}$$where$$\begin{array}{}\text{(20)}& B_{11}=\frac{\partial F_1}{\partial H}=r1-\frac{2H}{K}-\frac{k_{1}P1+Iq-1}{{1+k_{2}H}^2}-m_1+m_2-m_{1}I,\end{array}$$$$\begin{array}{}\text{(21)}& B_{12}=\frac{\partial F_1}{\partial I}=-\frac{k_{1}q-1HP}{1+k_{2}H}-m_2-m_{1}H,\end{array}$$$$\begin{array}{}\text{(22)}& B_{13}=\frac{\partial F_1}{\partial P}=-\frac{k_{1}H1+Iq-1}{1+k_{2}H},\end{array}$$$$\begin{array}{}\text{(23)}& B_{21}=\frac{\partial F_2}{\partial H}=\frac{\theta rI}{K}+\frac{k_1{k}_{2}q-1I1-IP}{{1+k_{2}H}^2},\end{array}$$$$\begin{array}{}\text{(24)}& B_{22}=\frac{\partial F_2}{\partial I}=\beta 1-2I-b+\lambda -\frac{\theta rH}{K}-\frac{k_{1}q-11-2IP}{1+k_{2}H}-m_2,\end{array}$$$$\begin{array}{}\text{(25)}& B_{23}=\frac{\partial F_2}{\partial P}=-\frac{k_{1}q-11-II}{1+k_{2}H},\end{array}$$$$\begin{array}{}\text{(26)}& B_{31}=\frac{\partial F_3}{\partial H}=-\frac{\omega 1+q-1I}{{1+k_{2}H}^2}\end{array}$$$$\begin{array}{}\text{(27)}& B_{32}=\frac{\partial F_3}{\partial I}=\frac{\omega qH}{1+k_{2}H},\end{array}$$$$\begin{array}{}\text{(28)}& B_{33}=\frac{\partial F_3}{\partial P}=\frac{\omega 1+q-1IH}{1+k_{2}H}-\gamma -m_3.\end{array}$$

(i) The Jacobian matrix of system (2) at $E_0=0,0,0$ is given by:$$\begin{array}{}\text{(29)}& J{E}_0=\begin{array}{ccc}r& 0& 0\\ \\ 0& \beta 1-\frac{1}{{\mathcal{R}}_0}& 0\\ \\ \omega & 0& -\gamma -m_3\end{array}.\end{array}$$

The eigenvalues are ${\eta }_1=r>0,{\eta }_2=\beta 1-1/{\mathcal{R}}_0$ and ${\eta }_3=-\gamma +m_3$. Thus, $E_0$ is unstable.

(ii) The Jacobian matrix of system (2) at $E_1=K,0,0$ is given by:$$\begin{array}{}\text{(30)}& J{E}_1=\begin{array}{ccc}-r& -m_2-m_{1}K& -\frac{k_{1}K}{1+k_{2}K}\\ \\ 0& \beta 1-\frac{1}{{\mathcal{R}}_1}& 0\\ \\ \frac{\omega }{{1+k_{2}K}^2}& \frac{\omega qK}{1+k_{2}K}& \gamma {\mathcal{R}}_{ec}-1-\frac{m_3}{\gamma }\end{array}.\end{array}$$

The characteristic polynomial is given by:$$\begin{array}{}\text{(31)}& \begin{array}{ccc}{P}_{E_1}\eta & =& \beta 1-\frac{1}{{\mathcal{R}}_1}-\eta {\eta }^2+r-\gamma {\mathcal{R}}_{ec}-1-\frac{m_3}{\gamma }\eta \hfill \\ & +& \beta 1-\frac{1}{{\mathcal{R}}_1}-\eta \frac{\gamma k_1{\mathcal{R}}_{ec}}{{1+k_{2}K}^2}-r\gamma {\mathcal{R}}_{ec}-1-\frac{m_3}{\gamma }=0.\end{array}\end{array}$$

Thus, the eigenvalues of $J{E}_1$ is ${\eta }_0=\beta 1-1/{\mathcal{R}}_1.$ According to Equation (16), we get ${\eta }_0<0.$ From Equation (16), we find$$\begin{array}{}\text{(32)}& \begin{array}{c}r-\gamma {\mathcal{R}}_{ec}-1-\frac{m_3}{\gamma }>0,\\ \frac{\gamma k_1{\mathcal{R}}_{ec}}{{1+k_{2}K}^2}-r\gamma {\mathcal{R}}_{ec}-1-\frac{m_3}{\gamma }>0.\end{array}\end{array}$$

In the light of the Routh–Hurwitz criteria, $E_1$ is locally asymptotically stable.

(iii) The Jacobian matrix of system (2) at $E_2=H^0,0,P^0$ is given by:$$\begin{array}{}\text{(33)}& J{E}_2=\begin{array}{ccc}{C}_{11}& C_{12}& C_{13}\\ \\ 0& C_{22}& 0\\ \\ C_{31}& C_{32}& 0\end{array},\end{array}$$where$$\begin{array}{}\text{(34)}& C_{11}=-\frac{m_{1}K-\gamma +m_3{m}_{1}K{k}_2+r}{K\omega -k_2{m}_3+\gamma },C_{12}=-\frac{m_1{m}_3+\gamma }{\omega -k_2{m}_3+\gamma },C_{13}=\frac{-k_1\gamma +m_3}{\omega },C_{21}=0,\end{array}$$$$\begin{array}{}\text{(35)}& \begin{array}{ccc}{C}_{22}& =& \beta 1-\frac{1}{{\mathcal{R}}_1}+\frac{K\gamma +m_3\theta r-m_{1}q-1k_2}{K\omega -k_2{m}_3+\gamma }\hfill \\ & -& \frac{q-11+k_{2}K\gamma {\mathcal{R}}_{ec}-1-m_3-\omega m_{1}K}{K\omega -k_2{m}_3+\gamma },\end{array}\end{array}$$$$\begin{array}{}\text{(36)}& C_{23}=0,C_{31}=\frac{{\omega -k_2{m}_3+\gamma }^2}{\omega },C_{32}=\gamma +m_3,C_{33}=0.\end{array}$$

The characteristic polynomial associated is as follows:$$\begin{array}{}\text{(37)}& P_{E_2}\eta =C_{22}-\eta {\eta }^2-C_{11}\eta -C_{13}{C}_{31}=0.\end{array}$$

Then, the eigenvalue of $J{E}_2$ is ${\eta }_0=C_{22}.$ From conditions (16) and (17), we get ${\eta }_0<0.$ According to Equation (18), we have $C_{11}<0.$ In the light of the Routh–Hurwitz criteria, $E_2$ is locally asymptotically stable.

3.3. Existence of Coexistence Equilibria

To compute the equilibrium solutions, we set the right-hand-side of system (2) to zero. Thus, we get:$$\begin{array}{}\text{(38)}& \begin{array}{cc}rH1-\frac{H}{K}-\frac{k_{1}PH1+Iq-1}{1+k_{2}H}-m_1+m_2-m_{1}I\hfill & H=0,\\ \beta 1-II-b+\lambda -\frac{\theta rH}{K}I-\frac{k_{1}1-Iq-1IP}{1+k_{2}H}-m_2\hfill & I=0,\\ \frac{\omega 1+Iq-1HP}{1+k_{2}H}-\gamma +m_3\hfill & P=0.\end{array}\end{array}$$

Assume that $0<H^{\ast }<K,0<I^{\ast }<1,$ and $P^{\ast }>0.$ By using the third equation of system (38) and dividing by $P^{\ast }$, we get:$$\begin{array}{}\text{(39)}& H^{\ast }=\frac{\gamma +m_3}{B_0},\end{array}$$$$\begin{array}{}\text{(40)}& \text{if}\omega >k_2\gamma +m_3,\end{array}$$where $B_0=\omega -k_2\gamma +m_3+\omega q-1I^{\ast }.$

By plugging $H^{\ast }$ in the second equation of system (38) we have:$$\begin{array}{}\text{(41)}& P^{\ast }=\frac{\beta K1-I^{\ast }{B}_0-B_{0}K{m}_2+\lambda +b-\theta r\gamma +m_3}{k_{1}K1-I^{\ast }q-1B_0},\end{array}$$$$\begin{array}{}\text{(42)}& \text{if}\frac{\beta K1-I^{\ast }{B}_0}{B_{0}K{m}_2+\lambda +b-\theta r\gamma +m_3}>1.\end{array}$$

By substituting $P^{\ast }$ and $H^{\ast }$ in the first equation of system (38), we get the following cubic equation verify by $I^{\ast }$:$$\begin{array}{}\text{(43)}& \begin{array}{c}{A}_3{I}^{{\ast }^3}+A_2{I}^{{\ast }^2}{A}_1{I}^{\ast }+A_0=0,\\ A_0=\frac{{\Gamma }_4-{\Gamma }_5}{\omega k_1{q-1}^{2}K\beta +m_1-m_2},\end{array}\end{array}$$$$\begin{array}{}\text{(44)}& \begin{array}{c}{A}_1=\frac{{\Gamma }_2-{\Gamma }_3}{\omega k_1{q-1}^{2}K\beta +m_1-m_2},\end{array}\end{array}$$$$\begin{array}{}\text{(45)}& \begin{array}{c}{A}_2=\frac{{\Gamma }_0-{\Gamma }_1}{\omega k_1{q-1}^{2}K\beta +m_1-m_2},\\ A_3=1,\end{array}\end{array}$$where$$\begin{array}{}\text{(46)}& \begin{array}{ccc}{\Gamma }_0& =& k_{1}Kq-11+m_1-m_2\omega -k_2\gamma +m_3+\omega q-1+\beta +m_{1}q-1,\hfill \\ {\Gamma }_1& =& \omega k_{1}K{q-1}^2\beta +r,\hfill \\ {\Gamma }_2& =& r{k}_{1}q-1q-1\omega K+\theta \gamma +m_3+k_{1}K\omega -k_2\gamma +m_{3}q-12m_1-m_2+\beta ,\hfill \\ {\Gamma }_3& =& k_{1}q-1Kr\omega -k_2\gamma +m_3+m_1\omega q-1+\beta \omega ,\hfill \\ {\Gamma }_4& =& r{k}_{1}q-1K\omega -k_2\gamma +m_3+\theta \gamma +m_3,\hfill \\ {\Gamma }_5& =& k_{1}K\omega -k_2\gamma +m_3\beta +m_{1}q-1.\hfill \end{array}\end{array}$$

Since $A_3>0$, thus, the number of possible positive real roots of Equation (43) depends on the signs of $A_2,A_1,A_0$. This can be analyzed using the Descartes Rules of Signs [25] of the polynomial: $LX=A_3{X}^3+A_2{X}^2+A_{1}X+A_0$ with $X=I^{\ast }$. The different cases for the positive real roots of $LX$ are presented in Table 1.

Table 1

Various possibilities for positive real roots of $L{X}^{\ast }=0$.

Cases	$A_3$	$A_2$	$A_1$	$A_0$	Number of sign changes	Number of positive roots
1	+	+	+	+	0	0
2	+	+	+	—	1	1
3	+	+	—	—	1	1
4	+	—	—	—	1	1
5	+	—	+	+	2	0,2
6	+	+	—	+	2	0,2
7	+	—	—	+	2	0,2
8	+	—	+	—	3	1,3

Then, the following result gives the existence of the coexistence equilibrium point [25].

Theorem 4.

System (2):

(i) has no feasible coexistence equilibria points if case 1 is satisfied;

(ii) has a unique feasible coexistence equilibrium points if the conditions (40) and (42) and cases 2, 3, and 4 are satisfied;

(iii) has two feasible coexistence equilibria or no feasible coexistence equilibria points if the conditions (40) and (42) and cases 5, 6, and 7 are satisfied;

(iv) has three feasible coexistence equilibria points or a unique feasible coexistence equilibrium point if the conditions (40) and (42) and case 8 are satisfied.

3.4. Local Stability Analysis of the Coexistence Equilibrium Point

Let us define ${\kappa }_2={\Delta }_1/{\Delta }_2\text{and}{\kappa }_3={\Delta }_3/{\Delta }_4$ where,$$\begin{array}{}\text{(47)}& \begin{array}{ccc}{\Delta }_1& =& {1+k_2{H}^{\ast }}^2\omega r{I}^{\ast }{H}^{{\ast }^2}\theta m_1-m_2+\beta 1+k_2{H}^{\ast }\hfill \\ & +& \omega k_1{k}_{2}K{H}^{{\ast }^2}{m}_1-m_{2}1-I^{\ast }+k_1{I}^{\ast }{P}^{\ast }+k_1{k}_2{I}^{\ast }q-1.\\ {\Delta }_2& =& {\omega }_{1}q-1H^{{\ast }^2}{P}^{\ast }r1+\theta I^{\ast }{1+k_2{H}^{\ast }}^2+k_1{k}_{2}Kq-1K1-I^{\ast }{P}^{\ast }\hfill \\ & +& \omega \beta k_1{k}_{2}K{H}^{{\ast }^2}{I}^{\ast }{1+k_2{H}^{\ast }}^2.\\ {\Delta }_3& =& \gamma +m_{3}1+k_2{H}^{\ast }\omega H^{\ast }K{m}_1-m_2{k}_{1}q-11-I^{\ast }+q\gamma +m_3+\beta I^{\ast }1+k_2{H}^{\ast },\hfill \\ {\Delta }_4& =& \omega H^{\ast }q-1k_1\omega r1-I^{\ast }q{H}^{{\ast }^2}1+k_2{H}^{\ast }+k_{1}q-1\gamma +m_{3}K{P}^{\ast }\hfill \\ & +& {\omega k_1}^2{k}_{2}Kqq-11-I^{\ast }{P}^{\ast }+k_1{k}_2{I}^{\ast }q-1H^{{\ast }^3}+k_1\gamma +m_{3}q-1I^{\ast }{P}^{\ast }1+k_2{H}^{\ast }.\end{array}\end{array}$$

Proposition 3.

If the conditions (ii) of Theorem 4 and the following conditions are satisfied:$$\begin{array}{}\text{(48)}& \frac{k_{1}Kq-11+k_2{H}^{\ast }{I}^{\ast }{P}^{\ast }}{{1+k_2{H}^{\ast }}^2\beta K{I}^{\ast }+r{H}^{\ast }+k_1{k}_{2}K{H}^{\ast }{P}^{\ast }+k_1{k}_2{I}^{\ast }q-1}<1,\end{array}$$$$\begin{array}{}\text{(49)}& {\kappa }_2<1,\end{array}$$$$\begin{array}{}\text{(50)}& {\kappa }_3<1,\end{array}$$$$\begin{array}{}\text{(51)}& \frac{k_{1}q-1P^{\ast }1+\omega q{I}^{\ast }{H}^{\ast }}{1+k_2{H}^{\ast }{m}_1-m_2\gamma +m_3+\beta \omega q{I}^{\ast }{H}^{\ast }}<1,\end{array}$$then, the coexistence equilibrium $E_3=H^{\ast },I^{\ast },P^{\ast }$ point is locally asymptotically stable.

Proof.

Indeed, by linearizing model (2) around the coexistence equilibrium point $E_3$ yields:$$\begin{array}{}\text{(52)}& J{E}_3=\begin{array}{ccc}{D}_{11}& D_{12}& D_{13}\\ \\ D_{21}& D_{22}& D_{23}\\ \\ D_{31}& D_{32}& 0\end{array},\end{array}$$where$$\begin{array}{}\text{(53)}& D_{11}=-\frac{k_1{k}_2{H}^{\ast }{P}^{\ast }+k_1{k}_2{I}^{\ast }q-1}{{1+k_2{H}^{\ast }}^2}-\frac{r{H}^{\ast }}{K},D_{12}=-\frac{k_{1}q-1H^{\ast }{P}^{\ast }}{1+k_2{H}^{\ast }}-m_2-m_1{H}^{\ast },\end{array}$$$$\begin{array}{}\text{(54)}& D_{13}=-\frac{\gamma +m_3}{\omega },D_{21}=\frac{\theta r{I}^{\ast }}{K}+\frac{k_1{k}_{2}q-11-I^{\ast }{P}^{\ast }}{{1+k_2{H}^{\ast }}^2},\end{array}$$$$\begin{array}{}\text{(55)}& D_{22}=\frac{k_{1}q-1I^{\ast }{P}^{\ast }-\beta I^{\ast }1+k_2{H}^{\ast }}{1+k_2{H}^{\ast }},D_{23}=-\frac{k_{1}q-11-I^{\ast }}{1+k_2{H}^{\ast }},D_{31}=-\frac{\gamma +m_3}{H^{\ast }1+k_2{H}^{\ast }},\end{array}$$$$\begin{array}{}\text{(56)}& D_{32}=\frac{\omega q{H}^{\ast }}{1+k_2{H}^{\ast }},D_{33}=0.\end{array}$$

Therefore, the characteristic polynomial is given by:$$\begin{array}{}\text{(57)}& P_{E_3}\eta =-{\eta }^3+a_1{\eta }^2+a_2\eta +a_3=0,\end{array}$$where$$\begin{array}{}\text{(58)}& a_1=D_{11}+D_{22},a_2=-D_{11}{D}_{22}+D_{23}{D}_{32}+D_{12}{D}_{21}+D_{13}{D}_{31},\text{and}\end{array}$$$$\begin{array}{}\text{(59)}& a_3=D_{13}{D}_{21}{D}_{32}-D_{11}{D}_{23}{D}_{32}+D_{12}{D}_{23}{D}_{31}-D_{13}{D}_{22}{D}_{31}.\end{array}$$

From expressions (48) and (50), we get, respectively, $-a_1>0$ and $-a_3>0.$ Moreover, we have:$$\begin{array}{}\text{(60)}& a_1{a}_2+a_3=D_{11}+D_{22}{D}_{13}{D}_{31}+D_{12}{D}_{21}-D_{11}{D}_{22}+D_{23}{D}_{12}{D}_{31}+D_{22}{D}_{32}.\end{array}$$

According to Equations (49) and (51), we get, respectively, $D_{13}{D}_{31}+D_{12}{D}_{21}-D_{11}{D}_{22}<0$ and $D_{12}{D}_{31}+D_{22}{D}_{32}<0.$ Therefore, $a_1{a}_2+a_3>0.$ Then, referring to the Routh–Hurwitz criterion, it then follows that $E_3$ is locally asymptotically stable.

3.5. Bifurcation Analysis of the Coexistence Equilibrium $E_3$

Here, we establish the conditions when Hopf-bifurcation occurs at $E_3$. The parameter $m_1$ play a major role in the topological changing of coexistence equilibrium. A Hopf-bifurcation occurs if there exist a certain parameter $m_1=m_{cr}$ verifying:$$\begin{array}{}\text{(61)}& \mathcal{F}{m}_{cr}=a_1{m}_{cr}{a}_2{m}_{cr}+a_3{m}_{cr}=0,\end{array}$$$$\begin{array}{}\text{(62)}& -a_2{m}_{cr}>0,\end{array}$$$$\begin{array}{}\text{(63)}& {\mathcal{R}}_e{\frac{d\eta m_1}{d{m}_1}}_{m_1=m_{cr}}\ne 0,\end{array}$$where $\eta$ is the solution of Equation (57). Let us consider:$$\begin{array}{}\text{(64)}& \begin{array}{ccc}{\beta }_0& =& k_{1}q-11-I^{\ast }\omega H^{\ast }{k}_{1}q-1I^{\ast }{P}^{\ast }-\beta I^{\ast }1+k_2{H}^{\ast }\omega K^2{1+k_2{H}^{\ast }}^3\hfill \\ {\beta }_1& =& \omega k_{1}q-1I^{\ast }{P}^{\ast }-\beta I^{\ast }1+k_2{H}^{\ast }{k_1{k}_{2}K{H}^{\ast }{P}^{\ast }+k_1{k}_2{I}^{\ast }q-1+r{H}^{\ast }{1+k_2{H}^{\ast }}^2}^2\hfill \\ {\beta }_2& =& {k_{1}q-1I^{\ast }{P}^{\ast }-\beta I^{\ast }1+k_2{H}^{\ast }}^2\omega K{k}_1{k}_{2}1+k_2{H}^{\ast }{P}^{\ast }+k_1{k}_2{I}^{\ast }q-1+r{H}^{\ast }{1+k_2{H}^{\ast }}^2\hfill \\ {\beta }_3& =& \omega K{1+k_2{H}^{\ast }}^{3}K{k}_1{k}_{2}1+k_2{H}^{\ast }{P}^{\ast }+k_1{k}_2{I}^{\ast }q-1+r{H}^{\ast }{1+k_2{H}^{\ast }}^2{\gamma +m_3}^2.\hfill \end{array}\end{array}$$

The following theorem gives the conditions that Hopf-bifurcation occurs [15, 26–28].

Theorem 5.

If the conditions (ii) of Theorem 4 and the following conditions are satisfied:$$\begin{array}{}\text{(65)}& \frac{m_{2}1+k_2{H}^{\ast }}{k_{1}q-1P^{\ast }}>1,\end{array}$$$$\begin{array}{}\text{(66)}& \frac{\omega K^{2}q{\gamma +m_3}^2{1+k_2{H}^{\ast }}^3-{\beta }_0-{\beta }_1+{\beta }_2}{{\beta }_3}>1,\end{array}$$$$\begin{array}{}\text{(67)}& \frac{{\beta }_4{\beta }_5{H}^{{\ast }^2}}{k_1\gamma +m_{3}q-11-I^{\ast }{K}^2{1+k_2{H}^{\ast }}^2}>1,\end{array}$$then a Hopf-bifurcation occurs for system (2) around the coexistence point $E_3=H^{\ast },I^{\ast },P^{\ast }$ when the prey migration rate $m_1$ crosses a critical value $m_1=m_{cr},$ where$$\begin{array}{}\text{(68)}& m_{cr}=m_2-\frac{k_{1}q-1P^{\ast }}{1+k_2{H}^{\ast }}+\frac{q{\gamma +m_3}^2{\omega }^2{1+k_2{H}^{\ast }}^3-{\beta }_0-{\beta }_1+{\beta }_2-{\beta }_3{H}^{{\ast }^2}}{\omega 1+k_2{H}^{\ast }{\beta }_4{\beta }_5{H}^{{\ast }^2}-k_1{K}^2\gamma +m_{3}q-11-I^{\ast }1+k_2{H}^{\ast }}.\end{array}$$

Proof.

Indeed, assuming that $H-H^{\ast }\simeq e^{\eta t},I-I^{\ast }\simeq e^{\eta t}$ and $P-P^{\ast }\simeq e^{\eta t}$, the characteristic polynomial defined by Equation (57) is given by:$$\begin{array}{}\text{(69)}& P_{E_3}\eta ;m_1=0\iff {\eta }^3-a_1{m}_1{\eta }^2-a_2{m}_1\eta -a_3{m}_1=0,\end{array}$$where$$\begin{array}{}\text{(70)}& \begin{array}{c}{a}_1{m}_1=D_{11}+D_{22}\hfill \\ a_2{m}_1=D_{23}{D}_{32}+D_{12}{D}_{21}+D_{13}{D}_{31}-D_{11}{D}_{22}\hfill \\ a_3{m}_1=D_{13}{D}_{31}{D}_{32}-D_{11}{D}_{23}{D}_{32}+D_{12}{D}_{23}{D}_{31}-D_{13}{D}_{22}{D}_{31}\hfill \end{array}.\end{array}$$

According to expressions (49), (65), and (66), we get $-a_2{m}_1>0$ and $m_{cr}>0.$

At $m_1=m_{cr}$, it is easy to see that:$$\begin{array}{}\text{(71)}& \mathcal{F}{m}_{cr}=0.\end{array}$$

According to $-a_3{m}_1=-a_1{m}_1{a}_2{m}_1$ and from Equation (69), we get:$$\begin{array}{}\text{(72)}& {\eta }^2-a_2\eta -a_1=0,\end{array}$$whose roots are$$\begin{array}{}\text{(73)}& {\eta }_0=i\sqrt{-a_2},{\eta }_1=-i\sqrt{-a_2},{\eta }_2=a_1,\end{array}$$where $a_1{m}_{cr}>0.$

These roots are functions of $m_1,$ where $m_1\in m_{cr}-ϵ,m_{cr}+ϵ,$ for $m_{cr}>ϵ>0$ and subsequently can be written as follows:$$\begin{array}{}\text{(74)}& {\eta }_0{m}_1=x_1{m}_1+i{x}_2{m}_1,{\eta }_1{m}_1=x_1{m}_1-i{x}_2{m}_1,{\eta }_2{m}_1=a_1{m}_1.\end{array}$$

We have:$$\begin{array}{}\text{(75)}& {\mathcal{R}}_e{\frac{d\eta }{d{m}_1}}_{m_1=m_{cr}}\ne 0.\end{array}$$

Differentiating Equation (69) with respect $m_1$ yields:$$\begin{array}{}\text{(76)}& \begin{array}{ccc}\frac{d\eta }{d{m}_1}& =& \frac{a_1^{\prime }{m}_1{\eta }^2+a_2^{\prime }{m}_1\eta +a_3^{\prime }{m}_1}{3{\eta }^2-2a_1{m}_1\eta -a_2{m}_1}{}_{\eta =i\sqrt{-a_2{m}_1}}\hfill \\ & =& \frac{a_1^{\prime }{m}_1{a}_2{m}_1+a_3^{\prime }{m}_1+i{a}_2^{\prime }{m}_1\sqrt{-a_2{m}_1}}{2a_2{m}_1-2i{a}_1{m}_1\sqrt{-a_2{m}_1}}\hfill \\ & =& \frac{a_1^{\prime }{m}_1{a}_2^2{m}_1+a_3^{\prime }{m}_1{a}_2{m}_1+a_1{m}_1{a}_2{m}_1{a}_2^{\prime }{m}_1}{2a_2^2{m}_1-a_1^2{m}_1{a}_2{m}_1}\hfill \\ & +& i\frac{\sqrt{-a_2{m}_1}{a}_1{a}_1^{\prime }{m}_1{a}_2{m}_1+a_3^{\prime }{m}_1{a}_1{m}_1+a_2^{\prime }{m}_1{a}_2{m}_1}{2a_2^2{m}_1-a_1^2{m}_1{a}_2{m}_1}\\ & =& \frac{a_2{a_1{a}_2+a_3}^{\prime }}{2a_2^2{m}_1-a_1^2{m}_1{a}_2{m}_1}\hfill \\ & +& i\sqrt{-a_2{m}_1}\frac{{a_2^2}^{\prime }{m}_1}{4a_2^2{m}_1-a_1^2{m}_1{a}_2{m}_1}+\frac{a_1{m}_1{a}_1^{\prime }{m}_1{a}_2{m}_1+a_3^{\prime }{m}_1}{2a_2^2{m}_1-a_1^2{m}_1{a}_2{m}_1}\\ & =& \frac{a_2\frac{dF{m}_1}{d{m}_1}}{2a_2^2{m}_1-a_1^2{m}_1{a}_2{m}_1}\hfill \\ & +& i\sqrt{-a_2{m}_1}\frac{{a_2^2}^{\prime }{m}_1}{4a_2^2{m}_1-a_1^2{m}_1{a}_2{m}_1}+\frac{a_1{m}_1{a}_1^{\prime }{m}_1{a}_2{m}_1+a_3^{\prime }{m}_1}{2a_2^2{m}_1-a_1^2{m}_1{a}_2{m}_1}.\end{array}\end{array}$$

By separating real and imaginary parts, one has:$$\begin{array}{}\text{(77)}& {\mathcal{R}}_e\frac{d\eta }{d{m}_1}=\frac{\frac{dF{m}_1}{d{m}_1}}{2a_2{m}_1-a_1^2{m}_1}.\end{array}$$

Verifying that ${\mathcal{R}}_e{d\eta /d{m}_1}_{m_1=m_{cr}}\ne 0,$ we have:$$\begin{array}{}\text{(78)}& \begin{array}{ccc}\frac{dF{m}_1}{d{m}_1}& =& a_1^{\prime }{m}_1{a}_2{m}_1+a_3^{\prime }{m}_1+a_1{m}_1{a}_2^{\prime }{m}_1\hfill \\ & =& a_3^{\prime }{m}_1+a_1{m}_1{a}_2^{\prime }{m}_1,\text{because}{a}_1^{\prime }{m}_1=0\hfill \\ & =& D_{23}{D}_{31}+D_{21}{D}_{11}+D_{22}{D}_{12}^{\prime }{m}_1\hfill \\ & =& \frac{k_{1}q-11-I^{\ast }\gamma +m_3{H}^{\ast }}{{1+k_2{H}^{\ast }}^2{H}^{\ast }}\hfill \\ & +& \frac{\theta r{I}^{\ast }{1+k_2{H}^{\ast }}^2+K{k}_1{k}_{2}q-11-I^{\ast }{P}^{\ast }{k}_{1}q-1I^{\ast }{P}^{\ast }-\beta I^{\ast }1+k_2{H}^{\ast }{H}^{\ast }}{K{1+k_2{H}^{\ast }}^3}\\ & -& \frac{k_1{k}_2{H}^{\ast }{P}^{\ast }+k_1{k}_2{I}^{\ast }q-1\theta r{I}^{\ast }{1+k_2{H}^{\ast }}^2+K{k}_1{k}_{2}q-11-I^{\ast }{P}^{\ast }{H}^{\ast }}{K{1+k_2{H}^{\ast }}^3}.\end{array}\end{array}$$

At $m_1=m_{cr},$ it is straightforward to see that $dF{m}_1/d{m}_1\ne 0,$ as a result:$$\begin{array}{}\text{(79)}& {\mathcal{R}}_e{\frac{d\eta }{d{m}_1}}_{m_1=m_{cr}}=\frac{\frac{dF{m}_1}{d{m}_1}}{2a_2{m}_1-a_1^2{m}_1}{}_{m_1=m_{cr}}\ne 0.\end{array}$$

Therefore, the transversality condition hold. Thus, Hopf-bifurcation occurs at $m_1=m_{cr}$ for system (2) around the coexistence equilibrium $E_3.$

Now, let us investigate the global stability of the coexistence equilibrium point [12, 26, 29–32].

3.6. Global Stability Analysis of the Coexistence Equilibrium Point

Let us define:$$\begin{array}{}\text{(80)}& M=D_1\mathrm{\setminus }0,0,0;K,1,0.\end{array}$$

Setting$$\begin{array}{}\text{(81)}& \begin{array}{ccc}{\Delta }_5& =& 2k_1{k}_{2}K{P}^{\ast }+q-1I^{\ast }+\omega K1+q-1I^{\ast }+k_1{k}_{2}Kq-1P^{\ast }1-I^{\ast }\hfill \\ & +& m_1-m_{2}K1+k_{2}K1+k_2{H}^{\ast }+\theta rK1+k_{2}K1+k_2{H}^{\ast },\\ {\Delta }_6& =& 1+k_2{H}^{\ast }2r+k_{1}Kq,\hfill \\ {\Delta }_7& =& k_{1}Kq-11+k_2{H}^{\ast }+2+k_2{P}^{\ast }1-I^{\ast }\hfill \\ & +& 1+k_{2}K1+k_2{H}^{\ast }\theta r+m_1-m_{2}K+\omega q-1K1+k_2{H}^{\ast },\\ {\Delta }_8& =& 2\beta K+k_{1}Kq-11+k_2{H}^{\ast }+k_{1}Kq-1P^{\ast }1-I^{\ast },\hfill \\ {\Delta }_9& =& \omega q-1I^{\ast }+1+k_2{H}^{\ast }1+H^{\ast },\text{and}\hfill \\ {\Delta }_{10}& =& k_{1}q-1P^{\ast }1-I^{\ast }+1+k_2{H}^{\ast }.\hfill \end{array}\end{array}$$

The following theorem give the global stability of coexistence equilibrium.

Theorem 6.

Setting ${\kappa }_4={\Delta }_5/{\Delta }_6,{\kappa }_5={\Delta }_7/{\Delta }_8\text{and}{\kappa }_6={\Delta }_9/{\Delta }_{10}.$ If the conditions (ii) of Theorem 4 and the following conditions are satisfied:$$\begin{array}{}\text{(82)}& {\kappa }_4<1,{\kappa }_5<1,\text{and}{\kappa }_5<1,\end{array}$$then the coexistence equilibrium $E_3=H^{\ast },I^{\ast },P^{\ast }$ is globally asymptotically stable in the region $M.$

Proof.

Indeed, let us consider the function $L$ defined by:$$\begin{array}{}\text{(83)}& LH,I,P=H-H^{\ast }-H^{\ast }\text{ln}\frac{H^{\ast }}{H}+I-I^{\ast }-I^{\ast }\text{ln}\frac{I^{\ast }}{I}+P-P^{\ast }-P^{\ast }\text{ln}\frac{P^{\ast }}{P}.\end{array}$$

It is straightforward to see that $L{H}^{\ast },I^{\ast },P^{\ast }=0$ and for all $H,I,P\in M$ such that $H,I,P\ne H^{\ast },I^{\ast },P^{\ast }$, $LH,I,P>0.$ Hence, $L$ is well-defined on $M$.

The time derivative of $L$ along the solutions of system (2) is given by:$$\begin{array}{}\text{(84)}& \begin{array}{ccc}\dot{L}H,I,P& =& \frac{H-H^{\ast }}{H}\dot{H}+\frac{I-I^{\ast }}{I}\dot{I}+\frac{P-P^{\ast }}{P}\dot{P}\hfill \\ & =& H-H^{\ast }r1-\frac{H}{K}-\frac{k_{1}P1+I{q}_1-1}{1+k_{2}H}-m_1+m_2-m_{1}I\hfill \\ & +& I-I^{\ast }\beta 1-I-b+\lambda -\frac{\theta rH}{K}-\frac{k_{1}1-Iq-1P}{1+k_{2}H}-m_2\\ & +& P-P^{\ast }\frac{\omega 1+Iq-1H}{1+k_{2}H}-\gamma -m_3.\end{array}\end{array}$$

Using the fact that $E_3$ is an equilibrium point of system (2), we have:$$\begin{array}{}\text{(85)}& \begin{array}{c}r=r\frac{H^{\ast }}{K}+\frac{k_1{P}^{\ast }1+I^{\ast }q-1}{1+k_{2}H}+m_1+m_2-m_1{I}^{\ast },\hfill \\ \beta =\beta I^{\ast }+b+\lambda -\frac{\theta r{H}^{\ast }}{K}+\frac{k_{1}1-I^{\ast }q-1P^{\ast }}{1+k_2{H}^{\ast }}+m_2,\hfill \\ \gamma =\frac{\omega 1+I^{\ast }q-1H^{\ast }}{1+k_2{H}^{\ast }}+m_3.\hfill \end{array}\end{array}$$

By substituting $r,\beta$, and $\gamma$ into Equation (84) yields:$$\begin{array}{}\text{(86)}& \dot{L}H,I,P=L_{1}H,I,P+L_{2}H,I,P+L_{3}H,I,P,\end{array}$$where$$\begin{array}{}\text{(87)}& \begin{array}{ccc}{L}_{1}H,I,P& =& H-H^{\ast }\frac{k_{1}1+I^{\ast }q-1P^{\ast }}{1+k_2{H}^{\ast }}-\frac{k_{1}1+Iq-1P}{1+k_{2}H}-\frac{r}{K}H-H^{\ast }\hfill \\ & +& m_1-m_{2}H-H^{\ast }I-I^{\ast },\end{array}\end{array}$$$$\begin{array}{}\text{(88)}& \begin{array}{ccc}{L}_{2}H,I,P& =& I-I^{\ast }\frac{\theta r}{K}H-H^{\ast }+\frac{k_{1}q-11-I^{\ast }{P}^{\ast }}{1+k_2{H}^{\ast }}-\frac{k_{1}q-11-IP}{1+k_{2}H}-\beta I-I^{\ast },\hfill \\ L_{3}H,I,P& =& P-P^{\ast }\frac{\omega 1+Iq-1H}{1+k_{2}H}-\frac{\omega 1+I^{\ast }q-1H^{\ast }}{1+k_2{H}^{\ast }}.\hfill \end{array}\end{array}$$

It is obvious to see that:$$\begin{array}{}\text{(89)}& \begin{array}{ccc}{k}_1\frac{1+I^{\ast }q-1P^{\ast }}{1+k_2{H}^{\ast }}-\frac{1+Iq-1P}{1+k_{2}H}& =& k_1\frac{k_2{P}^{\ast }H-H^{\ast }-P-P^{\ast }1+k_2{H}^{\ast }}{1+k_2{H}^{\ast }1+k_{2}H}\hfill \\ & +& k_{1}q-1\frac{k_2{I}^{\ast }H-H^{\ast }-I-I^{\ast }1+k_2{H}^{\ast }}{1+k_2{H}^{\ast }1+k_{2}H},\end{array}\end{array}$$$$\begin{array}{}\text{(90)}& \begin{array}{ccc}\frac{k_{1}q-11-I^{\ast }{P}^{\ast }}{1+k_2{H}^{\ast }}-\frac{k_{1}q-11-IP}{1+k_{2}H}& =& \frac{k_1{k}_{2}q-1P^{\ast }1-I^{\ast }H-H^{\ast }}{1+k_2{H}^{\ast }1+k_{2}H}\hfill \\ & -& \frac{k_{1}q-1P^{\ast }1-I^{\ast }P-P^{\ast }}{1+k_2{H}^{\ast }1+k_{2}H}\\ & +& \frac{k_{1}q-1P^{\ast }1-I^{\ast }I-I^{\ast }}{1+k_2{H}^{\ast }1+k_{2}H}\\ & +& \frac{k_{1}q-11+k_2{H}^{\ast }P-P^{\ast }I-I^{\ast }}{1+k_2{H}^{\ast }1+k_{2}H},\\ \omega \frac{1+Iq-1H}{1+k_{2}H}-\frac{1+I^{\ast }q-1H^{\ast }}{1+k_2{H}^{\ast }}& =& \frac{\omega 1+q-1I^{\ast }H-H^{\ast }}{1+k_2{H}^{\ast }1+k_{2}H}\hfill \\ & +& \frac{\omega q-1H^{\ast }1+k_2{H}^{\ast }I-I^{\ast }}{1+k_2{H}^{\ast }1+k_{2}H}\\ & +& \frac{\omega q-11+k_2{H}^{\ast }I-I^{\ast }H-H^{\ast }}{1+k_2{H}^{\ast }1+k_{2}H}.\end{array}\end{array}$$

Plugging the above expressions into Equation (86) and rearranging gives:$$\begin{array}{}\text{(91)}& \begin{array}{ccc}\dot{L}H,I,P& =& \frac{k_1{k}_2{P}^{\ast }+q-1I^{\ast }}{1+k_2{H}^{\ast }1+k_{2}H}-\frac{r}{K}{H-H^{\ast }}^2+\frac{k_1{P}^{\ast }q-11-I^{\ast }}{1+k_2{H}^{\ast }1+k_{2}H}-\beta {I-I^{\ast }}^2\hfill \\ & +& \frac{k_{1}q-11+k_2{H}^{\ast }}{1+k_2{H}^{\ast }1+k_{2}H}P-P^{\ast }{I-I^{\ast }}^2\\ & +& \frac{\omega 1+q-1I^{\ast }-k_{1}1+k_2{H}^{\ast }}{1+k_2{H}^{\ast }1+k_{2}H}H-H^{\ast }P-P^{\ast }\\ & +& \frac{k_{1}q-1k_2{P}^{\ast }1-I^{\ast }-1+k_2{H}^{\ast }}{1+k_2{H}^{\ast }1+k_{2}H}+m_1-m_2+\frac{\theta r}{K}H-H^{\ast }I-I^{\ast }\\ & +& \frac{\omega q-11+k_2{H}^{\ast }{H}^{\ast }-k_{1}q-1P^{\ast }1-I^{\ast }}{1+k_2{H}^{\ast }1+k_{2}H}P-P^{\ast }I-I^{\ast }\\ & +& \frac{\omega q-11+k_2{H}^{\ast }}{1+k_2{H}^{\ast }1+k_{2}H}H-H^{\ast }I-I^{\ast }P-P^{\ast }.\end{array}\end{array}$$

Noticing that,$$\begin{array}{}\text{(92)}& \frac{e1+q-1I^{\ast }}{1+k_2{H}^{\ast }}>1,\frac{k_2{P}^{\ast }1-I^{\ast }}{1+k_2{H}^{\ast }}>1,\frac{e{H}^{\ast }1+k_2{H}^{\ast }}{P^{\ast }1-I^{\ast }}>1,\text{and}1<1+k_{2}H<1+k_{2}K,\mathrm{\forall }H\in 0,K.\end{array}$$

Applying the following classical relation:$$\begin{array}{}\text{(93)}& xy\le \frac{x^2}{2}+\frac{y^2}{2},\mathrm{\forall }x,y\in \mathbb{R},\end{array}$$on$$\begin{array}{}\text{(94)}& P-P^{\ast }\times 1,H-H^{\ast }P-P^{\ast },H-H^{\ast }I-I^{\ast },H-H^{\ast }I-I^{\ast }P-P^{\ast },\end{array}$$we get:$$\begin{array}{}\text{(95)}& \frac{k_{1}q-11+k_2{H}^{\ast }}{1+k_2{H}^{\ast }1+k_{2}H}P-P^{\ast }{I-I^{\ast }}^2<\frac{k_{1}q-11+k_2{H}^{\ast }}{1+k_2{H}^{\ast }1+k_{2}H}{I-I^{\ast }}^2\frac{{P-P^{\ast }}^2}{2}+\frac{1}{2},\end{array}$$$$\begin{array}{}\text{(96)}& \begin{array}{ccc}\frac{\omega 1+q-1I^{\ast }-k_{1}1+k_2{H}^{\ast }}{1+k_2{H}^{\ast }1+k_{2}H}H-H^{\ast }P-P^{\ast }& <& \frac{\omega 1+q-1I^{\ast }-k_{1}1+k_2{H}^{\ast }}{1+k_2{H}^{\ast }1+k_{2}H}\\ & \times & \frac{{H-H^{\ast }}^2}{2}+\frac{{P-P^{\ast }}^2}{2},\end{array}\end{array}$$$$\begin{array}{}\text{(97)}& \begin{array}{c}\frac{k_{1}q-1k_2{P}^{\ast }1-I^{\ast }-1+k_2{H}^{\ast }}{1+k_2{H}^{\ast }1+k_{2}H}+m_1-m_2+\frac{\theta r}{K}H-H^{\ast }I-I^{\ast }<C_{1}H,I,\end{array}\end{array}$$$$\begin{array}{}\text{(98)}& \begin{array}{c}\frac{\omega q-11+k_2{H}^{\ast }{H}^{\ast }-k_{1}q-1P^{\ast }1-I^{\ast }}{1+k_2{H}^{\ast }1+k_{2}H}P-P^{\ast }I-I^{\ast }<C_{2}H,I,P\\ \frac{\omega q-11+k_2{H}^{\ast }}{1+k_2{H}^{\ast }1+k_{2}H}H-H^{\ast }I-I^{\ast }P-P^{\ast }<\frac{\omega q-11+k_2{H}^{\ast }}{1+k_2{H}^{\ast }1+k_{2}H}\\ \times \frac{{H-H^{\ast }}^2}{2}+\frac{{I-I^{\ast }}^2{P-P^{\ast }}^2}{2},\end{array}\end{array}$$where$$\begin{array}{}\text{(99)}& \begin{array}{ccc}{C}_{1}H,I& =& \frac{k_{1}q-1k_2{P}^{\ast }1-I^{\ast }-1+k_2{H}^{\ast }}{1+k_2{H}^{\ast }1+k_{2}H}+m_1-m_2+\frac{\theta r}{K}\frac{{H-H^{\ast }}^2}{2}+\frac{{I-I^{\ast }}^2}{2}.\hfill \\ C_{2}H,I,P& =& \frac{\omega q-11+k_2{H}^{\ast }{H}^{\ast }-k_{1}q-1P^{\ast }1-I^{\ast }}{1+k_2{H}^{\ast }1+k_{2}H}\frac{{P-P^{\ast }}^2}{2}+\frac{{I-I^{\ast }}^2}{2}.\hfill \end{array}\end{array}$$

By using the above relations, one has$$\begin{array}{}\text{(100)}& \begin{array}{ccc}\dot{L}H,I,P& <& \frac{2k_1{k}_{2}K{P}^{\ast }+q-1I^{\ast }+\omega K1+q-1I^{\ast }+k_1{k}_{2}Kq-1P^{\ast }1-I^{\ast }}{2K1+k_2{H}^{\ast }}{H-H^{\ast }}^2\\ & +& \frac{m_1-m_{2}K1+k_{2}K1+k_2{H}^{\ast }+\theta rK1+k_{2}K1+k_2{H}^{\ast }}{2K1+k_2{H}^{\ast }}{H-H^{\ast }}^2\\ & -& \frac{1+k_2{H}^{\ast }2r+k_{1}Kq}{2K1+k_2{H}^{\ast }}{H-H^{\ast }}^2+\frac{k_{1}Kq-11+k_2{H}^{\ast }+2+k_2{P}^{\ast }1-I^{\ast }}{2K1+k_2{H}^{\ast }}{I-I^{\ast }}^2\\ & +& \frac{1+k_{2}K1+k_2{H}^{\ast }\theta r+m_1-m_{2}K+\omega q-1K1+k_2{H}^{\ast }}{2K1+k_2{H}^{\ast }}{I-I^{\ast }}^2\\ & -& \frac{2\beta K+k_{1}Kq-11+k_2{H}^{\ast }+k_{1}Kq-1P^{\ast }1-I^{\ast }}{2K1+k_2{H}^{\ast }}{I-I^{\ast }}^2\\ & +& \frac{\omega q-1I^{\ast }+1+k_2{H}^{\ast }1+H^{\ast }-k_{1}q-1P^{\ast }1-I^{\ast }+1+k_2{H}^{\ast }}{2K1+k_2{H}^{\ast }}{P-P^{\ast }}^2.\end{array}\end{array}$$

However, by grouping term by term, we obtain:$$\begin{array}{}\text{(101)}& \dot{L}H,I,P<\frac{{\Delta }_5-{\Delta }_6}{2K1+k_2{H}^{\ast }}{H-H^{\ast }}^2+\frac{{\Delta }_7-{\Delta }_8}{2K1+k_2{H}^{\ast }}{I-I^{\ast }}^2+\frac{{\Delta }_9-{\Delta }_{10}}{2K1+k_2{H}^{\ast }}{P-P^{\ast }}^2.\end{array}$$

Consequently, $\dot{L}H,I,P<0$ whenever ${\kappa }_4<1,{\kappa }_5<1$ and ${\kappa }_6<1.$ In addition, $\dot{L}H,I,P=0$ if and only if $H,I,P=H^{\ast },I^{\ast },P^{\ast }.$ Therefore, thanks to LaSalle’s invariance principle, the coexistence equilibria $E_3$ is globally asymptomatically stable in the region $M.$

4. Numerical Simulation

In this section, numerical simulation has been performed to illustrate our theoretical findings by using the parameters given in Table 2 [1, 14]. Our goal here is to illustrate numerically the effect of migration on the dynamics of system (2). Through Figures 2–4 we present the evolution of the population of the prey and predators. Figure 2 shows the convergence of the solution of system (2) toward the coexistence equilibrium $E_3$ with decreasing of the infectious prey. By increasing the parameter $m_1=0.25$ and maintaining the other values fixed in Table 2, from Figure 3, we get oscillations of the preys and predators trajectories of system (2). Now, if we consider $m_1=0.3,$ a limit cycle occurs arising from Hopf-bifurcation with always a decreasing of the infectious prey, see Figure 4.

Table 2

Numerical values of the parameters of system (2) used in the simulation.

Parameters	Values	References
$r$	1.5	[8, 14]
$K$	366	Estimated
$m_1$	0.2	Estimated
$m_2$	0.03	Estimated
$m_3$	0.02	Estimated
$k_1$	0.1	Estimated
$k_2$	0.001	Estimated
$\omega$	0.015	[8, 14]
$\beta$	0.5	Estimated
$q$	2.5	Estimated
$\gamma$	1	[14]
$b$	0.15	[14]
$\theta$	0.02	Estimated
$\lambda$	0.5	Estimated

[figure(s) omitted; refer to PDF]

We continue our numerical analysis in order to study the effect of migration in the community of prey and predators. From Figure 5(e), we observe that the trajectories of system (2) approach asymptotically toward $E_3$ with the persistence of the disease in the prey population, see Figure 5(b) for $m_1=0.56$. Thus, Theorem 6 holds. By increasing the migration rate to $m_1=0.63$, from Figure 6, an instability take place around the equilibrium point $E_3$. This result confirms Proposition 3. Now, for the critical value $m_{cr}=0.73,$ Figure 7(a)–7(d) show an oscillations of the prey and predator trajectories around $E_3$ with always persistence of infectious disease. We also observe that system (2) undergoes a limit cycle arising from the Hopf-bifurcation (Figure 7(e)). Thus, Theorem 5 holds.

Remark 7.

The dynamics of system (2) present variations when the migration parameters are varied while keeping the other parameters fixed. Thus, when the migration rate of the susceptible prey population $m_1\in 0.2,0.3$, we observe the extinction of the infectious prey illustrated by Figures 2–4(b)). When $m_1\in 0.56,0.7$, it is straightforward to see that infectious prey population increases over time (Figure 5(b)), with decreasing of the density of predator population (Figure 5(c)) which will ultimately lead to an unstable system (Figures 6 and 7). With the regard to this observation, we can conclude that the migration rate $m_1$ is a key parameter that govern our model.

[figure(s) omitted; refer to PDF]

5. Conclusion

The process of migration in the dynamics of a prey–predator model in the presence of infectious disease play a major role in the natural mechanisms of regulation of species. It is in this line of thought that we are interested in this paper, to the study of the dynamics of prey–predator populations with infectious disease to describe the effects of migration in the dynamics of species. The formulation of the model derives from an ODE system by considering Holling function response of type II to represent the strategy of predation between the prey and the predator. The mathematical analysis allowed us first to establish that the model is ecologically and epidemiologically well-posedness. Thus, the existence, the positivity and the boundedness of the solutions are proved. Moreover, we established the conditions of existence of the coexistence equilibria. Under certain thresholds ${\mathcal{R}}_0,$${\mathcal{R}}_1$, ${\mathcal{R}}_{ec}$, and migration parameter, we were able to investigate the local stability of equilibria by using Routh–Hurwitz criterion. In order to show the long-term coexistence in the prey and predator communities, we established the global stability of the coexistence equilibrium with an appropriate Lyapunov function. In addition, we have described the conditions of existence of Hopf bifurcation around coexistence equilibrium in order to highlight periodic variations of the number of prey and predators due to the effect of migration.

We used different scenarios in the model numerical simulation in order to show the effect of migration on the dynamics of the prey and predator populations. Thus, for $m_1\in 0.2,0.3$, we observe a coexistence of the prey and predator species and extinction of disease in the infected prey population illustrated by Figure 2(b)–2(e). But if $m_1\in 0.56,0.7$, we observe that the trajectories approach asymptotically to $E_3,$ with the persistence of disease in the prey population, see Figure 5(b)-5(e). Also, we remark a decrease in the density of prey and predator populations see, respectively, Figures 5(a) and 5(c). These results are in perfect adequacy with the mathematical results established in Section $3$. In light of these observations, we are led to conclude that the migration rate is a key parameter that governs the dynamics of our model and is useful for understanding the dynamics of species of prey and predators.

In last, in this present work, the model is formulated by using the ordinary differential equations. However, it has been proven that the use of fractional derivatives gives a more realistic description of most biological issues. Therefore, for our future study, it will be interesting to consider fractional derivatives while formulating an eco-epidemiological model with migration that would give a better description of the ecological process [33, 34].

Disclosure

A preprint of this paper presenting the main findings is available online: https://papers.ssrn.com/sol3/papers.cfm?abstract_id=4347564 [31].

Authors’ Contributions

All the authors read and approved the final manuscript.

Acknowledgments

The authors are also very grateful to professor Pierre-Alexandre Bliman of the Sorbonne Université for his insightful and constructive comments and suggestions which helped to improve the original manuscript. This work was supported by the PDI-MSC program. The authors are very much thankful to Prof. Christophe Cambier of the Sorbonne Université (France), the coordinator of the program.

References

[1] A. Savadogo, B. Sangaré, H. Ouedraogo, "A mathematical analysis of Hopf-bifurcation in a prey–predator model with nonlinear functional response," Advances in Difference Equations, vol. 2021,DOI: 10.1186/s13662-021-03437-2, 2021.

[2] T. Chowdhury, S. Chakraborty, J. Chattopadhyay, "Migratory effect of middle predator in a tri-trophic food chain model," Mathematical Methods in the Applied Sciences, vol. 33 no. 14, pp. 1699-1711, DOI: 10.1002/mma.1286, 2010.

[3] W. O. Kermack, A. G. McKendrick, "A contribution to the mathematical theory of epidemics," Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character, vol. 115 no. 772, pp. 700-721, DOI: 10.1098/rspa.1927.0118, 1927.

[4] C. Arora, V. Kumar, "Dynamics of predator–prey system with migrating species and disease in prey population," Differential Equations and Dynamical Systems, vol. 29 no. 1, pp. 87-112, DOI: 10.1007/s12591-020-00529-5, 2021.

[5] S. Biswas, S. Samanta, J. Chattopadhyay, "A cannibalistic eco-epidemiological model with disease in predator population," Journal of Applied Mathematics and Computing, vol. 57 no. 1-2, pp. 161-197, DOI: 10.1007/s12190-017-1100-9, 2018.

[6] S. Djilali, B. Ghanbari, "The influence of an infectious disease on a prey–predator model equipped with a fractional-order derivative," Advances in Difference Equations, vol. 2021,DOI: 10.1186/s13662-020-03177-9, 2021.

[7] M. Haque, D. Greenhalgh, "When a predator avoids infected prey: a model-based theoretical study," Mathematical Medicine and Biology, vol. 27 no. 1, pp. 75-94, DOI: 10.1093/imammb/dqp007, 2010.

[8] A. Savadogo, B. Sangaré, H. Ouedraogo, "A mathematical analysis of prey–predator population dynamics in the presence of an SIS infectious disease," Research in Mathematics, vol. 9 no. 1,DOI: 10.1080/27658449.2021.2020399, 2022.

[9] A. Savadogo, H. Ouedraogo, B. Sangaré, W. Ouedraogo, "Mathematical analysis of a fish-plankton eco-epidemiological system," Nonlinear Studies, vol. 27 no. 1, 2020.

[10] S. Sarwardi, M. Haque, E. Venturino, "A Leslie-Gower Holling-type II ecoepidemic model," Journal of Applied Mathematics and Computing, vol. 35 no. 1-2, pp. 263-280, DOI: 10.1007/s12190-009-0355-1, 2011.

[11] S. Kant, V. Kumar, "Stability analysis of predator–prey system with migrating prey and disease infection in both species," Applied Mathematical Modelling, vol. 42, pp. 509-539, DOI: 10.1016/j.apm.2016.10.003, 2017.

[12] B. Traoré, O. Koutou, B. Sangaré, "A global mathematical model of malaria transmission dynamics with structured mosquito population and temperature variations," Nonlinear Analysis: Real World Applications, vol. 53,DOI: 10.1016/j.nonrwa.2019.103081, 2020.

[13] T. Chowdhury, "Effect of fast migration on prey disease in a patchy system," Turcomat, vol. 11 no. 3, pp. 964-977, 2020.

[14] J. J. Tewa, V. Y. Djeumen, S. Bowong, "Predator–prey model with Holling response function of type II and SIS infectious disease," Applied Mathematical Modelling, vol. 37 no. 7, pp. 4825-4841, DOI: 10.1016/j.apm.2012.10.003, 2013.

[15] A. S. Abdulghafour, R. K. Naji, "A study of a diseased prey–predator model with refuge in prey and harvesting from predator," Journal of Applied Mathematics, vol. 2018,DOI: 10.1155/2018/2952791, 2018.

[16] J. Chattopadhyay, O. Arino, "A predator–prey model with disease in the prey," Nonlinear Analysis: Theory, Methods & Applications, vol. 36 no. 6, pp. 747-766, DOI: 10.1016/S0362-546X(98)00126-6, 1999.

[17] 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.

[18] H. W. Hethcote, W. Wang, L. Han, Z. Ma, "A predator–prey model with infected prey," Theoretical Population Biology, vol. 66 no. 3, pp. 259-268, DOI: 10.1016/j.tpb.2004.06.010, 2004.

[19] R. Bhattacharyya, B. Mukhopadhyay, "Spatial dynamics of nonlinear prey–predator models with prey migration and predator switching," Ecological Complexity, vol. 3 no. 2, pp. 160-169, DOI: 10.1016/j.ecocom.2006.01.001, 2006.

[20] R. Kiran, M. Roy, S. Abbas, A. Taraphder, "Effect of population migration and punctuated lockdown on the spread of infectious disease," Nonautonomous Dynamical Systems, vol. 8 no. 1, pp. 251-266, DOI: 10.1515/msds-2020-0137, 2021.

[21] G.-Q. Sun, Z. Jin, Q.-X. Liu, L. Li, "Dynamical complexity of a spatial predator–prey model with migration," Ecological Modelling, vol. 219 no. 1-2, pp. 248-255, DOI: 10.1016/j.ecolmodel.2008.08.009, 2008.

[22] L. Q. Gao, H. W. Hethcote, "Disease transmission models with density-dependent demographics," Journal of Mathematical Biology, vol. 30 no. 7, pp. 717-731, DOI: 10.1007/BF00173265, 1992.

[23] H. A. Antosiewicz, "Ordinary differential equations (G. Birkhoff and G. C. Rota)," SIAM Review, vol. 5 no. 2, pp. 160-161, DOI: 10.1137/1005043, 1963.

[24] Y.-H. Hsieh, C.-K. Hsiao, "Predator-prey model with disease infection in both populations," Mathematical Medicine and Biology, vol. 25 no. 3, pp. 247-266, DOI: 10.1093/imammb/dqn017, 2008.

[25] X. Wang, "A simple proof of descartes’s rule of signs," The American Mathematical Monthly, vol. 111 no. 6, pp. 525-526, DOI: 10.1080/00029890.2004.11920108, 2004.

[26] A. B. Diabaté, B. Sangaré, O. Koutou, "Mathematical modeling of the dynamics of vector-borne diseases transmitted by mosquitoes: taking into account aquatic stages and gonotrophic cycle," Nonautonomous Dynamical Systems, vol. 9 no. 1, pp. 205-236, DOI: 10.1515/msds-2022-0155, 2022.

[27] D. Greenhalgh, "Hopf bifurcation in epidemic models with a latent period and nonpermanent immunity," Mathematical and Computer Modelling, vol. 25 no. 2, pp. 85-107, DOI: 10.1016/S0895-7177(97)00009-5, 1997.

[28] H. Ouedraogo, W. Ouedraogo, B. Sangaré, "Bifurcation and stability analysis in complex cross-diffusion mathematical model of phytoplankton-fish dynamics," Journal of Partial Differential Equations, vol. 32 no. 3, pp. 207-228, DOI: 10.4208/jpde, 2019.

[29] C.-H. Chiu, "Lyapunov functions for the global stability of competing predators," Journal of Mathematical Analysis and Applications, vol. 230 no. 1, pp. 232-241, DOI: 10.1006/jmaa.1998.6198, 1999.

[30] O. Koutou, B. Traoré, B. Sangaré, "Analysis of schistosomiasis global dynamics with general incidence functions and two delays," International Journal of Applied and Computational Mathematics, vol. 7,DOI: 10.1007/s40819-021-01188-y, 2021.

[31] A. Savadogo, B. Sangaré, W. Ouedraogo, "A Complex Dynamic of Eco-Epidemiological Mathematical Model with Migration," 2023. Preprint

[32] S. Sarwardi, M. Haque, E. Venturino, "Global stability and persistence in LG–Holling type II diseased predator ecosystems," Journal of Biological Physics, vol. 37 no. 1, pp. 91-106, DOI: 10.1007/s10867-010-9201-9, 2011.

[33] M. Al-Refai, "On weighted Atangana–Baleanu fractional operators," Advances in Difference Equations, vol. 2020,DOI: 10.1186/s13662-019-2471-z, 2020.

[34] Q. Richard, M. Choisy, T. Lefèvre, R. Djidjou-Demasse, "Human-vector malaria transmission model structured by age, time since infection and waning immunity," Nonlinear Analysis: Real World Applications, vol. 63,DOI: 10.1016/j.nonrwa.2021.103393, 2022.