key: cord-0839870-42pmlclt authors: Zhang, Yue; Li, Xue; Zhang, Xianghua; Yin, Guisheng title: Stability and Hopf Bifurcation Analysis of an Epidemic Model with Time Delay date: 2021-07-01 journal: Comput Math Methods Med DOI: 10.1155/2021/1895764 sha: 2c13f21612965e495f46717a4a8b6df3598ab4d1 doc_id: 839870 cord_uid: 42pmlclt Epidemic models are normally used to describe the spread of infectious diseases. In this paper, we will discuss an epidemic model with time delay. Firstly, the existence of the positive fixed point is proven; and then, the stability and Hopf bifurcation are investigated by analyzing the distribution of the roots of the associated characteristic equations. Thirdly, the theory of normal form and manifold is used to drive an explicit algorithm for determining the direction of Hopf bifurcation and the stability of the bifurcation periodic solutions. Finally, some simulation results are carried out to validate our theoretic analysis. Today, the serious epidemics, such as SARS and H1N1, are still threatening the life of people continually. Plenty of mathematical models have been proposed to analyze the spread and the control of these diseases [1] [2] [3] [4] [5] [6] [7] . However, many infectious diseases, for instance, gonorrhea and syphilis, occur and spread amongst the mature, while some epidemics, for example, chickenpox and FMD, only result in infection and death in immature. For this reason, stage structure should be taken into consideration in models. Aliello and Freedman [8] proposed a stagestructured model described by where xðtÞ is the immature population density and yðtÞ represents the density of the mature population. α, γ, τ, and β are all positive constants. α is the birth rate, and γ is the natural death rate; τ is the time from birth to maturity; β is the death rate of the mature because of the competition with each other. And then, many infectious diseases with sage structure have been built and investigated [9] [10] [11] [12] [13] [14] . Xiao and Chen [15] improved (1) by separating the population into mature and immature and supposing that only the immature were susceptible to the infection. Based on the model in [15] , supposing that only the mature were susceptible, Jia and Li [16] built a new one as follows: where xðtÞ, yðtÞ, and RðtÞ are the susceptible, infectious, and recovered mature population densities, respectively; zðtÞ denotes the immature population density. All the parameters are positive constants. α, β, γ, and τ are the same as those in (1) ; m is the transmission coefficient describing the infection between the susceptible and the infectious; c is the death rate because of the epidemic; g is the recovery rate; αe −γτ xðt − τÞ denotes the population who were born at t − τ and survive at t. In systems (1) and (2), the time delay was also taken into consideration. Indeed, time delay plays an important role in the epidemic system, making the models more accurate. In recent years, delays have been introduced in more and more epidemic and predator-prey systems [17] [18] [19] . In this paper, on the basis of (2), we further assume that (1) Both the susceptible and the infectious have fertility, while in (2) , only the susceptible is fertile (2) For the infectious, there is competition with all the susceptible and the infectious, while for the susceptible, there is only competition between generations Meanwhile, all the death of the susceptible, the same as that in (1), is only due to the competition. To simplify model (2), we denote γ + c + g = d, and let bðxðtÞ + yðtÞÞ present the transmission from immature to mature. As a consequence, the new epidemic model could be described as follows: where w is the death rate of the mature because of the competition. We can notice that zðtÞ depends on xðtÞ and yðtÞ and R ðtÞ depends on yðtÞ; however, xðtÞ and yðtÞ have nothing to do with zðtÞ and RðtÞ. According to Qu and Wei [20] , we will mainly focus on xðtÞ and yðtÞ, that is, The rest of the paper is organized as follows. In Section 2, we calculate the steady states of system (4) and prove the existence and uniqueness of the positive equilibrium in particular. And then, the stability of the two nonzero equilibria and the existence of the Hopf bifurcation are investigated in Sections 3 and 4, respectively. In Section 5, the direction and stability of the Hopf bifurcation at the positive equilibrium are studied by using the center manifold theorem and the normal form theory [21] . And in the last section, some numerical simulations are carried out to validate the theoretical analysis. In this section, we discuss the existence of the equilibria of (4) and the positive one in particular. The equilibria are the solutions of the equations (5), Clearly, E 1 ð0, 0Þ and E 2 ðb/w, 0Þ are two equilibria of (4). In the following, we will focus on the existence of the positive equilibrium. (4) has one positive equilibrium E 3 ðx * , y * Þ, where Proof. Positive equilibrium is the positive solution of the equations (7), From the second equation of (7), we have Taking (8) into the first equation of (7), we can obtain which leads to where Together with we can know that both of the two solutions of (9) are positive, where ð13Þ Computational and Mathematical Methods in Medicine then, So, Therefore, if bðm − wÞ > dw, (4) has the unique positive equilibrium E 3 ðx * , y * Þ. 3. Stability Analysis of the Equilibrium E 2 ðb/w, 0Þ In this section, we analyze the stability of the equilibrium E 2 ðb/w, 0Þ. For convenience, the new variables uðtÞ = xðtÞ − b/w and vðtÞ = yðtÞ are introduced, and then, around E 2 ðb/w, 0Þ, the system (4) could be linearized as (18) : whose characteristic equation is given by from which, we can get that or Obviously, if bðm − wÞ > dw, then which implies that the equilibrium E 2 is unstable. If bðm − wÞ < dw, then As a consequence, we will discuss the other roots of (19) , that is, the roots of (21), under the condition bðm − wÞ < dw. For τ = 0, equation (21) becomes whose root is For τ > 0, if iωðω > 0Þ is a root of (21), then (27) can be obtained by separating the real and the imaginary parts, which leads to 3 Computational and Mathematical Methods in Medicine from which, we can get the unique positive root Let Then, when τ = τ j , (21) has a pair of purely imaginary roots ±iω 0 . Suppose which is the root of (21) such that To investigate the distribution of λðτÞ, we will discuss the trend of λðτÞ at τ = τ j . Substituting λðτÞ into (21) and taking the derivative with respect to τ, we can get which yields, dλ dτ Together with (27), we have sign Re dλ dτ which means that when undergoes τ = τ j , λðτÞ will add a pair of roots with positive real parts. That is, with the increase of τ, the number of roots with positive real part is increasing, leading to the change of the stability of the system (4). Therefore, the distribution of the roots of (21) could be obtained. Let ω 0 and τ j ( j = 0, 1, 2, ⋯ ) be defined by (29) and (30), respectively. (1) If bðm − wÞ > dw, then (19) has at least one positive root (2) If bðm − wÞ < dw, and τ = 0, then both roots of (19) are negative (3) If bðm − wÞ < dw, and τ > 0, then (19) has a pair of simple imaginary roots ±iω 0 at τ = τ j ; furthermore, if τ < τ 0 , then all the roots of (19) have negative real parts; if τ ∈ ðτ j , τ j+1 Þ, (19) has 2ðj + 1Þ roots with positive real parts Together with condition (35), the Hopf bifurcation theorem [21] , and Lemma 2, the following theorem could be stated. In this section, we analyze the stability of the positive equilibrium E 3 ðx * , y * Þ. For convenience, the new variables uðtÞ = xðtÞ − x * and vðtÞ = yðtÞ − y * are introduced, and then, around E 3 ðx * , y * Þ, the system (4) could be linearized as (36): whose characteristic equation is given by where For τ = 0, equation (37) becomes Firstly, computing λ 1 λ 2 , we have where Δ is the same as that in (10). So, which implies that the real parts of λ 1 and λ 2 have the same signs. Then, (λ 1 + λ 2 ) is calculated: Together with (41), we can get that both the real parts of the two roots of (39) are negative. For τ > 0, equation (37) can be rewritten as where If iωðω > 0Þ is a root of (37), then Separating the real and the imaginary parts, we have which leads to Let z = ω 2 > 0, and then, (47) can be rewritten as Firstly, computing z 1 z 2 , we get where has been proved in (41). Then, we will calculate If H(4-1): and then, which implies that (48) has one unique positive solution z 0 = ω 2 0 , where If H(4-2): ðm − wÞ ffiffiffi ffi Δ p − wB > 0, then and then, Let If Δ 1 < 0, then (48) has no real roots. If Δ 1 > 0 and z 1 + z 2 = 2p − r 2 + s 2 < 0, then (48) has two negative roots, and there no positive ω for (47); If Δ 1 > 0 and z 1 + z 2 = 2p − r 2 + s 2 > 0, (1) If ðm − wÞ ffiffiffi ffi Δ p − wB < 0, then (47) has one positive root ω 0 In the following, we will discuss the expression of τ j . where a > 0, t > 0, then (4) If f ðaÞ < 0, gðaÞ < 0, then In conclusion, If ðm − wÞ ffiffiffi ffi Δ p − wB < 0, substituting ω 0 defined in (54) into (64), we can get f ðω 0 Þ and gðω 0 Þ. Together with Lemma 5, the expression of τ j could be obtained. If If That is, when τ = τ j , the characteristic equation (37) has a pair of purely imaginary roots ±iω 0 . Suppose λðτÞ = αðτÞ + iωðτÞ is the root of (37), and then, we have To investigate the distribution of the λðτÞ, we will discuss the trend of λðτÞ at τ = τ j . Substituting λðτÞ into (37) and taking the derivative with respect to τ, we can get Together with (46) and (54), we have sign Re dλ dτ which means that when undergoes τ = τ j , λðτÞ will add a pair of roots with positive real parts. That is, with the increase of τ, the number of roots with positive real part is increasing, leading to the change of the stability of the system (4). If ðm − wÞ ffiffiffi ffi Δ p − wB > 0, Δ 1 > 0 and z 1 + z 2 = 2p − r 2 + s 2 > 0, then f ðω 0 Þ and gðω 0 Þ could be calculated by substituting ω ± defined in (58) into (64). According to Lemma 5, we can get the expression of τ ± j : If If That is, when τ = τ ± j , the characteristic equation (37) has a pair of purely imaginary roots ±iω ± . Let λðτÞ = αðτÞ + iωðτÞ be the root of (37), satisfying To investigate the distribution of the λðτÞ, we will discuss the trend of λðτÞ at τ = τ ± j . Using the same method, we have sign Re dλ dτ This implies that and dλ dτ which means that when undergoes τ = τ + j , λðτÞ will add a pair of roots with positive real parts, while undergoes τ = τ − j , λðτÞ will lose a pair of roots with positive real parts; if τ − j > τ + j+1 , then the characteristic equation (37) must have roots with positive real parts for τ > τ + j+1 . In conclusion, the distribution of the roots of (37) could be obtained. Let ω 0 be defined by (54), and τ j ( j = 0, 1, 2, ⋯ ) be defined by (66) or (68), and ω ± be defined by (58), and τ ± j ( j = 0, 1, 2, ⋯ ) be defined by (74) or (76), respectively. (1) When ðm − wÞ ffiffiffi ffi Δ p − wB > 0, if Δ 1 < 0, or z 1 + z 2 = 2p − r 2 + s 2 < 0, then all the roots of (37) are with negative real parts (2) When ðm − wÞ ffiffiffi ffi Δ p − wB < 0, (37) has a pair of simple imaginary roots ±iω 0 at τ = τ j ; furthermore, if τ ∈ ½0, τ 0 Þ, then all the roots of (37) are with negative real parts; if τ ∈ ðτ j , τ j+1 Þ, then (37) has 2ðj + 1Þ roots with positive real parts (3) When ðm − wÞ ffiffiffi ffi Δ p − wB > 0, if Δ 1 > 0, and z 1 + z 2 = 2p − r 2 + s 2 > 0, then (37) has a pair of simple imaginary roots ±iω ± at τ = τ ± j ; if τ ∈ ½0, τ + 0 Þ or τ ∈ ðτ − j , τ + j+1 Þ, then all the roots of (37) are with negative real parts; if τ ∈ ðτ + j+1 , τ − j+1 Þ, then (37) has a pair of roots with positive real parts; if τ − h > τ + h+1 , for τ > τ + h+1 , (37) must have roots with positive real parts Together with conditions (72) and (78), the Hopf bifurcation theorem [21] , and Lemma 6, the following theorem could be stated. Let τ j ( j = 0, 1, 2, ⋯ ) be defined by (66) or (68), and τ ± j ( j = 0, 1, 2, ⋯ ) be defined by (74) or (76), respectively, and then we have (4) is asymptotically stable (2) When ðm − wÞ ffiffiffi ffi Δ p − wB < 0, then the equilibrium E 3 ðx * , y * Þ is asymptotically stable when τ ∈ ½0, τ 0 Þ, and it is unstable when τ > τ 0 . System (4) undergoes a Hopf bifurcation at the equilibrium E 3 ðx * , y * Þ for +∞Þ. We call that system (4) undergoes ðh + 1Þ switches System (4) undergoes a Hopf bifurcation at the equilibrium E 3 ðx * , y * Þ for τ = τ ± j . Bifurcation at E 3 ðx * , y * Þ In the previous section, we have already gotten some conditions making that the system (4) undergoes a Hopf bifurcation at the positive equilibrium E 3 ðx * , y * Þ when τ = τ ± j , j = 0, 1, 2, ⋯ . In this section, under the conditions in Theorem 7, the direction of Hopf bifurcation and stability of the periodic solutions from E 3 will be investigated by using the center manifold and normal form theories [21] . Without loss of the generality, let τ be the critical value of τ = τ + j ðτ − j Þ, j = 0, 1, 2, ⋯ . For convenience, let τ = τ + ρ, ρ ∈ R, uðtÞ = xðτtÞ − x * , vðtÞ = yðτtÞ − y * , and then system (4) undergoes Hopf bifurcation at ρ = 0; with τ normalized by the time scaling t ⟶ t/τ, (4) could be rewritten as Choose the space as C = Cð½−1, 0, R 2 Þ; for any Φ = ðΦ 1 , and where s and n are the same as those in (37). According to Riesz's representation theorem, there is a 2 × 2 matrix ηðθ, ρÞ (θ ∈ ½−1, 0), whose components are bounded variation functions such that In fact, ηðθ, ρÞ could be chosen as where For any Φ ∈ C 1 ð½−1, 0, R 2 Þ, we define and Computational and Mathematical Methods in Medicine And then, system (4) could be translated into where For and a bilinear form where η θ ð Þ = η θ, 0 ð Þ: ð93Þ According to (87) and (91), we can get that Að0Þ and A * are adjoint operators. From the analysis in Section 4, we know that ±iω τ are a pair of eigenvalues of Að0Þ and also eigenvalues of A * , where ω is ω + or ω − defined in (58). It is easy to verify that vectors qðθÞ and q * ðsÞ are the eigenvalues of Að0Þ and A * corresponding to the eigenvalues iω τ and −iω τ, where For convenience, let then Taking (5-18) into (5) (6) (7) (8) (9) (10) (11) (12) (13) (14) (15) (16) (17) , we obtain By (111) and (112), the following can be gotten _ W 20 θ ð Þ = 2iω τW 20 θ ð Þ + g 20 q θ ð Þ + g 20 q θ ð Þ, _ W 11 θ ð Þ = g 11 q θ ð Þ + g 11 q θ ð Þ: Computational and Mathematical Methods in Medicine Because by integration, we have and where E 1 and E 2 are unknown. From (112) and the definition of A, we obtain Then, together with (108), we get Substituting (115) and (119) That is from which, E 1 can be determined. Figure 4 : The dynamic behaviors of system (4) with τ = 23 ∈ ðτ + 1 , τ − 1 Þ. (a) and (b) are the waveform plot and phase for system (4), respectively. E 3 ðx * , y * Þ is unstable with the emergence of oscillation and the initial condition is ð0:2, 0:2Þ. Figure 3 : The dynamic behaviors of system (4) with τ = 11 ∈ ðτ − 0 , τ + 1 Þ. (a) and (b) are the waveform plot and phase for system (4), respectively. E 3 ðx * , y * Þ is asymptotically stable and the initial condition is ð0:1,0:06Þ. Computational and Mathematical Methods in Medicine By (115), (117), and (119), we have That is from which, E 2 can be determined. Figure 6 : The dynamic behaviors of system (4) with τ = 36 ∈ ðτ + 2 ,+∞Þ. (a) and (b) are the waveform plot and phase for system (4), respectively. E 3 ðx * , y * Þ is unstable with the emergence of oscillation and the initial condition is ð0:15, 0:14Þ. Figure 5 : The dynamic behaviors of system (4) with τ = 30 ∈ ðτ − 1 , τ + 2 Þ. (a) and (b) are the waveform plot and phase for system (4), respectively. E 3 ðx * , y * Þ is asymptotically stable and the initial condition is ð0:15, 0:02Þ. Short-term predictions and prevention strategies for COVID-19: a modelbased study Fractional model for the spread of COVID-19 subject to government intervention and public perception Dynamics of a novel nonlinear stochastic SIS epidemic model with double epidemic hypothesis Dynamics of an ultradiscrete SIR epidemic model with time delay A stochastic SIRS epidemic model with nonlinear incidence rate A contribution to the mathematical theory of epidemics Stability analysis of delayed SIR epidemic models with a class of nonlinear incidence rates A time-delay model of single-species growth with stage structure A delayed epidemic model with stage-structure and pulses for pest management strategy Stability of Hopf bifurcation of a delayed SIRS epidemic model with stage structure Stability and Hopf bifurcation analysis of an eco-epidemic model with a stage structure Dynamics of an epidemic model with delays and stage structure Global stability of an agestructure epidemic model with imperfect vaccination and relapse Dynamics of an SIR epidemic model with stage structure and pulse vaccination An SIS epidemic model with stage structure and a delay Qualitative analysis of an SIR epidemic model with stage structure Global Hopf bifurcation and permanence of a delayed SEIRS epidemic model Stability and Hopf bifurcation of fractional-order complex-valued single neuron model with time delay Hopf bifurcation analysis in a diffusive predator-prey system with delay and surplus killing effect Bifurcation analysis in a time-delay model for prey-predator growth with stage-structure Theory and Applications of Hopf Bifurcation The authors declare that they have no competing interests. All authors conceived the study, carried out the proofs, and approved the final manuscript. This work is supported by the Fundamental Research Funds for the Central Universities under grant no. 3072021CF0609. We chooseand then,Using the same method as that in [21] , the center manifold C 0 at ρ = 0 is first computed. Suppose that x t is the solution of (81) when ρ = 0, and defineThen, on the center manifold C 0 , we can get thatwhereIn the direction q * and q * , z and z are local coordinates for center manifold C 0 . It is not difficult to note that when x t is real, W is real. And the real solutions are considered only.For any solution x t ∈ C 0 , since ρ = 0, we can get thatwhereWe rewrite Together with (83), we haveComparing with (105), we getAll the numbers in the expressions of g 20 , g 11 , and g 02 are known; however, W 20 and W 11 in g 21 are unknown, which will be computed in the following.From (89) and (99),whereWhen θ ∈ ½−1, 0,Comparing the coefficients in the expansions of (109) and (110), we have (115) and (116), W 20 and W 11 could be obtained; furthermore, g 21 can be calculated.Then, the important parameters can be obtained:which determine the quantities of the bifurcation of system (4) at τ = τ, where μ 2 determines the direction of the bifurcation: Hopf bifurcation is subcritical (supercritical) if μ 2 < 0 (μ 2 > 0) and the bifurcating periodic solutions exist for τ < τ (τ > τ); β 2 determines the stability of the bifurcating periodic solutions, which are stable (unstable) if β 2 < 0 (β 2 > 0); T 2 determines the period of the periodic solutions, which decreases (increases) if T 2 < 0 (T 2 > 0).From (78) in Section 4, we know that Re λ ′ ð τÞ > 0 when τ = τ + j , and Re λ ′ ð τÞ < 0 when τ = τ − j .Theorem 8. If ðm − wÞ ffiffiffi ffi Δ p − wB < 0, then the Hopf bifurcation of system (4) at positive equilibrium E 3 ðx * , y * Þ when τ = τ j is supercritical (subcritical) and the bifurcating periodic solutions on the manifold are stable (unstable) if Re ½C 1 ð0Þ < 0 (Re ½C 1 ð0Þ > 0).− wB > 0, Δ 1 > 0, and z 1 + z 2 = 2 p − r 2 + s 2 > 0, then when τ = τ + j , the Hopf bifurcation of system (4) at the positive equilibrium E 3 ðx * , y * Þ is supercritical (subcritical) and the bifurcating periodic solutions on the manifold are stable (unstable) if Re ½C 1 ð0Þ < 0 (Re ½C 1 ð0Þ > 0); when τ = τ − j , the Hopf bifurcation of system (4) at the positive equilibrium E 3 ðx * , y * Þ is subcritical (supercritical) and the bifurcating periodic solutions on the manifold are stable (unstable) if Re ½C 1 ð0Þ < 0 (Re ½C 1 ð0Þ > 0). In this section, some numerical simulations are carried out to support our theoretical analysis.There are so many different cases that only the most particular one in (3) of Theorem 7 is considered in this section. The coefficients are chosen as follows: b = 0:4, d = 0:8, w = 1, and m = 6. Then, the conditions of (3) where τ − 2 > τ + 3 . From Theorem 7, we know that the positive equilibrium E 3 ð0:169906, 0:049528Þ should be asymptotically stable when τ ∈ ½0, τ + 0 Þ ∪ ðτ − 0 , τ + 1 Þ ∪ ðτ − 1 , τ + 2 Þ, and it is unstable when τ ∈ ðτ + 0 , τ − 0 Þ ∪ ðτ + 1 , τ − 1 Þ ∪ ðτ + 2 ,+∞Þ. The system (4) undergoes 3 switches.All the simulation results for τ in the six different intervals are in consistent with the theoretical analysis, which are shown in Figures 1-6 . The data used to support the findings of the study are included within the article.