A study on fractional order infectious chronic wasting disease model in deers

Abstract

In this paper, the fractional order analysis of the behaviour of the four-dimensional chronic wasting disease (CWD) communities model has been presented. The complexity and dynamical behaviour of the CWD model have been calculated using two fractional derivatives. CWD, a neurological disease that affects deer, has resulted in many deaths and infections among deer populations around the world. To better understand and tackle this eco-epidemiological issue, we used two numerical schemes utilizing the Caputo fractional operator and the Atangana-Baleanu (AB) fractional operator. We have investigated the stability of the eco-epidemiological CWD model. The fixed point theory is a volumetric role play in analyzing the existence and uniqueness of the solution. We use bifurcation diagrams, time series diagrams, and phase diagrams to analyze fractional-order eco-epedimological systems with derivative orders and parameters varying. We examine the approximate result of the eco-epidemiological CWD model with the Atangana-Baleanu (AB) operator and the Caputo operator, and we also compare both solutions. We do a brief analysis of the simulated results, which reveals that the suggested methodologies are novel, dependable, and remarkably easy to implement. In addition to determining the direction, stability, bifurcating, and numerical solutions, graphic depiction and graph bifurcation provide better information about the proposed model.

Keywords:

• Atangana-Baleanu (AB) operator
• bifurcation
• caputo operator
• chronic wasting disease (CWD) model
• existence and uniqueness
• stability
• Toufik-Atangana numerical scheme

1. Introduction

Fractional calculus has become increasingly important in various interdisciplinary fields such as Mathematics, Physics, Biology, and Engineering. Fractional calculus is associated with Mathematics, which is as old as calculus, and it deals with arbitrary orders of differentiation and integration (Mainardi, 2012). In recent times, a lot of biological models have been analyzed using fractional derivatives (Caputo, 1969; Kumar, Alshahrani, Yakout, Abdel-Aty, & Kumar, 2021). The Caputo fractional operator is the most commonly used operator for simulating real-world problems. The Riemann-Liouville (RL) operator is also closely related to the Caputo fractional operator, but the main difference between them is the singularity quality of the kernels. Using these operators may not always lead to better results when studying the dynamics of various models. Therefore, researchers suggest using fractional operators with non-singular kernels to better understand the dynamics of models. Atangana and Baleanu developed a non-singular kernel derivative operator using the Mittag-Leffer function, which has a non-local and non-singular kernel, making it particularly useful for individuals engaged in numerical modelling of real-world problems. The Atangana-Baleanu (AB) derivative has gained significant importance and popularity in recent decades in various fields, particularly biological models (Alkahtani & Atangana, 2016; Kumar, Alshahrani, et al., 2021). In fractional calculus, many derivatives are there but the Riemann Liouville (RL), Caputo, and Atangana-Baleanu (AB) operators play major roles in the analysis of biological models but still need more particular attention (Losada & Nieto, 2015; Uçar, Uçar, Özdemir, & Hammouch, 2019). Fractional order derivatives can be used to illustrate a variety of natural occurrences and facts that have non-local, intricate dynamical behaviour. When examining the system, such operators make it easier to solve the system. As is well known, biological models have memory effects and hereditary characteristics, and fractional calculus is a better way to describe them. Fractional-order systems, as opposed to integer-order systems, results that are more accurate when describing the complicated behaviour of epidemic diseases. In this paper, the authors introduce a new method for modelling the dynamics of influenza A disease and the Caputo-Fabrizio fractional derivative operator into the model, a more accurate representation of the disease’s behaviour (Evirgen, Esmehan, Sümeyra, & Özdemir, 2023). In this paper, Veeresha et al. have given information about Zika virus disease spread through the bite of an infected mosquito from the Aedes species. It can become a severe epidemic if not controlled during its early stages. The nonlinear partial differential model is analyzed and solved using the q-homotopy analysis transform method. The Atangana-Baleanu (AB) definition of fractional derivative is utilized in the process (Veeresha, Akinyemi, Oluwasegun, Şenol, & Oduro, 2022).

Chronic wasting disease (CWD) is a highly infectious and dangerous disease that affects deer (Williams & Young, 1980). The disease was first identified in captive mule deer at a wildlife research lab in northern Colorado, USA, and southeastern Wyoming in 1967 (Sigurdson, 2008). In 1996, CWD was discovered among free-ranging deer and elk in northeastern Colorado and southeastern Wyoming. Apart from North America, this disease has only been found in South Korea, where it is believed to have originated from Canada (Miller et al., 2000). From 1996 to 1999, the prevalence of this disease was estimated to be approximately 2% in white-tailed deer, less than 1% in elk, and 5% in mule deer, based on disease surveillance of harvested animals in this area (Miller et al., 2000). Captive deer and elk are highly susceptible to this disease and it can spread quickly among them. However, we still do not fully understand the source and transmission mechanisms of the prions that cause CWD.

In the past, Chronic Wasting Disease (CWD) used to be considered a rare and foreign disease. However, in recent years, this disease has become increasingly concerning as it has emerged as a transmissible branch of prion diseases (Barlow, 1996; Lafferty & Holt, 2003). The most effective ways to eliminate diseases are through vaccination and prevention. Unfortunately, there is still no vaccination strategy available for lethal illnesses such as CWD. This disease is quite peculiar as it affects both farmed and wild animals. Animals that are infected with CWD can directly and indirectly affect other animals. When these infected animals leave behind excrement and carcasses, they can be harmful to susceptible animals. To truly understand the disease and its incidence in wildlife, we need to conduct long-term investigations. The narration of the symptoms of a certain disease is derived from the scanning of infected animals (Thompson Hobbs, 2006). This disease can be transmitted naturally through bodily fluids such as sputum, urine, and faeces. Infected deer blood in platelets and B cells carries the infectiousness but plasma is not infected. Recent research suggests that CWD prions may also be present in the skeletal muscle of infected animals. Symptoms of this disease in animals include excessive salivation, physiological changes, increased thirst, loss of coordination, changes in behaviour, difficulty swallowing, and more. CWD prions can remain in the environment for at least two years, making the disease highly infectious for the cervine family, including bucks, moose, and elk (Hsieh & Hsiao, 2008; Mukherjee, 2003; 2003; Venturino, 2002). Unfortunately, there is currently no perfect medication or treatment for this disease. The primary approach to managing the disease is to hunt infected animals and remove them from their families. In this article, we aim to analyze the eco-epidemic model with indirect transmission and explore ways to control CWD transmission through predation. A reliable eco-epidemiological model can provide vital information about diseases, and further research in this field can be conducted with the assistance of such models.

Some numerical schemes have been developed as mathematical tools for solving biological models and fractional-order differential equations (Kumar, Kumar, & Jleli, 2020). Recently, a new fractional operator was created to address issues with the power-law kernel and the existing fractional operator. Literature is abundant on various numerical methods, such as the predictor-corrector method (Evirgen, Esmehan, et al., 2023; 2023; Joshi & Yavuz, 2023), Newton interpolation formula (Rahman, Arfan, & Baleanu, 2023), Adam-Bashforth method, Toufik-Atangana technique (Joshi, Yavuz, Townley, & Jha, 2023) and Runge-Kutta technique (Fatima, Yavuz, Ur Rahman, & Al-Duais, 2023). The Toufik-Atangana approach is highly regarded as an effective numerical technique for solving non-linear systems (Khan & Atangana, 2020; Toufik & Atangana, 2017). To overcome the shortcomings of the popular Adams-Bashforth approach, the Toufik-Atangana numerical system is introduced. The novel numerical approach has the combination of the two-step Lagrange polynomial and the fundamental theorem of fractional calculus. This approach converges to the solution rapidly and with great accuracy. We used this numerical technique for investigating the eco-epidemiological system with the help of the Atangana-Baleanu (AB) fractional operator (Toufik & Atangana, 2017). In this paper, we have used this scheme with different fractional operators in this biological model.

The work is organized as follows: Section 2 examines fractional calculus definitions and basic concepts. In section 3, we examine Caputo and Atangana-Baleanu’s (AB) non-integer eco-epidemiological model. In section 4, we examine equilibrium points and stability. In section 5, we address the uniqueness and existence of the proposed model solution. In section 6, we present the numerical technique for the Toufik-Atangana (TA) with the Caputo operator. TA numerical scheme using the ABC operator is presented in section 7. Graphical analysis and numerical results are presented in section 8. The study’s conclusions are discussed in section 9.

2. Definitions and basic concepts

In this portion, we have present some essential elucidation, theorems, and symbols of the fractional calculus theory that are applied in the task.

Definition 2.1

(Uçar et al., 2019) Suppose that $\mho$ $\in$ $H^1(a,b)$ and $a<b,$ so the Caputo fractional differential operator is defined as follows (1) $${}_a^C{D}_t^{\phi }[\mho (t)]=\frac{1}{\Gamma (m-\phi )}{\int }_a^t{\mho }^{(m)}(\mathfrak{z}){(t-\mathfrak{z})}^{(m-\phi -1)}d\mathfrak{z},\mathrm{\hspace{1em}}\mathrm{\hspace{1em}}m-1<\phi <m\in \mathbb{N},$$(1) $\Gamma (.)$ denotes a gamma function.

Definition 2.2.

The Caputo integral is defined as follows (2) $${}_a^C{I}_t^{\phi }[\mho (t)]=\frac{1}{\Gamma (\phi )}{\int }_a^t\mho (\mathfrak{z}){(t-\mathfrak{z})}^{\phi -1}d\mathfrak{z}.$$(2)

Definition 2.3

(Uçar et al., 2019). Suppose that $\mho$ $\in$ $H^1(a,b),$ $a<b,$ $\phi$ $\in$ [0,1]. The Atangana-Baleanu (AB) operator in Riemann-Liouville type is defined as follows (3) $${{}^{\mathit{ABRL}}}_a{D}_t^{\phi }[\mho (t)]=\frac{B(\phi )}{1-\phi }\frac{d}{dt}{\int }_a^t\mho (\mathfrak{z}){\mathbb{E}}_{\phi }\left(-\frac{\phi {(t-\mathfrak{z})}^{\phi }}{1-\phi }\right)d\mathfrak{z}.$$(3)

It is usually the case that B(0) = B(1) = 1, and $B(\phi )$ is a normalization function.

Definition 2.4

(Uçar et al., 2019) Suppose that $\mho$ $\in$ $H^1(a,b),$ $a<b,$ and $\phi$ $\in$ [0,1]. Then the Atangana-Baleanu (AB) operator in Caputo type is defined as follows (4) $${{}^{\mathit{ABC}}}_a{D}_t^{\phi }[\mho (t)]=\frac{B(\phi )}{1-\phi }{\int }_a^t{\mho }^{\prime }(\mathfrak{z}){\mathbb{E}}_{\phi }\left(-\frac{\phi {(t-\mathfrak{z})}^{\phi }}{1-\phi }\right)d\mathfrak{z}.$$(4)

Definition 2.5

(Atangana & Koca, 2016) Suppose that $\prime a\prime$ be the base point and $\phi$ is the order of operator in the Atangana-Baleanu (AB) integral operator is defined as follows (5) $${{}^{AB}}_a{I}_t^{\phi }[\mho (t)]=\frac{1-\phi }{B(\phi )}\mho (t)+\frac{\phi }{B(\phi )\Gamma (\phi )}{\int }_a^t\mho (\mathfrak{z}){(t-\mathfrak{z})}^{\phi -1}d\mathfrak{z}.$$(5)

Theorem 2.1

(Chen, Petras, & Xue, 2009; Irudayaraj et al., 2020; Santosh Kumar, 2014) Matignon’s stability theorem and its condition is (6) $$|\arg (\mathit{eig}(\mathbb{M}))|=|\arg (\overline{{\lambda }_{\mathbf{j}}})|>\phi \frac{\pi }{2}, \mathbf{j}=1,2,\dots ,m,$$(6) where 0 < $\phi$ < 1 and $\mathit{eig}(\mathbb{M})$ (i.e $\overline{{\lambda }_{\mathbf{j}}}$) represents the eigenvalues of matrix $\mathbb{M}.$

3. Caputo and Atangana-Baleanu fractional model of eco-epidemiological

Some years ago, CWD was called zombie deer disease. This deadly- infectious and neurological infection influences the deer family. Based on experiments, it has been seen that CWD can be transmitted to vulnerable animals through the residues of excreta that are left behind in the environment by infected animals and their carcasses (Packer, Holt, Hudson, Lafferty, & Dobson, 2003). This mode of transmission is considered older than the traditional models of direct contact between animals (Miller & Wild, 2004). It can be challenging to observe these diseases thoroughly since a prolonged epizootic may induce low, frequently undiagnosed infection morbidity or mortality. Long-term investigations will be beneficial to comprehending the incidence and geographical dynamics of chronic wildlife disease. It’s significant to understand the patterns and dynamics of chronic disease in wildlife, and one way to do that is through modelling the complexity of the situation. By doing so, we can gain a better understanding of how these diseases spread over time and space, which can help us develop effective strategies for managing and preventing them. Although CWD prions can remain contagious in the environment for years, relatively little is known about the possible consequences of indirect CWD dynamic transfer. The proposed model obeyed a hunt-predator method and split the hunt population into two parts, the first part infected animals and the second part non-infected animals. The susceptible-infected-recovered (SIR) model in infection transfer among humans is the foundation of the eco-epidemiological model (Allen, Brauer, Van den Driessche, & Wu, 2008). This model is expressed in terms of nonlinear equations. Hunter’s functional reaction is supposed to be the Holling type II. Based on the aforementioned, the integer-order CWD model (Maji, Mukherjee, & Kesh, 2018) is as follows: (7) $$\begin{array}{c}\frac{dS}{dt}=S\left[r\left(1-\frac{S+I}{g}\right)-\beta E\right],\\ \frac{dI}{dt}=\beta SE-\nu I-\frac{\mathit{aIP}}{m+I},\\ \frac{dE}{dt}=ϵI-\varrho E,\\ \frac{dP}{dt}=P\left(-d+\frac{bI}{m+I}\right),\end{array}$$(7) with initial conditions S(0) > 0, I(0) > 0, E(0) > 0, P(0) $\ge$ 0.

