# The Dynamics of a Delayed Ecoepidemiological Model with Nonlinear Incidence Rate

Authors: Reem Mudar Hussien; Raid Kamel Naji
Journal: Journal of Applied Mathematics
 (2023)

---

## Abstract

In this paper, the general framework for calculating the stability of equilibria, Hopf bifurcation of a delayed prey-predator system with an SI type of disease in the prey population, is investigated. The impact of the incubation period delay on disease transmission utilizing a nonlinear incidence rate was taken into account. For the purpose of explaining the predation process, a modified Holling type II functional response was used. First, the existence, uniform boundedness, and positivity of the solutions of the considered model system, along with the behavior of equilibria and the existence of Hopf bifurcation, are studied. The critical values of the delay parameter for which stability switches and the nature of the Hopf bifurcation by using normal form theory and center manifold theorem are identified. Additionally, using numerical simulations and a hypothetical dataset, various dynamic characteristics are discovered, including stability switches, chaos, and Hopf bifurcation scenarios.


---

## Body

## 1. Introduction

Due to its universality and importance in population dynamics, interspecies interaction is frequent in the real-world ecosystem. How effective an organism consumes or uses its food resource determines whether it will survive. As a result of these interactions, each predator directly impacts the ecology in question, including the number of prey. The population typically has a saturation effect in real life. Predator functional response to prey population, which defines the amount of prey consumed per predator per unit of time and plays a major role in the stability and bifurcation dynamics of the underlying system, is thus the key component in prey-predator interaction. Over the past few decades, several functional responses have been created [1]. In natural life, the constant environment is a rare occurrence. Several factors, like mating habits, food availability, seasonal weather effects, harvesting, death and birth rates, infectious illness rates, and other significant population rates, have an impact on the natural environment. Infectious disease is the main biotic hazard to the entire prey-predator system, both in terms of suffering and social and economic implications. As a result, the ecosystem is constrained by various abiotic and biotic elements. For many years, mathematical biologists have been attempting to combine the two main fields of study ecology and epidemiology. This merging of infection into the ecological model is known as ecoepidemiological [2]. There have been numerous studies into ecoepidemiological models, which are described as mathematical representations of ecological systems where the disease affects prey and/or predators. Disease in the predator, disease in the prey, and disease in both groups can all be taken into account when examining the impact of disease on an environment. When the disease is present in the prey, the predator may consume both susceptible and diseased animals. There are several studies on prey-predator dynamics, disease, and infection, as well as any potential ecological and biological effects, as stated in [2–12].A well-known truth is that time delays exist in every biological process; for example, many diseases usually pass through a number of stages during their life cycle. So, the infection process is not instantaneous and may sometimes be known as an incubation delay. Time delay has a tremendous impact on system dynamics. Several ecoepidemiological models involving delay in different factors have been proposed and studied (see [13–21]). Additionally, the dynamics of disease are significantly influenced by the force of infection. The majority of earlier ecoepidemiological models relied on the mass action law, frequency dependence, and saturating disease transmission in populations of prey or predators. This resulted from numerous interactions with infected people prior to the disease’s outbreak. As a result, a three-dimensional nonlinear prey-predator model with a nonlinear saturated incidence rate, time delay, and Holling type II functional response is taken into consideration in this article. Various diseases are conveyed across people through contact, contaminated food ingestion, or environmental disease vectors. The primary goals of the current study are to determine how disease dynamics in the prey population may be impacted by the delay time and how the presence of disease in the prey population impacts the model’s behavior.The following is how the paper is set up: Section2 develops a mathematical model and establishes the positivity and boundedness of the solutions. The presence of potential equilibrium points is established in Section 3. Section 4 examined the occurrence of Hopf bifurcation and the equilibrium points’ local asymptotic stability. The bifurcating periodic solution’s stability and direction are established in Section 5. The numerical system findings are provided and cross-checked with the analytical results in Section 6. In Section 7, a brief conclusion is made.


## 2. Model Construction

