Dynamics of a Fractional-Order Eco-Epidemiological Model with Two Disease Strains in a Predator Population Incorporating Harvesting

Abstract

In this paper, a fractional-order eco-epidemiological model with two disease strains in the predator population incorporating harvesting is formulated and analyzed. The model assumes that the population is divided into a prey population, a susceptible predator population, a predator population infected by the first disease, and a predator population infected by the second disease. A mathematical analysis and numerical simulations are performed to explain the dynamics and properties of the proposed fractional-order eco-epidemiological model. The positivity, boundedness, existence, and uniqueness of the solutions are examined. The basic reproduction number and some sufficient conditions for the existence of four equilibrium points are obtained. In addition, some sufficient conditions are proposed to ensure the local and global asymptotic stability of the equilibrium points. Theoretical results are illustrated through numerical simulations, which also highlight the effect of the fractional order.

1. Introduction

The relationship between predators and their prey is a fundamental topic in mathematical ecology due to its widespread occurrence and ecological significance [1]. The interactions between prey and predators were studied for the first time by the famous mathematicians Lotka and Volterra [2]. After that, many predator–prey models have been established and studied by mathematicians and ecologists, for example [3,4,5,6,7].

Investigating the spread of infections within populations is a crucial area of mathematical biology, offering insights into predicting the impacts of such infections [8,9,10]. An eco-epidemiological model studies the dynamics of predator–prey interactions in the context of infectious diseases, which may affect either the prey population only [11,12], the predator population only [13,14,15], or both simultaneously [16,17].

Fractional-order differential equations are a generalized form of classical ordinary differential equations, extending their order to non-integer values [18]. Fractional-order models offer advantages over integer-order models, including memory effects and hereditary dynamics, which better capture complex system behaviors [19]. Models based on fractional-order differential equations may offer a more accurate representation of complex systems and elucidate the interactions between prey and predator species, particularly in the context of infectious diseases affecting the predator population [20]. Fractional-order derivatives have found widespread application across various scientific and engineering disciplines [21]. For a more comprehensive understanding of fractional-order differential equations, one can refer to [22,23,24,25,26,27,28] and the references contained therein. These references explore the application of the fractional order in ecological, epidemiological, and biological–economic models, emphasizing the analysis of stability, bifurcation, and memory effects to understand complex dynamic systems.

In this paper, a fractional-order eco-epidemiological model incorporating two disease strains within the predator population and the effects of harvesting is proposed and studied. The population is assumed to consist of prey, susceptible predators, predators infected by the first disease, and predators infected by the second disease. For instance, the black-footed ferret relies exclusively on Prairie dogs as its primary food source. This black-footed ferret population can be infected by the Sylvatic plague and the Canine distemper virus [29]. The positivity, boundedness, existence, and uniqueness of the solutions for the fractional-order model are examined. Additionally, the basic reproduction number is derived, along with sufficient conditions for the existence of four equilibrium points. The main contribution of this paper is establishing sufficient conditions to guarantee the local and global asymptotic stability of the equilibrium points in the proposed model. The theoretical findings are further illustrated through numerical simulations.

The structure of this paper is as follows. In the next section, the model formulation, positivity, boundedness, existence, and uniqueness are proposed. In Section 3, the equilibrium points, basic reproduction number, local stability, and global stability of the proposed fractional-order model are investigated. Section 4 presents numerical simulations to illustrate the theoretical results obtained. Finally, the conclusions are provided in Section 5.

2. Model Formulation

Following [30], the eco-epidemiological model incorporating two disease strains affecting the predator population can be described as follows.

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\displaystyle \frac{dx}{dt}}& =\widehat{r}\left(1-{\displaystyle \frac{x}{\widehat{k}}}\right)x-\widehat{a}xy,x\left(0\right)=x_0,\hfill \\ \hfill {\displaystyle \frac{dy}{dt}}& =\widehat{e}\widehat{a}xy-\widehat{\lambda }yz-\widehat{\beta }yw+\widehat{\gamma }z+\widehat{\phi }w-\widehat{m}y,y\left(0\right)=y_0,\hfill \\ \hfill {\displaystyle \frac{dz}{dt}}& =\widehat{\lambda }yz-\widehat{\gamma }z-\widehat{d}z+\widehat{\theta }w,z\left(0\right)=z_0,\hfill \\ \hfill {\displaystyle \frac{dw}{dt}}& =\widehat{\beta }yw-\widehat{\phi }w-\widehat{\nu }w-\widehat{\theta }w,w\left(0\right)=w_0.\hfill \end{array}\end{array}$$

(1)

The model (1) categorizes the populations into four classes: the prey population $x\left(t\right)$, the susceptible predator population $y\left(t\right)$, the predator population infected with the first disease $z\left(t\right)$, and the predator population infected with the second disease $w\left(t\right)$. It is assumed that disease transmission occurs within the predator populations, while the susceptible predators feed on the prey. All parameters in model (1) are non-negative for $t\ge 0$ and are detailed in

Table 1. Parameter descriptions.

Symbol	Description
$\widehat{r}$	Intrinsic growth rate of prey
$\widehat{k}$	Prey carrying capacity
$\widehat{a}$	Predation rate of susceptible predator
$\widehat{H}$	Prey harvesting
$\widehat{e}$	susceptible predator conversion efficiency
$\widehat{\lambda }$	Transmission coefficient of the first disease in predator
$\widehat{\beta }$	Transmission coefficient of the second disease in predator
$\widehat{m}$	Natural mortality rate of susceptible predator
$\widehat{\gamma }$	First disease recovery rate
$\widehat{\phi }$	Second disease recovery rate
$\widehat{d}$	Natural plus first disease mortality rate
$\widehat{\nu }$	Natural plus second disease mortality rate
$\widehat{\theta }$	Mutation factor of diseases.

.

This paper seeks to investigate the dynamic properties of a generalization of the eco-epidemiological model described in (1) through the introduction of the Caputo fractional derivative of order q$\left({}^{c}D^q\right)$ and prey harvesting $\left(\widehat{H}\right)$ as follows.

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {}^{c}D^{q}x\left(t\right)& =\widehat{r}\left(1-{\displaystyle \frac{x}{\widehat{k}}}\right)x-\widehat{a}xy-\widehat{H}x,x\left(0\right)=x_0,\hfill \\ \hfill {}^{c}D^{q}y\left(t\right)& =\widehat{e}\widehat{a}xy-\widehat{\lambda }yz-\widehat{\beta }yw+\widehat{\gamma }z+\widehat{\phi }w-\widehat{m}y,y\left(0\right)=y_0,\hfill \\ \hfill {}^{c}D^{q}z\left(t\right)& =\widehat{\lambda }yz-\widehat{\gamma }z-\widehat{d}z+\widehat{\theta }w,z\left(0\right)=z_0,\hfill \\ \hfill {}^{c}D^{q}w\left(t\right)& =\widehat{\beta }yw-\widehat{\phi }w-\widehat{\nu }w-\widehat{\theta }w,w\left(0\right)=w_0,\hfill \end{array}\end{array}$$