In this scenario, S represents the density of the susceptible animals, while I denotes the infected ones. E represents the mass of infectious material in the environment, and P is the density of the predator. In this context, the various parameters are defined as follows:

• r represents the maximum per-capita growth rate,

• d represents the death rate of the predator,

• g represents the environmental carrying capacity of the prey,

• ϱ represents the mass-specific rate of loss of infectious material from the environment,

• ν represents the death rate of infected animals through CWD,

• ϵ represents the per-capita rate of excretion of infectious material by infected animals,

• a represents the predation coefficient,

• m represents a half-saturation constant,

• b represents the conversion coefficient of the predator,

• β represents the indirect transmission coefficient for the disease.

When animals that are vulnerable to infection come in contact with infectious material in their surroundings, they have contracted the disease. These infected animals can then transmit the infection to others through contact with the contaminated environment, thus facilitating the spread of the disease. This model can be applied to the spread of various other diseases such as tuberculosis in livestock (Brennan, Kemp, & Christley, 2008), viral hepatitis A (Ajelli, Iannelli, Manfredi, & Ciofi Degli Atti, 2008), Vibrio cholerae (Das & Mukherjee, 2012), avian influenza (Breban, Drake, Stallknecht, & Rohani, 2009), and so on.

The fractional Caputo derivative is applied in system (7)(7) $$\begin{array}{c}\frac{dS}{dt}=S\left[r\left(1-\frac{S+I}{g}\right)-\beta E\right],\\ \frac{dI}{dt}=\beta SE-\nu I-\frac{\mathit{aIP}}{m+I},\\ \frac{dE}{dt}=ϵI-\varrho E,\\ \frac{dP}{dt}=P\left(-d+\frac{bI}{m+I}\right),\end{array}$$(7) . Thus, the existing model is replaced by the following model: (8) $$\begin{array}{c}{{}^C}_0{D}_t^{\phi }S=S\left[r\left(1-\frac{S+I}{g}\right)-\beta E\right],\\ {{}^C}_0{D}_t^{\phi }I=\beta SE-\nu I-\frac{\mathit{aIP}}{m+I},\\ {{}^C}_0{D}_t^{\phi }E=ϵI-\varrho E,\\ {{}^C}_0{D}_t^{\phi }P=P\left(-d+\frac{bI}{m+I}\right),\end{array}$$(8) with initial conditions S(0) > 0, I(0) > 0, E(0) > 0, P(0) $\ge$ 0.

To gain a fractional derivative model, we operate the Atangana-Baleanu operator in system (7)(7) $$\begin{array}{c}\frac{dS}{dt}=S\left[r\left(1-\frac{S+I}{g}\right)-\beta E\right],\\ \frac{dI}{dt}=\beta SE-\nu I-\frac{\mathit{aIP}}{m+I},\\ \frac{dE}{dt}=ϵI-\varrho E,\\ \frac{dP}{dt}=P\left(-d+\frac{bI}{m+I}\right),\end{array}$$(7) . We found the new model is an expression of the Atangana-Baleanu (AB) operator: (9) $$\begin{array}{c}{{}^{\mathit{ABC}}}_0{D}_t^{\phi }S=S\left[r\left(1-\frac{S+I}{g}\right)-\beta E\right],\\ {{}^{\mathit{ABC}}}_0{D}_t^{\phi }I=\beta SE-\nu I-\frac{\mathit{aIP}}{m+I},\\ {{}^{\mathit{ABC}}}_0{D}_t^{\phi }E=ϵI-\varrho E,\\ {{}^{\mathit{ABC}}}_0{D}_t^{\phi }P=P\left(-d+\frac{bI}{m+I}\right),\end{array}$$(9) with initial conditions S(0) > 0, I(0) > 0, E(0) > 0, P(0) $\ge$ 0.

4. Equilibrium points and stability analysis

The fractional-order eco-epidemiological model is employed here to examine its stability. At first, our work is to find out the equilibrium points of the system (9)(9) $$\begin{array}{c}{{}^{\mathit{ABC}}}_0{D}_t^{\phi }S=S\left[r\left(1-\frac{S+I}{g}\right)-\beta E\right],\\ {{}^{\mathit{ABC}}}_0{D}_t^{\phi }I=\beta SE-\nu I-\frac{\mathit{aIP}}{m+I},\\ {{}^{\mathit{ABC}}}_0{D}_t^{\phi }E=ϵI-\varrho E,\\ {{}^{\mathit{ABC}}}_0{D}_t^{\phi }P=P\left(-d+\frac{bI}{m+I}\right),\end{array}$$(9) . (10) $$\begin{array}{c}S\left[r\left(1-\frac{S+I}{g}\right)-\beta E\right]=0,\\ \beta SE-\nu I-\frac{\mathit{aIP}}{m+I}=0,\\ ϵI-\varrho E=0,\\ P\left(-d+\frac{bI}{m+I}\right)=0.\end{array}$$(10)

We have found four equilibrium points:

Table

List of Equilibrium Points
Serial No	Equilibrium points	Conditions
1.	E0 = (0,0,0,0)	–
2.	E1 = (g,0,0,0)	$\nu \varrho$ > $\beta gϵ$
3.	E2 = $(\widehat{S},\widehat{I},\widehat{E},0)$	$\frac{r}{\varrho }\widehat{S}$ > g > $\widehat{S}$ and d < $\frac{b\widehat{I}}{m+\widehat{I}}$
4.	E3 = $({\mathbb{S}}^{*},{\mathbb{I}}^{*},{\mathbb{E}}^{*},{\mathbb{P}}^{*})$	b > $\frac{d(gr\varrho +\mathit{\text{gmϵβ}}+mr\varrho )}{gr\varrho }$ and ${\mathbb{S}}^{*}$ > $\frac{\varrho \nu }{ϵ\beta }$

Let (S, I, E, P) be any point, then the Jacobian matrix for system (9)(9) $$\begin{array}{c}{{}^{\mathit{ABC}}}_0{D}_t^{\phi }S=S\left[r\left(1-\frac{S+I}{g}\right)-\beta E\right],\\ {{}^{\mathit{ABC}}}_0{D}_t^{\phi }I=\beta SE-\nu I-\frac{\mathit{aIP}}{m+I},\\ {{}^{\mathit{ABC}}}_0{D}_t^{\phi }E=ϵI-\varrho E,\\ {{}^{\mathit{ABC}}}_0{D}_t^{\phi }P=P\left(-d+\frac{bI}{m+I}\right),\end{array}$$(9) is described as (11) $$J(S,I,E,P)=\left[\begin{array}{cccc}r-\frac{2rS}{g}-\frac{rI}{g}-\beta E& -\frac{rS}{g}& -\beta S& 0\\ \beta E& -\nu -\frac{\mathit{maP}}{{(m+I)}^2}& \beta S& -\frac{aI}{m+I}\\ 0& ϵ& -\varrho & 0\\ 0& \frac{\mathit{mbP}}{{(m+I)}^2}& 0& -d+\frac{bI}{m+I}\end{array}\right]$$(11)

In this portion, we have investigated the behaviour of the equilibrium points of the fractional epidemiological system (9)(9) $$\begin{array}{c}{{}^{\mathit{ABC}}}_0{D}_t^{\phi }S=S\left[r\left(1-\frac{S+I}{g}\right)-\beta E\right],\\ {{}^{\mathit{ABC}}}_0{D}_t^{\phi }I=\beta SE-\nu I-\frac{\mathit{aIP}}{m+I},\\ {{}^{\mathit{ABC}}}_0{D}_t^{\phi }E=ϵI-\varrho E,\\ {{}^{\mathit{ABC}}}_0{D}_t^{\phi }P=P\left(-d+\frac{bI}{m+I}\right),\end{array}$$(9) with the help of the Matignon condition and the Jacobian matrix.

Theorem 4.1.

The trivial equilibrium point $E_0=(0,0,0,0)$ is unstable.

Proof 4.1.

The Jacobian matrix (11)(11) $$J(S,I,E,P)=\left[\begin{array}{cccc}r-\frac{2rS}{g}-\frac{rI}{g}-\beta E& -\frac{rS}{g}& -\beta S& 0\\ \beta E& -\nu -\frac{\mathit{maP}}{{(m+I)}^2}& \beta S& -\frac{aI}{m+I}\\ 0& ϵ& -\varrho & 0\\ 0& \frac{\mathit{mbP}}{{(m+I)}^2}& 0& -d+\frac{bI}{m+I}\end{array}\right]$$(11) evaluated at the equilibrium point $E_0=(0,0,0,0).$ $$J(E_0)=\left[\begin{array}{cccc}r& 0& 0& 0\\ 0& -\nu & 0& 0\\ 0& ϵ& -\varrho & 0\\ 0& 0& 0& -d\end{array}\right]$$

The eigenvalues of $J(E_0)$ are ${\overline{\lambda }}_1$ = r, ${\overline{\lambda }}_2$ = $-\nu ,{\overline{\lambda }}_3$ = $-\varrho$ and ${\overline{\lambda }}_4$ =—d. With the help of Matignon’s condition, we observed that $|\mathit{arg}({\overline{\lambda }}_i)|$ = π > $\frac{\phi \pi }{2}$, i = 2, 3, 4 for all 0 < $\phi$ < 1. But ${\overline{\lambda }}_1$ does not satisfy the Matignon condition. Therefore equilibrium point E₀ is unstable.

Theorem 4.2.

The equilibrium point $E_1=(g,0,0,0)$ of system (9)(9) $$\begin{array}{c}{{}^{\mathit{ABC}}}_0{D}_t^{\phi }S=S\left[r\left(1-\frac{S+I}{g}\right)-\beta E\right],\\ {{}^{\mathit{ABC}}}_0{D}_t^{\phi }I=\beta SE-\nu I-\frac{\mathit{aIP}}{m+I},\\ {{}^{\mathit{ABC}}}_0{D}_t^{\phi }E=ϵI-\varrho E,\\ {{}^{\mathit{ABC}}}_0{D}_t^{\phi }P=P\left(-d+\frac{bI}{m+I}\right),\end{array}$$(9) is

1. If ν ϱ > $ϵ\beta g$ then the equilibrium point E₁ is locally asymptotically stable (LAS).

2. If ν ϱ < $ϵ\beta g$ then the equilibrium point E₁ is unstable.

Proof 4.2.

The Jacobian matrix (11)(11) $$J(S,I,E,P)=\left[\begin{array}{cccc}r-\frac{2rS}{g}-\frac{rI}{g}-\beta E& -\frac{rS}{g}& -\beta S& 0\\ \beta E& -\nu -\frac{\mathit{maP}}{{(m+I)}^2}& \beta S& -\frac{aI}{m+I}\\ 0& ϵ& -\varrho & 0\\ 0& \frac{\mathit{mbP}}{{(m+I)}^2}& 0& -d+\frac{bI}{m+I}\end{array}\right]$$(11) evaluated at the equilibrium point $E_1=(g,0,0,0).$ $$J(E_1)=\left[\begin{array}{cccc}-r& -r& -\beta g& 0\\ 0& -\nu & \beta g& 0\\ 0& ϵ& -\varrho & 0\\ 0& 0& 0& -d\end{array}\right]$$

The characteristic equation of $J(E_1)$ is following as $$(\overline{\lambda }+r)(\overline{\lambda }+d)({\overline{\lambda }}^2+\overline{\lambda }(\varrho +\nu )+\varrho \nu -\beta gϵ)=0.$$

Therefore, eigenvalues are ${\overline{\lambda }}_1$ =—r, ${\overline{\lambda }}_2$ = $-\nu ,{\overline{\lambda }}_3$ = −d and ${\overline{\lambda }}_4$ = $-\frac{(\nu \varrho -ϵ\beta g)}{\nu }$. All the eigenvalues have negative real parts, then $|\mathit{arg}({\overline{\lambda }}_i)|$ = π > $\frac{\phi \pi }{2},i=1,2,3,4$ for all 0 < $\phi$ < 1. So, we can say that Matiginon’s condition hold, if $\nu \varrho$ > $ϵ\beta g$. Then the equilibrium point E₁ is locally asymptotically stable (LAS). We assume that if $\nu \varrho$ < $ϵ\beta g$ then the Matiginon’s condition not holds. Hence, the equilibrium point E₁ is unstable in this condition.

Theorem 4.3.

The equilibrium point $E_2=(\widehat{S},\widehat{I},\widehat{E},0)$ of system (9)(9) $$\begin{array}{c}{{}^{\mathit{ABC}}}_0{D}_t^{\phi }S=S\left[r\left(1-\frac{S+I}{g}\right)-\beta E\right],\\ {{}^{\mathit{ABC}}}_0{D}_t^{\phi }I=\beta SE-\nu I-\frac{\mathit{aIP}}{m+I},\\ {{}^{\mathit{ABC}}}_0{D}_t^{\phi }E=ϵI-\varrho E,\\ {{}^{\mathit{ABC}}}_0{D}_t^{\phi }P=P\left(-d+\frac{bI}{m+I}\right),\end{array}$$(9) is