In this section, a mathematical model that describes an ecoepidemiological system incorporating time delay is proposed and studied. The conventional Lotka-Volterra prey-predator model is composed of two species and is written as(1)dxdt=rx1−xK−α1xy,dydt=α2xy−μ2y,where the positive parameters r, K, α1, α2, and μ2 represent the prey growth rate, the carrying capacity, the attack rate, the conversion rate of food from prey to predator, and the mortality rate. The dynamic of system (1) was studied by many authors (see, for example, [22, 23] and the references therein). However, the dynamics of a prey-predator model with Holling type II schemes is written in the form:
(2)dxdt=rx1−xK−α1xyβ+x,dydt=α2xyβ+x−μ2y,where β>0 represents the half-saturation constant, while the rest of the parameters have the same meaning as above. The literature has extensively examined the dynamics of system (2) (for instance, see [24, 25] and the references therein). The Holling type II prey-predator model with SI-type disease in the prey was further examined by Naji and Mustafa [4]. They used the nonlinear incidence rate to describe the transition of disease in their ecoepidemiological prey-predator model. They assumed that the predator eats both susceptible and infected prey by the modified Holling type II functional response, which led to the following model:
(3)dSdt=rS1−SK−δI1+γIS−α1SYβ+S+mI,dIdt=δI1+γIS−α2IYβ+S+mI−μ1I,dYdt=θ1S+θ2Iβ+S+mIY−μ2Y.Here, the prey population is divided into two compartments,St and It, which stand for the susceptible and infected populations at time t, respectively. As of time t, Yt represents the predator population. The sum St+It indicates the entire biomass of the prey population, which is obviously considered in the model above to be owing to disease infection. Only susceptible prey is capable of reproducing, per the logistic rule with intrinsic growth rate r and carrying capacity K, but the infected population perishes before developing the ability for reproduction. The disease does not spread among prey populations through genetic transmission, but rather through close contact with infected prey, as indicated by the incidence rate δI/1+γI [26], where δ is the contact rate and γ is the inhibition effect resulting from the contact between the infected individuals themselves.Model (3) simulates a real-world system that adheres to the following presumptions in light of the aforementioned:
(i)
Susceptible and diseased compartments are created within the prey. The nonlinear incidence rate indicates that the disease is spread through contact(ii)
There are no indications of immunological development or recovery in the infected population. Additionally, it degrades exponentially with the death rate (μ1), which takes into account both infectious diseases and natural causes(iii)
Because they only consume readily accessible food, predators are unable to distinguish between healthy and unhealthy populations. In order to consume both susceptible and infected prey, they use modified Holling type II functional responsesα1S/β+S+mI and α2I/β+S+mI, where α1 and α2 are the attack rates of the susceptible and infected prey, respectively, and m is the preferred rate for the infected due to their weakness. The predator, on the other hand, grows at rates of θ1 and θ2 and, in the absence of food, dies at a rate of μ2 (predator mortality rate)(iv)
Since most diseases do not actually spread immediately away whenSt and It come into contact, they typically take some time to touch-transfer to the newly transferred individuals. The transmission term uses a delay time as the most realistic factorKeeping the above in mind, model (3) is therefore adjusted to include the delay term in the incidence rate, and the resulting system with delay is as follows:
(4)dSdt=rS1−SK−δIt−τ1+γIt−τS−α1SYβ+S+mI,dIdt=δIt−τ1+γIt−τS−α2IYβ+S+mI−μ1I,dYdt=θ1S+θ2Iβ+S+mIY−μ2Y.Here, the incubation periodτ∈ℝ+ is when the infectious agent develops in the host, and only after that does the infected prey become infectious. Therefore, the number of actively infected prey at the time t is arising from the contacts of the actual population of susceptible and infected prey at time t−τ, where τ is the discrete-time delay.All the system parameters are assumed to be positive, and system (4) initial conditions are
(5)Sθ=φ1θ,Iθ=φ2θ,Yθ=φ3θ;θ∈−τ,0;φiθ>0,i=1,2,3,where φ=φ1,φ2,φ3∈C−τ,0,ℝ+3 with ℝ+3=S,I,Y∈ℝ3:S≥0,I≥0,Y≥0. Here, C represents the Banach space of all continuous functions φ with the norm φ=Sup−τ≤θ≤0φ1θ,φ2θ,φ3θ.Now the qualitative behavior of the solution of system (4) is investigated in the following theorem.Theorem 1.
All the solutions of system (4) with initial conditions (5) are positive and bounded in the interior of ℝ+3.Proof.
Since the functions in the right-hand side of system (4) are completely continuous and have continuous partial derivatives, then they are locally Lipschitzian on C. Hence, the solution St,It,Yt of system (4) with initial conditions (5) exists and is unique on 0,α where 0<α<∞ [27].
Now, from system (4), it is obtained that
(6)St=S0exp∫0tr1−SvK−δIv−τ1+γIv−τ−α1Yvβ+Sv+mIvdv,It=I0exp∫0tδIv−τIv1+γIv−τS−α2Yvβ+Sv+mIv−μ1dv,Yt=Y0exp∫0tθ1Sv+θ2Ivβ+Sv+mIv−μ2dv.
Therefore, due to the positivity of initial conditions (5) and the exponential function, all the solutions of system (4) are positive. Now, to prove the boundedness of the solutions of system (4) for all t≥0, consider St,It,Yt as any solution of system (4) with their given initial conditions. From the first equation, it is deduced directly that
(7)limt⟶∞SupSt≤K.
LetM=St+It+Yt; then, system (4) gives that
(8)dMdt=rS1−SK−α1−θ1SYβ+S+mI−α2−θ2IYβ+S+mI−μ1I−μ2Y.
Since the predator growth rate constants cannot be exceeding the predator attack rate constants, hence, from a biological point of view,θi≤αi,i=1,2 is always true. Then
(9)dMdt≤rS−μ1I−μ2Y≤θ−μM,where μ=minμ1,μ2 and θ=r+μK. By solving the above linear differential inequality, it is obtained that limt⟶∞Mt≤θ/μ, which completes the proof.


## 3. The Existence of the Equilibrium Points

In this section, the existence of various equilibrium points (EPs) is discussed, and then, the basic reproduction number is determined. System (2) has the following nonnegative equilibrium points.The trivial equilibrium point (TEP) and disease-predator-free equilibrium point (DPFEP) unconditionally exist and are denoted to them asP0=0,0,0 and P1=S0,0,0 with S0=K, respectively.The predator-free equilibrium point (PFEP) is represented byP2=S¯,I¯,0, where
(10)S¯=μ11+γI¯δ,while I¯ is determined as a positive root of the following equation:
(11)A1I¯2+A2I¯+A3=0,where A1=μ1rγ2>0, A2=2μ1γr+Kδ2−rγKδ, and A3=rμ1−Kδ.Applying Descartes’ rule of signs [28], there exists a unique positive real root I¯ of equation (11) that is given by
(12)I¯=−A2+A22−4A1A32A1,if and only if the following condition is met:
(13)μ1<Kδ.The infected prey-free equilibrium point (IPYFEP), which is denoted byP3=S~,0,Y~, is determined as
(14)S~=μ2βθ1−μ2,Y~=rβ+S~K−S~α1K.Clearly, the IPYFEP exists if the following conditions are met:(15)μ2<θ1,S~<K.The coexistence equilibrium point (CEP) that is represented byP4=S∗,I∗,Y∗ is obtained by solving the following set of algebraic equations:
(16)rα11−SK−δIα11+γI∗=Yβ+S+mI,δSα21+γI−μ1α2=Yβ+S+mI,θ1S+θ2Iβ+S+mI−μ2=0.Direct computation gives that(17)S=μ2β+μ2m−θ2Iθ1−μ2.Now, using the first two equations of system (16) with equation (17) gives the following quadratic equation:
(18)fI=B1I2+B2I+B3=0,where
(19)B1=rγα2μ2m−θ2,B2=rγα2μ2β+Kθ1−μ2α2δ−α2rγ−α1γμ1+μ2m−θ2rα2+α1Kδ,B3=rα2μ2β−Kθ1−μ2+α1Kδβμ2−μ1θ1−μ2.Applying Descartes’ rule of signs on equation (18) shows that it has a unique positive real root I∗ if and only if B1 and B3 have opposite signs. Then, using the obtained positive root in equation (17) gives that S∗=SI∗. Straightforward computation shows that S∗,I∗ is positive if and only if one of the following sets of conditions is met:
(20a)θ2m<μ2<θ1,β<Kθ1−μ2μ2.rα2+α1μ1rα2+α1Kδ,or
(20b)θ1<μ2<θ2m,μ2βθ2−μ2m<I∗.Finally, substituting the resulting positive value ofS∗,I∗ in the first equation of system (16) gives that
(21)Y∗=β+S∗+mI∗α1K1+γI∗rK−S∗1+γI∗−δKI∗.Consequently, it is easy to verify that CEP exists uniquely in the interior ofℝ+3; if and only if in addition to either (20a) or (20b), the following condition is met:
(22)0<δKI∗K−S∗1+γI∗<r.