(2)

For $q\in (0,1)$, the Caputo fractional derivative ${}^{c}D^q$ is employed [18]. In model (2), the right-hand-side terms have a dimension of ${\left(time\right)}^{-1}$, while the left-hand-side terms have a dimension of ${\left(time\right)}^{-q}$. To ensure dimensional consistency, model (2) is reformulated as follows:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {}^{c}D^{q}x\left(t\right)& ={\widehat{r}}^q\left(1-{\displaystyle \frac{x}{\widehat{k}}}\right)x-{\widehat{a}}^{q}xy-{\widehat{H}}^{q}x,x\left(0\right)=x_0,\hfill \\ \hfill {}^{c}D^{q}y\left(t\right)& =\widehat{e}{\widehat{a}}^{q}xy-{\widehat{\lambda }}^{q}yz-{\widehat{\beta }}^{q}yw+{\widehat{\gamma }}^{q}z+{\widehat{\phi }}^{q}w-{\widehat{m}}^{q}y,y\left(0\right)=y_0,\hfill \\ \hfill {}^{c}D^{q}z\left(t\right)& ={\widehat{\lambda }}^{q}yz-{\widehat{\gamma }}^{q}z-{\widehat{d}}^{q}z+{\widehat{\theta }}^{q}w,z\left(0\right)=z_0,\hfill \\ \hfill {}^{c}D^{q}w\left(t\right)& ={\widehat{\beta }}^{q}yw-{\widehat{\phi }}^{q}w-{\widehat{\nu }}^{q}w-{\widehat{\theta }}^{q}w,w\left(0\right)=w_0.\hfill \end{array}\end{array}$$

(3)

For simplicity, model (3) is redefined using new parameter representations [31]:

$$\begin{array}{cc}\hfill & {\widehat{r}}^q=r,\widehat{k}=k,{\widehat{a}}^q=a,{\widehat{H}}^q=H,\widehat{e}=e,{\widehat{\lambda }}^q=\lambda ,{\widehat{\beta }}^q=\beta ,{\widehat{\gamma }}^q=\gamma ,{\widehat{\phi }}^q=\phi ,{\widehat{m}}^q=m,{\widehat{d}}^q=d,\hfill \\ \hfill & {\widehat{\theta }}^q=\theta ,{\widehat{\nu }}^q=\nu .\hfill \end{array}$$

Then, the model (3) can be reformulated as follows:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {}^{c}D^{q}x\left(t\right)& =r\left(1-{\displaystyle \frac{x}{k}}\right)x-axy-Hx,x\left(0\right)=x_0,\hfill \\ \hfill {}^{c}D^{q}y\left(t\right)& =eaxy-\lambda yz-\beta yw+\gamma z+\phi w-my,y\left(0\right)=y_0,\hfill \\ \hfill {}^{c}D^{q}z\left(t\right)& =\lambda yz-\gamma z-dz+\theta w,z\left(0\right)=z_0,\hfill \\ \hfill {}^{c}D^{q}w\left(t\right)& =\beta yw-\phi w-\nu w-\theta w,w\left(0\right)=w_0,\hfill \end{array}\end{array}$$

(4)

It is to be noted that the integer-order model (1) given in [30] cannot be sustained at a stable coexistence equilibrium level. However, the fractional-order model (4) proposed in this paper can be sustained at the stable coexistence equilibrium level. To the best of our knowledge, no prior studies have explored the dynamics of a fractional-order eco-epidemiological model with two disease strains in the predator population that incorporates harvesting (4).

2.1. Positivity and Boundedness

This subsection investigates the positivity and boundedness of the solutions for the fractional-order eco-epidemiological model (4). The positivity of the solution of model (4) with positive initial conditions is now investigated. Following model (4), one has

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {}^{c}D^q{x\left(t\right)|}_{x=0}& =0,\hfill \\ \hfill {}^{c}D^q{y\left(t\right)|}_{y=0}& =\gamma z+\phi w\ge 0,\hfill \\ \hfill {}^{c}D^q{z\left(t\right)|}_{z=0}& =\theta w\ge 0,\hfill \\ \hfill {}^{c}D^q{w\left(t\right)|}_{w=0}& =0.\hfill \end{array}\end{array}$$

Furthermore, the model satisfies the Lipschitz condition, as established in Theorem 2. Based on the positivity property, Theorem 5 and Theorem 6 of [32], the solutions of the fractional-order model (4) remain non-negative for $t\ge 0$.

The boundedness of the solutions for model (4) is established in the following theorem:

Theorem 1. 

All the solutions of model (4) starting in ${\mathbb{R}}_{+}^4$ are uniformly bounded.

Proof. 

Let $\varphi \left(t\right)=x\left(t\right)+y\left(t\right)+z\left(t\right)+w\left(t\right)$; then,

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {}^{c}D^q\varphi \left(t\right)& ={}^{c}D^{q}x\left(t\right)+{}^{c}D^{q}y\left(t\right)+{}^{c}D^{q}z\left(t\right)+{}^{c}D^{q}w\left(t\right)\hfill \\ \hfill & =r\left(1-{\displaystyle \frac{x}{k}}\right)x+(e-1)axy-Hx-my-dz-\nu w\hfill \\ \hfill & \le r\left(1-{\displaystyle \frac{x}{k}}\right)x-Hx-my-dz-\nu w\hfill \\ \hfill & \le -{\displaystyle \frac{r{x}^2}{k}}+rx-\sigma (x+y+z+w),\hfill \end{array}\end{array}$$

where $\sigma =min\left\{H,m,d,\nu \right\}$; thus,

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {}^{c}D^q\varphi \left(t\right)+\sigma \varphi \left(t\right)& \le -{\displaystyle \frac{r{x}^2}{k}}+rx\hfill \\ \hfill & \le -\frac{r}{k}{\left(x-\frac{k}{2}\right)}^2+{\displaystyle \frac{rk}{4}}\hfill \\ \hfill & \le {\displaystyle \frac{rk}{4}}.\hfill \end{array}\end{array}$$