(i) If $\frac{r}{\varrho }$ $\widehat{S}$ > g > $\widehat{S}$ and d < $\frac{b\widehat{I}}{m+\widehat{I}}$ then the equilibrium point E₂ is locally asymptotically stable.

Proof 4.3.

The Jacobian matrix (11)(11) $$J(S,I,E,P)=\left[\begin{array}{cccc}r-\frac{2rS}{g}-\frac{rI}{g}-\beta E& -\frac{rS}{g}& -\beta S& 0\\ \beta E& -\nu -\frac{\mathit{maP}}{{(m+I)}^2}& \beta S& -\frac{aI}{m+I}\\ 0& ϵ& -\varrho & 0\\ 0& \frac{\mathit{mbP}}{{(m+I)}^2}& 0& -d+\frac{bI}{m+I}\end{array}\right]$$(11) around the equilibrium point E₂ is given as $$\begin{array}{c}J(E_2)=\left[\begin{array}{cccc}-\frac{r\widehat{S}}{g}& -\frac{r\widehat{S}}{g}& -\beta \widehat{S}& 0\\ \beta \widehat{E}& -\nu & \beta \widehat{S}& -\frac{a\widehat{I}}{m+\widehat{I}}\\ 0& ϵ& -\varrho & 0\\ 0& 0& 0& -d+\frac{b\widehat{I}}{m+\widehat{I}}\end{array}\right]\\ w\mathit{here} \widehat{S}=\frac{\varrho \nu }{ϵ\beta }, \widehat{I}=\frac{\varrho \widehat{E}}{ϵ}, \widehat{E}=\frac{r(g\beta ϵ-\varrho \nu )}{\beta (r\varrho +gϵ\beta )}.\end{array}$$

So, the characteristic equation of $J(E_2)$ is following as $$\begin{array}{c}(\overline{\lambda }+d-\frac{b\widehat{I}}{m+\widehat{I}})({\overline{\lambda }}^3+{\wp }_1{\overline{\lambda }}^2+{\wp }_2\overline{\lambda }+{\wp }_3)=0,\\ w\mathit{here} {\wp }_1=r\frac{\widehat{S}}{g}+\nu +\varrho , {\wp }_2=r\frac{\widehat{S}}{g}(\nu +\varrho +\beta \widehat{E}), {\wp }_3=\frac{r\beta \varrho \widehat{S}\widehat{E}}{g}+{\beta }^2ϵ\widehat{S}\widehat{E}.\end{array}$$ $\text{Now}, \mathfrak{P}(\overline{\lambda })={\overline{\lambda }}^3+{\wp }_1{\overline{\lambda }}^2+{\wp }_2\overline{\lambda }+{\wp }_3=0,$ and $D(\mathfrak{P})$ is discriminant of the polynomial $\mathfrak{P}(\overline{\lambda }).$

1. When $D(\mathfrak{P})>0,{\wp }_1$, and ${\wp }_3$ are greater than zero and ${\wp }_1{\wp }_2$ - ${\wp }_3$ > 0 and $\phi \in [0,1),$ in this case the equilibrium point E₂ is LAS.

2. When $D(\mathfrak{P})<0,{\wp }_1\ge 0,{\wp }_2\ge 0,$ ${\wp }_3>0$ and $\phi \in [0,\frac{2}{3})$ in this case the equilibrium point E₂ is LAS.

First, we assume that ${\wp }_1$, and ${\wp }_3$ are greater than zero and ${\wp }_1{\wp }_2$ - ${\wp }_3$ > 0. We see clearly one of the root of the characteristic equation is negative (i.e ${\overline{\lambda }}_1$ = $-d+\frac{b\widehat{I}}{m+\widehat{I}}$). With the use of given conditions and our assumption, we found that each of the eigenvalues has negative real parts. Consequently, the equilibrium point $E_2$ is satisfied Matignon’s condition (i.e $|\mathit{arg}({\overline{\lambda }}_i)|$ > $\frac{\phi \pi }{2}$) and the present point is locally asymptotically stable.

Theorem 4.4.

The equilibrium point $E_3=({\mathbb{S}}^{*},{\mathbb{I}}^{*},{\mathbb{E}}^{*},{\mathbb{P}}^{*})$ is locally asymptotically stable (LAS) or unstable.

Proof 4.4.

Now, the equilibrium point E₃ putting in the Jacobian matrix (11)(11) $$J(S,I,E,P)=\left[\begin{array}{cccc}r-\frac{2rS}{g}-\frac{rI}{g}-\beta E& -\frac{rS}{g}& -\beta S& 0\\ \beta E& -\nu -\frac{\mathit{maP}}{{(m+I)}^2}& \beta S& -\frac{aI}{m+I}\\ 0& ϵ& -\varrho & 0\\ 0& \frac{\mathit{mbP}}{{(m+I)}^2}& 0& -d+\frac{bI}{m+I}\end{array}\right]$$(11) , we have $$\begin{array}{l}J(E_3)=\left[\begin{array}{cccc}{\mathfrak{a}}_{11}& {\mathfrak{a}}_{12}& {\mathfrak{a}}_{13}& {\mathfrak{a}}_{14}\\ {\mathfrak{a}}_{21}& {\mathfrak{a}}_{22}& {\mathfrak{a}}_{23}& {\mathfrak{a}}_{24}\\ {\mathfrak{a}}_{31}& {\mathfrak{a}}_{32}& {\mathfrak{a}}_{33}& {\mathfrak{a}}_{34}\\ {\mathfrak{a}}_{41}& {\mathfrak{a}}_{42}& {\mathfrak{a}}_{43}& {\mathfrak{a}}_{44}\end{array}\right]\\ w\mathit{here} {\mathfrak{a}}_{11}=-\frac{r{\mathbb{S}}^{*}}{g}, {\mathfrak{a}}_{12}=-\frac{r{\mathbb{S}}^{*}}{g}, {\mathfrak{a}}_{13}=-\beta {\mathbb{S}}^{*}, \\ {\mathfrak{a}}_{14}=0,\\ {\mathfrak{a}}_{21}=-\beta {\mathbb{E}}^{*}, {\mathfrak{a}}_{22}=-\nu -\frac{ma{\mathbb{P}}^{*}}{{(m+{\mathbb{I}}^{*})}^2}, {\mathfrak{a}}_{23}=\beta {\mathbb{S}}^{*}, \\ {\mathfrak{a}}_{24}=-\frac{a{\mathbb{I}}^{*}}{m+{\mathbb{I}}^{*}},\\ {\mathfrak{a}}_{31}=0, {\mathfrak{a}}_{32}=ϵ, {\mathfrak{a}}_{33}=-\varrho , {\mathfrak{a}}_{34}=0,\\ {\mathfrak{a}}_{41}=0, {\mathfrak{a}}_{42}=\frac{mb{\mathbb{P}}^{*}}{{(m+{\mathbb{I}}^{*})}^2}, {\mathfrak{a}}_{43}=0, {\mathfrak{a}}_{44}=0,\\ {\mathbb{S}}^{*}=\frac{gr\varrho (b-d)-\mathit{dgm}\beta ϵ-\mathit{dmr}\varrho }{r\varrho (b-d)}, {\mathbb{I}}^{*}=\frac{dm}{b-d},\\ {\mathbb{E}}^{*}=\frac{dmϵ}{\varrho (b-d)}, {\mathbb{P}}^{*}=\frac{bm({\mathbb{S}}^{*}ϵ\beta -\varrho \nu )}{a\varrho (b-d)}.\end{array}$$

Therefore, $J(E_3)$ has the characteristic equation $$\begin{array}{l}\mathfrak{P}(\overline{\lambda })={\overline{\lambda }}^4+{\mathfrak{h}}_1{\overline{\lambda }}^3+{\mathfrak{h}}_2{\overline{\lambda }}^2+{\mathfrak{h}}_3\overline{\lambda }+{\mathfrak{h}}_4=0,\\ w\mathit{here} {\mathfrak{h}}_1=-{\mathfrak{a}}_{11}-{\mathfrak{a}}_{22}-{\mathfrak{a}}_{33},\\ {\mathfrak{h}}_2={\mathfrak{a}}_{11}{\mathfrak{a}}_{22}+{\mathfrak{a}}_{11}{\mathfrak{a}}_{33}-{\mathfrak{a}}_{23}{\mathfrak{a}}_{32}-{\mathfrak{a}}_{24}{\mathfrak{a}}_{42}+{\mathfrak{a}}_{22}{\mathfrak{a}}_{33}-{\mathfrak{a}}_{12}{\mathfrak{a}}_{21},\\ {\mathfrak{h}}_3={\mathfrak{a}}_{11}{\mathfrak{a}}_{23}{\mathfrak{a}}_{32}+{\mathfrak{a}}_{11}{\mathfrak{a}}_{24}{\mathfrak{a}}_{42}-{\mathfrak{a}}_{11}{\mathfrak{a}}_{22}{\mathfrak{a}}_{33}+{\mathfrak{a}}_{24}{\mathfrak{a}}_{42}{\mathfrak{a}}_{33}\\ +{\mathfrak{a}}_{12}{\mathfrak{a}}_{21}{\mathfrak{a}}_{33}-{\mathfrak{a}}_{13}{\mathfrak{a}}_{21}{\mathfrak{a}}_{32},\\ {\mathfrak{h}}_4=-{\mathfrak{a}}_{11}{\mathfrak{a}}_{24}{\mathfrak{a}}_{42}{\mathfrak{a}}_{33}.\end{array}$$

We analyze the behavior of equilibrium point E₃ and $D(\mathfrak{P})$ is discriminant of the polynomial $\mathfrak{P}(\overline{\lambda }).$ $$D(\mathfrak{P})=\left|\begin{array}{ccccccc}1& {\mathfrak{h}}_1& {\mathfrak{h}}_2& {\mathfrak{h}}_3& {\mathfrak{h}}_4& 0& 0\\ 0& 1& {\mathfrak{h}}_1& {\mathfrak{h}}_2& {\mathfrak{h}}_3& {\mathfrak{h}}_4& 0\\ 0& 0& 1& {\mathfrak{h}}_1& {\mathfrak{h}}_2& {\mathfrak{h}}_3& {\mathfrak{h}}_4\\ 4& 3{\mathfrak{h}}_1& 2{\mathfrak{h}}_2& {\mathfrak{h}}_3& 0& 0& 0\\ 0& 4& 3{\mathfrak{h}}_1& 2{\mathfrak{h}}_2& {\mathfrak{h}}_3& 0& 0\\ 0& 0& 4& 3{\mathfrak{h}}_1& 2{\mathfrak{h}}_2& {\mathfrak{h}}_3& 0\\ 0& 0& 0& 4& 3{\mathfrak{h}}_1& 2{\mathfrak{h}}_2& {\mathfrak{h}}_3\end{array}\right|.$$

(i) Let ${\Delta }_1,{\Delta }_2,{\Delta }_3$ and ${\Delta }_4$ are Routh-Hurwitz determinants $${\Delta }_1={\mathfrak{h}}_1,{\Delta }_2=\left|\begin{array}{cc}{\mathfrak{h}}_1& 1\\ {\mathfrak{h}}_3& {\mathfrak{h}}_2\end{array}\right|,{\Delta }_3=\left|\begin{array}{ccc}{\mathfrak{h}}_1& 1& 0\\ {\mathfrak{h}}_3& {\mathfrak{h}}_2& {\mathfrak{h}}_1\\ 0& {\mathfrak{h}}_4& {\mathfrak{h}}_3\end{array}\right|,{\Delta }_4=\left|\begin{array}{cccc}{\mathfrak{h}}_1& 1& 0& 0\\ {\mathfrak{h}}_3& {\mathfrak{h}}_2& {\mathfrak{h}}_1& 0\\ 0& {\mathfrak{h}}_4& {\mathfrak{h}}_3& {\mathfrak{h}}_2\\ 0& 0& 0& {\mathfrak{h}}_4\end{array}\right|,$$when $\mathit{Det}({\Delta }_i)>0,i=1,2,3,4$ equilibrium point E₃ to be LAS $\forall$ $\phi \in [0,1)$ and this is sufficient conditions for this point.

(ii) When $D(\mathfrak{P})>0,{\mathfrak{h}}_1>0,{\mathfrak{h}}_2<0$ and $\phi \in [\frac{2}{3},1]$, in this case the equilibrium point E₃ is unstable.

(iii) When $D(\mathfrak{P})<0,{\mathfrak{h}}_1>0,{\mathfrak{h}}_2>0,{\mathfrak{h}}_3>0,{\mathfrak{h}}_4>0$ and $\phi \in [0,\frac{1}{2})$, in this case point E₃ is LAS. Again, when $D(\mathfrak{P})<0,{\mathfrak{h}}_1<0,{\mathfrak{h}}_2>0,{\mathfrak{h}}_3<0,{\mathfrak{h}}_4>0$ in this case E₃ is unstable.

(iv) When $D(\mathfrak{P})<0,{\mathfrak{h}}_1>0,{\mathfrak{h}}_2>0,{\mathfrak{h}}_3>0,{\mathfrak{h}}_4>0$ and ${\mathfrak{h}}_2=\frac{{\mathfrak{h}}_1{\mathfrak{h}}_4}{{\mathfrak{h}}_3}+\frac{{\mathfrak{h}}_3}{{\mathfrak{h}}_1},$ in this case the point E₃ is LAS $\forall$ $\phi \in [0,1).$