## 4. Local Stability of Equilibrium Points

LetP=S¯,I¯,Y¯ be an arbitrary equilibrium point of system (4). Let u1t=St−S¯, u2t=It−I¯, and u3t=Yt−Y¯, respectively. As a result, after removing the bar, system (4) becomes equivalent to the following linear system below:
(23)du1dt=q11St+q12It+q13Yt+p12It−τ,du2dt=q21St+q22It+q23Yt+p22It−τ,du3dt=q31St+q32It+q33Yt,where
(24)q11=r−2rSK−δI1+γI−α1Yβ+mIβ+S+mI2,q12=α1mSYβ+S+mI2,q13=−α1Sβ+S+mI,q21=δI1+γI+α2IYβ+S+mI2,q22=−α2βY+α2SYβ+S+mI2−μ1,q23=−α2Iβ+S+mI,q31=Yβθ1+θ1m−θ2Iβ+S+mI2,q32=Yβθ2+θ2−θ1mSβ+S+mI2,q33=θ1S+θ2Iβ+S+mI−μ2,p12=−δS1+γI2,p22=δS1+γI2.Then, the characteristic equation of system (4) at P can be determined by
(25)q11−λq12+p12e−λτq13q21q22+p22e−λτ−λq23q31q32q33−λ=0.The derived local stability analysis of system (4) around its equilibrium points is presented in the subsequent theorems.Theorem 2.
The TEP of system (4) is a saddle point for all τ≥0.Proof.
SubstitutingP0=0,0,0 in general characteristic equation (25) gives that
(26)r−λ000−μ1−λ000−μ2−λ=0.
Therefore, the eigenvalues of system (4) at TEP are
(27)λ01=r>0,λ02=−μ1<0,λ03=−μ2<0.
The TEP becomes a saddle point for allτ≥0 when there is a positive eigenvalue that is independent of τ. This shows that the populations cannot go extinct at the same time because there is always a positive eigenvalue.Now before the next theorems are given, the following proposition that is given in [19] is presented.Proposition 3.
Suppose thatC1>0 and C2>0; then, the following is obtained:
(i)
IfC1<C2, then all the roots of ρ+C1−C2e−ρτ=0 have positive real parts for τ<1/C22−C12cos−1C1/C2(ii)
IfC1>C2, then all the roots of ρ+C1−C2e−ρτ=0 have negative real parts for any τTheorem 4.
The DPFEP of system (4) is locally asymptotically stable for all τ≥0 if the following conditions are met:
(28)δS0<μ1,(29)θ1S0β+S0<μ2.
However, it is an unstable point for all0≤τ<τc=1/δS02−μ12cos−1μ1/δS0, if condition (28) is reflected.Proof.
According to linear system (23) and their characteristic equation (25) at the DPFEP, the characteristic equation becomes as follows:
(30)−r−λδS0e−λτ−μ1−λθ1S0β+S0−μ2−λ=0.
Clearly, equation (30) has two explicit real eigenvalues given by λ11=−r<0 and λ12=θ1S0/β+S0−μ2. It is easy to verify that the second eigenvalue λ12 is negative due to condition (29).
On the other hand, the rest of roots of equation (30) follows from the transcendental equation:
(31)δS0e−λτ−μ1−λ=0.
Without a doubt, forτ=0, equation (31) has a singular root given by λ13=δS0−μ1, which is always negative due to condition (28) and positive when condition (28) is reflected. Therefore, when condition (28) is satisfied or reflected, the DPFEP is locally asymptotically stable for τ=0 or a saddle point (28).
Now forτ>0, due to a proposition (3), all the roots of equation (31) have negative real part roots or positive real part roots for τ<τc=1/δS02−μ12cos−1μ1/δS0 when condition (28) is satisfied or reflected, respectively. Consequently, the DPFEP is locally asymptotically stable for all τ>0 if conditions (28) and (29) are met, while it is an unstable point for 0<τ<τc when condition (28) is reflected. This completes the proof.Theorem 5.
The PFEP of system (4) is locally asymptotically stable for all τ≥0 if the following conditions are met:
(32)θ1S¯+θ2I¯<μ2β+S¯+mI¯,(33)r<2rS¯K+δI¯1+γI¯.Proof.
According to linear system (23) and their characteristic equation (25) at the PFEP, the characteristic equation becomes as follows:
(34)λ2−c11+c22λ+c11c22+e−λτ−d22λ+c11d22−c21d12c33−λ=0,where
(35)c11=r1−2S¯K−δI¯1+γI¯,c13=0,c13=−α1S¯β+S¯+mI¯,c21=δI¯1+γI¯,c22=−μ1,c23=−α2I¯β+S¯+mI¯,c31=c32=0,c33=θ1S¯+θ2I¯β+S¯+mI¯−μ2,d12=−δS¯1+γI¯2,d22=δS¯1+γI¯2.Clearly, equation (36) has one explicit real eigenvalue given by λ21=c33, which is negative under condition (32). However, the other eigenvalues are obtained from the equation:
(36)λ2−c11+c22λ+c11c22+e−λτ−d22λ+c11d22−c21d12=0.Forτ=0, it is easy to verify that equation (36) has two eigenvalues with negative real part if and only if condition (33) holds.On the other hand, forτ>0, the eigenvalues of system (23) at PFEP are determined from equation (36). So to look whether the system is locally asymptotically stable or not and undergoes a Hopf bifurcation near PFEP, it is assumed that λ=μ+iσ where μ=0, σ>0.Therefore, direct computation by substitutingλ=iσ in equation (36) and then separate the real and imaginary parts gives that
(37)−σ2+c11c22+cosστ.c11d22−c21d12−d22σsinστ=0,(38)−c11+c22σ−d22σcosστ−sinστ.c11d22−c21d12=0.By squaring equations (37) and (38) and adding to each other, it is obtained that
(39)σ4+c11+c222−2c11c22+d222c11d22−c21d122σ2+c11c222=0.Or equivalently, whenρ=σ2,
(40)ρ2+c112+c222+d222c11d22−c21d122ρ+c11c222=0.According to the discard rule of sign, equation (40) has no positive root. Hence, there is no σ>0, for which λ=iσ satisfies equation (36). Therefore, all the eigenvalues of equation (36) do not intersect the imaginary axis as τ>0 and still have negative real parts. Thus, the PFEP is asymptotically stable for all τ≥0 if the given conditions are met.Theorem 6.
The IPYFEP of system (4) is locally asymptotically stable for all τ≥0 if the following conditions are met:
(41)r<2rS~K+α1βY~β+S~2,(42)θ1S~β+S~<μ2,(43)δS~<α2Y~β+S~+μ1.
However, it is an unstable point for allτ<τcc, if condition (43) is reflected.Proof.
According to linear system (23) and their characteristic equation (25) at the IPYFEP, the characteristic equation becomes as follows:
(44)e22+h22e−λτ−λλ2−e11+e33λ+e11e33−e13e31=0,where
(45)e11=r−2rS~K−α1βY~β+S~2,e12=α1mS~Y~β+S~2,e13=−α1S~β+S~,e21=e23=0,e22=−α2Y~β+S~−μ1,e31=βθ1Y~β+S~2,e32=θ2β+S~−θ1mS~β+S~2Y~,e33=θ1S~β+S~−μ2,h12=−δS~,h22=δS~.According the Routh-Hurwitz criterion, it is easy to verify that equation (43) has two eigenvalues with negative real parts that are given by
(46)λ31=e11+e33+e11+e332−4e11e33−e13e312,λ32=e11+e33−e11+e332−4e11e33−e13e312,if and only if conditions (41) and (42) are met.However, the other eigenvalues of equation (44) are determined from the equation:
(47)λ−e22−h22e−λτ=0⇔λ+H1−H2e−λτ=0.Clearly,H1=−e22>0 and H2=h22>0. Not that, for τ=0, equation (47) has a unique negative root λ=e22+h22 provided that condition (43) holds. Therefore, the IPYFEP is locally asymptotically stable.Now, forτ>0, an application to the proposition (3) shows that all roots of the equation (47) have negative real parts under condition (43) and they are positive real parts for all τ<τcc if condition (43) is reflected, where
(48)τcc=1δS~2−α2Y~/β+S~+μ12cos−1α2Y~/β+S~+μ1δS~.Therefore, the IPYFEP is locally asymptotically stable for allτ≥0 if conditions (41)–(43) are met, and it is unstable for all τ<τcc if condition (3) is reflected.Now, to study the local stability of the CEP of delayed system (4), it is assumed that aij=qijS∗,I∗,Y∗ and bij=pijS∗,I∗,Y∗ in linear system (23). Then, the characteristic equation becomes as follows:
(49)Pλ+Qλe−λτ=0,where
(50)Pλ=λ3+b1λ2+b2λ+b3,Qλ=b4λ2+b5λ+b6,b1=−a11+a22+a33,b2=a22a11+a33−a12a21−a23a32+a11a33−a13a31,b3=−a11a22a33−a23a32−a12a21a33−a23a31+a13a21a32−a22a31,b4=−b22,b5=a11b22−b12a21+b22a33,b6=−a11a33b22−a13a31b22−b12a21a33−a23a31.Theorem 7.
The CEP of system (4) is locally asymptotically stable for any τ=0, if and only if the following conditions are satisfied:
(51)b1+b4>0,(52)b3+b6>0,(53)b1+b4b2+b5−b3+b6>0.Proof.
Forτ=0, equation (49) reduces to
(54)λ3+b1+b4λ2+b2+b5λ+b3+b6=0.
If the requirements (51)–(53) are satisfied, it is determined by the well-known Routh–Hurwitz criterion that system (4) with τ=0 is locally asymptotically stable at CEP. The CEP becomes unstable if these conditions do not meet this requirement.Theorem 8.
If conditions (51)–(53) hold, then there is τ0>0, such that the CEP of system (4) is locally asymptotically stable for all τ<τ0, and it is unstable for τ>τ0 provided that the following condition holds:
(55)b32<b62.
Finally, system (4) undergoes a Hopf bifurcation at τ=τ0 provided that
(56)3β022+2Φ1β02+Φ2≠0,where all the new symbols are defined in the proof.Proof.
From theorem (11), the CEP is locally asymptotically stable under conditions (51)–(53) at τ=0; hence, all the roots of the characteristic equation (49) at τ=0 have negative real parts. Moreover, by Butler’s lemma [29], P4 remains stable for τ<τ0 for a specific value τ0. A characteristic root of (49) must intersect the imaginary axis if instability occurs for a specific value of the delay, according to corollary 2.4 in Ruan and Wei [30]. Assume that (49) has a root that is entirely imaginary, iβ, with β>0. Then, in (49), by separating the real and imaginary components, it is obtained that
(57)b6−b4β2cosβτ+b5βsinβτ=b1β2−b3,(58)b4β2−b6sinβτ+b5βcosβτ=β3−b2β.
By squaring and adding (57) and (58) with eliminatingτ, the following algebraic equation ofβis derived:
(59)β6+Φ1β4+Φ2β2+Φ3=0,where Φ1=b12−b42−2b2,Φ2=2b4b6−b52−2b1b3+b22, and Φ3=b32−b62.
Substitutingβ2=Z in (59) gives the following third order equation:
(60)fZ=Z3+Φ1Z2+Φ2Z+Φ3=0.
Clearly, condition (55) guarantees that equation (60) will include at least one positive root denoted by β0=Z, which leads to instability for P4. In the following, the stability’s change of P4 in relation to τ is determined.
From equations (56) and (58), the following is reached:
(61)sinβτ=β3−b2βb4β2−b6−b5βb3−b1β2b4β2−b62+b5β2,cosβτ=b5ββ3−b2β+b3−b1β2b4β2−b6b4β2−b62+b5β2.
Then,τn corresponding to β0 can be obtained as
(62)τn=1β0arccosb5ββ3−b2β+b3−b1β2b4β2−b6b4β2−b62+b5β2+2nπ,n=0,1,2,⋯.
Define(63)τ0=minτn,n≥0.
Now, to complete the proof of the theorem, the occurrence of Hopf bifurcation atτ=τ0 will be proved.From the above, it is proved that the characteristic equation (49) of system (4) has a pair of complex conjugate roots, given by λ=ατ±iβτ, which are pure imaginary at τ=τ0 that is ατ0=0. Therefore, the proof follows if we can prove the transversality criterion [31].
(64)dReλτdττ=τ0≠0.Taking the derivative of equation (49) with respect to τ gives that
(65)dλdτ−1=3λ2+2b1λ+b2eλτλb4λ2+b5λ+b6+2b4λ+b5λb4λ2+b5λ+b6−τλ,substituting λ=iβ0 into equation (65) yields
(66)ddτReλλ=iβ0=Redλdτ−1λ=iβ0=Reb2−3β03+2ib1β0cosβ0τ+isinβ0τ−b5β02+iβ0b6−b4β02+Re2ib4β0+b5−b5β02+iβ0b6−b4β02=β023β04+2b12−4b2−2b42β02+b22−2b1b3−b52+2b4b6b5β022+b6−b4β022=f′β02b5β022+b6−b4β022,where f′β02=3β022+2Φ1β02+Φ2≠0 due to condition (56).So, we can get(67)SignddτReλτ=τ0=SignRedλdττ=τ0−1=Signf′β02.Thus,Signd/dτReλτ=τ0≠0, so the Hopf bifurcation takes place at τ=τ0.If we assume thatd/dτReλτ=τ0<0, then the characteristic equation has roots with positive real parts when τ=τ0. It is in contradiction with the local stability of the positive equilibrium point. Hence, d/dτReλτ=τ0>0.