By using the Lemma 9 in [33], then,

$$0\le \varphi \left(t\right)\le \varphi \left(0\right)E_q(-\sigma t^q)+{\displaystyle \frac{rk}{4}}{t}^q{E}_{q,q+1}(-\sigma t^q),$$

Here, $E_q$ denotes the Mittag–Leffler function. Using Lemma 5 and Corollary 6 from [33], it is derived that

$$0\le \varphi \left(t\right)\le {\displaystyle \frac{rk}{4\sigma }},\mathrm{as}t\to \infty .$$

As a result, all solutions of model (4) with initial conditions in ${\mathbb{R}}_{+}^4$ are uniformly bounded within the region S, where

$$S=\left\{(x,y,z,w)\in R_{+}^4:\varphi \left(t\right)\le {\displaystyle \frac{rk}{4\sigma }}+ϵ,ϵ>0\right\}.$$

(5)

□

2.2. Existence and Uniqueness

The existence and uniqueness of solutions for the fractional-order model (4) within the region $M\times \left(0,T\right],$ where where

$$M=\left\{(x,y,z,w)\in {\mathbb{R}}^4:max\left(\left|x\right|,\left|y\right|,\left|z\right|,\left|w\right|\right)\le h\right\},$$

are investigated as follows:

Theorem 2. 

For each $X_0=(x_0,y_0,z_0,w_0)\in M$, there exists a unique solution $X\left(t\right)\in M$ of model (4) with the initial condition $X_0$, which is defined for all $t\ge 0$.

Proof. 

Consider a mapping $L\left(X\right)=(L_1\left(X\right),L_2\left(X\right),L_3\left(X\right),L_4\left(X\right))$, where

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill L_1\left(X\right)& =r\left(1-{\displaystyle \frac{x}{k}}\right)x-axy-Hx,\hfill \\ \hfill L_2\left(X\right)& =eaxy-\lambda yz-\beta yw+\gamma z+\phi w-my,\hfill \\ \hfill L_3\left(X\right)& =\lambda yz-\eta z+\theta w,\hfill \\ \hfill L_4\left(X\right)& =\beta yw-\xi w,\hfill \end{array}\end{array}$$

(6)

For any $X,\overline{X}\in M,$ it follows from (6) that

$$\begin{array}{cc}\hfill ∥L\left(X\right)-L\left(\overline{X}\right)∥=& \left|L_1\left(X\right)-L_1\left(\overline{X}\right)\right|+\left|L_2\left(X\right)-L_2\left(\overline{X}\right)\right|+\left|L_3\left(X\right)-L_3\left(\overline{X}\right)\right|+\left|L_4\left(X\right)-L_4\left(\overline{X}\right)\right|\hfill \\ \hfill =& \left|r\left(1-{\displaystyle \frac{x}{k}}\right)x-axy-Hx-r\left(1-{\displaystyle \frac{\overline{x}}{k}}\right)\overline{x}+a\overline{x}\overline{y}+H\overline{x}\right|\hfill \\ \hfill & +\left|eaxy-\lambda yz-\beta yw+\gamma z+\phi w-my-ea\overline{x}\overline{y}+\lambda \overline{y}\overline{z}+\beta \overline{y}\overline{w}-\gamma \overline{z}-\phi \overline{w}+m\overline{y}\right|\hfill \\ \hfill & +\left|\lambda yz-\eta z+\theta w-\lambda \overline{y}\overline{z}+\eta \overline{z}-\theta \overline{w}\right|+\left|\beta yw-\xi w-\beta \overline{y}\overline{w}+\xi \overline{w}\right|\hfill \\ \hfill \le & r|x-\overline{x}|+\frac{r}{k}|x-\overline{x}\left|\right|x+\overline{x}|+a(1+e)|xy-\overline{x}y+\overline{x}y-\overline{x}\overline{y}|\hfill \\ \hfill & +H|x-\overline{x}|+2\lambda |yz-\overline{y}z+\overline{y}z-\overline{y}\overline{z}|\hfill \\ \hfill & +2\beta |yw-\overline{y}w+\overline{y}w-\overline{y}\overline{w}|+\gamma |z-\overline{z}|+\phi |w-\overline{w}|\hfill \\ \hfill & +m|y-\overline{y}|+\eta |z-\overline{z}|+\theta |w-\overline{w}|+\xi |w-\overline{w}|\hfill \\ \hfill \le & r|x-\overline{x}|+\frac{2rh}{k}|x-\overline{x}|+a(1+e)h|x-\overline{x}|\hfill \\ \hfill & +a(1+e)h|y-\overline{y}|+H|x-\overline{x}|+2\lambda h|y-\overline{y}|+2\lambda h|z-\overline{z}|\hfill \\ \hfill & +2\beta h|y-\overline{y}|+2\beta h|w-\overline{w}|+\gamma |z-\overline{z}|+\phi |w-\overline{w}|\hfill \\ \hfill & +m|y-\overline{y}|+\eta |z-\overline{z}|+\theta |w-\overline{w}|+\xi |w-\overline{w}|\hfill \\ \hfill \le & \left(r+{\displaystyle \frac{2rh}{k}}+a(1+e)h+H\right)|x-\overline{x}|\hfill \\ \hfill & +\left(a(1+e)h+2\lambda h+2\beta h+m\right)|y-\overline{y}|\hfill \\ \hfill & +\left(2\lambda h+\gamma +\eta \right)|z-\overline{z}|\hfill \\ \hfill & +\left(2\beta h+\phi +\theta +\xi \right)|w-\overline{w}|\hfill \\ \hfill \le & U∥X-\overline{X}∥,\hfill \end{array}$$

where

$$U=max\left\{r+{\displaystyle \frac{2rh}{k}}+a(1+e)h+H,a(1+e)h+2\lambda h+2\beta h+m,2\lambda h+\gamma +\eta ,2\beta h+\phi +\theta +\xi \right\}.$$

Thus, $L\left(X\right)$ satisfies the Lipschitz condition, proving the existence and uniqueness of solutions for model (4) under the given initial conditions. □

3. Model Analysis

This section investigates the equilibrium points, basic reproduction number, local stability, and global stability of the fractional-order eco-epidemiological model (4).

3.1. Equilibrium Points