(v) E₃ equilibrium point is LAS, necessary condition for this point is if ${\mathfrak{h}}_4>0.$

5. Essential uniqueness and existent for the solution for the model

Within this segment, we investigate the essential uniqueness and existent of the solution for the Atangana-Baleanu fractional model with system (9)(9) $$\begin{array}{c}{{}^{\mathit{ABC}}}_0{D}_t^{\phi }S=S\left[r\left(1-\frac{S+I}{g}\right)-\beta E\right],\\ {{}^{\mathit{ABC}}}_0{D}_t^{\phi }I=\beta SE-\nu I-\frac{\mathit{aIP}}{m+I},\\ {{}^{\mathit{ABC}}}_0{D}_t^{\phi }E=ϵI-\varrho E,\\ {{}^{\mathit{ABC}}}_0{D}_t^{\phi }P=P\left(-d+\frac{bI}{m+I}\right),\end{array}$$(9) . Applying the AB-fractional integral operator into system (9)(9) $$\begin{array}{c}{{}^{\mathit{ABC}}}_0{D}_t^{\phi }S=S\left[r\left(1-\frac{S+I}{g}\right)-\beta E\right],\\ {{}^{\mathit{ABC}}}_0{D}_t^{\phi }I=\beta SE-\nu I-\frac{\mathit{aIP}}{m+I},\\ {{}^{\mathit{ABC}}}_0{D}_t^{\phi }E=ϵI-\varrho E,\\ {{}^{\mathit{ABC}}}_0{D}_t^{\phi }P=P\left(-d+\frac{bI}{m+I}\right),\end{array}$$(9) , we obtain (12) $$\begin{array}{c}S(t)-S(0)=\frac{1-\phi }{B(\phi )}\left[S(t)\left\{r(1-\frac{S(t)+I(t)}{g})-\beta E(t)\right\}\right]\\ +\frac{\phi }{B(\phi )\Gamma (\phi )}{\int }_0^t{(t-\omega )}^{\phi -1}\left[S(\omega )\left\{r(1-\frac{S(\omega )+I(\omega )}{g})-\beta E(\omega )\right\}\right]d\omega ,\\ I(t)-I(0)=\frac{1-\phi }{B(k)}\left[\beta S(t)E(t)-\nu I(t)-\frac{aI(t)P(t)}{m+I(t)}\right]\\ +\frac{\phi }{B(\phi )\Gamma (\phi )}{\int }_0^t{(t-\omega )}^{\phi -1}\left[\beta S(\omega )E(\omega )-\nu I(\omega )-\frac{aI(\omega )P(\omega )}{m+I(\omega )}\right]d\omega ,\\ E(t)-E(0)=\frac{1-\phi }{B(\phi )}\left[ϵI(t)-\varrho E(t)\right]+\frac{\phi }{B(\phi )\Gamma (\phi )}{\int }_0^t{(t-\omega )}^{k-1}\left[ϵI(\omega )-\varrho E(\omega )\right]d\omega ,\\ P(t)-P(0)=\frac{1-\phi }{B(\phi )}\left[P(t)\left\{-d+\frac{bI(t)}{m+I(t)}\right\}\right]\\ +\frac{\phi }{B(\phi )\Gamma (\phi )}{\int }_0^t{(t-\omega )}^{\phi -1}\left[P(\omega )\left\{-d+\frac{bI(\omega )}{m+I(\omega )}\right\}\right]d\omega .\end{array}$$(12)

We define the following kernels: (13) $$\begin{array}{c}{N}_1(t,S)=S(t)\left\{r(1-\frac{S(t)+I(t)}{g})-\beta E(t)\right\},\\ N_2(t,I)=\beta S(t)E(t)-\nu I(t)-\frac{aI(t)P(t)}{m+I(t)},\\ N_3(t,E)=ϵI(t)-\varrho E(t),\\ N_4(t,P)=P(t)\left\{-d+\frac{bI(t)}{m+I(t)}\right\}.\end{array}$$(13)

Theorem 5.1.

If only the undermentioned dissimilitude holds: $$0\le r\left(1+\frac{2{\theta }_1+{\theta }_2}{g}\right)+\beta {\theta }_3<1.$$

Then the kernel N₁ convinces the Lipschitz condition and contraction mapping.

Proof 5.1.

We assume that S and S₁ are any two functions, so we get $$\begin{array}{c}\Vert N_1(t,S)-N_1(t,S_1)\Vert =\Vert S(t)\left[r\right(1-\frac{S(t)+I(t)}{g})-\beta E(t)]\\ -S_1(t)\left[r\right(1-\frac{S_1(t)+I(t)}{g})-\beta E(t)]\Vert ,\\ \le \left[r\right(1+\frac{\Vert I(t)\Vert +\Vert S(t)\Vert +\Vert S_1(t)\Vert }{g})+\beta \Vert E(t)\Vert ]\Vert S(t)-S_1(t)\Vert ,\\ ={\Upsilon }_1\Vert S(t)-S_1(t)\Vert .\end{array}$$

Let ${\Upsilon }_1=0\le r\left(1+\frac{2{\theta }_1+{\theta }_2}{g}\right)+\beta {\theta }_3<1$, ourself do suppose that S, I, E, and P are positive and bounded functions, i.e. $\Vert S(t)\Vert \le {\theta }_1,$ $\Vert I(t)\Vert \le {\theta }_2,\Vert E(t)\Vert \le {\theta }_3$ and $\Vert P(t)\Vert \le {\theta }_4,$ where θ₁, θ₂, θ₃ and θ₄ are some positive constants. (14) $$\left|\left|N_1(t,S)-N_1(t,S_1)\right|\right|\le {\Upsilon }_1\left|\left|S(t)-S_1(t)\right|\right|.$$(14)

Hence, kernel N₁ is satisfies the Lipschitz condition and since $0\le r\left(1+\frac{2{\theta }_1+{\theta }_2}{g}\right)+\beta {\theta }_3<1,N_1$ is also a contraction.

In this, similar way for the kernels N₂, N₃ and N₄ can be obtained using $(I,I_1),(E,E_1)$, and $(P,P_1)$ respectively, as follow: (15) $$\begin{array}{c}\Vert N_2(t,I)-N_2(t,I_1)\Vert \le {\Upsilon }_2\Vert I(t)-I_1(t)\Vert ,\\ \Vert N_3(t,E)-N_3(t,E_1)\Vert \le {\Upsilon }_3\Vert E(t)-E_1(t)\Vert ,\\ \Vert N_4(t,P)-N_4(t,P_1)\Vert \le {\Upsilon }_4\Vert P(t)-P_1(t)\Vert .\end{array}$$(15)

With help of equations (14)(14) $$\left|\left|N_1(t,S)-N_1(t,S_1)\right|\right|\le {\Upsilon }_1\left|\left|S(t)-S_1(t)\right|\right|.$$(14) and (15)(15) $$\begin{array}{c}\Vert N_2(t,I)-N_2(t,I_1)\Vert \le {\Upsilon }_2\Vert I(t)-I_1(t)\Vert ,\\ \Vert N_3(t,E)-N_3(t,E_1)\Vert \le {\Upsilon }_3\Vert E(t)-E_1(t)\Vert ,\\ \Vert N_4(t,P)-N_4(t,P_1)\Vert \le {\Upsilon }_4\Vert P(t)-P_1(t)\Vert .\end{array}$$(15) , we apply the above mentioned kernels, then system (12)(12) $$\begin{array}{c}S(t)-S(0)=\frac{1-\phi }{B(\phi )}\left[S(t)\left\{r(1-\frac{S(t)+I(t)}{g})-\beta E(t)\right\}\right]\\ +\frac{\phi }{B(\phi )\Gamma (\phi )}{\int }_0^t{(t-\omega )}^{\phi -1}\left[S(\omega )\left\{r(1-\frac{S(\omega )+I(\omega )}{g})-\beta E(\omega )\right\}\right]d\omega ,\\ I(t)-I(0)=\frac{1-\phi }{B(k)}\left[\beta S(t)E(t)-\nu I(t)-\frac{aI(t)P(t)}{m+I(t)}\right]\\ +\frac{\phi }{B(\phi )\Gamma (\phi )}{\int }_0^t{(t-\omega )}^{\phi -1}\left[\beta S(\omega )E(\omega )-\nu I(\omega )-\frac{aI(\omega )P(\omega )}{m+I(\omega )}\right]d\omega ,\\ E(t)-E(0)=\frac{1-\phi }{B(\phi )}\left[ϵI(t)-\varrho E(t)\right]+\frac{\phi }{B(\phi )\Gamma (\phi )}{\int }_0^t{(t-\omega )}^{k-1}\left[ϵI(\omega )-\varrho E(\omega )\right]d\omega ,\\ P(t)-P(0)=\frac{1-\phi }{B(\phi )}\left[P(t)\left\{-d+\frac{bI(t)}{m+I(t)}\right\}\right]\\ +\frac{\phi }{B(\phi )\Gamma (\phi )}{\int }_0^t{(t-\omega )}^{\phi -1}\left[P(\omega )\left\{-d+\frac{bI(\omega )}{m+I(\omega )}\right\}\right]d\omega .\end{array}$$(12) becomes (16) $$\begin{array}{c}S(t)=S(0)+\frac{1-\phi }{B(\phi )}{N}_1(t,S)+\frac{\phi }{B(\phi )\Gamma (\phi )}{\int }_0^t{(t-\omega )}^{\phi -1}{N}_1(\omega ,S)d\omega ,\\ I(t)=I(0)+\frac{1-\phi }{B(\phi )}{N}_2(t,I)+\frac{\phi }{B(\phi )\Gamma (\phi )}{\int }_0^t{(t-\omega )}^{\phi -1}{N}_2(\omega ,I)d\omega ,\\ E(t)=E(0)+\frac{1-\phi }{B(\phi )}{N}_3(t,E)+\frac{\phi }{B(\phi )\Gamma (\phi )}{\int }_0^t{(t-\omega )}^{\phi -1}{N}_3(\omega ,E)d\omega ,\\ P(t)=P(0)+\frac{1-\phi }{B(\phi )}{N}_4(t,P)+\frac{\phi }{B(\phi )\Gamma (\phi )}{\int }_0^t{(t-\omega )}^{\phi -1}{N}_4(\omega ,P)d\omega .\end{array}$$(16)

We now introduce the following recursive formulas: (17) $$\begin{array}{c}{S}_e(t)=\frac{1-\phi }{B(\phi )}{N}_1(t,S_{e-1})+\frac{\phi }{B(\phi )\Gamma (\phi )}{\int }_0^t{(t-\omega )}^{\phi -1}{N}_1(\omega ,S_{e-1})d\omega ,\\ I_e(t)=\frac{1-\phi }{B(\phi )}{N}_2(t,I_{e-1})+\frac{\phi }{B(\phi )\Gamma (\phi )}{\int }_0^t{(t-\omega )}^{\phi -1}{N}_2(\omega ,I_{e-1})d\omega ,\\ E_e(t)=\frac{1-\phi }{B(\phi )}{N}_3(t,E_{e-1})+\frac{\phi }{B(\phi )\Gamma (\phi )}{\int }_0^t{(t-\omega )}^{\phi -1}{N}_3(\omega ,E_{e-1})d\omega ,\\ P_e(t)=\frac{1-\phi }{B(\phi )}{N}_4(t,P_{e-1})+\frac{\phi }{B(\phi )\Gamma (\phi )}{\int }_0^t{(t-\omega )}^{\phi -1}{N}_4(\omega ,P_{e-1})d\omega ,\end{array}$$(17) where the initial condition are $$\begin{array}{c}{S}_0(t)=S(0),\\ I_0(t)=I(0),\\ E_0(t)=E(0),\\ P_0(t)=P(0).\end{array}$$

The differences between the consecutive terms for the recursive formulas can be written as (18) $$\begin{array}{c}{\chi }_e(t)=S_e(t)-S_{e-1}(t)\\ =\frac{1-\phi }{B(\phi )}\left[N_1(t,S_{e-1})-N_1(t,S_{e-2})\right]\\ +\frac{\phi }{B(\phi )\Gamma (\phi )}{\int }_0^t{(t-\omega )}^{\phi -1}\left[N_1(\omega ,S_{e-1})-N_1(\omega ,S_{e-2})\right]d\omega ,\end{array}$$(18) (19) $$\begin{array}{c}{\varphi }_e(t)=I_e(t)-I_{e-1}(t)\\ =\frac{1-\phi }{B(\phi )}[N_2(t,I_{e-1})-N_2(t,I_{e-2})]\\ +\frac{\phi }{B(\phi )\Gamma (\phi )}{\int }_0^t{(t-\omega )}^{\phi -1}[N_2(\omega ,I_{e-1})-N_2(\omega ,I_{e-2})]d\omega ,\end{array}$$(19) (20) $$\begin{array}{c}{\psi }_e(t)=E_e(t)-E_{e-1}(t)\\ =\frac{1-\phi }{B(\phi )}[N_3(t,E_{e-1})-N_3(t,E_{e-2})]\\ +\frac{\phi }{B(\phi )\Gamma (\phi )}{\int }_0^t{(t-\omega )}^{\phi -1}[N_3(\omega ,E_{e-1})-N_3(\omega ,E_{e-2})]d\omega ,\end{array}$$(20) (21) $$\begin{array}{c}{\xi }_e(t)=P_e(t)-P_{e-1}(t)\\ =\frac{1-\phi }{B(\phi )}[N_4(t,P_{e-1})-N_4(t,P_{e-2})]\\ +\frac{\phi }{B(\phi )\Gamma (\phi )}{\int }_0^t{(t-\omega )}^{\phi -1}{N}_4[(\omega ,P_{e-1})-N_4(\omega ,P_{e-2})]d\omega .\end{array}$$(21)