## 5. Stability and Direction of the Hopf Bifurcation

The investigation above revealed the circumstances under which system (4) experiences a Hopf bifurcation at CEP as the delay parameter moves through the value τ0. This section examines the direction of Hopf bifurcation and necessary conditions for the stability of the bifurcating periodic solution of system (4) using Hassard’s center manifold theorem and normal form theory [31]. The following results are attained for system (4).Theorem 9.
Under the following determined quantities, the stability and direction of bifurcating periodic solution can be specified as follows:(1)
IfM2>0M2<0, then the Hopf bifurcation is supercritical (subcritical) and the bifurcating periodic solutions exist for τ>τ0τ<τ0(2)
Ifν2<0ν2>0, then the Hopf bifurcation is stable (unstable)(3)
IfT2>0T2<0, then the period of the bifurcating periodic solutions is increased (decreased), where M2, ν2, and T2 are given as follows:(68)C10=i2β0τ0g11g20−2g112−g0223+g212,M2=−ReC10Redλ/dττ0,ν2=2ReC10,T2=−ImC10+M2Imdλ/dττ0β0τ0.Proof.
Letu1t=St−S∗, u2t=It−I∗, u3t=Yt−Y∗, and τ=τ0+a, where τ0 is define by equation (63) and a∈ℝ, so a=0 is the Hopf bifurcation value of system (4). Rescaling the time delay t⟶t/τ, then system (4) can be transformed into a functional differential equation in C=C−1,0,ℝ3 as follows:
(69)u′t=Laut+fa,ut,where ut=u1t,u2t,u3tT∈C=C−1,0,ℝ3 and La:C⟶ℝ3, f:ℝ×C⟶ℝ3 with
(70)Laϕ=τ0+aMϕ0+Nϕ−1,where
(71)M=f10001f01001f00101f10002f01002f00102f1003f0103f0013=qij3×3,N=0f0001100f000120000=0p12e−λτ00p22e−λτ0000,with qij and p12,p22 are given in equation (23), while the nonlinear term is given by
(72)fa,ϕ=τ0+aF1F2F3,where
(73)F1=∑i+l+j+k≥21i!l!j!k!filjk1ϕ1i0ϕ2l0ϕ3j0ϕ~2k−1,F2=∑i+l+j+k≥21i!l!j!k!filjk2ϕ1i0ϕ2l0ϕ3j0ϕ~2k−1,F3=∑i+l+j≥21i!l!j!filj3ϕ1i0ϕ2l0ϕ3j0,ϕθ=ϕ1θ,ϕ2θ,ϕ3θ∈C,−1≤θ≤0,with
(74)filjk1ϕ1i0ϕ2l0ϕ3j0ϕ~2k−1=∂i+l+j+kf1∂ϕ1iϕ2lϕ3jϕ~2kϕ1,ϕ2,ϕ3,ϕ~2=0,0,0,−1,filjk2ϕ1i0ϕ2l0ϕ3j0ϕ~2k−1=∂i+l+j+kf2∂ϕ1iϕ2lϕ3jϕ~2kϕ1,ϕ2,ϕ3,ϕ~2=0,0,0,−1,filj3ϕ1i0ϕ2l0ϕ3j0=∂i+l+jf3∂ϕ1iϕ2lϕ3jϕ1,ϕ2,ϕ3=0,0,0.
Moreover, direct computation gives the following higher derivatives:(75)f11001=α1mY∗β+mI∗+3α1mY∗S∗β+S∗+mI∗3,f10101=−α1β+mI∗β+S∗+mI∗2,f10011=−δ1+γI∗2,f01101=α1S∗mβ+S∗+mI∗2,f01011=0,f00111=0,f20001=−rK−α1Y∗β+S∗+mI∗−α1Y∗β+mI∗β+S∗+mI∗3,f02001=α1m2S∗Y∗β+S∗+mI∗3,f00201=0,f00021=12δS∗−γ2δS∗I∗γI∗+21+γI∗4,f11002=−α2Y∗β+S∗−α2mY∗I∗β+S∗+mI∗3,f10102=−α2I∗β+S∗+mI∗2,f10012=δ1+γI∗2,f01102=−α2β+S∗β+S∗+mI∗2,f01012=0,f00112=0,f20002=−α2I∗Y∗β+S∗+mI∗3,f02002=α2Y∗+2α2I∗Y∗m2β+S∗+mI∗3,f00202=0,f00022=12γδS∗γI∗−21+γI∗4,f1103=mθ2Y∗β+S∗+mI∗2−θ2Y∗+2mY∗θ1S∗+θ2I∗β+S∗+mI∗β+S∗+mI∗4,f1013=θ1β+θ1mI∗−θ2I∗β+S∗+mI∗2,f0113=θ2β+S∗−θ1mS∗β+S∗+mI∗2,f0203=12−mθ2Y∗+2m2Y∗θ1S∗+θ2I∗β+S∗+mI∗3,f0023=0,f2003=θ1Y∗β+S∗+mI∗2−θ1Y∗−2θ1Y∗β+S∗+mI∗+2θ1S∗Y∗β+S∗+mI∗−2θ2I∗Y∗β+S∗+mI∗β+S∗+mI∗4.
According to the Riesz representation theorem, there is a3×3 matrix function ψv,a,v∈−1,0, such that
(76)Laϕ=∫−10dψv,aϕv,ϕ∈C.
In fact, it can choose(77)ψv,a=τ0+aMδv−Nδv+1,where δ is the Dirac delta function.
Forϕ∈C−1,0,ℝ3, it is defined that
(78)Aaϕv=dϕvdv,v∈−1,0,∫−10dψσ,aϕσ,v=0,Raϕv=0,v∈−1,0,fa,ϕ,v=0.
Thus, system (70) is equivalent to the abstract operator differential equation
(79)u′t=Aaut+Raut,where ut=ut+v,v∈−1,0.
Now, forφ∈C1−1,0,ℝ3∗, the adjoint operater A∗ of A0 is defined as
(80)A∗φσ=−dφσdσ,σ∈0,1,∫−10φ−tdϕTt,0,σ=0,and a bilinear form
(81)φσ,ϕv=φ¯0ϕ0−∫v=−10∫∈=0vφ¯∈−vdψvϕ∈d∈,where ψv=ψv,0. Clearly, A0 and A∗ are adjoint operators. Thus, for a=0, by a simple computation, we can calculate qv=1,q1,q2Teiβ0τ0v be the eigenvector of A0 belonging to the eigenvalue iβ0τ0, and q∗σ=D1,q1∗,q2∗Te−iβ0τ0σ is the eigenvector of A∗ that associated with the eigenvalue −iβ0τ0, where
(82)q1=−f0013−iβ0q2−f1003f0103,q2=f10001f0103−iβ0f0103−f01001f1003−f00011f1003e−iτ0β0f01001f0013−iβ0f01001+f00011f0013e−iτ0β0−iβ0f00011e−iτ0β0−f00101,q1∗=−f10001−iβ0−f1003q2∗f10002,q2∗=f00102f10001+f00102iβ0+f1003f00102−f00101f10002f0013f10002+iβ0f10002.
From bilinear inner product (81), we get
(83)q∗σ,qv=D¯1+q¯1∗q1+q¯2∗q2+q1τ0e−iβ0τ0f00011+q1q¯1∗τ0e−iβ0τ0f00012.
LetD=1+q1∗q¯1+q2∗q¯2+τ0e−iβ0τ0q¯1f00011+q1∗f00012−1, where D¯ is the conjugate complex number of D, then q∗,q=1 and q∗,q¯=0.
In the following, using similar arguments as in [5], it can be determined the properties of the Hopf bifurcation:
(84)g20=2τ0D¯P1+P5q¯1∗+P9q¯2∗,g11=τ0D¯P2+P6q¯1∗+P10q¯2∗,g02=2τ0D¯P3+P7q¯1∗+P11q¯2∗,g21=2τ0D¯P4+P8q¯1∗+P12q¯2∗,where
(85)P1=q1f11001+q2f10101+q¯1eiβ0τ0f10011+q1q2f01101+12f20001+12q12f02001+12q12e−2iβ0τ0f00021,P2=q1+q¯1f11001+f10011+q2+q¯2f10101+q1q¯2+q¯1q2f01101+f20001+q1q¯1f02001+q1q¯1f00021,P3=q¯1f11001+q¯2f10101+q¯1eiβ0τ0f10011+q1q¯2f01101+12f20001+12q12f02001+12q12e2iβ0τ0f00021,P4=12q¯1w2010+q1w1110+w1120+12w2020f11001+q2w1110+12q¯2w2010+12w2030+w1130f10101+12q¯1eiβ0τ0w2010+q1e−iβ0τ0w1110+w112−1+12w202−1f10011+q2w1120+12q¯2w2020+12q¯1w2030+q1w1130f01101+12w2010+2w1110f20001+q¯1w2010+2q1w1110f02001+12q¯1eiβ0τ0w202−1+q1e−iβ0τ0w112−1f00021,P5=q1f11002+q2f10102+q¯1eiβ0τ0f10012+q1q2f01102+12f20002+12q12f02002+12q12e−2iβ0τ0f00022,P6=q1+q¯1f11002+f10012+q2+q¯2f10102+q1q¯2+q¯1q2f01102+f20002+q1q¯1f02002+q1q¯1f00022,P7=q¯1f11002+q¯2f10102+q¯1eiβ0τ0f10012+q1q¯2f01102+12f20002+12q12f02002+12q12e2iβ0τ0f00022,P8=12q¯1w2010+q1w1110+w1120+12w2020f11002+q2w1110+12q¯2w2010+12w2030+w1130f10102+12q¯1eiβ0τ0w2010+q1e−iβ0τ0w1110+w112−1+12w202−1f10012+q2w1120+12q¯2w2020+12q¯1w2030+q1w1130f01102+12w2010+2w1110f20002+q¯1w2010+2q1w1110f02002+12q¯1eiβ0τ0w202−1+q1e−iβ0τ0w112−1f00022,P9=q1f1103+q2f1013+q1q2f0113+12f2003+12q12f0203,P10=q1+q¯1f1103+q2+q¯2f1013+q1q¯2+q¯1q2f0113+f2003+q¯1q1f0203,P11=q¯1f1103+q¯2f1013+q¯1q¯2f0113+12f2003+q¯12f0203,P12=12q¯1w2010+q1w1110+w1120+12w2020f1103+q2w1110+12q¯2w2010+12w2030+w1130f1013+q2w1120+12q¯2w2020+12q¯1w2030+q1w1130f0113+w2010+2w1110f2003+q¯1w2010+2q1w1110f0203,with
(86)w20v=ig20β0τ0q0eiβ0τ0v+ig¯023β0τ0q¯0e−iβ0τ0v+U1e2iβ0τ0v,(87)w11v=−ig11β0τ0q0eiβ0τ0v+ig¯11β0τ0q¯0e−iβ0τ0v+U2.
Notice thatU1=U11,U12,U13T∈ℝ3 and U2=U21,U22,U23T∈ℝ3 are constant vectors that can be determined from the following equations:
(88)J1∗U1=2J1,J2∗U2=−J2,where
(89)J1∗=2iβ0τ0I−∫−10dψve2iβ0τ0v,J2∗=∫−10dψv,J1=P1P5P9T,J2=P2P6P10T.
Accordingly, it is obtained that(90)J1∗=2iβ0−f10001−f01001−f00011e−2iβ0τ0v−f00101−f100022iβ0−f01002−f00012e−2iβ0τ0v−f00102−f1003−f01032iβ0−f0013,J2∗=−f10001−f01001−f00011−f00101−f10002−f01002−f00012−f00102−f1003−f0103−f0013.
Thus,U1i=2△i/△,i=1,2,3, where △=DetJ1∗ and △i is the value of the determinant Xi, where Xi is formed by replacing the ith column vector of J1∗ by J1 for i=1,2,3. Similarly, U2i=2△¯i/△¯,i=1,2, and 3, where △¯=DetJ2∗ and △¯i is the value of the determinant Yi, where Yi is formed by replacing the ith column vector of J2∗ by J2 for i=1,2,3.
Consequently,w20v and w11v can be computed using equations (86)–(87). Then, the expressions given in equation (68) can be determined depending on those given in equation (84), and the proof is done.