This subsection and the next will utilize the basic reproduction number $\left({\Re }_0\right)$ of model (4) to determine the existence and stability of its equilibrium points. The basic reproduction number $\left({\Re }_0\right)$ can be obtained by using the next-generation method [34]. One can rewrite the fractional-order model (4) as follows

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {}^{c}D^{q}w\left(t\right)& =\beta yw-\xi w,\hfill \\ \hfill {}^{c}D^{q}z\left(t\right)& =\lambda yz-\eta z+\theta w,\hfill \\ \hfill {}^{c}D^{q}y\left(t\right)& =eaxy-\lambda yz-\beta yw+\gamma z+\phi w-my,\hfill \\ \hfill {}^{c}D^{q}x\left(t\right)& =r\left(1-{\displaystyle \frac{x}{k}}\right)x-axy-Hx,\hfill \end{array}\end{array}$$

(7)

where $\xi =\phi +\nu +\theta$ and $\eta =\gamma +d$. The model (7) can subsequently be expressed as follows:

$$D^{q}X\left(t\right)=f\left(X\right)-v\left(X\right),$$

where

$$f\left(X\right)=\left[\begin{array}{c}{f}_1\\ f_2\\ f_3\\ f_4\end{array}\right]=\left[\begin{array}{c}\beta yw\\ \lambda yz\\ eaxy\\ 0\end{array}\right],v\left(X\right)=\left[\begin{array}{c}{v}_1\\ v_2\\ v_3\\ v_4\end{array}\right]=\left[\begin{array}{c}\xi w\\ \eta z-\theta w\\ \lambda yz+\beta yw-\gamma z-\phi w+my\\ -r\left(1-{\displaystyle \frac{x}{k}}\right)x+axy+Hx\end{array}\right].$$

The matrices $F\left(X\right)$ and $V\left(X\right)$ are defined as

$$F\left(X\right)=\left[\begin{array}{cccc}{\displaystyle \frac{\partial f_1}{\partial w}}& {\displaystyle \frac{\partial f_1}{\partial z}}& {\displaystyle \frac{\partial f_1}{\partial y}}& {\displaystyle \frac{\partial f_1}{\partial x}}\\ {\displaystyle \frac{\partial f_2}{\partial w}}& {\displaystyle \frac{\partial f_2}{\partial z}}& {\displaystyle \frac{\partial f_2}{\partial y}}& {\displaystyle \frac{\partial f_2}{\partial x}}\\ {\displaystyle \frac{\partial f_3}{\partial w}}& {\displaystyle \frac{\partial f_3}{\partial z}}& {\displaystyle \frac{\partial f_3}{\partial y}}& {\displaystyle \frac{\partial f_3}{\partial x}}\\ {\displaystyle \frac{\partial f_4}{\partial w}}& {\displaystyle \frac{\partial f_4}{\partial z}}& {\displaystyle \frac{\partial f_4}{\partial y}}& {\displaystyle \frac{\partial f_4}{\partial x}}\end{array}\right],V\left(X\right)=\left[\begin{array}{cccc}{\displaystyle \frac{\partial v_1}{\partial w}}& {\displaystyle \frac{\partial v_1}{\partial z}}& {\displaystyle \frac{\partial v_1}{\partial y}}& {\displaystyle \frac{\partial v_1}{\partial x}}\\ {\displaystyle \frac{\partial v_2}{\partial w}}& {\displaystyle \frac{\partial v_2}{\partial z}}& {\displaystyle \frac{\partial v_2}{\partial y}}& {\displaystyle \frac{\partial v_2}{\partial x}}\\ {\displaystyle \frac{\partial v_3}{\partial w}}& {\displaystyle \frac{\partial v_3}{\partial z}}& {\displaystyle \frac{\partial v_3}{\partial y}}& {\displaystyle \frac{\partial v_3}{\partial x}}\\ {\displaystyle \frac{\partial v_4}{\partial w}}& {\displaystyle \frac{\partial v_4}{\partial z}}& {\displaystyle \frac{\partial v_4}{\partial y}}& {\displaystyle \frac{\partial v_4}{\partial x}}\end{array}\right].$$

Thus,

$$F\left(X\right)=\left[\begin{array}{cccc}\beta y& 0& \beta w& 0\\ 0& \lambda y& \lambda z& 0\\ 0& 0& aex& aey\\ 0& 0& 0& 0\end{array}\right],$$

$$V\left(X\right)=\left[\begin{array}{cccc}\xi & 0& 0& 0\\ -\theta & \eta & 0& 0\\ \beta y-\phi & \lambda y-\gamma & \beta w+\lambda z+m& 0\\ 0& 0& ax& r\left({\displaystyle \frac{2x}{k}}-1\right)+ay+H\end{array}\right].$$

To obtain the eigenvalues of $F·V^{-1},$ at equilibrium point $E_1\left({\displaystyle \frac{k(r-H)}{r}},0,0,0\right),$ the equation

$$\left|F·V^{-1}-\mu I\right|=0,$$

can be solved. Based on Theorem 2 in [34], the basic reproduction number of Model 2 is given as ${\Re }_0={\displaystyle \frac{ae(r-H)k}{rm}}.$ Additionally, the threshold parameters will be utilized to establish the conditions for the existence and stability of the equilibrium points of model (4):

$${\Re }_2={\displaystyle \frac{\beta \left(ae\right(r-H)k-rm)}{e{a}^{2}k\xi }},{\Re }_{22}={\displaystyle \frac{\lambda \left(ae\right(r-H)k-rm)}{e{a}^{2}k\eta }},{\Re }_3={\displaystyle \frac{\beta \eta }{\xi \lambda }}.$$

It is to be noted that the basic reproduction number $\left({\Re }_0\right)$ and the threshold parameters (${\Re }_2$ and ${\Re }_{22}$) depend on the prey harvesting $\left(H\right)$. This means that the prey harvesting $\left(H\right)$ has a crucial effect on the existence and stability conditions of equilibrium points of model (4).

The fractional-order eco-epidemiological model (4) has four equilibrium points:

• $E_0=(0,0,0,0)$, which always exists.
• $E_1=\left(\frac{k}{r}(r-H),0,0,0\right),$ which exists if $r>H$.
• $E_2=(x_2,y_2,0,0)$ where

  $$x_2={\displaystyle \frac{m}{ae}},y_2={\displaystyle \frac{ae(r-H)k-rm}{e{a}^{2}k}}={\displaystyle \frac{rm}{e{a}^{2}k}}({\Re }_0-1).$$
  Therefore, $E_2$ exists if ${\Re }_0>1.$