Note that: (22) $$\begin{array}{c}{S}_e(t)=\sum _{j=1}^e{\chi }_j,\\ I_e(t)=\sum _{j=1}^e{\varphi }_j,\\ E_e(t)=\sum _{j=1}^e{\psi }_j,\\ P_e(t)=\sum _{j=1}^e{\xi }_j.\end{array}$$(22)

We solve the equation (18)(18) $$\begin{array}{c}{\chi }_e(t)=S_e(t)-S_{e-1}(t)\\ =\frac{1-\phi }{B(\phi )}\left[N_1(t,S_{e-1})-N_1(t,S_{e-2})\right]\\ +\frac{\phi }{B(\phi )\Gamma (\phi )}{\int }_0^t{(t-\omega )}^{\phi -1}\left[N_1(\omega ,S_{e-1})-N_1(\omega ,S_{e-2})\right]d\omega ,\end{array}$$(18) , and the norm is applied on both sides of the equation (18)(18) $$\begin{array}{c}{\chi }_e(t)=S_e(t)-S_{e-1}(t)\\ =\frac{1-\phi }{B(\phi )}\left[N_1(t,S_{e-1})-N_1(t,S_{e-2})\right]\\ +\frac{\phi }{B(\phi )\Gamma (\phi )}{\int }_0^t{(t-\omega )}^{\phi -1}\left[N_1(\omega ,S_{e-1})-N_1(\omega ,S_{e-2})\right]d\omega ,\end{array}$$(18) , then we obtain (23) $$\begin{array}{c}\Vert {\chi }_e(t)\Vert =\Vert S_e(t)-S_{e-1}(t)\Vert \\ =\Vert \frac{1-\phi }{B(\phi )}[N_1(t,S_{e-1})-N_1(t,S_{e-2})]\\ +\frac{\phi }{B(\phi )\Gamma (\phi )}{\int }_0^t{(t-\omega )}^{\phi -1}[N_1(\omega ,S_{e-1})-N_1(\omega ,S_{e-2})]d\omega \Vert .\end{array}$$(23)

Applying the triangle inequality on equation (23)(23) $$\begin{array}{c}\Vert {\chi }_e(t)\Vert =\Vert S_e(t)-S_{e-1}(t)\Vert \\ =\Vert \frac{1-\phi }{B(\phi )}[N_1(t,S_{e-1})-N_1(t,S_{e-2})]\\ +\frac{\phi }{B(\phi )\Gamma (\phi )}{\int }_0^t{(t-\omega )}^{\phi -1}[N_1(\omega ,S_{e-1})-N_1(\omega ,S_{e-2})]d\omega \Vert .\end{array}$$(23) , we obtain (24) $$\begin{array}{c}\Vert {\chi }_e(t)\Vert =\Vert S_e(t)-S_{e-1}(t)\Vert \\ \le \frac{1-\phi }{B(\phi )}\Vert [N_1(t,S_{e-1})-N_1(t,S_{e-2})]\Vert \\ +\frac{\phi }{B(\phi )\Gamma (\phi )}\Vert {\int }_0^t{(t-\omega )}^{\phi -1}[N_1(\omega ,S_{e-1})-N_1(\omega ,S_{e-2})]d\omega \Vert .\end{array}$$(24)

So the kernel N₁ propitiate the Lipschitz condition with Lipschitz constant ${\Upsilon }_1$, we find (25) $$\begin{array}{c}\Vert {\chi }_e(t)\Vert =\Vert S_e(t)-S_{e-1}(t)\Vert \\ \le \frac{1-\phi }{B(\phi )}{\Upsilon }_1\Vert S_{e-1}-S_{e-2}\Vert \\ +\frac{\phi }{B(\phi )\Gamma (\phi )}{\Upsilon }_1{\int }_0^t{(t-\omega )}^{\phi -1}\Vert S_{e-1}-S_{e-2}\Vert d\omega ,\end{array}$$(25) thus, we obtain (26) $$\Vert {\chi }_e(t)\Vert \le \frac{1-\phi }{B(\phi )}{\Upsilon }_1\Vert {\chi }_{e-1}(t)\Vert +\frac{\phi }{B(\phi )\Gamma (\phi )}{\Upsilon }_1{\int }_0^t{(t-\omega )}^{\phi -1}\Vert {\chi }_{e-1}(\omega )\Vert d\omega .$$(26)

In a similar manner, we get (27) $$\begin{array}{c}\Vert {\varphi }_e(t)\Vert \le \frac{1-\phi }{B(\phi )}{\Upsilon }_2\Vert {\varphi }_{e-1}(t)\Vert +\frac{\phi }{B(\phi )\Gamma (\phi )}{\Upsilon }_2{\int }_0^t{(t-\omega )}^{\phi -1}\Vert {\varphi }_{e-1}(\omega )\Vert d\omega ,\\ \Vert {\psi }_e(t)\Vert \le \frac{1-\phi }{B(\phi )}{\Upsilon }_3\Vert {\psi }_{e-1}(t)\Vert +\frac{\phi }{B(\phi )\Gamma (\phi )}{\Upsilon }_3{\int }_0^t{(t-\omega )}^{\phi -1}\Vert {\psi }_{e-1}(\omega )\Vert d\omega ,\\ \Vert {\xi }_e(t)\Vert \le \frac{1-\phi }{B(\phi )}{\Upsilon }_4\Vert {\xi }_{e-1}(t)\Vert +\frac{\phi }{B(\phi )\Gamma (\phi )}{\Upsilon }_4{\int }_0^t{(t-\omega )}^{\phi -1}\Vert {\xi }_{e-1}(\omega )\Vert d\omega .\end{array}$$(27)

Theorem 5.2.

If there exists a time $t_0>0$, then following inequalities hold: (28) $$\frac{1-\phi }{B(\phi )}{\Upsilon }_1+\frac{t_0^{\phi }}{B(\phi )\Gamma (\phi )}{\Upsilon }_1<1,$$(28) then proposed model of solution exist.

Proof 5.2.

We know that $S(t),I(t),E(t)$ and P(t) are bounded and the kernels satisfy Lipschitz condition.Using equations (26)(26) $$\Vert {\chi }_e(t)\Vert \le \frac{1-\phi }{B(\phi )}{\Upsilon }_1\Vert {\chi }_{e-1}(t)\Vert +\frac{\phi }{B(\phi )\Gamma (\phi )}{\Upsilon }_1{\int }_0^t{(t-\omega )}^{\phi -1}\Vert {\chi }_{e-1}(\omega )\Vert d\omega .$$(26) and (27)(27) $$\begin{array}{c}\Vert {\varphi }_e(t)\Vert \le \frac{1-\phi }{B(\phi )}{\Upsilon }_2\Vert {\varphi }_{e-1}(t)\Vert +\frac{\phi }{B(\phi )\Gamma (\phi )}{\Upsilon }_2{\int }_0^t{(t-\omega )}^{\phi -1}\Vert {\varphi }_{e-1}(\omega )\Vert d\omega ,\\ \Vert {\psi }_e(t)\Vert \le \frac{1-\phi }{B(\phi )}{\Upsilon }_3\Vert {\psi }_{e-1}(t)\Vert +\frac{\phi }{B(\phi )\Gamma (\phi )}{\Upsilon }_3{\int }_0^t{(t-\omega )}^{\phi -1}\Vert {\psi }_{e-1}(\omega )\Vert d\omega ,\\ \Vert {\xi }_e(t)\Vert \le \frac{1-\phi }{B(\phi )}{\Upsilon }_4\Vert {\xi }_{e-1}(t)\Vert +\frac{\phi }{B(\phi )\Gamma (\phi )}{\Upsilon }_4{\int }_0^t{(t-\omega )}^{\phi -1}\Vert {\xi }_{e-1}(\omega )\Vert d\omega .\end{array}$$(27) , and applying the recursive method, we find (29) $$\begin{array}{c}\Vert {\chi }_e(t)\Vert \le \Vert S(0)\Vert {[\frac{1-\phi }{B(\phi )}{\Upsilon }_1+\frac{t^{\phi }}{B(\phi )\Gamma (\phi )}{\Upsilon }_1]}^e,\\ \Vert {\varphi }_e(t)\Vert \le \Vert I(0)\Vert {[\frac{1-\phi }{B(\phi )}{\Upsilon }_2+\frac{t^{\phi }}{B(\phi )\Gamma (\phi )}{\Upsilon }_2]}^e,\\ \Vert {\psi }_e(t)\Vert \le \Vert E(0)\Vert {[\frac{1-\phi }{B(\phi )}{\Upsilon }_3+\frac{t^{\phi }}{B(\phi )\Gamma (\phi )}{\Upsilon }_3]}^e,\\ \Vert {\xi }_e(t)\Vert \le \Vert P(0)\Vert {[\frac{1-\phi }{B(\phi )}{\Upsilon }_4+\frac{t^{\phi }}{B(\phi )\Gamma (\phi )}{\Upsilon }_4]}^e.\end{array}$$(29)

We will prove the above functions are solutions of system (9)(9) $$\begin{array}{c}{{}^{\mathit{ABC}}}_0{D}_t^{\phi }S=S\left[r\left(1-\frac{S+I}{g}\right)-\beta E\right],\\ {{}^{\mathit{ABC}}}_0{D}_t^{\phi }I=\beta SE-\nu I-\frac{\mathit{aIP}}{m+I},\\ {{}^{\mathit{ABC}}}_0{D}_t^{\phi }E=ϵI-\varrho E,\\ {{}^{\mathit{ABC}}}_0{D}_t^{\phi }P=P\left(-d+\frac{bI}{m+I}\right),\end{array}$$(9) . Let assume that (30) $$\begin{array}{c}S(t)-S(0)=S_e(t)-{\mathfrak{A}}_e(t),\\ I(t)-I(0)=I_e(t)-{\mathfrak{B}}_e(t),\\ E(t)-E(0)=E_e(t)-{\mathfrak{C}}_e(t),\\ P(t)-P(0)=P_e(t)-{\mathfrak{D}}_e(t).\end{array}$$(30)

Then, we have (31) $$\begin{array}{c}\Vert {\mathfrak{A}}_e(t)\Vert =\Vert \frac{1-\phi }{B(\phi )}[N_1(t,S)-N_1(t,S_{e-1})]\\ +\frac{\phi }{B(\phi )\Gamma (\phi )}{\int }_0^t{(t-\omega )}^{\phi -1}[N_1(\omega ,S)-N_1(\omega ,S_{e-1})]d\omega \Vert ,\\ \le \frac{1-\phi }{B(\phi )}\Vert [N_1(t,S)-N_1(t,S_{e-1})]\Vert \\ +\frac{\phi }{B(\phi )\Gamma (\phi )}{\int }_0^t{(t-\omega )}^{\phi -1}\Vert [N_1(\omega ,S)-N_1(\omega ,S_{e-1})]d\omega \Vert ,\\ \le \frac{1-\phi }{B(\phi )}{\Upsilon }_1\Vert S-S_{e-1}\Vert +\frac{t^{\phi }}{B(\phi )\Gamma (\phi )}{\Upsilon }_1\Vert S-S_{e-1}\Vert .\end{array}$$(31)

Repeating this process recursively, it becomes (32) $$\Vert {\mathfrak{A}}_e(t)\Vert \le {[\frac{1-\phi }{B(\phi )}+\frac{t^{\phi }}{B(\phi )\Gamma (\phi )}]}^{e+1}{\Upsilon }_1^{e+1}\widehat{a}.$$(32)

At the point t₀, we get (33) $$\Vert {\mathfrak{A}}_e(t)\Vert \le {[\frac{1-\phi }{B(\phi )}+\frac{t_0^{\phi }}{B(\phi )\Gamma (\phi )}]}^{e+1}{\Upsilon }_1^{e+1}\widehat{a}.$$(33)

As e tends to infinity,then limit in the equation (33)(33) $$\Vert {\mathfrak{A}}_e(t)\Vert \le {[\frac{1-\phi }{B(\phi )}+\frac{t_0^{\phi }}{B(\phi )\Gamma (\phi )}]}^{e+1}{\Upsilon }_1^{e+1}\widehat{a}.$$(33) , we get $$\Vert {\mathfrak{A}}_e(t)\Vert \to 0.$$

Similarly, we get $$\Vert {\mathfrak{B}}_e(t)\Vert \to 0,\Vert {\mathfrak{C}}_e(t)\Vert \to 0,\Vert {\mathfrak{D}}_e(t)\Vert \to 0.$$

As a result, the existence of proposed model is proven.

Theorem 5.3.

The system (9)(9) $$\begin{array}{c}{{}^{\mathit{ABC}}}_0{D}_t^{\phi }S=S\left[r\left(1-\frac{S+I}{g}\right)-\beta E\right],\\ {{}^{\mathit{ABC}}}_0{D}_t^{\phi }I=\beta SE-\nu I-\frac{\mathit{aIP}}{m+I},\\ {{}^{\mathit{ABC}}}_0{D}_t^{\phi }E=ϵI-\varrho E,\\ {{}^{\mathit{ABC}}}_0{D}_t^{\phi }P=P\left(-d+\frac{bI}{m+I}\right),\end{array}$$(9) has a unique solution if satisfies the following conditions: (34) $$\left(1-\frac{1-\phi }{B(\phi )}{\Upsilon }_j-\frac{t^{\phi }}{B(\phi )\Gamma (\phi )}{\Upsilon }_j\right)>0. \mathit{for} j=1,2,3,4.$$(34)

