Adv Syst Sci Appl 2025; 1:55–85 Published online at https://ijassa.ipu.ru. Stability Analysis of a Delayed Epidemic Model with General Incidence and Treatment Functions Amine Bernoussi1*, Khalid Hattaf2,3, Brahim El Boukari4 1Laboratory: Équations aux dérivées partielles, Algèbre et Géométrie spectrales, Faculty of Science, Ibn Tofail University, BP 133, 14000 Kenitra, Morocco 2Equipe de Recherche en Modélisation et Enseignement des Mathématiques (ERMEM), Centre Régional des Métiers de l’Education et de la Formation (CRMEF), 20340 Derb Ghalef, Casablanca, Morocco 3Laboratory of Analysis, Modeling and Simulation (LAMS), Faculty of Sciences Ben M’Sick, Hassan II University of Casablanca, P.O Box 7955 Sidi Othman, Casablanca, Morocco 4Laboratory of Applied Mathematics and Scientific Calculus (LMACS), Higher school of technology, Sultan Moulay Slimane University, 23000 Béni Mellal, Morocco Abstract: In this paper, we study an SIR epidemic model with general nonlinear incidence function, general function of treatment and two discrete time delays, the first described the time delay due to the latent period of the disease and the second is the time delay due to the period that the infected individuals use to move into the recovered class. Lyapunov’s method is used to show the global stability of the disease-free equilibrium if the basic reproduction number R0 ≤ 1, while if R0 > 1 and under some conditions of delays, the existence of Hopf bifurcation appears for the endemic equilibrium. Keywords: Treatment, nonlinear incidence, time delay, global stability, Hopf bifurcation. 1. INTRODUCTION Mathematical modeling has been used in several areas, it is considered as a decision support tool on any situation. In the epidemiological area, it has contributed to the epidemiological surveillance of the disease because it makes to predict the health consequences of actions as varied as vaccination, quarantine or the distribution of screening tests. Among the models used in epidemiology, there are the SIR models which are compartment models, where S stands for susceptible subpopulation, I is infected subpopulation, and R is recovered subpopulation. The incidence rate of a disease measures how fast the disease is spreading and it plays an important role in the research of mathematical epidemiology. In many previous epidemic models, the bilinear incidence rate βSI was frequently used [13, 21, 22, 26, 27, 29, 34, 40, 43, 44]. However, there are some advantages for adopting more general forms of incidence rates. For instance, Capasso et al. [14–16] observed that the incidence rate may increase more slowly as I increases. So, they proposed a saturated incidence rate βSI 1+αI which was used in [2, 9, 33, 41], where α is the saturation factor that measures the inhibitory effect, βI measures the infection force of the disease and 1 1+αI measures the inhibition effect from the behavioral change of the susceptible individuals when their number increases or from the crowding effect of the infected individuals. This incidence rate seems more reasonable than the bilinear incidence rate βSI , because it includes the behavioral change and crowding ∗Corresponding author: amine.bernoussi@yahoo.fr 56 A. BERNOUSSI, K. HATTAF, B. EL BOUKARI effect of the infected individuals and prevents the unboundedness of the contact rate by choosing suitable parameters. In the last years, many forms of incidence function have been considered by the researchers in mathematical epidemiology. For example, the first one is the saturated incidence βSI d+S+I [3], where β and d are the positive constants. The second one is the Beddington-Deangelis incidence βSI 1+α1S+α2I [6], where α1 and α2 are the positive constants. The effect of saturation factor (refers to α1 and α2) stems from epidemic control and the protection mesures. The third one is the standard incidence βSI N [18, 26]. In addition, a recent Hattaf-Yousfi incidence βSI α0+α1S+α2I+α3SI that includes the three above functions was introduced in [24] and used in [36]. Models with incidence functions of the form g(I)h(S) have been studied in [30, 37]. The most general incidence function f(S, I)I [25] which generalizes the previous incidence functions has been studied by the many authors [4,10–12]. It is well known that treatment is an important and effective method to prevent and control the spread of various infectious diseases. Therefore, it very important to adopt a suitable treatment function. For instance, Wang and Ruan [40] introduced a constant treatment in SIR model as follows: { r, I > 0, 0, I = 0, which simulated a limited capacity for treatment. Further, Wang [39] considered the following piecewise linear treatment function: { rI, 0 ≤ I ≤ I0, rI0, I > I0, where I0 is the infective level at which the health care system reaches capacity. Based on this, Zhang et al. [42] proposed the following saturated treatment function rI 1+kI , where r is the maximal medical resources supplied per unit time and k is the saturation factor that measures the effect of the infected being delayed treatment. A very general from of treatment function T (I) was considered by Elazzouzi et al. [19]. Epidemiological models can contain a delay that is either discrete [8, 17] or continuous [7,10,35] because the delays appear in differential equations to describe the time lag between the action on the system and the system’s response to this action, or because a some threshold must be reached before the system is activated. In [33], Liu proposed an SEIR epidemic model with saturated incidence and saturated treatment function with two discrete delays: the time delay due to the latent period of the disease and the time delay due to the period that the infected individuals use to move into the recovered class. Motivated by the above works, we propose a mathematical SIR model that generalizes the above models and incorporates the general nonlinear incidence function f(S, I)I , general function treatment T (I) and two discrete time delays, the first τ1 described the time delay due to the latent period of the disease and the second τ2 is the time delay due to the period that the infected individuals use to move into the recovered class. This model is given by the following nonlinear system: dS(t) dt = A− µS(t) + γ1I(t)− f(S(t), I(t))I(t), dI(t) dt = f(S(t− τ1), I(t− τ1))I(t− τ1)− η1I(t)− γ2I(t− τ2)− T (I(t− τ2)), dR(t) dt = T (I(t− τ2)) + γ2I(t− τ2)− µR(t), (1.1) where η1 = µ+ γ1 + α. Here A is the constant recruitment rate into the population, S(t) represents the number of individuals who are susceptible to disease, that is, who are not yet infected at time t, I(t) represents the number of infected individuals who are infectious and are able to spread the disease by contact with susceptible individuals at time t, R(t) is the Copyright © 2025 ASSA. Adv Syst Sci Appl (2025) STABILITY ANALYSIS OF A DELAYED EPIDEMIC MODEL... 57 number of individuals who have been infected and temporarily recovered at time t, µ is the natural death rate of the population, τ1 is the units of time after infection expressing latent period, τ2, is the time delay due to the period that the infected individuals use to move into the recovered class, γ1 is the transfer rate from the infected class to the susceptible class, γ2 is the transfer rate from the infected class to the recovered class, α is the disease-induced death rate, T (I) is the general treatment function and f(S, I)I is the incidence function, i.e., the number of susceptible individuals infected through their contacts with the infectious individuals. Model (1.1) generalizes several special cases existing in the literature. For example, Abta et al. [1] studied the model (1.1) with f(S, I)I = βSI 1+α1S+α2I , γ1 = 0, τ2 = 0 and T (I) = 0. In [25], Hattaf et al. studied the model (1.1) with γ1 = 0, τ2 = 0 and T (I) = 0. In [40], Wang and Ruan studied the model (1.1) with f(S, I)I = βSI, γ1 = 0, τ1 = τ2 = 0 and T (I) = r. In [39], Wang studied the model (1.1) with f(S, I)I = βSI, γ1 = 0, τ1 = τ2 = 0 and T (I) = { rI, 0 ≤ I ≤ I0, rI0, I > I0. In [5], Balamuralitharan et al. studied the model (1.1) with f(S, I)I = βSI, γ1 = 0, τ1 = τ2 = 0 and T (I) = rI . In this paper, the initial condition for the system (1.1) is: S(θ) = φ1(θ), I(θ) = φ2(θ), R(θ) = φ3(θ), θ ∈ [−τ, 0], (1.2) where τ = max{τ1, τ2} and φ = (φ1, φ2, φ3) T ∈ C such that φi(θ) ≥ 0 for −τ ≤ θ ≤ 0 and i = 1, 2, 3. The space C denotes the Banach space C([−τ, 0],R3 +0) of continuous functions mapping the interval [−τ, 0] into R3 +0 with the supremum norm, where R+0 = {x ∈ R | x ≥ 0}. Also, we assume that φi(0) > 0 for i = 1, 2, 3. On the other hand, the first two equation in system (1.1) do not depend on the third equation, and therefore this equation can be omitted without loss of generality. System (1.1) can be rewritten as{ dS(t) dt = A− µS(t) + γ1I(t)− f(S(t), I(t))I(t), dI(t) dt = f(S(t− τ1), I(t− τ1))I(t− τ1)− η1I(t)− γ2I(t− τ2)− T (I(t− τ2)), (1.3) the incidence function f(S, I)I : R+ × R+ → R+ is a continuously differentiable and locally Lipschitz function on R+ × R+ and satisfying the following hypotheses (see [11, 25]): (H0) f(0, I) = 0 for I ≥ 0; (H1) f(S, I) is a strictly monotone increasing function of S > 0, for any fixed I ≥ 0; (H2) f(S, I) is a monotone decreasing function of I ≥ 0, for any fixed S ≥ 0; (H3) ϕ(S, I) = f(S, I))I is a monotone increasing function of I ≥ 0, for any fixed S ≥ 0. Moreover, the treatment function T : R+ → R+ is a locally Lipschitz continuous differentiable function on R+ satisfying the following hypotheses (see [19]) (T0) T (0) = 0; (T1) T (I) I is a monotone increasing function of I > 0. The rest of the paper is organized as follows. In Section 2, we carry out mathematical analysis about the basic reproduction number and the existence of equilibria of the model (1.1). In Section 3, we study the global stability of the disease-free equilibrium. Section 4 is devoted to the stability and Hopf bifurcation of the endemic equilibrium. An application of our results and some numerical simulations are presented in Section 5. At the end, we present some concluding remarks. Copyright © 2025 ASSA. Adv Syst Sci Appl (2025) 58 A. BERNOUSSI, K. HATTAF, B. EL BOUKARI 2. MATHEMATICAL ANALYSIS In this section, we prove the basic results which guarantee the positivity of solutions as well as the existence and uniqueness of the endemic equilibrium for system (1.3) under initial condition (1.2). Lemma 2.1: The closed set: Ω := { (S, I) ∈ (R+)2|S + I ≤ A µ } is positively invariant with respect to system (1.3). Proof Let N(t) = S(t) + I(t), and (S, I) ∈ (R+)2. Then it follows from system (1.3) that dN(t) dt = A− µS(t)− µI(t)− αI(t)− γ2I(t− τ2)− f(S(t), I(t))I(t) + f(S(t− τ1), I(t− τ1))I(t− τ1)− T (I(t− τ2)) = A− µN(t)− ∫ t t−τ1 f(S(σ), I(σ))I(σ)dσ − T (I(t− τ2))− αI(t)− γ2I(t− τ2). By the hypotheses, we have T (I(t− τ2)) ≥ 0 and f(S(t), I(t))I(t) ≥ 0 for all t ≥ 0 and τ2 ≥ 0. Then dN(t) dt ≤ A− µN(t), which implies that N(t) ≤ A µ when N(0) ≤ A µ . This completes the proof. Hence, we discuss system (1.3) in the closed set Ω. Next, we discuss the existence of equilibria for system (1.3). The basic reproduction number of system (1.3) can be defined by R0 = f(A µ , 0) η1 + η2 , where η2 = γ2 + T ′ (0). Note that the system (1.3) always has a disease-free equilibrium P0 = (A µ , 0). On the other hand, to prove the existence and uniqueness of an endemic equilibrium, we need the following Lemma as in [19]. Lemma 2.2: Assume that the assumptions (T0) and (T1), are satisfied. Then b2 − a2u− T (u) = 0, for a2 > 0 and b2 > 0, has a unique positive solution u1, and b2 − a2u− T (u) > 0, for all u ∈ [0, u1), and b2 − a2u− T (u) ≤ 0 for all u ∈ [u1, b2 a2 ]. Proof Let L be the function defined on R+ by L(u) := b2 − a2u− T (u), we have: L(0) = b2 > 0 and L( b2 a2 ) = −T ( b2 a2 ) < 0 and L ′ = −a2 − T ′ < 0. Since L is continuous and is strictly monotone decreasing function, then the equation L(u) = 0, has a unique positive solution in the interval (0, b2 a2 ). This completes the proof. Copyright © 2025 ASSA. Adv Syst Sci Appl (2025) STABILITY ANALYSIS OF A DELAYED EPIDEMIC MODEL... 59 Theorem 2.1: Under the hypotheses (H1), (H2), (T0) and (T1) if R0 > 1, then system (1.3) admits a unique endemic equilibrium P ∗ = (S∗, I∗), with S∗ = A µ − (µ+ γ2 + α)I∗ + T (I∗) µ , and I∗ is the unique solution of the following equation: f( A µ − (µ+ γ2 + α)I + T (I) µ , I) = η1 + γ2 + T (I) I . Proof For simplicity, we put K(I) = T (I) I . At an equilibrium point (S, I) of system (1.3), the following equations hold.{ A− µS + γ1I − f(S, I)I = 0, f(S, I)I − (η1 + γ2)I − T (I) = 0. (2.4) Substituting the second equation into the first equation of (2.4), we obtain the following system: S = A µ − (µ+ γ2 + α)I + T (I) µ , f(S, I)I = (η1 + γ2)I + T (I). (2.5) If I = 0, we obtain the disease-free equilibrium point P0 = (A µ , 0). If I ̸= 0, then using (2.5) we get the following equation f( A µ − (µ+ γ2 + α)I + T (I) µ , I) = η1 + γ2 +K(I). We have S = A µ − [(µ+γ2+α)I+T (I)] µ , which implies that S ≤ A µ . By Lemma 2.2, we have A µ − (µ+γ2+α)I+T (I) µ > 0 if and only if I ∈ [0, I1), where I1 is a unique positive solution of the equation A µ − [(µ+γ2+α)I+T (I)] µ = 0. Hence, there is no positive equilibrium point if S > A µ or I ≥ I1. Now, we consider the function g defined on the interval [0, I1] as follows g(I) := f( A µ − (µ+ γ2 + α)I + T (I) µ , I)− (η1 + γ2 +K(I)). We have lim I→0 g(I) = f( A µ , 0)− (η1 + η2) = (η1 + η2) (f(A µ , 0) η1 + η2 − 1 ) = (η1 + η2)(R0 − 1) > 0 for R0 > 1, and g(I1) = −(η1 + γ2 +K(I1)) < 0. Copyright © 2025 ASSA. Adv Syst Sci Appl (2025) 60 A. BERNOUSSI, K. HATTAF, B. EL BOUKARI Furthermore, g ′ (I) = −µ+ γ2 + α + T ′ µ ∂f ∂S + ∂f ∂I −K ′ . According to the hypotheses (H1), (H2) and (T1), we have g ′ (I) < 0. Hence, there exists a unique endemic equilibrium P ∗ = (S∗, I∗) with I∗ ∈ (0, I1) and S∗ > 0 satisfies the equations S∗ = A µ − (µ+γ2+α)I∗+T (I∗) µ . This completes the proof. 3. GLOBAL STABILITY OF THE DISEASE-FREE EQUILIBRIUM In this section, we establish the global stability of the disease-free equilibrium P0 of system (1.3). We can defined Rc = f(A µ , 0) η1 . Note that R0 < Rc. Theorem 3.1: Suppose the hypotheses (H1), (H2), (T0) and (T1) hold. • If R0 > 1, then the disease-free equilibrium P0 of system (1.3) is unstable. • If R0 ≤ 1 then we distinguish three cases: (i) The disease-free equilibrium P0 of system (1.3) is globally asymptotically stable whenever τ1 ≥ τ2. (ii) When Rc ≤ 1, the disease-free equilibrium P0 of system (1.3) is globally asymptotically stable for all τ1 ≥ 0 and τ2 ≥ 0, then in particularly for τ2 ≥ τ1. (iii) When Rc > 1, the equation (3.9) has a purely imaginary root. If more h ′ 0(ω ∗ 22,0) ̸= 0, then there exists a positive τ ∗22,0 where system (1.3) undergoes Hopf bifurcation at P0 when τ2 = τ ∗22,0. However, the steady state P0 is locally asymptotically stable when τ2 ∈ [0, τ ∗22,0) and unstable when τ2 > τ ∗22,0, where τ ∗22,0 is given by τ ∗22,0 = 1 ω∗ 22,0 arccos ( f(A µ , 0) cos (ω∗ 22,0τ1)− η1 η2 ), and ω∗ 22,0 is the positive root of equation (3.9), and h0 are defined in equation (3.13). Proof • If R0 > 1, then the characteristic equation at P0 is given by (λ+ µ)(λ− f( A µ , 0)e−λτ1 + η1 + η2e −λτ2) = 0. (3.6) Obviously, λ = −µ is eigenvalue for (3.6), and hence, the stability of P0 is determined by the distribution of the roots of equation λ− f( A µ , 0)e−λτ1 + η1 + η2e −λτ2 = 0. (3.7) Copyright © 2025 ASSA. Adv Syst Sci Appl (2025) STABILITY ANALYSIS OF A DELAYED EPIDEMIC MODEL... 61 We put Ψ(λ) = λ− f( A µ , 0)e−λτ1 + η1 + η2e −λτ2 . We have lim λ→+∞ Ψ(λ) = +∞, and Ψ(0) = −f( A µ , 0) + η1 + η2 = (η1 + η2)(1− f(A µ , 0) (η1 + η2) ) = (η1 + η2)(1−R0). Since the function Ψ(λ) is continuous on the interval [0,+∞), we conclude that the equation Ψ(λ) = 0 has a positive real root and the disease-free equilibrium is unstable when R0 > 1. • If R0 ≤ 1, then there are three cases: - First case τ1 > τ2 Consider the following Lyapunov functional V0(t) = ∫ S(t−τ1) A µ (1− f(A µ , 0) f(σ, 0) )dσ + I + η1 ∫ t t−τ2 I(ξ)dξ + ∫ t−τ2 t−τ1 f(A µ , 0) f(S(σ), 0) f(S(σ), I(σ))I(σ)dσ. We will show that dV0(t) dt ≤ 0 for all t ≥ 0. We have dV0(t) dt = (1− f(A µ , 0) f(S(t− τ1), 0) )(A− µS(t− τ1) + γ1I(t− τ1) − f(S(t− τ1), I(t− τ1))I(t− τ1)) + f(S(t− τ1), I(t− τ1))I(t− τ1)− η1I(t)− γ2I(t− τ2)− T (I(t− τ2)) + η1I(t)− η1I(t− τ2) + f(A µ , 0) f(S(t− τ2), 0) f(S(t− τ2), I(t− τ2))I(t− τ2)) − f(A µ , 0) f(S(t− τ1), 0) f(S(t− τ1), I(t− τ1))I(t− τ1)) = µ(1− f(A µ , 0) f(S(t− τ1), 0) )( A µ − S(t− τ1)) + γ1I(t− τ1)(1− f(A µ , 0) f(S(t− τ1), 0) ) + f( A (µ) , 0) f(S(t− τ2), 0) f(S(t− τ2), I(t− τ2))I(t− τ2)) − (η1 + γ2)I(t− τ2)− T (I(t− τ2)). Copyright © 2025 ASSA. Adv Syst Sci Appl (2025) 62 A. BERNOUSSI, K. HATTAF, B. EL BOUKARI Then we have dV0(t) dt = µ(1− f(A µ , 0) f(S(t− τ1), 0) )( A µ − S(t− τ1)) + γ1I(t− τ1)(1− f(A µ , 0) f(S(t− τ1), 0) ) + ( (η1 + γ2)I(t− τ2) + T (I(t− τ2)) ) × ( f(A µ , 0)f(S(t− τ2), I(t− τ2))I(t− τ2) f(S(t− τ2), 0)[(η1 + γ2)I(t− τ2) + T (I(t− τ2))] − 1 ) = µ(1− f(A µ , 0) f(S(t− τ1), 0) )( A µ − S(t− τ1)) + γ1I(t− τ1)(1− f(A µ , 0) f(S(t− τ1), 0) ) + ( (η1 + γ2)I(t− τ2) + T (I(t− τ2)) ) × ( R0 η1 + η2 η1 + γ2 +K(I(t− τ2)) f(S(t− τ2), I(t− τ2)) f(S(t− τ2), 0) − 1 ) . Furthermore, it follows from the hypotheses (H2), (T0) and (T1) that f(S(t− τ2), I(t− τ2)) f(S(t− τ2), 0) ≤ 1, and η1 + η2 η1 + γ2 +K(I(t− τ2)) = η1 + γ2 + limI→0 T (I(t−τ2)) I(t−τ2) η1 + γ2 + T (I(t−τ2)) I(t−τ2) ≤ 1. Hence, dV0(t) dt ≤ µ(1− f(A µ , 0) f(S(t− τ1), 0) )( A µ − S(t− τ1)) + γ1I(t− τ1)(1− f(A µ , 0) f(S(t− τ1), 0) ) + [(η1 + γ2)I(t− τ2) + T (I(t− τ2))](R0 − 1). By the hypothesis (H1), we obtain that (1− f(A µ , 0) f(S(t− τ1), 0) )( A µ − S(t− τ1)) ≤ 0, and I(t− τ1)(1− f(A µ , 0) f(S(t− τ1), 0) ) ≤ 0, where equality holds if and only if S = A µ . Since R0 ≤ 1, ensures that dV0(t) dt ≤ 0, for all t ≥ 0, Thus, the disease-free equilibrium P0 is stable and dV0(t) dt = 0 holds if and only if S = A µ and [(η1 + γ2)I(t− τ2) + T (I(t− τ2))] ( R0 η1+η2 η1+γ2+K(I(t−τ2)) f(S(t−τ2),I(t−τ2)) f(S(t−τ2),0) − 1 ) = 0. We dis- cuss two cases: Copyright © 2025 ASSA. Adv Syst Sci Appl (2025) STABILITY ANALYSIS OF A DELAYED EPIDEMIC MODEL... 63 • If R0 < 1, then it follows from S = A µ and Lemma 2.1 that I = 0. • If R0 = 1, then it follows from S = A µ and Lemma 2.1 that I = 0. By the above discussions, we deduce that {P0} is the largest invariant set in{ (S, I)|dV0(t) dt = 0 } . From the Lyapunov-LaSalle theorem [31, 32], we conclude that P0 is globally asymptotically stable. - Second case τ1 = τ2 = τ Consider the following Lyapunov functional W0(t) = ∫ S(t−τ) A µ (1− f(A µ , 0) f(σ, 0) )dσ + I + η1 ∫ t t−τ I(ξ)dξ. We will show that dW0(t) dt ≤ 0 for all t ≥ 0. We have dW0(t) dt = (1− f(A µ , 0) f(S(t− τ), 0) )(A− µS(t− τ) + γ1I(t− τ) − f(S(t− τ), I(t− τ))I(t− τ)) + f(S(t− τ), I(t− τ))I(t− τ)− η1I(t)− T (I(t− τ)) + η1I(t) − (η1 + γ2)I(t− τ) = µ(1− f(A µ , 0) f(S(t− τ), 0) )( A µ − S(t− τ)) + γ1I(t− τ)(1− f(A µ , 0) f(S(t− τ), 0) ) + f(A µ , 0) f(S(t− τ), 0) f(S(t− τ), I(t− τ))I(t− τ)) − (η1 + γ2)I(t− τ)− T (I(t− τ)). Then we have dW0(t) dt = µ(1− f(A µ , 0) f(S(t− τ), 0) )( A µ − S(t− τ)) + γ1I(t− τ)(1− f(A µ , 0) f(S(t− τ), 0) ) + ( (η1 + γ2)I(t− τ) + T (I(t− τ)) ) × ( f(A µ , 0)f(S(t− τ), I(t− τ)) f(S(t− τ), 0)(η1 + γ2 +K(I(t− τ))) − 1 ) = µ(1− f(A µ , 0) f(S(t− τ), 0) )( A µ − S(t− τ)) + γ1I(t− τ)(1− f(A µ , 0) f(S(t− τ), 0) ) + ( (η1 + γ2)I(t− τ) + T (I(t− τ)) ) × ( R0 η1 + η2 η1 + γ2 +K(I(t− τ)) f(S(t− τ), I(t− τ)) f(S(t− τ), 0) − 1 ) . Furthermore, it follows from the hypotheses (H2), (T0) and (T1), we have f(S(t− τ), I(t− τ)) f(S(t− τ), 0) ≤ 1, Copyright © 2025 ASSA. Adv Syst Sci Appl (2025) 64 A. BERNOUSSI, K. HATTAF, B. EL BOUKARI and η1 + η2 (η1 + γ2) +K(I(t− τ)) = η1 + γ2 + limI→0 T (I(t−τ)) I(t−τ) η1 + γ2 + T (I(t−τ)) I(t−τ) ≤ 1. Then dW0(t) dt ≤ µ(1− f(A µ , 0) f(S(t− τ), 0) )( A µ − S(t− τ)) + γ1I(t− τ)(1− f(A µ , 0) f(S(t− τ), 0) ) + [(η1 + γ2)I(t− τ) + T (I(t− τ))](R0 − 1). By the hypothesis (H1), we get that (1− f(A µ , 0) f(S(t− τ), 0) )( A µ − S(t− τ)) ≤ 0, and I(t− τ)(1− f(A µ , 0) f(S(t− τ), 0) ) ≤ 0, where equality holds if and only if S = A µ . Since R0 ≤ 1, we deduce that dW0(t) dt ≤ 0 for all t ≥ 0. Thus, the disease-free equilibrium P0 is stable and dW0(t) dt = 0, holds if and only if S = A µ , and [(η1 + γ2)I(t− τ) + T (I(t− τ))] ( R0 η1+η2 η1+γ2+K(I(t−τ)) f(S(t−τ),I(t−τ)) f(S(t−τ),0) − 1 ) = 0. We dis- cuss two cases: • If R0 < 1, then it follows from S = A µ and Lemma 2.1 that I = 0. • If R0 = 1, then it follows from S = A µ and Lemma 2.1 that I = 0. Therefore, {P0} is the largest invariant set in { (S, I)|dW0(t) dt = 0 } . From the Lyapunov- LaSalle theorem [31, 32], we conclude that P0 is globally asymptotically stable. - Third case τ2 > τ1 Suppose that Rc ≤ 1. Consider the following Lyapunov functional U0(t) = ∫ S(t) A µ (1− f(A µ , 0) f(σ, 0) )dσ + I + ∫ t t−τ1 f(S(ξ), I(ξ))I(ξ)dξ. Copyright © 2025 ASSA. Adv Syst Sci Appl (2025) STABILITY ANALYSIS OF A DELAYED EPIDEMIC MODEL... 65 We will show that dU0(t) dt ≤ 0 for all t ≥ 0. We have dU0(t) dt = (1− f(A µ , 0) f(S(t), 0) )(A− µS(t) + γ1I(t)− f(S(t), I(t))I(t)) + f(S(t− τ), I(t− τ))I(t− τ)− η1I(t)− γ2I(t− τ2)− T (I(t− τ2)) + f(S(t), I(t))I(t)− f(S(t− τ), I(t− τ))I(t− τ) = µ(1− f(A µ , 0) f(S(t), 0) )( A µ − S(t)) + γ1I(t)(1− f(A µ , 0) f(S(t), 0) ) + f(A µ , 0) f(S(t), 0) f(S(t), I(t))I(t))− η1I(t)− γ2I(t− τ2)− T (I(t− τ2)). Then we have dU0(t) dt = µ(1− f(A µ , 0) f(S(t), 0) )( A µ − S(t)) + γ1I(t)(1− f(A µ , 0) f(S(t), 0) ) + η1I(t) (f(A µ , 0)f(S(t), I(t)) f(S(t), 0)η1 − 1 ) − γ2I(t− τ2)− T (I(t− τ2)). Hence, dU0(t) dt = µ(1− f(A µ , 0) f(S(t), 0) )( A µ − S(t)) + γ1I(t)(1− f(A µ , 0) f(S(t), 0) ) + η1I(t) ( Rc f(S(t), I(t)) f(S(t), 0) − 1 ) − γ2I(t− τ2)− T (I(t− τ2)). (3.8) Furthermore, it follows from the hypotheses (H2) that f(S(t), I(t)) f(S(t), 0) ≤ 1, and by the hypotheses, we have I(t− τ2) ≥ 0 and T (I(t− τ2)) ≥ 0. Then dU0(t) dt ≤ µ(1− f(A µ , 0) f(S(t), 0) )( A µ − S(t)) + γ1I(t)(1− f(A µ , 0) f(S(t), 0) ) + η1I(t) ( Rc − 1 ) . By the hypothesis (H1), we obtain that (1− f(A µ , 0) f(S(t), 0) )( A µ − S(t)) ≤ 0, Copyright © 2025 ASSA. Adv Syst Sci Appl (2025) 66 A. BERNOUSSI, K. HATTAF, B. EL BOUKARI and I(t− τ)(1− f(A µ , 0) f(S(t), 0) ) ≤ 0, where equality holds if and only if S = A µ . Since Rc ≤ 1, ensures that dU0(t) dt ≤ 0 for all t ≥ 0, Thus, the disease-free equilibrium P0 is stable and by the equation (3.8) dU0(t) dt = 0, holds if and only if S = A µ and I = 0. By the above discussions, we deduce that {P0} is the largest invariant set in{ (S, I)|dU0(t) dt = 0 } . From the Lyapunov-LaSalle theorem [31, 32], we conclude that P0 is globally asymptotically stable for all τ1 ≥ 0 and τ2 ≥ 0. Then P0 is globally asymptotically stable for all τ2 > τ1. Now, we suppose that R0 < 1 and Rc > 1 The characteristic equation at P0 is given by the equation (3.6), obviously, λ = −µ, is eigenvalue for (3.6), and hence, the stability of P0 is determined by the distribution of the roots of equation λ− f( A µ , 0)e−λτ1 + η1 + η2e −λτ2 = 0. (3.9) Equation (3.9) has a purely imaginary root iω22, with ω22 > 0 if and only if f( A µ , 0) cos (ω22τ1)− η1 = η2 cos (ω22τ2), (3.10) ω22 + f( A µ , 0) sin (ω22τ1) = η2 sin (ω22τ2). (3.11) Squaring and adding the squares together, we obtain ω2 22 + η21 + f( A µ , 0)2 − η22 − 2η1f( A µ , 0) cos (ω22τ1) + 2ω22f( A µ , 0) sin (ω22τ1) = 0. (3.12) We put h0(ω22) = ω2 22 + η21 + f( A µ , 0)2 − η22 − 2η1f( A µ , 0) cos (ω22τ1) + 2ω22f( A µ , 0) sin (ω22τ1). (3.13) On the other hand, we have −1 ≤ cos (ω22τ1) ≤ 1 and − 1 ≤ sin (ω22τ1) ≤ 1. Then we have g01(ω22) ≤ h0(ω22) ≤ g02(ω22), where g01(ω22) = ω2 22 + η21 + f( A µ , 0)2 − η22 − 2η1f( A µ , 0)− 2ω22f( A µ , 0), and g02(ω22) = ω2 22 + η21 + f( A µ , 0)2 − η22 + 2η1f( A µ , 0) + 2ω22f( A µ , 0). Copyright © 2025 ASSA. Adv Syst Sci Appl (2025) STABILITY ANALYSIS OF A DELAYED EPIDEMIC MODEL... 67 We have lim ω22→+∞ g01(ω22) = lim ω22→+∞ g02(ω22) = +∞. Then lim ω22→+∞ h0(ω22) = +∞, and h0(0) = η21 + f( A µ , 0)2 − η22 − 2η1f( A µ , 0) = (f( A µ , 0)− η1) 2 − η22 = ( f( A µ , 0)− η1 − η2 )( f( A µ , 0)− η1 + η2 ) = ( η1 + η2 )( R0 − 1 )( η1(Rc − 1) + η2 ) ≤ 0. Since h0(ω22) is continuous in [0,+∞), then the equation (3.9) has at least one positive root. We assume that equation (3.9) admits a finite family of solution ω∗ 22,o, with o = 1, ...,m and m ∈ N. By the equation (3.10), we have τ f2,o = 1 ω∗ 22,o ( arccos ( f(A µ , 0) cos (ω∗ 22,oτ1)− η1 η2 ) + 2Πf ) , f = 0, 1, ...; o = 1, 2, ...,m. Then ±iω∗ 22,o is a pair of purely imaginary root of equation (3.9), with τ2 = τ f2,o, f = 1, 2, ....; o = 1, 2, ...,m. Clearly, lim f→∞ τ f2,o = ∞, o = 1, 2, ...,m. Thus, we can define τ ∗22,0 = τ f02,o0 = min f=0,1...,o=1,2,...,m (τ f2,o), ω ∗ 22,0 = ω∗ 2,o0 . Lets show that iω∗ 22,0 is simple, consider the branche of characteristic roots λ(τ2) = x1(τ2) + iy1(τ2), of equation (3.9) bifurcating from iω∗ 22,0 at τ2 = τ ∗22,0. By derivation (3.9) with respect to the delay τ2, we obtain dλ dτ2 ( 1 + τ1f( A µ , 0)e−λτ1 − τ2η2e −λτ2 ) = λη2e −λτ2 . (3.14) If we suppose, by contradiction, that iω∗ 22,0 is not simple, the right hand side (3.14) gives iω∗ 22,0η2 = 0, and leads a contradiction with the fact that ω∗ 22,0 > 0, and η2 > 0. Next, we need to grantee the transversality condition of the Hopf bifurcation theorem (see [23]). Clearly, λ(τ2) = x1(τ2) + iy1(τ2), is a root of equation (3.9) if and only if x1 + η1 − e−x1τ1f( A µ , 0) cos (y1τ1) = −e−x1τ2η2 cos (y1τ2), (3.15) Copyright © 2025 ASSA. Adv Syst Sci Appl (2025) 68 A. BERNOUSSI, K. HATTAF, B. EL BOUKARI y1 + e−x1τ1f( A µ , 0) sin (y1τ1) = e−x1τ2η2 sin (y1τ2). (3.16) Let x1(τ ∗ 22,0) and y1(τ ∗ 22,0) satisfying x1(τ ∗ 22,0) = 0, and y1(τ ∗ 22,0) = ω∗ 22,0. By differentiating equations (3.15) and (3.16) with respect to τ2 and then set τ2 = τ ∗22,0. Doing this, we get G1,0 dx1(τ ∗ 22,0) dτ2 +G2,0 dy1(τ ∗ 22,0) dτ2 = h1,0, −G2,0 dx1(τ ∗ 22,0) dτ2 +G1,0 dy1(τ ∗ 22,0) dτ2 = h2,0, (3.17) where G1,0 = 1 + τ1f( A µ , 0) cos (ω∗ 22,0τ1)− τ2η2 cos (ω ∗ 22,0τ2), G2,0 = τ1f( A µ , 0) sin (ω∗ 22,0τ1)− τ2η2 sin (ω ∗ 22,0τ2), h1,0 = ω∗ 22,0η2 sin (ω ∗ 22,0τ2), and h2,0 = ω∗ 22,0η2 cos (ω ∗ 22,0τ2). Calculating dx1(τ∗22,0) dτ2 , we get dx1(τ ∗ 22,0) dτ2 = G1,0h1,0 −G2,0h2,0 G2 1,0 +G2 2,0 . Therefore, according to the equation (3.10) and (3.11) we have dx1(τ ∗ 22,0) dτ2 = ω∗ 22,0h ′ (ω∗ 22,0) 2(G2 1,0 +G2 2,0) ̸= 0. This proves the Theorem. 4. STABILITY AND HOPF BIFURCATION OF THE ENDEMIC EQUILIBRIUM In the next, we will study the local stability of the positive equilibrium P ∗ with respect to the time delay. We set x = S − S∗, and y = I − I∗. Then the linearized system of equations around the equilibrium point P ∗ is given as follows: dx dt = (−µ− ∂f(S∗, I∗) ∂S I∗)x(t)− ( ∂f(S∗, I∗) ∂I I∗ + f(S∗, I∗)− γ1)y(t), dy dt = ∂f(S∗, I∗) ∂S I∗x(t− τ1) + ( ∂f(S∗, I∗) ∂I I∗ + f(S∗, I∗))y(t− τ1)− η1y(t) −(γ2 + T ′ (I∗))y(t− τ2). For simplicity, we put η3 = γ2 + T ′ (I∗), and η4 = µ+ I∗ ∂f(S ∗,I∗) ∂S . Hence, the characteristic equation at P ∗ is given by ∆(λ, τ1, τ2) = ∣∣∣∣ λ+ η4 I∗ ∂f(S ∗,I∗) ∂I + f(S∗, I∗)− γ1 −I∗ ∂f(S ∗,I∗) ∂S e−λτ1 λ+ η1 − (I∗ ∂f(S ∗,I∗) ∂I + f(S∗, I∗))e−λτ1 + η3e −λτ2 ∣∣∣∣ = 0. Copyright © 2025 ASSA. Adv Syst Sci Appl (2025) STABILITY ANALYSIS OF A DELAYED EPIDEMIC MODEL... 69 By a simple computation, we get ∆(λ, τ1, τ2) = P (λ) +Q(λ)e−λτ1 + J(λ)e−λτ2 = 0, (4.18) where P (λ) = λ2 + A1λ+B, Q(λ) = Cλ+D, J(λ) = Eλ+ F, A1 = η1 + η4, B = η1η4, C = −( ∂f(S∗, I∗) ∂I I∗ + f(S∗, I∗)), D = −η4( ∂f(S∗, I∗) ∂I I∗ + f(S∗, I∗)) + ∂f(S∗, I∗) ∂S I∗( ∂f(S∗, I∗) ∂I I∗ + f(S∗, I∗))− γ1 ∂f(S∗, I∗) ∂S I∗ = −µ( ∂f(S∗, I∗) ∂I I∗ + f(S∗, I∗))− γ1 ∂f(S∗, I∗) ∂S I∗, E = η3, and F = η4η3. By the hypotheses (H2) and (T1), we easily deduce that (see [20, 28]) T ′ (I∗)− T (I∗) I∗ ≥ 0, (4.19) and f(S∗, I∗)− ∂ϕ(S∗, I∗) ∂I ≥ 0. (4.20) Several cases arise. 4.1. Case τ1 = τ2 = 0 Theorem 4.1: If R0 > 1 and τ1 = τ2 = 0, then the endemic equilibrium P ∗ is locally asymptotically stable. Proof When τ1 = τ2 = 0, the characteristic equation (4.18) reads as λ2 + (A1 + C + E)λ+ (B +D + F ) = 0 By using the second equation of system (1.3), we have f(S∗, I∗) = (η1 + γ2) + T (I∗) I∗ . (4.21) Copyright © 2025 ASSA. Adv Syst Sci Appl (2025) 70 A. BERNOUSSI, K. HATTAF, B. EL BOUKARI Then A1 + C + E = η1 + η4 − ( ∂f(S∗, I∗) ∂I I∗ + f(S∗, I∗)) + η3 = η4 − ∂f(S∗, I∗) ∂I I∗ + (T ′ (I∗)− T (I∗) I∗ ), and B +D + F = η1η4 − µ( ∂f(S∗, I∗) ∂I I∗ + f(S∗, I∗))− γ1 ∂f(S∗, I∗) ∂S I∗ + η4η3 = µ(T ′ (I∗)− T (I∗) I∗ ) + ∂f(S∗, I∗) ∂S I∗[(µ+ α) + η3]− µ ∂f(S∗, I∗) ∂I I∗. By the hypotheses (H1), (H2), (T1), and equation (4.19), we have A1 + C + E > 0 and B +D + F > 0. Hence, according to the Routh-Hurwitz criterion, the endemic equilibrium P ∗ is locally asymptotically stable. This completes the proof. 4.2. Case τ1 > 0 and τ2 = 0 Theorem 4.2: If R0 > 1, τ1 > 0 and τ2 = 0, then the endemic equilibrium P ∗ is locally asymptotically stable. Proof When τ1 > 0 and τ2 = 0 the characteristic equation (4.18) becomes ∆(λ, τ1) = λ2 + (A1 + E)λ+B + F + (λC +D)e−λτ1 = 0. (4.22) Equation (4.22) has a purely imaginary root iω1, with ω1 > 0 if and only if ω2 1 − (B + F ) = D cos (ω1τ1) + ω1C sin (ω1τ1), and ω1(A1 + E) = D sin(ω1τ1)− ω1C cos (ω1τ1). Squaring and adding the squares together, we obtain ω4 1 + cω2 1 + d = 0, (4.23) with c = (A1 + E)2 − 2(B + F )− C2 and d = (B + F )2 −D2. Letting z1 = ω2 1, then equation (4.23) becomes the following equation z21 + cz1 + d = 0. On the other hand, we have d = (B + F )2 −D2 = (B + F −D)(B + F +D). Using the case τ1 = τ2 = 0, we have (B + F +D) > 0 and by using the hypotheses (H1), (H3) and (T1), we have (B + F −D) > 0 then we have d > 0. Now, we will prove c > 0. Copyright © 2025 ASSA. Adv Syst Sci Appl (2025) STABILITY ANALYSIS OF A DELAYED EPIDEMIC MODEL... 71 Indeed c = (A1 + E)2 − 2(B + F )− C2, we have (A1 + E)2 = (η1 + η3 + η4) 2. By using the equation (4.21), find that (A1 + E)2 = ( f(S∗, I∗) + η4 + T ′ (I∗)− T (I∗) I∗ )2 = f(S∗, I∗)2 + η24 + 2f(S∗, I∗)η4 + (T ′ (I∗)− T (I∗) I∗ )2 + 2(T ′ (I∗)− T (I∗) I∗ )f(S∗, I∗) + 2(T ′ (I∗)− T (I∗) I∗ )η4, and −2(B + F ) = −2η4(η1 + η3) = −2f(S∗, I∗)η4 − 2η4(T ′ (I∗)− T (I∗) I∗ ), and −C2 = −( ∂f(S∗, I∗) ∂I I∗ + f(S∗, I∗))2 = −( ∂ϕ(S∗, I∗) ∂I )2. Then we have c = ( f(S∗, I∗)− ∂ϕ(S(t), I(t)) ∂I )( f(S∗, I∗) + ∂ϕ(S(t), I(t)) ∂I ) + η24 + (T ′ (I∗)− T (I∗) I∗ )2 + 2(T ′ (I∗)− T (I∗) I∗ )f(S∗, I∗). By the hypothesis (H3) and equations (4.19) and (4.20), we have c > 0. Then the equation (4.23) has non positive solution. Consequently, using the case τ1 = τ2 = 0, the endemic equilibrium P ∗ is locally asymptotically stable . This completes the proof. 4.3. Case τ1 = τ2 = τ Consider the assumption: η1 − γ2 + µ∂ϕ(S∗,I∗) ∂I + γ1 ∂f(S∗,I∗) ∂S I∗ η4 < T ′ (I∗). (4.24) Theorem 4.3: If R0 > 1 and (4.24) holds, then there exists a positive τ0 where system (1.3) undergoes Hopf bifurcation at P ∗ when τ = τ0. However, the steady state P ∗ is locally asymptotically stable Copyright © 2025 ASSA. Adv Syst Sci Appl (2025) 72 A. BERNOUSSI, K. HATTAF, B. EL BOUKARI when τ ∈ [0, τ0) and unstable when τ > τ0. Here τ0 is given by τ0 = 1 ω0 arccos ( (D + F )(ω2 0 −B)− ω2 0A1(C + E) ω2 0(C + E)2 + (D + F )2 ), and ω0 = √ −x+ √ (x2 − 4y) 2 , where x and y are defined in equation (4.28). Proof When τ1 = τ2 = τ , the characteristic equation (4.18) becones ∆(λ, τ) = λ2 + A1λ+B + [λ(C + E) + (D + F )]e−λτ = 0. (4.25) Equation (4.25) has a purely imaginary root iω, with ω > 0 ∆(iω, τ) = 0, if and only if ω2 −B = ω(C + E) sin (ωτ) + (D + F ) cos (ωτ), (4.26) ωA1 = (D + F ) sin (ωτ)− ω(C + E) cos (ωτ). (4.27) Squaring and adding the squares together, we obtain ω4 + xω2 + y = 0, (4.28) with x = A2 1 − 2B − (C + E)2, and y = B2 − (D + F )2. Letting z = ω2, equation (4.28), becomes the following equation z2 + xz + y = 0. (4.29) On the other hand, we have y = B2 − (D + F )2 = ( B +D + F )( B −D − F ) . Using the case τ1 = τ2 = 0, we have B + F +D > 0, and the assumption (4.24) implies that B −D − F < 0, then we have y < 0. Consequently, the equation (4.29) has a unique solution positive z0 = −x+ √ (x2 − 4y) 2 . Then the equation (4.28) has a unique positive solution ω0 = √ −x+ √ (x2 − 4y) 2 . Copyright © 2025 ASSA. Adv Syst Sci Appl (2025) STABILITY ANALYSIS OF A DELAYED EPIDEMIC MODEL... 73 On the other hand, the equations (4.26) and (4.27) imply that τ0 = 1 ω0 arccos ( (D + F )(ω2 0 −B)− ω2 0A1(C + E) ω2 0(C + E)2 + (D + F )2 ). Lets show that iω0 is simple, consider the branche of characteristic roots λ(τ) = k(τ) + ij(τ), of equation (4.25) bifurcating from iω0 at τ = τ0. By derivation (4.25) with respect to the delay τ, we obtain dλ dτ ( 2λ+ A1 + e−λτ ([C + E]− [λ(C + E) + (D + F )]τ) ) = λ[λ(C + E) + (D + F )]e−λτ . (4.30) If we suppose, by contradiction, that iω0 is not simple, the right hand side (4.30) gives iω0(C + E) + (D + F ) = 0. On the other hand, the equation (4.24) implies that B −D − F < 0, and by the hypotheses (H1), (H3) and (T1), we have B > 0,−D > 0 and F > 0 then we have −D < F. This contradicts −D = F. Next, we need to guarantee the transversality condition of the Hopf bifurcation theorem (see [23]). Clearly, λ(τ) = k(τ) + ij(τ) is a root of equation (4.25) if and only if k2 − j2 + A1k +B = −e−kτ ( [k(C + E) + (D + F )] cos (jτ) + j(C + E) sin (jτ) ) , (4.31) and 2kj + jA1 = −e−kτ ( j(C + E) cos (jτ)− [k(C + E) + (D + F )] sin (jτ) ) . (4.32) Let k(τ0) and j(τ0), satisfying k(τ0) = 0, and j(τ0) = ω0. By differentiating equations (4.31) and (4.32) with respect to τ and then set τ = τ0. Doing this, we get G1 dk(τ0) dτ +G2 dj(τ0) dτ = h1, −G2 dk(τ0) dτ +G1 dj(τ0) dτ = h2, where G1 = A1 + [(C + E)− τ0(D + F )] cos (ω0τ0)− τ0ω0(C + E) sin (ω0τ0), G2 = −2ω0 + τ0ω0(C + E) cos (ω0τ0) + [(C + E)− τ0(D + F )] sin (ω0τ0), h1 = ω0(D + F ) sin (ω0τ0)− ω2 0(C + E) cos (ω0τ0), and h2 = ω2 0(C + E) sin (ω0τ0) + ω0(D + F ) cos (ω0τ0). Also, we have dk(τ0) dτ = G1h1 −G2h2 G2 1 +G2 2 , Therefore, according to the equations (4.26) and (4.27) we have dk(τ0) dτ = ω2 0 √ x2 − 4y G2 1 +G2 2 . This completes the proof. Copyright © 2025 ASSA. Adv Syst Sci Appl (2025) 74 A. BERNOUSSI, K. HATTAF, B. EL BOUKARI 4.4. Case τ2 > 0 and τ1 = 0 Theorem 4.4: If R0 > 1 and (4.24) holds, then there exists a positive τ2,0 where system (1.3) undergoes Hopf bifurcation at P ∗ when τ2 = τ2,0. However, the steady state P ∗ is locally asymptotically stable when τ2 ∈ [0, τ2,0) and unstable when τ2 > τ2,0. Here τ2,0 is given by τ2,0 = 1 ω2,0 arccos [ ω2 2,0(F − E[A1 + C])− (B +D)F ω2 2,0E 2 + F 2 ], and ω2,0 = √ −p+ √ (p2 − 4q) 2 , where p and q are defined in equation (4.36). Proof When τ2 > 0 and τ1 = 0, the characteristic equation (4.18) becomes ∆(λ, τ2) = λ2 + (A1 + C)λ+B +D + (λE + F )e−λτ2 = 0. (4.33) Equation (4.33), has a purely imaginary root iω2, with ω2 > 0, ∆(iω2, τ2) = 0, if and only if ω2 2 − (B +D) = ω2E sin (ω2τ2) + F cos (ω2τ2), (4.34) ω2(A1 + C) = F sin (ω2τ2)− ω2E cos (ω2τ2). (4.35) Squaring and adding the squares together, we get ω4 2 + pω2 2 + q = 0, (4.36) with p = (A1 + C)2 − 2(B +D)− E2, and q = (B +D)2 − F 2. Letting z2 = ω2 2, equation (4.36) becomes the following equation z22 + pz2 + q = 0. (4.37) On the other hand, we have q = (B +D)2 − F 2 = ( B +D + F )( B +D − F ) . Using the case τ1 = τ2 = 0, we have B +D + F > 0, and the equation (4.24) implies that B −D < F, and by the hypotheses (H1), (H2), and (T1), we have B > 0, D < 0 and F > 0, then we have B +D < F, then, we have q < 0. Consequently, the equation (4.37) has a unique solution positive z2,0 = −p+ √ (p2 − 4q) 2 . Copyright © 2025 ASSA. Adv Syst Sci Appl (2025) STABILITY ANALYSIS OF A DELAYED EPIDEMIC MODEL... 75 Then the equation (4.36) has a unique positive solution ω2,0 = √ −p+ √ (p2 − 4q) 2 . On the other hand, the equations (4.34) and (4.35) imply that τ2,0 = 1 ω2,0 arccos [ ω2 2,0(F − E[A1 + C])− (B +D)F ω2 2,0E 2 + F 2 ]. Lets show that iω2,0 is simple, consider the branche of characteristic roots λ(τ2) = u(τ2) + iv(τ2), of equation (4.33) bifurcating from iω2,0 at τ2 = τ2,0. By derivation (4.33) with respect to the delay τ2, we obtain dλ dτ2 ( 2λ+ A1 + C + e−λτ2(E − [λE + F ]τ2) ) = λ(λE + F )e−λτ2 . (4.38) If we suppose, by contradiction, that iω2,0 is not simple, the right hand side (4.38) gives iω2,0E + F = 0, and leads a contradiction with the fact that ω2,0 > 0, E > 0 and F > 0. Next we need to guarantee the transversality condition of the Hopf bifurcation theorem (see [23]). Clearly, λ(τ2) = u(τ2) + iv(τ2), is a root of equation (4.33), if and only if u2 − v2 + (A1 + C)u+B +D = −e−uτ2 ( [uE + F ] cos (vτ2) + vE sin (vτ2) ) , (4.39) 2uv + v(A1 + C) = −e−uτ2 ( vE cos (vτ2)− (uE + F ) sin (vτ2) ) . (4.40) Let u(τ2,0) and v(τ2,0), satisfying u(τ2,0) = 0, and v(τ2,0) = ω2,0. By differentiating equations (4.39) and (4.40) with respect to τ2 and then set τ2 = τ2,0. Doing this, we get G3 du(τ2,0) dτ2 +G4 dv(τ2,0) dτ2 = h3, −G4 du(τ2,0) dτ2 +G3 dv(τ2,0) dτ2 = h4, where G3 = A1 + C + (E − τ2,0F ) cos (ω2,0τ2,0)− τ2,0ω2,0E sin (ω2,0τ2,0), G4 = −2ω2,0 + τ2,0ω2,0E cos (ω2,0τ2,0) + (E − τ2,0F ) sin (ω2,0τ2,0), h3 = ω2,0F (0) sin (ω2,0τ2,0)− ω2 2,0E cos (ω2,0τ2,0), and h4 = ω2 2,0E sin (ω2,0τ2,0) + ω2,0F (0) cos (ω2,0τ2,0). Further, we have du(τ2,0) dτ2 = G3h3 −G4h4 G2 3 +G2 4 . Therefore, according to the equations (4.34) and (4.35) we have du(τ2,0) dτ2 = ω2 2,0 √ p2 − 4q G2 3 +G2 4 . This completes the proof. Copyright © 2025 ASSA. Adv Syst Sci Appl (2025) 76 A. BERNOUSSI, K. HATTAF, B. EL BOUKARI 4.5. Case τ1 > 0 and τ2 > 0 Theorem 4.5: If R0 > 1, G5h5 −G6h6 ̸= 0 and (4.24) holds, then there exists a positive τ ∗2,0 where system (1.3) undergoes Hopf bifurcation at P ∗ when τ2 = τ ∗2,0. However, the steady state P ∗ is locally asymptotically stable when τ2 ∈ [0, τ ∗2,0) and unstable when τ2 > τ ∗2,0. Here τ ∗2,0 is given by τ ∗2,0 = 1 ω∗ 2,0 arccos[ (F − EA1)(ω ∗ 2,0) 2 − FB + ω∗ 2,0(ED − CF ) sin (ω∗ 2,0τ1) (ω∗ 2,0E)2 + F 2 − (DF + EC(ω∗ 2,0) 2) cos (ω∗ 2,0τ1) (ω∗ 2,0E)2 + F 2 ]. where ω∗ 2,0 is the positive root of the equation (4.43), and G5, G6, h5, h6 as defined in (4.48). Proof Equation (4.18) has a purely imaginary root iω∗ 2, with ω∗ 2 > 0, ∆(iω∗ 2, τ1, τ2) = 0, if and only if (ω∗ 2) 2 −B − ω∗ 2C sin (ω∗ 2τ1)−D cos (ω∗ 2τ1) = Eω∗ 2 sin (ω ∗ 2τ2) + F cos (ω∗ 2τ2), (4.41) and A1ω ∗ 2 + ω∗ 2C cos (ω∗ 2τ1)−D sin (ω∗ 2τ1) = F sin (ω∗ 2τ2)− Eω∗ 2 cos (ω ∗ 2τ2). (4.42) Squaring and adding the squares together, we obtain (ω∗ 2) 4 + (ω∗ 2) 2(C2 + A2 1 − E2) +B2 +D2 − F 2 + 2 ( [−(ω∗ 2) 3C + ω∗ 2(BC − A1D)] sin (ω∗ 2τ1) + [(ω∗ 2) 2(A1C −D) +BD] cos (ω∗ 2τ1) ) = 0. (4.43) We define h(ω∗ 2) = (ω∗ 2) 4 + (ω∗ 2) 2(C2 + A2 1 − E2) +B2 +D2 − F 2 + 2 ( [−(ω∗ 2) 3C + ω∗ 2(BC − A1D)] sin (ω∗ 2τ1) + [(ω∗ 2) 2(A1C −D) +BD] cos (ω∗ 2τ1) ) . (4.44) On the other hand, we have −1 ≤ sin (ω∗ 2τ1) ≤ 1 and − 1 ≤ cos (ω∗ 2τ1) ≤ 1, and by the hypotheses (H1), (H3) and (T1), we have BD < 0, A1C < 0,−D > 0,−C > 0, BC < 0 and − A1D > 0. Then g2(ω ∗ 2) ≤ h(ω∗ 2) ≤ g1(ω ∗ 2), where g1(ω ∗ 2) = (ω∗ 2) 4 − 2(ω∗ 2) 3C + (ω∗ 2) 2[C2 + A2 1 − E2 − 2A1C − 2D] − ω∗ 2(2BD + 2A1D) +B2 +D2 − F 2 − 2BD, Copyright © 2025 ASSA. Adv Syst Sci Appl (2025) STABILITY ANALYSIS OF A DELAYED EPIDEMIC MODEL... 77 and g2(ω ∗ 2) = (ω∗ 2) 4 + 2(ω∗ 2) 3C + (ω∗ 2) 2[C2 + A2 1 − E2 + 2A1C + 2D] + ω∗ 2(2BD + 2A1D) +B2 +D2 − F 2 + 2BD. We have lim ω∗ 2→+∞ g1(ω ∗ 2) = lim ω∗ 2→+∞ g2(ω ∗ 2) = +∞. Then we have lim ω∗ 2→+∞ h(ω∗ 2) = +∞, and h(0) = B2 +D2 + 2BD − F 2 = (B +D)2 − F 2 = ( B +D + F )( B +D − F ) . By the case τ1 = τ2 = 0, we have (B +D + F ) > 0, and by the case (τ2 > 0 and τ1 = 0) if the equation (4.24) is satisfied then we have (B +D − F ) < 0. Therefore, h(0) < 0. Since h is continuous in [0,+∞), then the equation (4.43) has at least one positive root. We assume that equation (4.43) admits a finite family of solution ω∗ 2,i with i = 1, 2, ..., n. n ∈ N. By the equations (4.41) and (4.42), we have τ li = 1 ω∗ 2,i ( arccos [ (F − EA1)(ω ∗ 2,i) 2 − FB + ω∗ 2,i(ED − CF ) sin (ω∗ 2,iτ1) (ω∗ 2,0E)2 + F 2 − (DF + EC(ω∗ 2,i) 2) cos (ω∗ 2,iτ1) (ω∗ 2,0E)2 + F 2 ] + 2Πl ) , l = 0, 1, ....; i = 1, 2, ..., n. Then ±iω∗ 2,i is a pair of purely imaginary roots of equations (4.18), with τ2 = τ li , l = 0, 1, ...; i = 1, 2, ..., n. Clearly, lim l→∞ τ li = ∞, i = 1, 2, ..., n. Thus, we can define τ ∗2,0 = τ l0i0 = min l=0,1...,i=1,2,...,n (τ li ), ω ∗ 2,0 = ω∗ 2,i0 . Lets show that iω∗ 2,0 is simple, consider the branche of characteristic roots λ(τ2) = u1(τ2) + iv1(τ2), of equation (4.18) bifurcating from iω∗ 2,0 at τ2 = τ ∗2,0. By derivation (4.18) with respect to the delay τ2, we obtain dλ dτ2 ( 2λ+ A1 + [C − τ1(λC +D)]e−λτ1 + [E − τ2(λE + F )]e−λτ2 ) = λ[λE + F ]e−λτ2 . (4.45) If we suppose, by contradiction, that iω∗ 2,0, is not simple, the right hand side (4.45) gives iω∗ 2,0E + F = 0, Copyright © 2025 ASSA. Adv Syst Sci Appl (2025) 78 A. BERNOUSSI, K. HATTAF, B. EL BOUKARI and leads a contradiction with the fact that ω∗ 2,0 > 0, E > 0 and F > 0. Next we need to guarantee the transversality condition of the Hopf bifurcation theorem (see [23]). Clearly λ(τ2) = u1(τ2) + iv1(τ2), is a root of equation (4.18) if and only if u2 1 − v21 + A1u1 + e−u1τ1 ( Cu1 cos (v1τ1) + Cv1 sin (v1τ1) +D cos (v1τ1) ) = −e−u1τ2 ( Eu1 cos (v1τ2) + Ev1 sin (v1τ2) + F cos (v1τ2) ) , (4.46) and 2u1v1 + A1v1 + e−u1τ1 ( Cv1 cos (v1τ1)− Cu1 sin (v1τ1)−D sin (v1τ1) ) = −e−u1τ2 ( Ev1 cos (v1τ2)− Eu1 sin (v1τ2)− F sin (v1τ2) ) . (4.47) Let u1(τ ∗ 2,0) and v1(τ ∗ 2,0), satisfying u1(τ ∗ 2,0) = 0, and v1(τ ∗ 2,0) = ω∗ 2,0. By differentiating equations (4.46) and (4.47) with respect to τ2, and then set τ2 = τ ∗2,0. Doing this, we get G5 du1(τ ∗ 2,0) dτ2 +G6 dv1(τ ∗ 2,0) dτ2 = h5, −G6 du1(τ ∗ 2,0) dτ2 +G5 dv1(τ ∗ 2,0) dτ2 = h6, (4.48) where G5 = A1 − τ1Cω∗ 2,0 sin (ω ∗ 2,0τ1) + (C − τ1D) cos (ω∗ 2,0τ1)− τ ∗2,0Eω∗ 2,0 sin (ω ∗ 2,0τ ∗ 2,0) + (E − τ ∗2,0F ) cos (ω∗ 2,0τ ∗ 2,0), G6 = −2ω∗ 2,0 + (C − τ1D) sin (ω∗ 2,0τ1) + τ1Cω∗ 2,0 cos (ω ∗ 2,0τ1) + (E − τ ∗2,0F ) sin (ω∗ 2,0τ ∗ 2,0) + τ ∗2,0Eω∗ 2,0 cos (ω ∗ 2,0τ ∗ 2,0), h5 = Fω∗ 2,0 sin (ω ∗ 2,0τ ∗ 2,0)− E(ω∗ 2,0) 2 cos (ω∗ 2,0τ ∗ 2,0), and h6 = E(ω∗ 2,0) 2 sin (ω∗ 2,0τ ∗ 2,0) + Fω∗ 2,0 cos (ω ∗ 2,0τ ∗ 2,0). Calculating du1(τ∗2,0) dτ2 , we get du1(τ ∗ 2,0) dτ2 = G5h5 −G6h6 G2 5 +G2 6 . This completes the proof. 5. NUMERICAL SIMULATION In this section, we shall give some simulations to illustrate the previous results. For the following, let’s consider the saturated incidence rate function: f(S, I)I = βSI 1 + α1I , Copyright © 2025 ASSA. Adv Syst Sci Appl (2025) STABILITY ANALYSIS OF A DELAYED EPIDEMIC MODEL... 79 and the treatment function T (I) = rI. We take the parameters of the system (1.3) as follows: A = 8, β = 0.955, µ = 0.5, γ2 = 0.935, α = 0.03 γ1 = 0.3, r = 0.9, and α1 = 0.5. Some parameters are for the spread of tuberculosis disease in Turkey from 2005 to 2015 (see, [38]). The rest of the parameters are the hypothetical set of parameter values. 5.1. Case τ1 = τ2 = 0 By applying Theorem 4.1, the endemic equilibrium P ∗ is locally asymptotically stable when R0 > 1 and τ1 = τ2 = 0. Figure 5.1 illustrates this result. Fig. 5.1. Dynamics of system (1.3) when τ1 = τ2 = 0 and R0 = 5.7336 > 1. 5.2. Case τ1 > 0 and τ2 = 0 According to Theorem 4.2, the endemic equilibrium P ∗ is locally asymptotically stable if R0 > 1, τ1 > 0 and τ2 = 0. Figure 5.2 demonstrates this result. Fig. 5.2. Dynamics of system (1.3) when τ1 = 1.7678, τ2 = 0 and R0 = 5.7336 > 1. Copyright © 2025 ASSA. Adv Syst Sci Appl (2025) 80 A. BERNOUSSI, K. HATTAF, B. EL BOUKARI 5.3. Case τ1 = τ2 = τ From Theorem 4.3, the following figures show that if the delay is below τ0 = 1.7678, the endemic equilibrium P ∗ is locally asymptotically stable (see, Figure 5.3a). When the value of the delay τ increases we lose the stability of the endemic equilibrium (see, Figure 5.3b) and P ∗ becomes unstable for τ > τ0 = 1.7678 (see, Figure 5.3c), and vice versa when the delay is decreasing the model converges rapidly to P ∗. (a) Dynamics of system (1.3) when τ1 = τ2 = 1.4178 and R0 = 5.7336 > 1. (b) For τ1 = τ2 = 1.7678, Hopf bifurcation occurs and periodic solutions appear for model (1.3) with R0 = 5.7336 > 1. (c) For τ1 = τ2 = 2.1178, all solutions (S, I) of model (1.3) are unstable with R0 = 5.7336 > 1. Fig. 5.3. Case τ1 = τ2 = τ Copyright © 2025 ASSA. Adv Syst Sci Appl (2025) STABILITY ANALYSIS OF A DELAYED EPIDEMIC MODEL... 81 5.4. Case τ1 = 0 and τ2 > 0 As in Theorem 4.4, the following figures show that if the delay is below τ2,0 = 0.8357, the endemic equilibrium P ∗ is locally asymptotically stable (see, Figure 5.4a). When the value of the delay τ2 increases we lose the stability of the endemic equilibrium (see, Figure 5.4b) and P ∗ becomes unstable for τ2 > τ2,0 = 0.8357 (see, Figure 5.4c), and vice versa when the delay is decreasing the model converges rapidly to P ∗. (a) Dynamics of system (1.3) when τ1 = 0, τ2 = 0.7857 and R0 = 5.7336 > 1. (b) For τ1 = 0 and τ2 = 0.8357, Hopf bifurcation occurs and periodic solutions appear for model (1.3) with R0 = 5.7336 > 1. (c) For τ1 = 0 and τ2 = 0.8651, all solutions (S, I) of model (1.3) are unstable with R0 = 5.7336 > 1. Fig. 5.4. Case τ1 = 0 and τ2 > 0 Copyright © 2025 ASSA. Adv Syst Sci Appl (2025) 82 A. BERNOUSSI, K. HATTAF, B. EL BOUKARI 5.5. Case τ1 > 0 and τ2 > 0 Based on Theorem 4.5, the following figures show that if the delay τ2 is below τ ∗2,0 = 0.9249 with τ1 fixed by taking τ1 = 3.6688, the endemic equilibrium P ∗ is locally asymptotically stable (see, Figure 5.5a). When the value of the delay τ2 increases we lose the stability of the endemic equilibrium (see, Figure 5.5b) and P ∗ becomes unstable for τ2 > τ ∗2,0 = 0.9249 (see ,Figure 5.5c), and vice versa when the delay is decreasing the model converges rapidly to P ∗. (a) Dynamics of system (1.3) when τ1 = 3.6688, τ2 = 0.8228 and R0 = 5.7336 > 1. (b) For τ1 = 3.6688 and τ2 = 0.9249, Hopf bifurcation occurs and periodic solutions appear for model (1.3) with R0 = 5.7336 > 1. (c) For τ1 = 3.6688 and τ2 = 0.9976, all solutions (S, I) of model (1.3) are unstable with R0 = 5.7336 > 1. Fig. 5.5. Case τ1 > 0 and τ2 > 0 6. CONCLUDING REMARKS In this work, we have proposed and analyzed a delayed SIR model with generalized incidence and treatment functions. Our analysis proved that the two delays τ1 and τ2 have a very big influence on the stability of the equilibrium points. In fact, • When R0 ≤ 1, we have: - If τ1 ≥ τ2, then the disease-free equilibrium is globally asymptotically stable for all Rc, this shows that the disease disappears. This biologically mean that when the latent period is greater than the healing period and the basic reproduction number is Copyright © 2025 ASSA. Adv Syst Sci Appl (2025) STABILITY ANALYSIS OF A DELAYED EPIDEMIC MODEL... 83 less than or equal to one, the infected individuals cured before becoming infectious and the disease will disappear. - If Rc ≤ 1, then the disease-free equilibrium is globally asymptotically stable independently of delays. - If Rc > 1 and τ1 is held fixed, then there exists a positive constant τ ∗22,0 such that when τ2 ∈ [0, τ ∗22,0), the disease-free equilibrium is locally asymptotically stable. When τ2 = τ ∗22,0 a Hopf bifurcation occurs and when τ2 > τ ∗22,0 the disease-free equilibrium is unstable. • When R0 > 1, the disease-free equilibrium is unstable and the proposed model admits an endemic equilibrium P ∗. In this case, we have: - If τ2 = 0 and τ1 ≥ 0, then the endemic equilibrium P ∗ is locally asymptotically stable. - If τ2 > 0, then there exists a positive constant τ̃ such that when τ2 ∈ [0, τ̃), the endemic equilibrium is locally asymptotically stable. When τ2 = τ̃ a Hopf bifurcation occurs and the host populations and the disease coexists within an oscillatory mode, and when τ2 > τ̃ the endemic equilibrium is unstable and the disease persists. From the aforementioned, we conclude that the longer the recovery period is for the infected individuals, the higher the spreading risk becomes and vice versa. ACKNOWLEDGMENTS The authors would like to thank the editor and the anonymous referees for their valuable comments and suggestions that have led to improvement of the quality of this work. REFERENCES 1. Abta, A., Kaddar, A. & Alaoui, H. T. (2012). Global stability for delay SIR and SEIR epidemic models with saturated incidence rates, Electronic Journal of Differential Equations. 2012(23),1-13, 2012. 2. Alshammari, F. S. & Akyildiz, F. T. (2021). Global Stability for Novel Complicated SIR Epidemic Models with the Nonlinear Recovery Rate and Transfer from Being Infectious to Being Susceptible to Analyze the Transmission of COVID-19, Journal of Function Spaces. Article ID 5207152, doi.org/10.1155/2021/5207152. 3. Anderson, R. M. & May, R. M. (1978). Regulation and stability of host-parasite population interactions: I. Regulatory processes, The Journal of Animal Ecology. 47(1), 219-267. 4. Avila-Vales, E.J. & Cervantes-Perez, A.G. (2019). Global stability for SIRS epidemic models with general incidence rate and transfer from infectious to susceptible, Boletin de la Sociedad Matematica Mexicana. 25, 637-658. 5. Balamuralitharan, S. & Radha, M. (2018). Bifurcation analysis in SIR epidemic model with treatment, National Conference on Mathematical Techniques and its Applications (NCMTA 18). IOP Conf. Series: Journal of Physics: Conf. Series 1000, (2018), 012169 doi :10.1088/1742-6596/1000/1/012169. 6. Beddington, J. R. (1975). Mutual interference between parasites or predators and its effect on searching efficiency, J. Animal Ecol. 44, 331-341. 7. Beretta, E., Hara, T., Ma, W. & Takeuchi, Y. (2002). Permanence of an SIR epidemic model with distributed time delays, Tohoku Math. J. 54, 581–591. 8. Beretta, E. & Takeuchi, Y. (1995). Global stability of an SIR epidemic model with time delays, J. Math. Biol. 33, 250–260. Copyright © 2025 ASSA. Adv Syst Sci Appl (2025) 84 A. BERNOUSSI, K. HATTAF, B. EL BOUKARI 9. Bernoussi, A. (2021). Bifurcation of periodic solution of a delayed SEIR epidemic model with nonlinear incidence rate, Journals of Applied Nonlinear Dynamics. 10(3), 351-367. 10. Bernoussi, A. (2023). Stability analysis of an SIR epidemic model with homestead-isolation on the susceptible and infectious, immunity, relapse and general incidence rate, International journal of biomathematics. 16(5), https: doi.org/10.1142/S1793524522501029. 11. Bernoussi, A. & Hattaf, K. (2021). Global dynamics of an SIRSI epidemic model with discrete delay and general incidence rate, Journals of Applied Nonlinear Dynamics. 10(3), 545-560. 12. Bernoussi, A. (2019). Global stability analysis of an SEIR epidemic model with relapse and general incidence rates, Journals of Applied Sciences (APPS). 21, 54-68. 13. Bjørnstad, N., Shea, K., Krzywinski, M. & Naomi, A. (2020). The SEIRS model for infectious disease dynamics, journal of Nature Methods. 17, 557–558. 14. Capasso, V. (1978). Global solution for a dikusive nonlinear deterministic epidemic model, In: SIAM Journal on Applied Mathematics. 35(2), 274-284. 15. Capasso, V. & Serio, G. (1978). A generalization of the Kermack-McKendrick deterministic epidemic model, In: Mathematical Biosciences. 42(1), 43-61. 16. Capasso, V., Grosso, E. & Serio, G. (1978). I modelli matematici nella indagine epidemiologica. Applicazione allepidemia di colera verificatasi in Bari nel 1973, In: Annali sclavo. 19, 193-208. 17. Cooke, K. L. (1979). Stability analysis for a vector disease model, Rocky Mountain J. Math. 9(1), 31–42. 18. de Jong, M.C.M., Diekmann, O. & Heesterbeek, H. (1996). How does transmission of infection depend on population size, In: Epidemic models: their structure and relation to data, Mollison D. (Ed.), Cambridge University Press, Cambridge., 84-94. 19. Elazzouzi, A., Alaoui, A. L., Tilioua, M. & Tridane, A. (2019). Global stability analysis for a generalized delayed SIR model with vaccination and treatment, journal of Advances in Difference Equations. 2019(1), 1-19. 20. Elkhaiar, S. & Kaddar, A. (2016). Stability Analysis of an SEIR Model with Treatment,Research in Applied Mathematics. 1, DOI: 10.11131/2017/101266. 21. Gabriela, M., Gomes, M., White, L. J. & Medley, G. F. (2005), The reinfection threshold, Journal of Theoretical Biology. 236, 111-113. 22. Guo, H., Li, M. Y. & Shuai, Z. (2006). Global stability of the endemic equilibrium of multigroup SIR epidemic models, In: Canadian Applied Mathematics Quarterly. 14, 259-284. 23. Hale, J.K. & Verduyn Lunel, S.M. (1993). Introduction to Functional Differential Equations, Springer- Ver- lag, New York. 24. Hattaf, K. & Yousfi, N. (2016). A class of delayed viral infection models with general incidence rate and adaptive immune response, Int. J. Dynam. Control. 4, 254-265. 25. Hattaf, K., Lashari, A. A., Louartassi, Y. & Yousfi, N. (2013). A delayed SIR epidemic model with a generalized incidence rate , Electronic Journal of Qualitative Theory of Differential Equations. 2013(3), 1-9. 26. Hethcote, H. W. (2000). The mathematics of infectious diseases, In: SIAM review. 42, 599-653. 27. Jiang, Z. & Wei, J. (2008). Stability and bifurcation analysis in a delayed SIR model, Chaos Solitons and Fractals. 35, 609-619. 28. Kaddar, A., Elkhaiar, S. & Eladnani, F. (2016). Global Asymptotic Stability of a Generalized SEIRS Epidemic Model, Differ Equ Dyn Syst. DOI 10.1007/s12591-016- 0313-y. 29. Kermack, W. O. & McKendrick, A. G. (1927). A contribution to the mathematical theory of epidemics, In: Proceedings of the Royal Society of London A: mathematical, physical and engineering sciences. 115, 772. The Royal Society. 700-721. Copyright © 2025 ASSA. Adv Syst Sci Appl (2025) STABILITY ANALYSIS OF A DELAYED EPIDEMIC MODEL... 85 30. Korobeinikov, A. & Maini, P. K. (2005). Non-linear incidence and stability of infectious disease models, In: Mathematical Medicine and Biology. 22(2) 113-128. 31. Kuang, Y. (1993). Delay Differential Equations with Applications in Population Dynamics, Academic Press, San Diego. 32. LaSalle, J.P. (1976). The stability of dynamical systems, Regional Conference Series in Applied Mathematics, SIAM, Philadelphia. 33. Liu, J. (2019), Bifurcation analysis for a delayed SEIR epidemic model with saturated incidence and saturated treatment function, Journal of Biological Dynamics. 13(1), 461- 480. 34. Mena-Lorcat, J. & Hethcote, H. W. (1992). Dynamic models of infectious diseases as regulators of population sizes, In: Journal of mathematical biology. 30, 693-716. 35. Nakata, Y., Enatsu, Y. & Muroya, Y. (2011). On the global stability of an SIRS epidemic model with distributed delays, Disc. Contin. Dyn. Syst. Suppl. 2011, 1119–1128. 36. Riad, D., Hattaf, K. & Yousfi, N. (2016). Dynamics of Capital-labour Model with Hattaf-Yousfi Functional Response, Journal of Advances in Mathematics and Computer Science. 18, 1-7. 37. Tang, Q., Teng, Z. & Abdurahman, X. (2017). A New Lyapunov Function for SIRS Epidemic Models, In: Bulletin of the Malaysian Mathematical Sciences Society. 40(1), 237-258. 38. Ucakan, S., Gulen, S. & Koklu, K. (2021). Analysing of Tuberculosis in Turkey through SIR, SEIR and BSEIR Mathematical Models, journal of Mathematical and Computer Modelling of Dynamical Systems. 27(1), 179-202. 39. Wang, W. (2006). Backward bifurcation of an epidemic model with treatment, Mathematical Biosciences. 201(1-2), 58–71. 40. Wang, W. & Ruan, S. (2004). Bifurcation in epidemic model with constant removal rate infectives, Journal of Mathematical Analysis and Applications. 291, 775-793. 41. Xu, R. & Ma, Z. (2009). Stability of a delayed SIRS epidemic model with a nonlinear incidence rate, Chaos, Soliton. Fract. 41(5), 2319-2325. 42. Zhang, J., Jia, J. & Song, X. (2014), Analysis of an SEIR epidemic model with saturated incidence and saturated treatment function. Scientific World Journal. 2014;2014:910421. doi: 10.1155-2014-910421. 43. Zhang, F., Li, Z. Z. & Zhang, F. (2008). Global stability of an SIR epidemic model with constant infectious period, Applied Mathematics and Computation. 199, 285-291. 44. Zhou, Y. & Liu, H. (2003). Stability of periodic solutions for an SIS model with pulse vaccination, Mathematical and Computer Modelling. 38, 299-308. Copyright © 2025 ASSA. Adv Syst Sci Appl (2025) Introduction Mathematical analysis Global stability of the disease-free equilibrium Stability and Hopf bifurcation of the endemic equilibrium Case 1 = 2 = 0 Case 1 > 0 and 2 = 0 Case 1 = 2 = Case 2 > 0 and 1 = 0 Case 1 > 0 and 2 > 0 Numerical simulation Case 1 = 2 = 0 Case 1 > 0 and 2 = 0 Case 1 = 2 = Case 1 = 0 and 2 > 0 Case 1 > 0 and 2 > 0 Concluding remarks