• $E_3=(x_3,y_3,z_3,0)$ where

  $$x_3={\displaystyle \frac{k(r-H)}{r}}-{\displaystyle \frac{a\eta k}{r\lambda }}=x_1-{\displaystyle \frac{a\eta k}{r\lambda }},y_3={\displaystyle \frac{\eta }{\lambda }},z_3={\displaystyle \frac{({\Re }_{22}-1)ek{a}^2{\eta }^2}{r(\eta -\gamma ){\lambda }^2}}.$$
  Then, $E_3$ exists if $x_1>{\displaystyle \frac{a\eta k}{r\lambda }}$ and ${\Re }_{22}>1$.
• $E_4=(x_4,y_4,z_4,w_4)$ where

  $$x_4={\displaystyle \frac{k\left(\beta \right(r-H)-a\xi )}{r\beta }},y_4={\displaystyle \frac{\xi }{\beta }},z_4={\displaystyle \frac{e{a}^{2}k{\xi }^2({\Re }_2-1)}{C_1}},w_4={\displaystyle \frac{e{a}^{2}k\lambda {\xi }^3}{\beta C_1}}({\Re }_2-1)({\Re }_3-1),$$

  where $C_1=r\beta (\lambda \xi (\theta +\phi -\xi )+\beta (\xi \eta -\gamma \theta -\eta \phi ))$. Therefore $E_4$ exists if $\beta >{\displaystyle \frac{a\xi }{r-H}},$ ${\Re }_2>1,$ ${\Re }_3>1$ and $C_1>0.$

3.2. Local Stability Analysis

In the following, the asymptotic stability of equilibrium points of model (4) is studied. The Jacobian matrix $\left(J\right(x,y,z,w\left)\right)$ of model (4) is as follows:

$$J(x,y,z,w)=\left(\begin{array}{cccc}r-{\displaystyle \frac{2rx}{k}}-ay-H& -ax& 0& 0\\ aey& aex-\lambda z-\beta w-m& \gamma -\lambda y& \phi -\beta y\\ 0& \lambda z& \lambda y-\eta & \theta \\ 0& \beta w& 0& \beta y-\xi \end{array}\right).$$

(8)

The stability analysis of the equilibrium point $E_0$ is not considered, as this state signifies the extinction of all populations. The $E_0$ is unstable.

Lemma 1. 

If ${\Re }_0<1$, then $E_1$ is locally asymptotically stable.

Proof. 

The $J\left(E_1\right)$ is

$$J\left(E_1\right)=\left(\begin{array}{cccc}H-r& {\displaystyle \frac{a(H-r)k}{r}}& 0& 0\\ 0& ({\Re }_0-1)m& \gamma & \phi \\ 0& 0& -\eta & \theta \\ 0& 0& 0& -\xi \end{array}\right).$$

(9)

The eigenvalues of $J\left(E_1\right)$ are ${\mu }_1=H-r,{\mu }_2=({\Re }_0-1)m,{\mu }_3=-\xi$ and ${\mu }_4=-\eta$. Thus $|arg\left({\mu }_{1,3,4}\right)|=\pi >{\displaystyle \frac{q\pi }{2}}$. If ${\Re }_0<1$, then $|arg\left({\mu }_2\right)|=\pi >{\displaystyle \frac{q\pi }{2}}$ for all $q\in (0,1)$. As demonstrated in [35,36], the proof is thus complete. □

Lemma 2. 

If $y_2<min\left\{{\displaystyle \frac{\xi }{\beta }},{\displaystyle \frac{\eta }{\lambda }}\right\}$, then $E_2$ is locally asymptotically stable.

Proof. 

The $J\left(E_2\right)$ is

$$J\left(E_2\right)=\left(\begin{array}{cccc}-{\displaystyle \frac{r{x}_2}{k}}& -a{x}_2& 0& 0\\ ae{y}_2& 0& \gamma -\lambda y_2& \phi -\beta y_2\\ 0& 0& \lambda y_2-\eta & \theta \\ 0& 0& 0& \beta y_2-\xi \end{array}\right).$$

(10)

The eigenvalues of $J\left(E_2\right)$ are ${\mu }_1=\beta y_2-\xi ,{\mu }_2=\lambda y_2-\eta$, and ${\mu }_{3,4}$ are the solutions of:

$${\mu }^2+\frac{r{x}_2}{k}\mu +e{a}^2{x}_2{y}_2=0,$$

(11)

since ${\displaystyle \frac{r{x}_2}{k}}>0$ and $e{a}^2{x}_2{y}_2>0$, the eigenvalues of Equation (11), exhibit negative real parts. If $y_2<min\left\{{\displaystyle \frac{\xi }{\beta }},{\displaystyle \frac{\eta }{\lambda }}\right\}$, then $|arg\left({\mu }_{1,2}\right)|=\pi >{\displaystyle \frac{q\pi }{2}}$ for all $q\in (0,1)$. As demonstrated in [35,36], the proof is thus complete. □

Lemma 3. 

If ${\displaystyle \frac{\gamma }{\lambda }}<y_3<{\displaystyle \frac{\xi }{\beta }}$, then $E_3$ is locally asymptotically stable.

Proof. 

The $J\left(E_3\right)$ is

$$J\left(E_3\right)=\left(\begin{array}{cccc}-{\displaystyle \frac{r{x}_3}{k}}& -a{x}_3& 0& 0\\ ae{y}_3& -{\displaystyle \frac{\gamma z_3}{y_3}}& \gamma -\lambda y_3& \phi -\beta y_3\\ 0& \lambda z_3& 0& \theta \\ 0& 0& 0& \beta y_3-\xi \end{array}\right).$$

The eigenvalues of $J\left(E_3\right)$ are ${\mu }_1=\beta y_3-\xi ,$ while the other three eigenvalues ${\mu }_{2,3,4}$ are the solutions of

$${\mu }^3+B_1{\mu }^2+B_2\mu +B_3=0,$$

(12)

where

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill B_1& ={\displaystyle \frac{r{x}_3}{k}}+{\displaystyle \frac{\gamma z_3}{y_3}},\hfill \\ \hfill B_2& =\left({\displaystyle \frac{r\gamma z_3}{k{y}_3}}+e{a}^2{y}_3\right)x_3+\lambda (\lambda y_3-\gamma )z_3,\hfill \\ \hfill B_3& =\frac{1}{k}\left(r\lambda (\lambda y_3-\gamma )x_3{z}_3\right).\hfill \end{array}\end{array}$$