## 6. Numerical Simulation

In this section, to study the effects of different parameters on the model dynamics and visualize the complex dynamics of system (4), some numerical simulations are carried out. MATLAB 2021A was used to run all simulations. Additionally, to run our simulations, the following set of biologically feasible hypothetical parameter values is used:
(91)r=0.5,K=10,δ=0.75,γ=1,α1=2,α2=1.5,θ1=0.4,θ2=0.4,m=2,β=4,μ1=0.15,μ2=0.1,τ=0.5.For dataset (91), the trajectories of system (4) starting from different initial values approach asymptotically to the CEP that is given by P4=0.58,1.12,0.25 as shown in Figure 1.Figure 1
The trajectories of system (4) utilizing dataset (91) with different sets of initial values. (a) Phase portrait approaches P4. (b) The time series.
(a)(b)According to Figure1, dataset (91) satisfies Theorems 7 and 8 and that τ=0.5<τ0 at which the system undergoes a Hopf bifurcation. Therefore, in Figure 2, the bifurcation point τ0 is shown.Figure 2
The trajectories of system (4) utilizing dataset (91) with different values of delay parameter τ. (a, b) The trajectory at τ=2.57 approaches to P4. (c, d) The trajectory at τ0=2.578 undergoes Hopf bifurcation. (e, f) The trajectory at τ=2.75 approaches to 3D periodic with large period. (g, h) The trajectory at τ=15 approaches to 2D periodic in the SY-plane. (i, j) The trajectory at τ=45 approaches to chaos.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)According to Figure2, system (4) undergoes a Hopf bifurcation at τ0=2.578; after that, the periodic dynamics became larger in the range τ∈2.578,5.5. It is observed that for the range τ∈5.6,30, system (4) approaches to periodic dynamic in the SY-plane; however, it returns to periodic for the range τ∈31,40, while it has a chaotic dynamic in the range τ∈41,54 and then return to periodic and so on. Therefore, the delay has a destabilizing effect on the dynamic of system (4).Now, the effects of different parameters on the model dynamics are investigated and presented in the following figures. Figure3 shows the influence of varying the parameter r, while Figure 4 shows the influence of varying the parameter δ.Figure 3
The trajectories of system (4) utilizing dataset (91) with different values of the parameter r. (a, b) The trajectory at r=0.25 approaches to P2=0.29,0.47,0. (c, d) The trajectory at r=0.78 approaches to 3D periodic. (e, f) The trajectory at r=1.5 approaches to 2D periodic in the SY-plane.
(a)(b)(c)(d)(e)(f)Figure 4
The trajectories of system (4) utilizing dataset (91) with different values of the parameter δ. (a, b) The trajectory at δ=0.1 approaches to 2D periodic in the SY-plane. (c, d) The trajectory at δ=0.3 approaches to 3D periodic. (e, f) The trajectory at δ=0.6 approaches to P4=0.86,0.69,0.65. (g, h) The trajectory at δ=0.9 approaches to P2=0.35,1.15,0.
(a)(b)(c)(d)(e)(f)(g)(h)Obviously, the parameterr has a destablizing effect on the system’s dynamics up to a critical value, and then, system (4) loses its persistence and approaches to periodic attractor in the SY-plane. However, the parameter δ has a stabilizing effect on system’s (4) dynamics up to a vital value, and then, system (4) loses its persistence and approaches to P2 in the SI-plane.In Figures5 and 6, the influence of the parameters γ and α1 is investigated, respectively.Figure 5
The trajectories of system (4) utilizing dataset (91) with different values of the parameter γ. (a, b) The trajectory at γ=0.15 approaches to 2D periodic in the SI-plane. (c, d) The trajectory at γ=0.5 approaches to P2=0.29,0.95,0.
(a)(b)(c)(d)Figure 6
The trajectories of system (4) utilizing dataset (91) with different values of the parameter α1. (a, b) The trajectory at α1=0.5 approaches to 2D periodic in the SY-plane. (c, d) The trajectory at α1=0.77 approaches to 3D periodic.
(a)(b)(c)(d)From Figures5 and 6, it is observed that both the parameters γ and α1 have positive effect on the persistence and stability of system (4), as the increasing of these parameters the system dynamic transfer from periodic in the boundary planes to the asymptotic at a coexistence point. Moreover, it is observed that the parameter α2 has opposite influence on the dynamic of system (4) as that obtained for α1, while the parameter θ1 has similar influence on the dynamic of system (4) as that obtained for r.Now, Figures7–9 explain the influence of the parameters θ2,μ1, and μ2, respectively.Figure 7
The trajectories of system (4) utilizing dataset (91) with different values of the parameter θ2. (a, b) The trajectory at θ2=0.3 approaches to P2=0.54,1.7,0. (c, d) The trajectory at θ2=0.6 approaches to P4=0.62,0.53,0.59. (e, f) The trajectory at θ2=0.65 approaches to 3D periodic.
(a)(b)(c)(d)(e)(f)Figure 8
The trajectories of system (4) utilizing dataset (91) with different values of the parameter μ1. (a, b) The trajectory at μ1=0.5 approaches to P4=1.18,0.22,0.85. (c, d) The trajectory at μ1=0.55 approaches to 3D periodic. (e, f) The trajectory at μ1=0.75 approaches to 2D periodic in the SY-plane.
(a)(b)(c)(d)(e)(f)Figure 9
The trajectories of system (4) utilizing dataset (91) with different values of the parameter μ2. (a, b) The trajectory at μ2=0.05 approaches to 2D periodic in the SY-plane. (c, d) The trajectory at μ2=0.07 approaches to 3D periodic. (e, f) The trajectory at μ2=0.12 approaches to P2=0.54,1.7,0.
(a)(b)(c)(d)(e)(f)From Figures7–9, it is observed that the parameter θ2 has a destablizing effect on the dynamic of system (4). However, the parameter μ1 destablizes system (4) first, and then after a critical value, the system loses its persistence and approaches to a periodic dynamic in the SY-plane, while the parameter μ2 stablizes system (4) first and then loses its persistence and approaches to P2 in the SI-plane.It is observed that the parametersm and β have similar influence on the dynamic of system (4) as that obtained for μ2. Finally, the influence of the parameter K on the dynamic of system (4) is shown in Figure 10.Figure 10
The trajectories of system (4) utilizing dataset (91) with different values of K and utilizing different sets of initial values. (a, b) The trajectory approaches to P1=S0,0,0 with S0=0.15 when K=0.15. (c, d) The trajectory approaches to P2=0.42,1.1,0 when K=0.2.
(a)(b)(c)(d)According to Figure10, dataset (91) satisfies the Theorem 4 in the range K∈0,0.2, while system (4) with set (91) satisfies Theorem 5 in the range K∈0.21,3.1. Otherwise, it approaches to the P4. So, the parameter K has a stablizing influence on system (4).


