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.
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: Section 2 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 where the positive parameters , , , , and 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: where 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:
Here, the prey population is divided into two compartments, and , which stand for the susceptible and infected populations at time , respectively. As of time , represents the predator population. The sum 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 and carrying capacity , 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 [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 (), 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 and , where and are the attack rates of the susceptible and infected prey, respectively, and is the preferred rate for the infected due to their weakness. The predator, on the other hand, grows at rates of and and, in the absence of food, dies at a rate of (predator mortality rate)(iv)Since most diseases do not actually spread immediately away when and 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 factor
Keeping 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:
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 is arising from the contacts of the actual population of susceptible and infected prey at time , where is the discrete-time delay.
All the system parameters are assumed to be positive, and system (4) initial conditions are where with . Here, represents the Banach space of all continuous functions with the norm .
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
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 . Hence, the solution of system (4) with initial conditions (5) exists and is unique on where [27].
Now, from system (4), it is obtained that
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 , consider as any solution of system (4) with their given initial conditions. From the first equation, it is deduced directly that
Let ; then, system (4) gives that
Since the predator growth rate constants cannot be exceeding the predator attack rate constants, hence, from a biological point of view, is always true. Then
where and . By solving the above linear differential inequality, it is obtained that , 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 as and with , respectively.
The predator-free equilibrium point (PFEP) is represented by , where while is determined as a positive root of the following equation: where , , and .
Applying Descartes’ rule of signs [28], there exists a unique positive real root of equation (11) that is given by if and only if the following condition is met:
The infected prey-free equilibrium point (IPYFEP), which is denoted by , is determined as
Clearly, the IPYFEP exists if the following conditions are met:
The coexistence equilibrium point (CEP) that is represented by is obtained by solving the following set of algebraic equations:
Direct computation gives that
Now, using the first two equations of system (16) with equation (17) gives the following quadratic equation: where
Applying Descartes’ rule of signs on equation (18) shows that it has a unique positive real root if and only if and have opposite signs. Then, using the obtained positive root in equation (17) gives that . Straightforward computation shows that is positive if and only if one of the following sets of conditions is met: or
Finally, substituting the resulting positive value of in the first equation of system (16) gives that
Consequently, it is easy to verify that CEP exists uniquely in the interior of ; if and only if in addition to either (20a) or (20b), the following condition is met:
4. Local Stability of Equilibrium Points
Let be an arbitrary equilibrium point of system (4). Let , , and , respectively. As a result, after removing the bar, system (4) becomes equivalent to the following linear system below: where
Then, the characteristic equation of system (4) at can be determined by
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 .
Proof. Substituting in general characteristic equation (25) gives that Therefore, the eigenvalues of system (4) at TEP are The TEP becomes a saddle point for all 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 that and ; then, the following is obtained: (i)If , then all the roots of have positive real parts for (ii)If , then all the roots of have negative real parts for any
Theorem 4. The DPFEP of system (4) is locally asymptotically stable for all if the following conditions are met: However, it is an unstable point for all , 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:
Clearly, equation (30) has two explicit real eigenvalues given by and . It is easy to verify that the second eigenvalue is negative due to condition (29).
On the other hand, the rest of roots of equation (30) follows from the transcendental equation:
Without a doubt, for , equation (31) has a singular root given by , 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 or a saddle point (28).
Now for , due to a proposition (3), all the roots of equation (31) have negative real part roots or positive real part roots for when condition (28) is satisfied or reflected, respectively. Consequently, the DPFEP is locally asymptotically stable for all if conditions (28) and (29) are met, while it is an unstable point for when condition (28) is reflected. This completes the proof.
Theorem 5. The PFEP of system (4) is locally asymptotically stable for all if the following conditions are met:
Proof. According to linear system (23) and their characteristic equation (25) at the PFEP, the characteristic equation becomes as follows: where
Clearly, equation (36) has one explicit real eigenvalue given by , which is negative under condition (32). However, the other eigenvalues are obtained from the equation:
For , 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 , 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 where , .
Therefore, direct computation by substituting in equation (36) and then separate the real and imaginary parts gives that
By squaring equations (37) and (38) and adding to each other, it is obtained that
Or equivalently, when ,
According to the discard rule of sign, equation (40) has no positive root. Hence, there is no , for which satisfies equation (36). Therefore, all the eigenvalues of equation (36) do not intersect the imaginary axis as and still have negative real parts. Thus, the PFEP is asymptotically stable for all if the given conditions are met.
Theorem 6. The IPYFEP of system (4) is locally asymptotically stable for all if the following conditions are met: However, it is an unstable point for all , 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: where
According the Routh-Hurwitz criterion, it is easy to verify that equation (43) has two eigenvalues with negative real parts that are given by if and only if conditions (41) and (42) are met.
However, the other eigenvalues of equation (44) are determined from the equation:
Clearly, and . Not that, for , equation (47) has a unique negative root provided that condition (43) holds. Therefore, the IPYFEP is locally asymptotically stable.
Now, for , 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 if condition (43) is reflected, where
Therefore, the IPYFEP is locally asymptotically stable for all if conditions (41)–(43) are met, and it is unstable for all if condition (3) is reflected.
Now, to study the local stability of the CEP of delayed system (4), it is assumed that and in linear system (23). Then, the characteristic equation becomes as follows: where
Theorem 7. The CEP of system (4) is locally asymptotically stable for any , if and only if the following conditions are satisfied:
Proof. For , equation (49) reduces to If the requirements (51)–(53) are satisfied, it is determined by the well-known Routh–Hurwitz criterion that system (4) with 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 , such that the CEP of system (4) is locally asymptotically stable for all , and it is unstable for provided that the following condition holds: Finally, system (4) undergoes a Hopf bifurcation at provided that 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 ; hence, all the roots of the characteristic equation (49) at have negative real parts. Moreover, by Butler’s lemma [29], remains stable for for a specific value . 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, , with . Then, in (49), by separating the real and imaginary components, it is obtained that
By squaring and adding (57) and (58) with eliminating, the following algebraic equation ofis derived:
where and
Substituting in (59) gives the following third order equation:
Clearly, condition (55) guarantees that equation (60) will include at least one positive root denoted by , which leads to instability for . In the following, the stability’s change of in relation to is determined.
From equations (56) and (58), the following is reached:
Then, corresponding to can be obtained as
Define
Now, to complete the proof of the theorem, the occurrence of Hopf bifurcation at 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 , which are pure imaginary at that is . Therefore, the proof follows if we can prove the transversality criterion [31].
Taking the derivative of equation (49) with respect to gives that substituting into equation (65) yields where due to condition (56).
So, we can get
Thus, , so the Hopf bifurcation takes place at .
If we assume that then the characteristic equation has roots with positive real parts when It is in contradiction with the local stability of the positive equilibrium point. Hence, .
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 . 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)If , then the Hopf bifurcation is supercritical (subcritical) and the bifurcating periodic solutions exist for (2)If , then the Hopf bifurcation is stable (unstable)(3)If , then the period of the bifurcating periodic solutions is increased (decreased), where , , and are given as follows:
Proof. Let , , , and , where is define by equation (63) and , so is the Hopf bifurcation value of system (4). Rescaling the time delay , then system (4) can be transformed into a functional differential equation in as follows:
where and , with
where
with and are given in equation (23), while the nonlinear term is given by
where
with
Moreover, direct computation gives the following higher derivatives:
According to the Riesz representation theorem, there is a matrix function , such that
In fact, it can choose
where is the Dirac delta function.
For , it is defined that
Thus, system (70) is equivalent to the abstract operator differential equation
where .
Now, for , the adjoint operater of is defined as
and a bilinear form
where . Clearly, and are adjoint operators. Thus, for , by a simple computation, we can calculate be the eigenvector of belonging to the eigenvalue , and is the eigenvector of that associated with the eigenvalue , where
From bilinear inner product (81), we get
Let , where is the conjugate complex number of , then and
In the following, using similar arguments as in [5], it can be determined the properties of the Hopf bifurcation:
where
with
Notice that and are constant vectors that can be determined from the following equations:
where
Accordingly, it is obtained that
Thus, , where and is the value of the determinant , where is formed by replacing the column vector of by for . Similarly, , and , where and is the value of the determinant , where is formed by replacing the column vector of by for .
Consequently, and 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:
For dataset (91), the trajectories of system (4) starting from different initial values approach asymptotically to the CEP that is given by as shown in Figure 1.
(a)
(b)
According to Figure 1, dataset (91) satisfies Theorems 7 and 8 and that at which the system undergoes a Hopf bifurcation. Therefore, in Figure 2, the bifurcation point is shown.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
According to Figure 2, system (4) undergoes a Hopf bifurcation at ; after that, the periodic dynamics became larger in the range . It is observed that for the range , system (4) approaches to periodic dynamic in the -plane; however, it returns to periodic for the range , while it has a chaotic dynamic in the range 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. Figure 3 shows the influence of varying the parameter , while Figure 4 shows the influence of varying the parameter .
(a)
(b)
(c)
(d)
(e)
(f)
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
Obviously, the parameter 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 -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 in the -plane.
In Figures 5 and 6, the influence of the parameters and is investigated, respectively.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
From Figures 5 and 6, it is observed that both the parameters and 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 has opposite influence on the dynamic of system (4) as that obtained for , while the parameter has similar influence on the dynamic of system (4) as that obtained for .
Now, Figures 7–9 explain the influence of the parameters , and , respectively.
(a)
(b)
(c)
(d)
(e)
(f)
(a)
(b)
(c)
(d)
(e)
(f)
(a)
(b)
(c)
(d)
(e)
(f)
From Figures 7–9, it is observed that the parameter has a destablizing effect on the dynamic of system (4). However, the parameter destablizes system (4) first, and then after a critical value, the system loses its persistence and approaches to a periodic dynamic in the -plane, while the parameter stablizes system (4) first and then loses its persistence and approaches to in the -plane.
It is observed that the parameters and have similar influence on the dynamic of system (4) as that obtained for . Finally, the influence of the parameter on the dynamic of system (4) is shown in Figure 10.
(a)
(b)
(c)
(d)
According to Figure 10, dataset (91) satisfies the Theorem 4 in the range , while system (4) with set (91) satisfies Theorem 5 in the range . Otherwise, it approaches to the . So, the parameter 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 , while it is an unstable otherwise and a Hopf bifurcation occurs at . 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 ; it has a 3D periodic dynamics, 2D periodic dynamics, and even chaotic dynamics for all values . 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 () leads to extinction in the predator species, and the solution is settled at the PFEP. However, increasing this parameter over a certain value (), 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 -plane.
Decreasing the contact rate below the destabilizes the CEP, and system (4) approaches asymptotically to 3D periodic dynamics (Hopf bifurcation occurs) while when the contact rate in the range , the infected population faces extinction and the system approaches asymptotically to 2D periodic dynamics in the -plane. On the other hand, increasing the contact rate above the vital value , 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 , the predator goes extinction and the system approaches asymptotically PFEP. However, for , the PFEP loses its stability and a Hopf bifurcation occurs in the interior of -plane. Moreover, the system approaches asymptotically PFEP for . 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 , and then, for , system (4) approaches asymptotically to a periodic dynamic in the interior of -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 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 -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 -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.
Data Availability
No data were used in this study.
Conflicts of Interest
The authors declare that there are no competing interests regarding the publication of this paper.