It is obvious that $B_1>0,B_2>0,B_3>0$ and $B_1{B}_2>B_3$ as long as $\lambda y_3>\gamma .$ By applying the Routh–Hurwitz criterion, it is established that all solutions of Equation (12) have negative real parts. Consequently, the equilibrium point $E_3$ is locally asymptotically stable when ${\displaystyle \frac{\gamma }{\lambda }}<y_3<{\displaystyle \frac{\xi }{\beta }}$. □

The stability of the equilibrium point $E_4$ is now investigated. The $J\left(E_4\right)$ is

$$J\left(E_4\right)=\left(\begin{array}{cccc}-{\displaystyle \frac{r{x}_4}{k}}& -a{x}_4& 0& 0\\ ae{y}_4& -{\displaystyle \frac{\gamma z_4+\phi w_4}{y_4}}& \gamma -\lambda y_4& \phi -\beta y_4\\ 0& \lambda z_4& -{\displaystyle \frac{\theta w_4}{z_4}}& \theta \\ 0& \beta w_4& 0& 0\end{array}\right).$$

The eigenvalues of $J\left(E_4\right)$ are the solutions of

$${\mu }^4+A_1{\mu }^3+A_2{\mu }^2+A_3\mu +A_4=0,$$

(13)

where

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill A_1& ={\displaystyle \frac{r{x}_4}{k}}+{\displaystyle \frac{C_2}{y_4}}+{\displaystyle \frac{\theta w_4}{z_4}},\hfill \\ \hfill A_2& ={\displaystyle \frac{\theta (k{C}_2+r{x}_4{y}_4)w_4}{k{y}_4{z}_4}}+{\displaystyle \frac{r{C}_2{x}_4}{k{y}_4}}+e{a}^2{x}_4{y}_4-\lambda C_3{z}_4-\beta C_4{w}_4,\hfill \\ \hfill A_3& ={\displaystyle \frac{1}{k{y}_4{z}_4}}(r\theta C_2{w}_4{x}_4+y_4(k\theta w_4(e{a}^2{x}_4{y}_4-\beta C_4{w}_4)-\beta w_4(k\theta C_3+r{C}_4{x}_4)z_4-r\lambda C_3{x}_4{z}_4^2)),\hfill \\ \hfill A_4& =-{\displaystyle \frac{r\beta \theta (C_3{z}_4+C_4{w}_4)x_4{w}_4}{k{z}_4}},\hfill \\ \hfill C_2& =\gamma z_4+\phi w_4,C_3=\gamma -\lambda y_4,C_4=\phi -\beta y_4.\hfill \end{array}\end{array}$$

The conditions for stability at $E_4$ can be derived using the proposition outlined in [37].

3.3. Global Stability Analysis

The global asymptotic stability of all four equilibrium points of the fractional-order model (4) is investigated as follows.

Theorem 3. 

The equilibrium point $E_1$ is globally asymptotically stable if ${\displaystyle \frac{ak(r-H)}{rm}}<1$.

Proof. 

A suitable positive definite Lyapunov function is considered as follows:

$$V=x-x_1-x_{1}ln\left({\displaystyle \frac{x}{x_1}}\right)+y+z+w.$$

By calculating the q-order derivative of V throughout the solution of (4) and applying Lemma 3.1 in [38],

$$\begin{array}{cc}\hfill {}^{c}D^{q}V\le & (x-x_1)\left(r-{\displaystyle \frac{rx}{k}}-ay-H\right)+eaxy-my-dz-\nu w\hfill \\ \hfill \le & (x-x_1)\left({\displaystyle \frac{r{x}_1}{k}}-{\displaystyle \frac{rx}{k}}-ay\right)+eaxy-my-dz-\nu w\hfill \\ \hfill \le & -\frac{r}{k}{(x-x_1)}^2+a(e-1)xy+(a{x}_1-m)y-dz-\nu w.\hfill \end{array}$$

Thus, ${}^{c}D^{q}V\le 0$ if ${\displaystyle \frac{a{x}_1}{m}}<1$ which is equivalent to ${\displaystyle \frac{ak(r-H)}{rm}}<1$. By applying Lemma 4.6 in [39], the equilibrium point $E_1$ is proven to be globally asymptotically stable. □

Theorem 4. 

The equilibrium point $E_2$ is globally asymptotically stable if $r\lambda k<4\sigma \gamma$.

Proof. 

A suitable positive definite Lyapunov function is considered as follows:

$$V=C_5\left(x-x_2-x_{2}ln\left({\displaystyle \frac{x}{x_2}}\right)\right)+y-y_2-y_{2}ln\left({\displaystyle \frac{y}{y_2}}\right)+z+w.$$

By calculating the q-order derivative of V throughout the solution of (4) and applying Lemma 3.1 in [38].

$$\begin{array}{cc}\hfill {}^{c}D^{q}V\le & C_5(x-x_2)\left(r-{\displaystyle \frac{rx}{k}}-ay-H\right)\hfill \\ \hfill & +(y-y_2)(aex-\lambda z-m)+\left({\displaystyle \frac{y-y_2}{y}}\right)\left(\gamma z-\beta yw+\phi w\right)\hfill \\ \hfill & +\lambda yz-\gamma z-dz+\theta w+\beta yw-\phi w-\nu w-\theta w\hfill \\ \hfill \le & -\frac{r{C}_5}{k}{(x-x_2)}^2-a{C}_5(x-x_2)(y-y_2)\hfill \\ \hfill & +ae(x-x_2)(y-y_2)+y_2\left(\lambda -{\displaystyle \frac{\gamma }{y}}\right)z\hfill \\ \hfill & +\left(-{\displaystyle \frac{\phi y_2}{y}}+\beta y_2-\nu \right)w-dz\hfill \\ \hfill \le & -\frac{r{C}_5}{k}{(x-x_2)}^2+a(e-C_5)(x-x_2)(y-y_2)\hfill \\ \hfill & +y_2\left(\lambda -{\displaystyle \frac{\gamma }{y_{max}}}\right)z.\hfill \end{array}$$

Suppose $C_5=e$. Thus, ${}^{c}D^q\le 0$ when $\lambda <{\displaystyle \frac{\gamma }{y_{max}}}$ which is equivalent to $r\lambda k<4\sigma \gamma$. Hence, the proof is established. □