## 7. Conclusion

In this work, a three-species ecoepidemiological model that describes the delay in the transmission of the disease within the prey population has been proposed and studied. The prey population was divided into two groups: susceptible and infected. Modified Holling type II functional response was used for describing the predation process, while the nonlinear type of incidence rate was used to describe disease transmission. All the properties of the solution including the positivity and boundedness solution were discussed. It was obtained that system (4) has at most five feasible equilibrium points. The existence and stability conditions for boundary equilibrium points were established. The analysis near the CEP showed that it is locally asymptotically stable for all 0<τ<τ0, while it is an unstable otherwise and a Hopf bifurcation occurs at τ=τ0. Furthermore, an explicit algorithm and a necessary criteria for the stability and direction of bifurcating periodic solutions were constructed utilizing the normal form method and center manifold theorem.According to the simulation study using dataset (91), the effects of different parameters on the model dynamics were investigated. The obtained outcomes can be summarized as follows.System (4) undergoes a Hopf bifurcation at τ0=2.578; it has a 3D periodic dynamics, 2D periodic dynamics, and even chaotic dynamics for all values τ0>2.578. Accordingly, the delay factor has a destabizing influence on the dynamics of system (4) as the delay factor increases.Decreasing the intrinsic growth rate of the susceptible prey below a certain value (r=0.45) leads to extinction in the predator species, and the solution is settled at the PFEP. However, increasing this parameter over a certain value (r=0.773), the CEP loses its stability and system (4) approaches 3D periodic (Hopf bifurcation occurs) first, and then, as the intrinsic growth rate values increase further, system (4) loses their persistence and the solution approaches 2D periodic dynamics in the SY-plane.Decreasing the contact rate below theδ=0.5 destabilizes the CEP, and system (4) approaches asymptotically to 3D periodic dynamics (Hopf bifurcation occurs) while when the contact rate in the range δ<0.3, the infected population faces extinction and the system approaches asymptotically to 2D periodic dynamics in the SY-plane. On the other hand, increasing the contact rate above the vital value δ=0.82, the predator population faces extinction and system (4) approaches asymptotically to the PFEP. Consequently, the parameter δ has a stabilizing effect on system’s (4) dynamics up to a vital value.Now, decreasing the inhibition rate so that it belongs to the rangeγ∈02,0.83, the predator goes extinction and the system approaches asymptotically PFEP. However, for γ<0.2, the PFEP loses its stability and a Hopf bifurcation occurs in the interior of SI-plane. Moreover, the system approaches asymptotically PFEP for γ>0.83. This indicates the inhibition rate’s stabilizing effect on system (4). On the other hand, decreasing the attack rate of the predator to susceptible prey causes destabilization of system (4), so that the CEP becomes unstable and a Hopf bifurcation occurs at α1=0.77, and then, for α1≤0.76, system (4) approaches asymptotically to a periodic dynamic in the interior of SY-plane. Moreover, the system still persists at the CEP otherwise.It has been noted that the attack rate of the predator on infected prey has an effect on the dynamics of system (4) that is opposite to that of the attack rate on susceptible prey. But just as with the intrinsic growth rate of the susceptible prey, the conversion rate of the biomass of the susceptible prey into the predator has an impact on the dynamics of system (4).Finally, the simulation results show that decreasing the conversion rate of the biomass of the infected prey into the predator causes extinction in the predator population, and system (4) approaches asymptotically to PFEP while increasing this rate destabilizing the CEP and Hopf bifurcation occurs in the interior of the first octant. However, increasing the death rate of the infected population leads to destabilizing the CEP first, and a Hopf bifurcation appears in the range μ1∈0.55,0.65 in the interior of the first octant. Then, extinction in the infected population takes place, and the solution approaches asymptotically to a periodic dynamic in the SY-plane. On the other hand, decreasing the death rate of the predator population destabilizes the CEP, and a Hopf bifurcation takes place first in the interior of the first octant, and then, an extinction in the infected population occurs, and then, the solution goes to the periodic dynamics in the interior of the SY-plane; however, increasing this rate leads to extinction in the predator population, and the solution approaches PFEP.Finally, from a different angle, Maiti et al. [16] examined a delayed ecoepidemiological model for infected prey and predator, when the predator only ingested the diseased prey, with a nonlinear incidence rate and Crowley-Martin functional response. They discovered that the model becomes unstable when the time delay rises, the contact rate rises, or the inhibitory effect declines. The model is found to become unstable as the time delay grows or/and the contact rate decreases or/and the inhibitory effect decreases in this study, where it is believed that the predator does not distinguish between healthy and nonhealthy prey and consumes both.



---

*Source: 1366763-2023-02-10.xml*