Proof 5.3.

Let assume that $S_1(t),I_1(t),E_1(t)$ and $P_1(t)$ be different solutions of the proposed model. We have (35) $$\begin{array}{c}S(t)-S_1(t)=\frac{1-\phi }{B(\phi )}[N_1(t,S)-N_1(t,S_1)]\\ +\frac{\phi }{B(\phi )\Gamma (\phi )}{\int }_0^t{(t-\omega )}^{\phi -1}\left[N_1(\omega ,S)-N_1(\omega ,S_1)\right]d\omega ,\end{array}$$(35) taking norm on both sides, (36) $$\begin{array}{c}\Vert S(t)-S_1(t)\Vert \le \frac{1-\phi }{B(\phi )}\Vert [N_1(t,S)-N_1(t,S_1)]\Vert \\ +\frac{\phi }{B(\phi )\Gamma (\phi )}{\int }_0^t{(t-\omega )}^{\phi -1}\Vert N_1(\omega ,S)-N_1(\omega ,S_1)\Vert d\omega ,\end{array}$$(36) the kernel satisfies Lipschitz condition, so $$\Vert S(t)-S_1(t)\Vert \le \frac{1-\phi }{B(\phi )}{\Upsilon }_1\Vert S(t)-S_1(t)\Vert +\frac{t^{\phi }}{B(\phi )\Gamma (\phi )}{\Upsilon }_1\Vert S(t)-S_1(t)\Vert ,$$then $$\Vert S(t)-S_1(t)\Vert (1-\frac{1-\phi }{B(\phi )}{\Upsilon }_1-\frac{t_{\phi }}{B(\phi )\Gamma (\phi )}{\Upsilon }_1)\le 0.$$

So $$\Vert S(t)-S_1(t)\Vert =0,$$we have (37) $$S(t)=S_1(t).$$(37)

Like this way, we obtain (38) $$\begin{array}{c}I(t)=I_1(t),\hfill \\ E(t)=E_1(t),\hfill \\ P(t)=P_1(t).\hfill \end{array}$$(38)

Therefore, the proof is completed for a unique solution of the fractional eco-epidemiological system.

6. Toufik-Atangana (TA) numerical scheme for eco-epidemiological model with Caputo operator

In this paper, authors have achieved various numerical simulations and results through the Toufik-Atangana scheme for host-parasitoid model (Kumar, Alshahrani, et al., 2021). We operate the Caputo fractional derivative in system (7)(7) $$\begin{array}{c}\frac{dS}{dt}=S\left[r\left(1-\frac{S+I}{g}\right)-\beta E\right],\\ \frac{dI}{dt}=\beta SE-\nu I-\frac{\mathit{aIP}}{m+I},\\ \frac{dE}{dt}=ϵI-\varrho E,\\ \frac{dP}{dt}=P\left(-d+\frac{bI}{m+I}\right),\end{array}$$(7) . We obtain new model in terms of Caputo sense: (39) $$\begin{array}{c}{{}^C}_0{D}_t^{\phi }S=S\left[r(1-\frac{S+I}{g})-\beta E\right],\\ {{}^C}_0{D}_t^{\phi }I=\beta SE-\nu I-\frac{\mathit{aIP}}{m+I},\\ {{}^C}_0{D}_t^{\phi }E=ϵI-\varrho E,\\ {{}^C}_0{D}_t^{\phi }P=P\left(-d+\frac{bI}{m+I}\right).\end{array}$$(39)