Theorem 5. 

The equilibrium point $E_3$ is globally asymptotically stable if $ae{x}_3<\lambda z_3+m$, $\lambda <{\displaystyle \frac{\gamma }{y_{max}}}+{\displaystyle \frac{d}{y_3}}$, $\beta y_3<\nu$, and $m{y}_3+\gamma z_3+d{z}_3<ae{x}_3{y}_3$.

Proof. 

A suitable positive definite Lyapunov function is considered as follows:

$$V=e\left(x-x_3-x_{3}ln\left({\displaystyle \frac{x}{x_3}}\right)\right)+y-y_3-y_{3}ln\left({\displaystyle \frac{y}{y_3}}\right)+z-z_3-z_{3}ln\left({\displaystyle \frac{z}{z_3}}\right)+w.$$

By calculating the q-order derivative of V throughout the solution of (4) and applying Lemma 3.1 in [38],

$$\begin{array}{cc}\hfill {}^{c}D^{q}V\le & e(x-x_3)\left(r-{\displaystyle \frac{rx}{k}}-ay-H\right)+\left(1-{\displaystyle \frac{y_3}{y}}\right)(aexy-\lambda yz-\beta yw+\gamma z+\phi w-my)\hfill \\ \hfill & +\left(1-{\displaystyle \frac{z_3}{z}}\right)(\lambda yz-\gamma z-dz+\theta w)+\beta yw-\phi w-\nu w-\theta w\hfill \\ \hfill \le & -\frac{re}{k}{(x-x_3)}^2+(ae{x}_3-\lambda z_3-m)y+d(z_3-z)\hfill \\ \hfill & +\left(\lambda y_3-{\displaystyle \frac{\gamma y_3}{y}}\right)z+(\beta y_3-\nu )w+(m{y}_3+\gamma z_3-ae{x}_3{y}_3)\hfill \\ \hfill \le & -\frac{re}{k}{(x-x_3)}^2+(ae{x}_3-\lambda z_3-m)y\hfill \\ \hfill & +\left(\lambda y_3-{\displaystyle \frac{\gamma y_3}{y_{max}}}-d\right)z+(\beta y_3-\nu )w+(m{y}_3+\gamma z_3+d{z}_3-ae{x}_3{y}_3).\hfill \end{array}$$

Thus, ${}^{c}D^q\le 0$ when $ae{x}_3<\lambda z_3+m$, $\lambda <{\displaystyle \frac{\gamma }{y_{max}}}+{\displaystyle \frac{d}{y_3}}$, $\beta y_3<\nu$, and $m{y}_3+\gamma z_3+d{z}_3<ae{x}_3{y}_3$. Consequently, the theorem is proven. □

Theorem 6. 

The equilibrium point $E_4$ is globally asymptotically stable if $ae{x}_4<\lambda z_4+\beta w_4+m$, $\lambda <{\displaystyle \frac{\gamma }{y_{max}}}+{\displaystyle \frac{d}{y_4}}$, $\beta y_4<{\displaystyle \frac{\phi y_4}{y_{max}}}+{\displaystyle \frac{\theta z_4}{z_{max}}}+\nu$, and $m{y}_4+\gamma z_4+\xi w_4+d{z}_4<ae{x}_4{y}_4$.

Proof. 

A suitable positive definite Lyapunov function is considered as follows:

$$V=e\left(x-x_4-x_{4}ln\left({\displaystyle \frac{x}{x_4}}\right)\right)+y-y_4-y_{4}ln\left({\displaystyle \frac{y}{y_4}}\right)+z-z_4-z_{4}ln\left({\displaystyle \frac{z}{z_4}}\right)+w-w_4-w_{4}ln\left({\displaystyle \frac{w}{w_4}}\right).$$

By calculating the q-order derivative of V throughout the solution of (4) and applying Lemma 3.1 in [38],

$$\begin{array}{cc}\hfill {}^{c}D^{q}V\le & e(x-x_4)\left(r-{\displaystyle \frac{rx}{k}}-ay-H\right)+\left(1-{\displaystyle \frac{y_4}{y}}\right)(aexy-\lambda yz-\beta yw+\gamma z+\phi w-my)\hfill \\ \hfill & +\left(1-{\displaystyle \frac{z_4}{z}}\right)(\lambda yz-\gamma z-dz+\theta w)+(w-w_4)(\beta y-\phi -\nu -\theta )\hfill \\ \hfill \le & -\frac{re}{k}{(x-x_4)}^2+(ae{x}_4-\lambda z_4-\beta w_4-m)y\hfill \\ \hfill & +\left(\lambda y_4-{\displaystyle \frac{\gamma y_4}{y}}\right)z+\left(\beta y_4-\nu -{\displaystyle \frac{\phi y_4}{y}}-{\displaystyle \frac{\theta z_4}{z}}\right)w\hfill \\ \hfill & +(m{y}_4+\gamma z_4+\xi w_4-ae{x}_4{y}_4)+d(z_4-z)\hfill \\ \hfill \le & -\frac{re}{k}{(x-x_4)}^2+(ae{x}_4-\lambda z_4-\beta w_4-m)y\hfill \\ \hfill & +\left(\lambda y_4-{\displaystyle \frac{\gamma y_4}{y_{max}}}-d\right)z+\left(\beta y_4-\nu -{\displaystyle \frac{\phi y_4}{y_{max}}}-{\displaystyle \frac{\theta z_4}{z_{max}}}\right)w\hfill \\ \hfill & +(m{y}_4+\gamma z_4+\xi w_4+d{z}_4-ae{x}_4{y}_4).\hfill \end{array}$$

Thus, ${}^{c}D^{q}V(x,y,z)\le 0,$ when $ae{x}_4<\lambda z_4+\beta w_4+m$, $\lambda <{\displaystyle \frac{\gamma }{y_{max}}}+{\displaystyle \frac{d}{y_4}}$, $\beta y_4<{\displaystyle \frac{\phi y_4}{y_{max}}}+{\displaystyle \frac{\theta z_4}{z_{max}}}+\nu$ and $m{y}_4+\gamma z_4+\xi w_4+d{z}_4<ae{x}_4{y}_4$. By applying Lemma 4.6 in [39], the equilibrium point $E_4$ is proven to be globally asymptotically stable. □

4. Numerical Simulations

This section presents numerical simulations performed using the numerical method described in [40,41]. The numerical simulations are conducted to illustrate the theoretical findings regarding the fractional order $\left(q\right)$ and stability of model (4). The parameter values used in the simulations are detailed in

Table 2. Parameter values for model (4).

Case	r	m	e	a	k	$\mathit{\gamma }$	$\mathit{\phi }$	$\mathit{\beta }$	$\mathit{\lambda }$	d	$\mathit{\nu }$	H	$\mathit{\theta }$	Figures
1	1	$0.5$	$0.07$	$0.02$	100	$0.9$	$0.3$	$0.2$	$0.4$	$0.25$	$0.4$	$0.01$	$0.01$	Figure 1
2	1	$0.5$	$0.7$	$0.5$	1000	$0.25$	$0.47$	$0.6$	$0.48$	$0.39$	$0.33$	$0.35$	$0.1$	Figure 2
3	1	$0.05$	$0.7$	$0.2$	5000	$0.9$	$0.3$	$0.2$	$0.4$	$0.25$	$0.4$	$0.01$	$0.01$	Figure 3
4	1	$0.05$	$0.6$	$0.5$	1000	$0.25$	$0.3$	$0.6$	$0.48$	$0.39$	$0.33$	$0.1$	$0.1$	Figure 4

, and most of them are given in [30].

In case 1 of Table 2, the fractional-order model (4) shows the equilibrium point $E_1=(99,0,0,0)$, where all populations are healthy, and no infections exist. In this case, ${\Re }_0=0.2772<1$, which indicates that $E_1$ is locally asymptotically stable. This coincides with Lemma 1 and is indicated in Figure 1. Figure 1 demonstrates that the populations remain stable across various values of the fractional order $\left(q\right)$, with the solutions reaching the equilibrium point $E_1=(99,0,0,0)$.

In case 2 of Table 2, the fractional-order model (4) shows the equilibrium point $E_2=(1.428,1.297,0,0)$. In this case, $y_2=1.29714<min\left\{{\displaystyle \frac{\xi }{\beta }}=1.5,{\displaystyle \frac{\eta }{\lambda }}=1.3\right\}$, which means that $E_2$ is locally asymptotically stable. This coincides with the result of Lemma 2 and is shown in Figure 2. It can be observed from Figure 2 that the oscillation of fractional-order model (4) decreases with decreasing the value of the fractional order $\left(q\right).$ Figure 2 illustrates that the populations maintain stability for various values of the fractional order $(q\in (0,1\left)\right)$, with the solutions reaching the equilibrium point $E_2=(1.428,1.297,0,0)$.

In case 3 of Table 2, the fractional-order model (4) shows the equilibrium point $E_3=(2075,2.875,3340.175,0)$. In this case, ${\displaystyle \frac{\gamma }{\lambda }}=2.25<y_3=2.875<{\displaystyle \frac{\xi }{\beta }}=3.55$, which indicates that $E_3$ is locally asymptotically stable. This coincides with the result of Lemma 3 and is indicated in Figure 3. In order to verify the Routh–Hurwitz criteria of Lemma 3, one has $B_1=1046.04>0$, $B_2=934.987>0$, $B_3=138.617>0$ and $B_1{B}_2-B_3=977891>0$. Therefore, the fractional-order model (4) exhibits local asymptotic stability around $E_3$, as demonstrated in Figure 3. Figure 3 shows that the populations remain stable for different values of fractional order $(q\in (0,1\left)\right),$ with the solutions reaching the equilibrium point $E_3=(2075,2.875,3340.175,0)$.

In case 4 of Table 2 the fractional-order model (4) shows the coexistence equilibrium point $E_4=(291.667,1.217,185.104,103.658)$, where all the populations in the ecosystem coexist: the prey $\left(x\right)$, susceptible predator $\left(y\right)$, predator infected by the first disease $\left(z\right)$, and predator infected by the second disease $\left(w\right)$ reach constant levels over time. In this case, the $E_4$ is locally asymptotically stable as shown in Figure 4. This means that the two infectious diseases will persist in the predator population. Figure 4 indicates that the populations remain stable for different values of fractional order $(q\in (0,1\left)\right)$, with the solutions reaching the equilibrium point $E_4=(291.667,1.217,185.104,103.658)$. It is to be noted that the integer-order model (1) given in [30] cannot be sustained at a stable coexistence equilibrium level. However, the newly proposed fractional-order model (4) can be sustained at the stable coexistence equilibrium level as illustrated in Figure 4 and coincides with Theorem 6. Therefore, the fractional order has a stabilization effect.

Figure 5 shows the 3D plot of the basic reproduction number ${\Re }_0$ when the predation rate of susceptible predator (a) and prey harvesting (H) varies. It is observed that as the predation rate of susceptible predator (a) increases, ${\Re }_0$ will increase and cross the threshold ${\Re }_0=1$, thus leading to the outbreak of the diseases. Moreover, when the prey harvesting (H) increases, ${\Re }_0$ will increase. Therefore, one can control the reproduction ${\Re }_0$ by reducing the predation rate of susceptible predator (a) and prey harvesting (H).

5. Conclusions

This paper proposed and analyzed a fractional-order eco-epidemiological model incorporating two disease strains in the predator population and harvesting. The model categorizes the populations into four groups: prey $\left(x\right)$, susceptible predators $\left(y\right)$, predators infected by the first disease $\left(z\right)$, and predators infected by the second disease $\left(w\right)$. The proposed model (4) has been analyzed to investigate its dynamical behavior. The model’s dynamics, including positivity, boundedness, and the existence and uniqueness of solutions, have been studied. The proposed eco-epidemiological model exhibits four non-negative equilibrium points, and the threshold parameters have been utilized to determine equilibrium existence and stability conditions. Furthermore, sufficient conditions for the locally asymptotic stability of the four equilibrium points have been derived. The global properties of the equilibrium points $E_1$, $E_2$, $E_3$ and $E_4$ have been investigated by constructing suitable Lyapunov functions. Numerical simulations have been performed to illustrate the theoretical findings, demonstrating the influence of the fractional order $\left(q\right)$ on the stability of the equilibrium points. It has been shown that the populations remain stable for different values of fractional order $(q\in (0,1\left)\right)$, though the solutions reach the obtained equilibrium points. It has been observed that the integer-order model (1) given in [30] cannot be sustained at a stable coexistence equilibrium level. However, it has been shown that the fractional-order model (4) can be sustained at the stable coexistence equilibrium level. Therefore, the fractional order has a stabilization effect. Future research will explore the inclusion of time delays in the system and analyze their potential effect.