First, assume we have the non-integer differential equation below: (40) $$\{\begin{array}{c}{{}^C}_0{D}_t^{\phi }\mathfrak{X}(t)=F(t,\mathfrak{X}(t)),\hfill \\ \mathfrak{X}(0)={\mathfrak{X}}_0.\hfill \end{array}$$(40)

Integrate equation (40)(40) $$\{\begin{array}{c}{{}^C}_0{D}_t^{\phi }\mathfrak{X}(t)=F(t,\mathfrak{X}(t)),\hfill \\ \mathfrak{X}(0)={\mathfrak{X}}_0.\hfill \end{array}$$(40) using fractional integrals (41) $$\mathfrak{X}(t)-\mathfrak{X}(0)=\frac{1}{\Gamma (\phi )}{\int }_0^{t}F(u,\mathfrak{X}(u)){(t-u)}^{\phi -1}du,$$(41) when $t=t_{e+1},e=0,1,2,\dots$ in equation (41)(41) $$\mathfrak{X}(t)-\mathfrak{X}(0)=\frac{1}{\Gamma (\phi )}{\int }_0^{t}F(u,\mathfrak{X}(u)){(t-u)}^{\phi -1}du,$$(41) , then we get (42) $$\begin{array}{c}\mathfrak{X}(t_{e+1})-\mathfrak{X}(0)=\frac{1}{\Gamma (\phi )}{\int }_0^{t_{e+1}}F(u,\mathfrak{X}(u)){(t_{e+1}-u)}^{\phi -1}du,\\ =\frac{1}{\Gamma (\phi )}\sum _{q=0}^e{\int }_0^{t_{e+1}}F(u,\mathfrak{X}(u)){(t_{e+1}-u)}^{\phi -1}du.\end{array}$$(42)

We are estimated the function $F(u,\mathfrak{X}(u))$ with the help of the two-step Lagrange polynomial interpolation over the closed interval $[t_q,t_{q+1}]$ and we also use the polynomial interpolation concept. Thus (43) $$L_q(u)\cong \frac{F(t_q,{\mathfrak{X}}_q)}{h}(u-t_{q-1})-\frac{F(t_{q-1},{\mathfrak{X}}_{q-1})}{h}(u-t_q).$$(43)

With help of equations (42)(42) $$\begin{array}{c}\mathfrak{X}(t_{e+1})-\mathfrak{X}(0)=\frac{1}{\Gamma (\phi )}{\int }_0^{t_{e+1}}F(u,\mathfrak{X}(u)){(t_{e+1}-u)}^{\phi -1}du,\\ =\frac{1}{\Gamma (\phi )}\sum _{q=0}^e{\int }_0^{t_{e+1}}F(u,\mathfrak{X}(u)){(t_{e+1}-u)}^{\phi -1}du.\end{array}$$(42) and (43)(43) $$L_q(u)\cong \frac{F(t_q,{\mathfrak{X}}_q)}{h}(u-t_{q-1})-\frac{F(t_{q-1},{\mathfrak{X}}_{q-1})}{h}(u-t_q).$$(43) , we get (44) $$\begin{array}{c}{\mathfrak{X}}_{e+1}(t)={\mathfrak{X}}_0+\frac{1}{\Gamma (\phi )}\sum _{q=0}^e(\frac{F(t_q,{\mathfrak{X}}_q)}{h}{\int }_{t_q}^{t_{q+1}}(t-t_{q-1}){(t_{q+1}-t)}^{\phi -1}dt\\ -\frac{F(t_{q-1},{\mathfrak{X}}_{q-1})}{h}{\int }_{t_q}^{t_{q+1}}(t-t_q){(t_{q+1}-t)}^{\phi -1}dt).\end{array}$$(44)

Toward simplifying the equation (44)(44) $$\begin{array}{c}{\mathfrak{X}}_{e+1}(t)={\mathfrak{X}}_0+\frac{1}{\Gamma (\phi )}\sum _{q=0}^e(\frac{F(t_q,{\mathfrak{X}}_q)}{h}{\int }_{t_q}^{t_{q+1}}(t-t_{q-1}){(t_{q+1}-t)}^{\phi -1}dt\\ -\frac{F(t_{q-1},{\mathfrak{X}}_{q-1})}{h}{\int }_{t_q}^{t_{q+1}}(t-t_q){(t_{q+1}-t)}^{\phi -1}dt).\end{array}$$(44) , we defined the undermentioned equations (45) $$A_{\phi ,q,1}=h^{\phi +1}\frac{{(e+1-q)}^{\phi }(e-q+2+\phi )-{(e-q)}^{\phi }(e-q+2+2\phi )}{\phi (\phi +1)},$$(45) (46) $$A_{\phi ,q,2}=h^{\phi +1}\frac{{(e+1-q)}^{\phi +1}-{(e-q)}^{\phi }(e-q+1+\phi )}{\phi (\phi +1)}.$$(46)

The approximate solution of equation (44)(44) $$\begin{array}{c}{\mathfrak{X}}_{e+1}(t)={\mathfrak{X}}_0+\frac{1}{\Gamma (\phi )}\sum _{q=0}^e(\frac{F(t_q,{\mathfrak{X}}_q)}{h}{\int }_{t_q}^{t_{q+1}}(t-t_{q-1}){(t_{q+1}-t)}^{\phi -1}dt\\ -\frac{F(t_{q-1},{\mathfrak{X}}_{q-1})}{h}{\int }_{t_q}^{t_{q+1}}(t-t_q){(t_{q+1}-t)}^{\phi -1}dt).\end{array}$$(44) is obtain with help of equations (45)(45) $$A_{\phi ,q,1}=h^{\phi +1}\frac{{(e+1-q)}^{\phi }(e-q+2+\phi )-{(e-q)}^{\phi }(e-q+2+2\phi )}{\phi (\phi +1)},$$(45) and (46)(46) $$A_{\phi ,q,2}=h^{\phi +1}\frac{{(e+1-q)}^{\phi +1}-{(e-q)}^{\phi }(e-q+1+\phi )}{\phi (\phi +1)}.$$(46) , then we get (47) $$\begin{array}{c}{\mathfrak{X}}_{e+1}(t)={\mathfrak{X}}_0+\frac{1}{\Gamma (\phi )}\sum _{q=0}^e(\frac{h^{\phi }F(t_q,{\mathfrak{X}}_q)}{\phi (\phi +1)}\\ \times \left\{{(e+1-q)}^{\phi }(e-q+2+\phi )-{(e-q)}^{\phi }(e-q+2+2\phi )\right\}\\ -\frac{h^{\phi }F(t_{q-1},{\mathfrak{X}}_{q-1})}{\phi (\phi +1)}\left\{{(e+1-q)}^{\phi +1}-{(e-q)}^{\phi }(e-q+1+\phi )\right\}),\end{array}$$(47) the numerical solution for proposed model are given as (48) $$\begin{array}{c}{S}_{e+1}(t)=S(0)+\frac{1}{\Gamma (\phi )}\sum _{q=0}^e(\frac{h^{\phi }{F}_1(t_q,S_q,I_q,E_q,P_q)}{\phi (\phi +1)}\{{(e+1-q)}^{\phi }(e-q+2+\phi )\\ -{(e-q)}^{\phi }(e-q+2+2\phi )\}\\ -\frac{h^{\phi }{F}_1(t_{q-1},S_{q-1},I_{q-1},E_{q-1},P_{q-1})}{\phi (\phi +1)}\left\{{(e+1-q)}^{\phi +1}-{(e-q)}^{\phi }(e-q+1+\phi )\right\}),\end{array}$$(48) (49) $$\begin{array}{c}{I}_{e+1}(t)=I(0)+\frac{1}{\Gamma (\phi )}\sum _{q=0}^e(\frac{h^{\phi }{F}_2(t_q,S_q,I_q,E_q,P_q)}{\phi (\phi +1)}\{{(e+1-q)}^{\phi }(e-q+2+\phi )\\ -{(e-q)}^{\phi }(e-q+2+2\phi )\}\\ -\frac{h^{\phi }{F}_2(t_{q-1},S_{q-1},I_{q-1},E_{q-1},P_{q-1})}{\phi (\phi +1)}\left\{{(e+1-q)}^{\phi +1}-{(e-q)}^{\phi }(e-q+1+\phi )\right\}),\end{array}$$(49) (50) $$\begin{array}{c}{E}_{e+1}(t)=E(0)+\frac{1}{\Gamma (\phi )}\sum _{q=0}^e(\frac{h^{\phi }{F}_3(t_q,S_q,I_q,E_q,P_q)}{\phi (\phi +1)}\{{(e+1-q)}^{\phi }(e-q+2+\phi )\\ -{(e-q)}^{\phi }(e-q+2+2\phi )\}\\ -\frac{h^{\phi }{F}_3(t_{q-1},S_{q-1},I_{q-1},E_{q-1},P_{q-1})}{\phi (\phi +1)}\left\{{(e+1-q)}^{\phi +1}-{(e-q)}^{\phi }(e-q+1+\phi )\right\}),\end{array}$$(50) (51) $$\begin{array}{c}{P}_{e+1}(t)=P(0)+\frac{1}{\Gamma (\phi )}\sum _{q=0}^e(\frac{h^{\phi }{F}_4(t_q,S_q,I_q,E_q,P_q)}{\phi (\phi +1)}\{{(e+1-q)}^{\phi }(e-q+2+\phi )\\ -{(e-q)}^{\phi }(e-q+2+2\phi )\}\\ -\frac{h^{\phi }{F}_4(t_{q-1},S_{q-1},I_{q-1},E_{q-1},P_{q-1})}{\phi (\phi +1)}\left\{{(e+1-q)}^{\phi +1}-{(e-q)}^{\phi }(e-q+1+\phi )\right\}),\end{array}$$(51) where $$\begin{array}{c}{F}_1(t,S,I,E,P)=S\left[r(1-\frac{S+I}{g})-\beta E\right],\\ F_2(t,S,I,E,P)=\beta SE-\nu I-\frac{\mathit{aIP}}{m+I},\\ F_3(t,S,I,E,P)=ϵI-\varrho E,\\ F_4(t,S,I,E,P)=P\left(-d+\frac{bI}{m+I}\right).\end{array}$$

7. Toufik-Atangana (TA) numerical scheme for eco-epidemiological model with ABC-operator

Fractional derivative with the biological model is very complicated work to find solutions because of their nonlinearity (Atangana & Baleanu, 2016). So many numerical methods are used for the solution of the fractional biological model. We present the complete details of this method how to developed the method with the help of the Lagrange interpolation polynomial of two-step (Kumar, Kumar, & Jleli, 2020; Toufik & Atangana, 2017). We apply this method and find the approximate solution of the eco-epidemiological model. First, we assume the non-integer differential equation as follows: (52) $$\{\begin{array}{c}{{}^{\mathit{ABC}}}_0{D}_t^{\phi }\mathfrak{X}(t)=F(t,\mathfrak{X}(t)),\hfill \\ \mathfrak{X}(0)={\mathfrak{X}}_0.\hfill \end{array}$$(52)

Apply the Atangana-Baleanu integral operator on equation (52)(52) $$\{\begin{array}{c}{{}^{\mathit{ABC}}}_0{D}_t^{\phi }\mathfrak{X}(t)=F(t,\mathfrak{X}(t)),\hfill \\ \mathfrak{X}(0)={\mathfrak{X}}_0.\hfill \end{array}$$(52) , then we have: (53) $$\mathfrak{X}(t)-\mathfrak{X}(0)=\frac{1-\phi }{B(\phi )}F(t,\mathfrak{X}(t))+\frac{\phi }{B(\phi )\Gamma (\phi )}{\int }_0^{t}F(u,\mathfrak{X}(u)){(t-u)}^{\phi -1}du,$$(53) when $t=t_{e+1}$ = $(e+1)h,$ we have (54) $$\begin{array}{c}\mathfrak{X}(t_{e+1})-\mathfrak{X}(0)=\frac{1-\phi }{B(\phi )}F(t_e,\mathfrak{X}(t_e))+\frac{\phi }{B(\phi )\Gamma (\phi )}{\int }_0^{t_{e+1}}F(u,\mathfrak{X}(u)){(t_{e+1}-u)}^{\phi -1}du,\\ =\frac{1-\phi }{B(\phi )}F(t_e,\mathfrak{X}(t_e))+\frac{\phi }{B(\phi )\Gamma (\phi )}\sum _{q=0}^e{\int }_{t_q}^{t_{q+1}}F(u,\mathfrak{X}(u)){(t_{e+1}-u)}^{\phi -1}du.\end{array}$$(54)

We are estimated the function $F(u,\mathfrak{X}(u))$ with the help of the Lagrange polynomial over the closed interval $[t_q,t_{q+1}]$ and we also use the polynomial interpolation concept (Alzahrani, El-Dessoky, & Baleanu, 2021; Toufik & Atangana, 2017). Thus, (55) $$H(u,F(u))\cong \widehat{L}(u)=\frac{F(t_q,\mathfrak{X}(t_q))}{h}(u-t_{q-1})-\frac{F(t_{q-1},\mathfrak{X}(t_{q-1}))}{h}(u-t_q).$$(55)

Now, we simplify equation (54)(54) $$\begin{array}{c}\mathfrak{X}(t_{e+1})-\mathfrak{X}(0)=\frac{1-\phi }{B(\phi )}F(t_e,\mathfrak{X}(t_e))+\frac{\phi }{B(\phi )\Gamma (\phi )}{\int }_0^{t_{e+1}}F(u,\mathfrak{X}(u)){(t_{e+1}-u)}^{\phi -1}du,\\ =\frac{1-\phi }{B(\phi )}F(t_e,\mathfrak{X}(t_e))+\frac{\phi }{B(\phi )\Gamma (\phi )}\sum _{q=0}^e{\int }_{t_q}^{t_{q+1}}F(u,\mathfrak{X}(u)){(t_{e+1}-u)}^{\phi -1}du.\end{array}$$(54) with the help of equation (55)(55) $$H(u,F(u))\cong \widehat{L}(u)=\frac{F(t_q,\mathfrak{X}(t_q))}{h}(u-t_{q-1})-\frac{F(t_{q-1},\mathfrak{X}(t_{q-1}))}{h}(u-t_q).$$(55) , we get (56) $$\begin{array}{c}\mathfrak{X}(t_{e+1})=\mathfrak{X}(0)+\frac{1-\phi }{B(\phi )}F(t_e,\mathfrak{X}(t_e))\\ +\frac{\phi }{B(\phi )\Gamma (\phi )}\sum _{q=0}^e(\frac{F(t_q,\mathfrak{X}(t_q))}{h}{\int }_{t_q}^{t_{q+1}}(u-t_{q-1}){(t_{e+1}-u)}^{\phi -1}du\\ -\frac{F(t_{q-1},\mathfrak{X}(t_{q-1}))}{h}{\int }_{t_q}^{t_{q+1}}(u-t_q){(t_{e+1}-u)}^{\phi -1}du).\end{array}$$(56)

We simplifying the equation (56)(56) $$\begin{array}{c}\mathfrak{X}(t_{e+1})=\mathfrak{X}(0)+\frac{1-\phi }{B(\phi )}F(t_e,\mathfrak{X}(t_e))\\ +\frac{\phi }{B(\phi )\Gamma (\phi )}\sum _{q=0}^e(\frac{F(t_q,\mathfrak{X}(t_q))}{h}{\int }_{t_q}^{t_{q+1}}(u-t_{q-1}){(t_{e+1}-u)}^{\phi -1}du\\ -\frac{F(t_{q-1},\mathfrak{X}(t_{q-1}))}{h}{\int }_{t_q}^{t_{q+1}}(u-t_q){(t_{e+1}-u)}^{\phi -1}du).\end{array}$$(56) , we get (57) $$\begin{array}{c}\mathfrak{X}(t_{e+1})=\mathfrak{X}(t_0)+\frac{1-\phi }{B(\phi )}F(t_e,\mathfrak{X}(t_e))+\frac{\phi }{B(\phi )}\sum _{q=0}^e(\frac{h^{\phi }F(t_q,\mathfrak{X}(t_q))}{\Gamma (\phi +2)}\\ \times \left\{{(e-q+1)}^{\phi }(e-q+2+\phi )-{(e-q)}^{\phi }(e-q+2+2\phi )\right\}\\ -\frac{h^{\phi }F(t_{q-1},\mathfrak{X}(t_{q-1}))}{\Gamma (\phi +2)}\left\{{(e-q+1)}^{\phi +1}-{(e-q)}^{\phi }(e-q+1+\phi )\right\}).\end{array}$$(57)

We apply the above scheme (Alzahrani et al., 2021) on the system (9)(9) $$\begin{array}{c}{{}^{\mathit{ABC}}}_0{D}_t^{\phi }S=S\left[r\left(1-\frac{S+I}{g}\right)-\beta E\right],\\ {{}^{\mathit{ABC}}}_0{D}_t^{\phi }I=\beta SE-\nu I-\frac{\mathit{aIP}}{m+I},\\ {{}^{\mathit{ABC}}}_0{D}_t^{\phi }E=ϵI-\varrho E,\\ {{}^{\mathit{ABC}}}_0{D}_t^{\phi }P=P\left(-d+\frac{bI}{m+I}\right),\end{array}$$(9) , then we find $$\begin{array}{c}S(t_{e+1})=S(t_0)+\frac{1-\phi }{B(\phi )}{F}_1(t_e,\mathfrak{X}(t_e))+\frac{\phi }{B(\phi )}\sum _{q=0}^e(\frac{h^{\phi }{F}_1(t_q,\mathfrak{X}(t_q))}{\Gamma (\phi +2)}\\ \times \left\{{(e-q+1)}^{\phi }(e-q+2+\phi )-{(e-q)}^{\phi }(e-q+2+2\phi )\right\}\\ -\frac{h^{\phi }{F}_1(t_{q-1},\mathfrak{X}(t_{q-1}))}{\Gamma (\phi +2)}\left\{{(e-q+1)}^{\phi +1}-{(e-q)}^{\phi }(e-q+1+\phi )\right\}),\\ I(t_{e+1})=I(t_0)+\frac{1-\phi }{B(\phi )}{F}_2(t_e,\mathfrak{X}(t_e))+\frac{\phi }{B(\phi )}\sum _{q=0}^e(\frac{h^{\phi }{F}_2(t_q,\mathfrak{X}(t_q))}{\Gamma (\phi +2)}\\ \times \left\{{(e-q+1)}^{\phi }(e-q+2+\phi )-{(e-q)}^{\phi }(e-q+2+2\phi )\right\}\\ -\frac{h^{\phi }{F}_2(t_{q-1},\mathfrak{X}(t_{q-1}))}{\Gamma (\phi +2)}\left\{{(e-q+1)}^{\phi +1}-{(e-q)}^{\phi }(e-q+1+\phi )\right\}),\end{array}$$ (58) $$\begin{array}{c}E(t_{e+1})=E(t_0)+\frac{1-\phi }{B(\phi )}{F}_3(t_e,\mathfrak{X}(t_e))+\frac{\phi }{B(\phi )}\sum _{q=0}^e(\frac{h^{\phi }{F}_3(t_q,\mathfrak{X}(t_q))}{\Gamma (\phi +2)}\\ \times \left\{{(e-q+1)}^{\phi }(e-q+2+\phi )-{(e-q)}^{\phi }(e-q+2+2\phi )\right\}\\ -\frac{h^{\phi }{F}_3(t_{q-1},\mathfrak{X}(t_{q-1}))}{\Gamma (\phi +2)}\left\{{(e-q+1)}^{\phi +1}-{(e-q)}^{\phi }(e-q+1+\phi )\right\}),\\ P(t_{e+1})=P(t_0)+\frac{1-\phi }{B(\phi )}{F}_4(t_e,\mathfrak{X}(t_e))+\frac{\phi }{B(\phi )}\sum _{q=0}^e(\frac{h^{\phi }{F}_4(t_q,\mathfrak{X}(t_q))}{\Gamma (\phi +2)}\\ \times \left\{{(e-q+1)}^{\phi }(e-q+2+\phi )-{(e-q)}^{\phi }(e-q+2+2\phi )\right\}\\ -\frac{h^{\phi }{F}_4(t_{q-1},\mathfrak{X}(t_{q-1}))}{\Gamma (\phi +2)}\left\{{(e-q+1)}^{\phi +1}-{(e-q)}^{\phi }(e-q+1+\phi )\right\}).\end{array}$$(58)

8. Numerical results and graphical analysis

Here, we present a new numerical technique with Caputo and ABC non-integer operators to simulate the non-linear fractional eco-epidemiological model. Throughout this simulation, we display graphic results for different fractional orders and infection rates. Monitoring and controlling the infection rate is crucial for determining the system’s behavior and operation. I followed the instructions and accomplished the task of the numerical scheme detailed in section 7. Simulations have been conducted using version R2016a (MATLAB 9.0) of MATLAB (The MathWorks Inc, 2016) to demonstrate the efficiency of the proposed approach. We used MATLAB-9 to analyse numerical results for different fractional orders and time durations. Through the simulations, we can observe how changes in the parameters and initial conditions impact the model’s predictions. It gave me a better understanding of the model’s dynamics and allowed me to conduct a more comprehensive analysis. The graphical results from the fractional order analysis were more informative and generalized than other relative works. Model parameters used for the simulation are as follows $$r=4.6667,\nu =1.5,m=1,a=1,g=7,ϵ=1,d=0.5,\varrho =1,b=1,$$

S(0) =1/3, I(0) =1, E(0) =1 and P(0) =1/3 are our initial values of proposed system.

In epidemiological modeling, infection rates provide a crucial insight into the dynamic of the system. Thus, we have taken the different infection rates and order of the derivative to observe different behaviors of the PP model. Based on classical models $\phi =1$ and β = 2, the system exhibits stable behavior in Figures 1 and 2. Figures 1 and 2 depict the three-dimensional and two-dimensional stable phase portrait representation of Caputo eco-epidemiological model (8)(8) $$\begin{array}{c}{{}^C}_0{D}_t^{\phi }S=S\left[r\left(1-\frac{S+I}{g}\right)-\beta E\right],\\ {{}^C}_0{D}_t^{\phi }I=\beta SE-\nu I-\frac{\mathit{aIP}}{m+I},\\ {{}^C}_0{D}_t^{\phi }E=ϵI-\varrho E,\\ {{}^C}_0{D}_t^{\phi }P=P\left(-d+\frac{bI}{m+I}\right),\end{array}$$(8) when the order of the derivative is $\phi =1$ and β = 2. Furthermore, in Figures 3 and 4, if we change the value of the fractional order $(\phi =0.99)$ of the epidemiological model, it is quickly stable compared to the previous classical model. Also, Figures 3 and 4 depict the three-dimensional and two-dimensional stable phase portrait representation of Atangana-Baleanu (AB) eco-epidemiological model (9)(9) $$\begin{array}{c}{{}^{\mathit{ABC}}}_0{D}_t^{\phi }S=S\left[r\left(1-\frac{S+I}{g}\right)-\beta E\right],\\ {{}^{\mathit{ABC}}}_0{D}_t^{\phi }I=\beta SE-\nu I-\frac{\mathit{aIP}}{m+I},\\ {{}^{\mathit{ABC}}}_0{D}_t^{\phi }E=ϵI-\varrho E,\\ {{}^{\mathit{ABC}}}_0{D}_t^{\phi }P=P\left(-d+\frac{bI}{m+I}\right),\end{array}$$(9) when the order of the derivative is $\phi =0.99$ and β = 2. When the orders of fractional derivatives are different, $\phi =1,0.99,0.97,0.95,0.93.$ Figures 5–7 illustrate the dynamic nature and complexity of the eco-epidemiological (8)(8) $$\begin{array}{c}{{}^C}_0{D}_t^{\phi }S=S\left[r\left(1-\frac{S+I}{g}\right)-\beta E\right],\\ {{}^C}_0{D}_t^{\phi }I=\beta SE-\nu I-\frac{\mathit{aIP}}{m+I},\\ {{}^C}_0{D}_t^{\phi }E=ϵI-\varrho E,\\ {{}^C}_0{D}_t^{\phi }P=P\left(-d+\frac{bI}{m+I}\right),\end{array}$$(8) based on the Caputo Toufic-Atangana (TA) Scheme and we observe that if we decrease the order of the derivative, the model quickly goes to the stable stage. During the simulation, we noticed that the susceptible animal population appeared to be fluctuating at different rates for different values of $\phi ,$ and eventually reached a stable equilibrium state. The infected population, on the other hand, showed a rapid increase and then a decrease, followed by a gradual decline for varying fractional orders. The amount of infectious material in the environment also followed a similar pattern, rising quickly and then gradually falling over time, with different paths for different fractional orders. Lastly, the predator density showed a rapid decrease and then a gradual decline, eventually stabilizing towards the end of the simulation period. Figures 8–10 show the comparative results relating to the eco-epidemiological model of Caputo and Atangana-Baleanu (AB) with $\beta =2.5$ and $\phi =1.$ We have examined the graphs and observed that when $\phi =1,$ the Caputo and AB derivatives models yield the same results with identical parameter values and initial conditions. However, for fractional orders, the proposed fractional models display different paths. We have also noticed significant deviations in the efficacy of both the Caputo and AB models when using the same parameter values, which is caused by the memory properties of the kernels involved in describing the fractional operators. Figures 11–13 depict the phase diagram for eco-epidemiological model (8)(8) $$\begin{array}{c}{{}^C}_0{D}_t^{\phi }S=S\left[r\left(1-\frac{S+I}{g}\right)-\beta E\right],\\ {{}^C}_0{D}_t^{\phi }I=\beta SE-\nu I-\frac{\mathit{aIP}}{m+I},\\ {{}^C}_0{D}_t^{\phi }E=ϵI-\varrho E,\\ {{}^C}_0{D}_t^{\phi }P=P\left(-d+\frac{bI}{m+I}\right),\end{array}$$(8) with $\phi =1$ and various $\beta =2,2.07,2.5$ respectively by Caputo Toufik-Atangana (TA) scheme. In the proposed model, it is essential to note that β is a sensitive parameter. Even a slight change in this parameter has a drastic impact on the behaviour of the system. Further, we have analysed the variation of parameter values with the state variable of the proposed model. Figures 14–16 represent the bifurcation diagram of the considered eco-epidemiological model for analysed the complex behaviour and stability analysis. Based on the findings, it has become apparent that the Caputo and AB approach is the most significant and reliable method for explaining physical processes, compared to the regular fractional and classical-order cases.

Figure 1. Three-dimensional numerical results for Caputo eco-epidemiological model (8)(8) $$\begin{array}{c}{{}^C}_0{D}_t^{\phi }S=S\left[r\left(1-\frac{S+I}{g}\right)-\beta E\right],\\ {{}^C}_0{D}_t^{\phi }I=\beta SE-\nu I-\frac{\mathit{aIP}}{m+I},\\ {{}^C}_0{D}_t^{\phi }E=ϵI-\varrho E,\\ {{}^C}_0{D}_t^{\phi }P=P\left(-d+\frac{bI}{m+I}\right),\end{array}$$(8) with β = 2 and $\phi =1.$

Figure 2. Two-dimensional numerical results for Caputo eco-epidemiological model (8)(8) $$\begin{array}{c}{{}^C}_0{D}_t^{\phi }S=S\left[r\left(1-\frac{S+I}{g}\right)-\beta E\right],\\ {{}^C}_0{D}_t^{\phi }I=\beta SE-\nu I-\frac{\mathit{aIP}}{m+I},\\ {{}^C}_0{D}_t^{\phi }E=ϵI-\varrho E,\\ {{}^C}_0{D}_t^{\phi }P=P\left(-d+\frac{bI}{m+I}\right),\end{array}$$(8) with β = 2 and $\phi =1.$

Figure 3. Three-dimensional numerical results for Atangana-Baleanu (AB) eco-epidemiological model (9)(9) $$\begin{array}{c}{{}^{\mathit{ABC}}}_0{D}_t^{\phi }S=S\left[r\left(1-\frac{S+I}{g}\right)-\beta E\right],\\ {{}^{\mathit{ABC}}}_0{D}_t^{\phi }I=\beta SE-\nu I-\frac{\mathit{aIP}}{m+I},\\ {{}^{\mathit{ABC}}}_0{D}_t^{\phi }E=ϵI-\varrho E,\\ {{}^{\mathit{ABC}}}_0{D}_t^{\phi }P=P\left(-d+\frac{bI}{m+I}\right),\end{array}$$(9) with β = 2 and $\phi =0.99.$

Figure 4. Two-dimensional numerical results for Atangana-Baleanu (AB) eco-epidemiological model (9)(9) $$\begin{array}{c}{{}^{\mathit{ABC}}}_0{D}_t^{\phi }S=S\left[r\left(1-\frac{S+I}{g}\right)-\beta E\right],\\ {{}^{\mathit{ABC}}}_0{D}_t^{\phi }I=\beta SE-\nu I-\frac{\mathit{aIP}}{m+I},\\ {{}^{\mathit{ABC}}}_0{D}_t^{\phi }E=ϵI-\varrho E,\\ {{}^{\mathit{ABC}}}_0{D}_t^{\phi }P=P\left(-d+\frac{bI}{m+I}\right),\end{array}$$(9) with β = 2 and $\phi =0.99.$

Figure 5. Three-dimensional numerical results for eco-epidemiological model (8)(8) $$\begin{array}{c}{{}^C}_0{D}_t^{\phi }S=S\left[r\left(1-\frac{S+I}{g}\right)-\beta E\right],\\ {{}^C}_0{D}_t^{\phi }I=\beta SE-\nu I-\frac{\mathit{aIP}}{m+I},\\ {{}^C}_0{D}_t^{\phi }E=ϵI-\varrho E,\\ {{}^C}_0{D}_t^{\phi }P=P\left(-d+\frac{bI}{m+I}\right),\end{array}$$(8) with $\beta =2.5$ and various fractional orders by Caputo Toufik-Atangana (TA) scheme.

Figure 6. Two-dimensional numerical results for eco-epidemiological model (8)(8) $$\begin{array}{c}{{}^C}_0{D}_t^{\phi }S=S\left[r\left(1-\frac{S+I}{g}\right)-\beta E\right],\\ {{}^C}_0{D}_t^{\phi }I=\beta SE-\nu I-\frac{\mathit{aIP}}{m+I},\\ {{}^C}_0{D}_t^{\phi }E=ϵI-\varrho E,\\ {{}^C}_0{D}_t^{\phi }P=P\left(-d+\frac{bI}{m+I}\right),\end{array}$$(8) with $\beta =2.5$ and various fractional orders by Caputo Toufik-Atangana (TA) scheme.

Figure 7. Time series numerical results for eco-epidemiological model (8)(8) $$\begin{array}{c}{{}^C}_0{D}_t^{\phi }S=S\left[r\left(1-\frac{S+I}{g}\right)-\beta E\right],\\ {{}^C}_0{D}_t^{\phi }I=\beta SE-\nu I-\frac{\mathit{aIP}}{m+I},\\ {{}^C}_0{D}_t^{\phi }E=ϵI-\varrho E,\\ {{}^C}_0{D}_t^{\phi }P=P\left(-d+\frac{bI}{m+I}\right),\end{array}$$(8) with $\beta =2.5$ and various fractional orders by Caputo Toufik-Atangana (TA) scheme.

Figure 8. Comparison of Three-dimensional numerical results for Caputo and Atangana-Baleanu (AB) eco-epidemiological model with $\beta =2.5$ and $\phi =1.$

Figure 9. Comparison of two-dimensional numerical results for Caputo and Atangana-Baleanu (AB) eco-epidemiological model with $\beta =2.5$ and $\phi =1.$

Figure 10. Comparison of time series numerical results for Caputo and Atangana-Baleanu (AB) eco-epidemiological model with $\beta =2.5$ and $\phi =1.$

Figure 11. Three-dimensional phase diagram for eco-epidemiological model (8)(8) $$\begin{array}{c}{{}^C}_0{D}_t^{\phi }S=S\left[r\left(1-\frac{S+I}{g}\right)-\beta E\right],\\ {{}^C}_0{D}_t^{\phi }I=\beta SE-\nu I-\frac{\mathit{aIP}}{m+I},\\ {{}^C}_0{D}_t^{\phi }E=ϵI-\varrho E,\\ {{}^C}_0{D}_t^{\phi }P=P\left(-d+\frac{bI}{m+I}\right),\end{array}$$(8) with $\phi =1$ and various $\beta =2,2.07,2.5$ respectively by Caputo Toufik-Atangana (TA) scheme.

Figure 12. Two-dimensional phase diagram for eco-epidemiological model (8)(8) $$\begin{array}{c}{{}^C}_0{D}_t^{\phi }S=S\left[r\left(1-\frac{S+I}{g}\right)-\beta E\right],\\ {{}^C}_0{D}_t^{\phi }I=\beta SE-\nu I-\frac{\mathit{aIP}}{m+I},\\ {{}^C}_0{D}_t^{\phi }E=ϵI-\varrho E,\\ {{}^C}_0{D}_t^{\phi }P=P\left(-d+\frac{bI}{m+I}\right),\end{array}$$(8) with $\phi =1$ and various $\beta =2,2.07,2.5$ respectively by Caputo Toufik-Atangana (TA) scheme.

Figure 13. Time series diagram for eco-epidemiological model (8)(8) $$\begin{array}{c}{{}^C}_0{D}_t^{\phi }S=S\left[r\left(1-\frac{S+I}{g}\right)-\beta E\right],\\ {{}^C}_0{D}_t^{\phi }I=\beta SE-\nu I-\frac{\mathit{aIP}}{m+I},\\ {{}^C}_0{D}_t^{\phi }E=ϵI-\varrho E,\\ {{}^C}_0{D}_t^{\phi }P=P\left(-d+\frac{bI}{m+I}\right),\end{array}$$(8) with $\phi =1$ and various $\beta =2,2.07,2.5$ respectively by Caputo Toufik-Atangana (TA) scheme.

Figure 14. Bifurcation plots of the proposed eco-epedimological model (7)(7) $$\begin{array}{c}\frac{dS}{dt}=S\left[r\left(1-\frac{S+I}{g}\right)-\beta E\right],\\ \frac{dI}{dt}=\beta SE-\nu I-\frac{\mathit{aIP}}{m+I},\\ \frac{dE}{dt}=ϵI-\varrho E,\\ \frac{dP}{dt}=P\left(-d+\frac{bI}{m+I}\right),\end{array}$$(7) : (a) β versus S (b) ε versus S.

Figure 15. Bifurcation plots of the proposed eco-epedimological model (7)(7) $$\begin{array}{c}\frac{dS}{dt}=S\left[r\left(1-\frac{S+I}{g}\right)-\beta E\right],\\ \frac{dI}{dt}=\beta SE-\nu I-\frac{\mathit{aIP}}{m+I},\\ \frac{dE}{dt}=ϵI-\varrho E,\\ \frac{dP}{dt}=P\left(-d+\frac{bI}{m+I}\right),\end{array}$$(7) : (a) r versus S (b) g versus S.

Figure 16. Bifurcation plots of the proposed eco-epedimological model (7)(7) $$\begin{array}{c}\frac{dS}{dt}=S\left[r\left(1-\frac{S+I}{g}\right)-\beta E\right],\\ \frac{dI}{dt}=\beta SE-\nu I-\frac{\mathit{aIP}}{m+I},\\ \frac{dE}{dt}=ϵI-\varrho E,\\ \frac{dP}{dt}=P\left(-d+\frac{bI}{m+I}\right),\end{array}$$(7) : (a) ν versus S (b) a versus S.

9. Conclusion

In this article, we have analyzed the eco-epidemic model that explains the spread of this disease within an animal community and how the environment influences this model. We have converted the eco-epidemiological CWD integer order model presented into a fractional order model. With the aid of fixed-point theory, we have proven the existence and uniqueness of solutions to the proposed model under certain conditions. We have investigated the stability of all equilibrium points of the non-integer eco-epidemiological CWD model with the help of Matignon’s condition. With the support of the Lagrange interpolation polynomial of two steps, we have used two novel numerical methods for locating the approximate solution. We have examined the behaviour of the proposed model by bifurcation plots for various parameter values. Finally, we have compared the Atangana-Baleanu and Caputo solutions of the eco-epidemiological CWD model, the graphical results give better information about the model and its solutions. Our study’s findings play a crucial role in improving the accuracy of the CWD model and developing effective strategies. In future work, we can employ other non-integer order operators and more effective numerical techniques for this model to learn more about the behaviour of the CWD model. These efforts will undoubtedly advance our understanding of this disease and help us combat it more effectively.

Disclosure statement

The authors declare that they have no conflicts of interests.

Data availability statement

Data will be made available on reasonable request.