key: cord-0845627-0h9e422r authors: Al-Darabsah, Isam title: A Time-delayed SVEIR Model for Imperfect Vaccine with a Generalized Nonmonotone Incidence and Application to Measles date: 2020-10-01 journal: Appl Math Model DOI: 10.1016/j.apm.2020.08.084 sha: d5d2610d9221881c9c08205e03b53dfede2babd6 doc_id: 845627 cord_uid: 0h9e422r In this paper, we investigate the effects of the latent period on the dynamics of infectious disease with an imperfect vaccine. We assume a general incidence rate function with a non-monotonicity property to interpret the psychological effect in the susceptible population when the number of infectious individuals increases. After we propose the model, we provide the well-posedness property by verifying the non-negativity and boundedness of the models solutions. Then, we calculate the effective reproduction number [Formula: see text]. The threshold dynamics of the system is obtained with respect to [Formula: see text]. We discuss the global stability of the disease-free equilibrium when [Formula: see text] and explore the system persistence when [Formula: see text]. Moreover, we prove the coexistence of an endemic equilibrium when the system persists. Then, we discuss the critical vaccination coverage rate that is required to eliminate the disease. Numerical simulations are provided to: (i) implement a case study regarding the measles disease transmission in the United States from 1963 to 2016; (ii) study the local and global sensitivity of [Formula: see text] with respect to the model parameters; (iii) discuss the stability of endemic equilibrium; and (iv) explore the sensitivity of the proposed model solutions with respect to the main parameters. imperfect is proposed in [6] . Due to the inherent complexity of epidemiological transmission, other works studied epidemiological models with random and stochastic processes [7] [8] [9] [10] [11] [12] . For example, in [12] , the author wrote a primer to introduce the topic of stochastic epidemic modeling from the perspective of an ODE epidemic framework where some mathematical methods for formulation and numerical simulation of stochastic epidemic models were presented. In [9] and [10] , the authors provided a full probabilistic description of the random SIS-type and SI-type epidemiological model, respectively. In [13] , the authors studied stochastic SIS epidemic model with vaccination. To integrate the effect of behavioral changes on the disease spreading dynamics, the authors in [14] introduce a nonlinear incidence rate of the form where S and I represent the numbers of susceptible and infectious individuals, respectively, β is the probability of transmission per contact per unit time, and the constant q 1 is positive while the constants q 2 , d are nonnegative. Here, the constant d measures the inhibitory effect. In (1) , βI q 1 measures the infection force of the disease and 1 (1 + dI q 2 ) represents the inhibition effect from the behavioral change of the susceptible individuals when the number of infectious individuals increases [15] . It also can be used to describe the crowding effect of infectious individuals [16] . There are three types of the infection force function g(I) = βI q 1 (1 + dI q 2 ) based on the values of q 1 and q 2 [15, 17] : (i) unbounded incidence function: q 1 > q 2 ; (ii) saturated incidence function: q 1 = q 2 ; and (iii) nonmonotone incidence function: q 1 < q 2 . In the literature, these types have been used in different scenarios. For example, the authors in [15] consider a nonmonotone incidence rate to represent the psychological effect with q 1 = 1 and q 2 = 2. In [16] , a saturated incidence rate is considered with q 1 = q 2 = 2. For more details and examples we refer the reader to [17] and references within. In this paper, we consider a general form the infection force function g(I), in the sense that g(I) is increasing when the population of infectious I is small and decreasing when I is large. From a biological point of view, this can be interpreted as the "psychological effect", that is, when a disease is spreading among a population and the number of infective individuals becomes large, the behavior of the population may tend towards reducing the number of contacts among individuals per unit time [15] . In 1979, Cooke introduced a "time delay" to represent the disease incubation period in studying the spread of an infectious disease transmitted by a vector in [18] . Since then, many authors have incorporated time delays in epidemic models in different scenarios, such as incubation period or latent period [19] [20] [21] [22] [23] , vaccination period [24] , asymptomatic carriage period [19] and immune period [21, 22] . More precisely, in [21] , the authors study an SEIRS epidemic model with a constant latent and immune periods. In [20] , a latent period and relapse are considered in a general mathematical model for disease transmission. In [22] , the authors proposed an SEIRS disease transmission model with latent and temporary immune periods as a distributed delays. In [19] , a disease transmission model with two delays in incubation and asymptomatic carriage periods is investigated. In addition, many authors studied time delayed epidemic models with vaccination [2, [25] [26] [27] . For example, the authors in [2] study a vaccination model with a time delay to represent the time that an unaware susceptible individual takes to become aware of the infection. In the real world, the latent period may vary from days, as in influenza A H1N1 [28] , to years, as in AIDS [29] . The latent period has a profound effect on the generation time, which is defined as the time period between a case becoming infected and its subsequent infection of another case [30] . Thus, the latent period has an influence on the epidemic growth [31] . The goal in this paper is to study the impact of the latent period on the epidemic in the presence of imperfect vaccine when the infection incidence rate is in a nonlinear general form. The paper is organized as follows. In the next section, we present the mathematical model and the study its well-posedness. Then, in Section 3, we calculate the effective reproduction number R E . In Section 4, we use various methods, such as the method of fluctuations, the theory of limiting systems and chain transitive sets, to establish the threshold dynamics with respect to R E . We show that the disease-free equilibrium is global stable when when R E < 1 and the coexistence of endemic equilibrium and system persistence when R E > 1. Then, we discuss the critical vaccination coverage in Section 5. In the numerical simulations section, Section 6, we consider an application to the measles disease transmission. Furthermore, we study the local and global sensitivities of R E , discuss the stability of the endemic equilibrium and examine the sensitivity of the model solutions with respect to the main parameters. Finally, we discuss our results in Section 7. We provide the details proofs in Section 8. Motivated by [5, 6] , we let N (t) be the total population at time t and divide it into five classes: S(t) is the susceptible population, V (t) is the population of vaccinated individuals, E(t) is the population of individuals who are infected but not yet infectious (exposed class), I(t) is the infectious population and R(t) is the recovered population. We make the following assumptions to describe the interaction among the five classes: • The individuals in S move to V when they are vaccinated. The infected individuals in S move to the exposed class E with transmission rate β, and remain there for a certain latent period τ (see e.g., [32] ); • The individuals in V are not in a fully protective level and not perfect. When the vaccinated individuals become infected, they move into E for τ period. We assume that this happens at a reduced transmission rate ξβ, where ξ ∈ [0, 1] is the reduction coefficient [32] . When individuals gain immunity, they move into the class R [5]; • After the latent period τ , individuals become infectious and move from E to I. Then after recovery, they move to R; • There is a recruitment rate into the population and the natural death rate is the same for all the compartments with an additional disease-related death rate in I. Theses assumptions lead to the following system of delay differential equations: The term e −µτ represents the probability for individuals to survive the latent period [t − τ, t], see [19] . The population flux among the five compartments is given in Figure 1 and a description of the parameters is given in Table 1 . In the system (2), the term βI f (I) is the infection force function which is the probability that a susceptible is infected in unit time t. To ensure a non-monotonic infection force, we consider the following assumptions [4, 33] : There is a ζ > 0 such that d dI > 0 for 0 < I < ζ and d dI Here, ζ is the highest infection level, that is, the incidence rate is increasing when I ∈ (0, ζ] and it is decreasing when I > ζ. The critical level ξ can be calculated from the first derivative of I f (I) . In the numerical simulations section 6, we present different examples and calculate the corresponding ζ, see Figure 5 . From an epidemiological point of view, the function I f (I) in assumption (Q 1 ) represents the effect of intervention strategies on the reduction. The assumption (Q 2 ) means the infection force function decreases when the number of infectious individuals increases since the individuals tend to reduce the number of contacts among them per unit time when they are under intervention polices. When ζ is relatively high, the behavior of susceptible individuals changes after the appearance of a high number of infection cases. Consequently, the nonlinear incidence rate is S βI f (I) + V ⋅ ξ ⋅ βI f (I) which is the total number of newly infected individuals resulting per unit time at time t [34] . Denote the Banach space Recruitment rate of susceptible humans β Effective contact rate ψ Vaccination coverage rate µ Natural mortality rate ξ The vaccine efficacy is = 1 − ξ 1 η The immunity development period τ Latent period γ Recovery rate δ Additional death rate due to the infectious disease is a normal cone of C 1 and the interior of C 1 is not empty. Let r > 0 and u = (u 1 , u 2 , u 3 , u 4 , u 5 ) ∶ [−τ, r) → R 5 be a continuous function. For t ≥ 0, define u t ∈ C 1 by u t (θ) = (u 1 (t + θ), u 2 (t + θ), u 3 (t + θ), u 4 (t + θ), u 5 (t + θ)) for all θ ∈ [−τ, 0]. The initial data set for system (2) is in form: The form of φ 3 (0) comes from the implicit solution of E(t) in system (2) which has form: The following theorem shows the nonnegativity and boundedness of model (2). Theorem 2.1. Let φ ∈ X. Then, system (2) has an ultimately bounded unique non-negative solution (S(t), V (t) , E(t), I(t), R(t)) for t ≥ 0 in C 1 . Furthermore, the region Proof of Theorem 2.1 is provided in Section 8.1. In system (2) , when E(t) = I(t) ≡ 0, the disease-free equilibrium is , 0, 0, ηψΛ µ(µ + η)(µ + ψ) always exists. The equations for the diseased classes E and I in the linearized system of (2) about E 0 can be rewritten as d dt where Since the matrix B is non-singular, the next infection operator AB −1 [35] . Hence, reproduction number R E for system (2) is the spectral radius of the matrix AB −1 , which is Since we introduce a vaccination program in model (2), R E is called the effective reproduction number which gives the actual number of secondary infections per infectious person at any time [4, 36] . Biologically, 1 γ+δ+µ is the time spent as an infectious individual and e −µτ is the survival rate of infected individual in latent period. Hence, first second term of R E gives the number of secondary infections of susceptible vaccinated individuals that one infected individual can produce in a disease-free population S 0 V 0 . From the epidemiological point of view, when R E is less than unity, the epidemic cannot develop and eventually the disease dies out. On the other hand, when R E > 1, the population of infected host grows and an outbreak occurs. In this section, we study the threshold dynamics of system (2) with respect to the sign of R E − 1 . Global stability of disease-free equilibrium. As we mentioned in Section 3, If R E > 1, the disease-free equilibrium E 0 is unstable for system (2) , and it is locally asymptotically stable if R E < 1. Since the equations of S, V and I are decoupled in (2), it is enough to study the system: . To obtain the global stability of E 0 in (2), first, we prove the following results, Theorems 4. with The following result shows the global stability of E 0 for system (2) . Proof of Theorem 4.4 is obtained by using the reverse Fatou lemma and the limiting systems and chain transitive sets theory, see Section 8.4. Theorem 4.4. When R E < 1, the disease-free equilibrium E 0 is globally asymptotically stable for system (2) in X . In this section, we prove the persistence of system (2) when Here ∂X 1 is the set of states without disease presence. The following results demonstrate the uniform persistence of the disease state in (2) with respect to X 1 . Furthermore, we can prove the uniform persistence of (2) with respect to X 1 . Theorem 4.6. If R E > 1, then the system (2) is uniformly persistent in (2), i.e., there is a positive number κ 2 > 0 such that every solution in system (2) with φ ∈ X 1 , satisfies Proofs of Theorems 4.5 and 4.6 are given in Sections 8.5 and 8.6, respectively. Let R E > 1, then from Theorem 4.6 and (6), we have . and hence, there is no intersection between Z 1 and Z 2 . When R E > 1, then there exists I * > 0 such that Z 1 (I * ) = Z 2 (I * ), see Figure 2 . We cannot guarantee the uniqueness of I * because Z 1 may intersect Z 2 in two intervals (0.ζ) and (ζ, ∞). Consequently, we have the following result. Theorem 4.7. If R E < 1, then (2) has no endemic equilibrium. If R E > 1, then an endemic equilibrium E 1 = (S * , V * , E * , I * , R * ) exists such that All functions are defined in (8) . Moreover S * ∈ (0, S 0 ) and V * ∈ (0, V 0 ). Regarding the linear stability of E 1 , the analysis seems to have a tedious calculations due to the third order transcendental characteristic equation. However, in the numerical simulations, we discuss the stability in Section 6.4. In this section, we discuss the critical vaccination coverage rate that eliminates the disease. Let R E (ψ) ∶= R E . When the vaccination is absent, i.e. ψ = 0, the effective reproduction number In fact R 0 is so-called the basic reproduction number which is the average number of secondary cases arising from one infectious individual in a totally susceptible population [37, 38] . Notice that, R E (ψ) can be written as Thus, Since we have ξµ µ+ψ R 0 ≤ R E ≤ R 0 , and hence, R 0 < 1 implies R E < 1, but the reverse is not true. Now, we find the critical level of vaccination to eradicate of the disease when R 0 > 1. Assume R 0 > 1, then Obviously * increases as R 0 increases. Hence, when (the vaccine efficacy) is not large enough and R 0 is high, the disease may not be eradicated even if everybody gets the vaccine. That is, R E (ψ) cannot become below 1 when ψ becomes high, see Figure 3 . Biologically, Lemma 5.1 indicates that ψ * is the vaccination coverage rate to eradicate of the disease. From V 0 and Theorem 2.1, we have The equality holds in (10) when δ = 0 in the model (2) . Hence, in a disease-free population, the proportion of vaccinated individuals is at least Thus, when R 0 > 1, the critical proportion of the population that should be vaccinated when the vaccination is imperfect and ψ = ψ * is given by Figures 4a and 4b show the contour plot of ρ ε in (R 0 , )and (τ, )-plane, respectively. We notice that when is fixed, ρ ε increases as R 0 increases. While τ does not have a noticeable effect on the value of ρ ε . Also, the figures are consistent with that fact that ∂ρε ∂ε ≤ 0. In this section, firstly, we do a case study regarding measles disease transmission in the United States from 1963 to 2016. Secondly, we study the local and global sensitivity of R E with respect to the parameters of model (2) . Thirdly, we discuss the stability of the endemic equilibrium. Finally, we investigate the sensitivity system of the model (2) with respect to main parameters. All numerical simulations, except the PRCC calculation, are obtained by using Wolfram Mathematica. The PRCC is calculated via the Matlab implementation 1 of the PRCC algorithm developed in [39] . 6.1. Case study. We use the model (2) to simulate the data of measles cases in the United States after the measles vaccine was licensed in 1963 until 2016. [44] . By estimating the natural mortality rate µ as we find an approximate baseline for µ, see Table 2 . Based on the fact that the duration of infection is much shorter than life expectancy, the authors in [45] estimatedR 0 , the expected number of secondary infections engendered by an infective introduced into a fully susceptible population, bŷ R 0 ≈ β × (Infectious period). By using the values ofR 0 and infectious period in [45, Table 2 ], we estimate the transmission rate β. Moreover, we estimate the recovery rate γ as We provide a range for β and γ in Table 2 corresponding to the minimum and maximum for these values. The data set [46] provides the number of deaths by Measles and number reported annual measles incidence in the United States from 1963 to 2016. By taking the ratio of these numbers at each year, we estimate the disease-induced mortality δ (see [47] ): Number of deaths during a specified period Number of persons at risk of dying during the period . A range of δ is given in Table 2 . Consider the nonlinear function I f (I) in (Q 2 ) with In [4, 33] , n was chosen as 2. Figure 5 shows the plot of I f (I) for n = 2, 3, 4, 5 in (13) . We notice that the value of the highest infection level ζ decreases as n increases. The value of ζ decreases from ζ = 31.6 when n = 2 to ζ = 3 when n = 5. Consequently, for small n, more infection cases need to appear to inhibit the behavior of susceptible and vaccinated individuals. In Figure 6a , we plot the numerical solution of (2) with n = 2 and n = 5 in (13) besides the real data of measles cases in the United States from 1963 to 2016. We also plot in Figure 6b the normalized error where I(t i ) is numerical solution of model (2) andÎ(t i ) is the reported measles data from United States. Note that the quantity E N is the normalized error with respect to value of I(t i ) in (2) . The concept of E N is adopted from [48] . We notice that the numerical solution with n = 2 in (13) gives a better fit. Through this section, we take n = 2 in (13). the proportional change in R E responding to a small variations of a single parameter p at one time. The normalized forward sensitivity index of R E (elasticity of R E ) measures such relative change in R E , denoted by Υ R E p , and defined as [49, 50] : From the explicit formula of R E in (5), we derive an analytical expression for Υ R E p to each parameter described in Table 1 . From (5) and (15), we have First, we notice that the sensitivity indices Υ R E Λ and Υ R E β for the parameters Λ and β are independent of any other parameters, hence, they are locally and globally valid. Also, both parameters are Table 2 to calculate numeric values for Υ R E p , see Figure 7 . For example, when µ increases by 100%, R E decreases by 45% while when ξ increases by 100%, R E decreases by 41%. Thus, due to absolute value of the sensitivity index we have that µ is more important for R E . From Figure 7 , we notice that the order of parameters from the highest importance to the lowest is Λ(and β), ψ, γ, µ, δ, ξ, η and τ . values [51] . We use the Latin Hypercube Sampling (LHS) design and Partial Rank Correlation Coefficient (PRCC) analysis technique to provide good insight on global sensitivity of R E on the uncertainties in its parameters [39] . The PRCC value varies in the interval [−1, 1] such that there is a perfect negative positive correlation when the value is −1 1. In addition, the PRCC value is statistically significant when PRCC value > 0.5. Figure 8 shows PRCC values of R E where the parameters sampling are produced from LHS with uniform distribution over the parameters values in Table 2 with 100, 000 samples. Form Figure 8 , we have that the most influential parameters on R E , ordered from highest to lowest, are β, γ, ξ, δ and Λ. Thus, R E = 3.88 > 1. We find that I * = 12.78 satisfies Z 1 (I) = Z 2 (I) (see Figure 9 -a), hence, E 1 = (1.59, 8.313, 6.99, 12.78, 19.18 ) exists uniquely besides E 0 = (5.26, 31.58, 0, 0, 63.16) . The stabilities of of E 0 and E 1 are determined by the sign of the real part of the roots of respectively. The distribution of ∆ 0 and ∆ 1 roots is shown in Figure 9 -b. It clear that E 0 is unstable and E 1 is locally stable, see Figure 9 -c. Figure 9 -d shows a phase diagram with different initial points. We can notice that the unique endemic equilibrium E 1 is globally asymptotically stable when R E > 1. 6.5. Sensitivity of model solutions. The sensitivity system of the model (2) with respect to a parameter p ∈ {ψ, µ, ξ, η, τ, γ, δ} is given by the partial derivative of X = (S, V, E, I, R) T with respect to p, denoted by X p = ∂X ∂p . Define f (S, V, E, I, R) ∶= dX dt , then by the Chain Rule and Clairaut's Theorem, we have The semi-relative sensitivity for X is represented by pX p , while the logarithmic sensitivity describes the percentage changes in the solutions induced by positive perturbations of the parameters p Xp X . The details method is given in [52] . Figures 10-11 show the semi-relative and logarithmic sensitivity curves for S, V and I with respect to the parameters β, ψ, ξ and τ . From the semi-relative sensitivity solution curves in Figure 10 , we can interpret that the perturbation of β or ψ has a big influence over S, V and I. For all t ∈ (0, 50) the effects of the parameters ψ and ξ are the same, while they are opposite to the effects of β and τ on the variables V and I. However, on the variable S, τ has an opposite effect to the parameter p ∈ {β, ψ, ξ}. More precisely, the effect of β is positive on I and negative on S and V . Moreover, at t ≈ 15, we notice that all parameters have the largest effects on V and solution curves of I start to stabilize at an equilibrium level. From the logarithmic sensitivity solution curves in Figure 11 , we observe that a perturbation of τ has the biggest positive impact upon solutions of S, V and I. For instance, At time t = 12, τ causes a roughly 42% change in the solution of V . On the other hand, β and ψ have negative influences on S and opposite influences over {V, I}. β has negative impact upon solution of V and positive upon I. The affects of ξ perturbation is not noticeable on S but small negative (positive) impact on V (I). Finally, since the function f (I) is independent of the value of R E , we consider f (I) = 1 + dI q and study the influence of the parameters d and q on the solution of the model (2). In Figure 12 , the semi-relative sensitivity solution curves show that d and q have a positive influences on S and V and a negative affects on I. The influence of q is relatively higher than that of d. On the other hand, the logarithmic sensitivity solution curves in Figure 11 show that q has a positive influence on {S, V, I}. There are small affects (around zero) of d on the three variables. Nowadays, epidemiological modeling plays a key role in studying communicable disease transmission and providing strategies for prevention and control of them. To date, vaccination is considered to be one of the most popular and effective methods of mitigation and prevention of epidemics. In the paper, we have studied the transmission dynamics of an infectious disease in the presence of imperfect vaccine by a stage-structured mathematical model. To interpret the "psychological effect" when an infectious disease is spread in a population, we have considered a general nonmonotone nonlinear incidence rate function. We have proved that the model is biologically well-posed by showing that the solutions of the proposed model exist (uniquely determined) and they are nonnegative and bounded. Since we have introduce a vaccination program, we have used the method in [35] to obtain an explicit expression for the effective reproduction number (R E ) which gives the actual number of secondary infections per infectious person at any time [4, 36] . Then we have obtained the threshold dynamics of the system with respect to R E : (i) we have shown that the disease-free equilibrium is globally stable when R E < 1; and (ii) we have discussed the system persistence and the coexistence of endemic equilibrium when R E > 1. Then, we have discussed the critical vaccination coverage rate that is required to eliminate the disease and the critical proportion of the population (ρ ) that should be vaccinated when the vaccination is imperfect. We have not notice any influence of the latent time τ on ρ when the vaccine efficacy is fixed. However, we have noticed that ρ increases as R 0 increases when is fixed. R 0 is the value of R E when the vaccination rate ψ is zero. Furthermore, through the theoretical analysis, we have found that when is not large enough and R 0 is high, the disease may not be eradicated even if everybody gets the vaccine. That is, R E cannot become below the unity even when ψ becomes high. As a case study, we have studied the measles disease transmission in the United States from 1963 to 2016. Our numerical simulation has indicated a significant decrease in the case number through this period which consistent with the fact that the risk of the measles virus transmission is very low in the United States since the measles vaccine was licensed in 1963. Moreover, through the numerical simulations of the model, we have following: • We have carried out the global and local sensitivities analysis for R E . We have found that the latent time τ does not have a noticeable affect on R E . On the other hand, we have seen and the transmission rate β has a high influence on R E . Regarding to the vaccination parameters, ψ and the reduction coefficient (ξ), they both have an opposite effect on the value of R E , ψ has a negative influence while ξ has a positive one; • When there is a unique endemic equilibrium E * , the system has shown local and global stability behaviour about this point; • The semi-relative sensitivity curves have shown that the perturbation of β and ψ has a big influence over the variables {S, V.I}. In addition, the logarithmic sensitivity curves have indicated that the latent period τ has exhibited the biggest positive impact upon solutions of {S, V.I}. This work suggests that τ affects the solution behaviour although it does not affect the value of R E ; • Although R E is independent of the function f , the sensitivity of model solutions has shown the influence of f on the solutions behaviour. For further work, it will be interesting to consider the "asymptomatic carriers". These carriers are individuals who have been infected and are able to transmit their illness without showing any symptoms [19] . For certain infectious diseases, asymptomatic carriers are a potential source for transmission such as Typhoid Fever, HIV and, most recently, the 2019-nCoV. In this section, we provide a proof for each result in the manuscript. f (φ 4 (−τ )) e −µτ − (µ + γ + δ)φ 4 (0), is continuous and Lipschitz in φ in each compact set in R × X because X is closed in C and for any φ ∈ X. Thus, there is a unique solution u(t, φ) of system (2) through (0, φ) on its maximal interval [0, r) of existence [53, Theorem 2.2.3] . Let φ ∈ X, if φ 4 (0) = 0, then F 4 (φ) ≥ 0. Consequently, F i (φ) ≥ 0 when φ i (0) = 0 for i = 1, 2, 5. Hence, it follows from [54, Theorem 5.2.1] that for i = 1, 2, 4, 5, the solutions S(t), V (t), I(t) and R(t) are non-negative for all t ∈ [0, r). Consequently, from (3) we obtain E(t) ≥ 0. In the later equation (replace ≤ by =), N (t) = Λ µ is globally asymptotically stable. Hence, by the comparison arguments [55, Lemma 1.2], we have (7), we know that S 0 is a unique solution and globally attractive. Since lim t→∞ S(t) = S 0 in (7), we have the following limiting system: Now, we lift the dynamics to the main system (7) from the limiting system (17) . Define the omega limit set of (S(0), V (0)) for the solution semi-flow Ψ t of (7) byω =ω(S(0), V (0)). Thus,ω is an internally chain transitive set for Ψ t [35, Lemma 1.2.1]. Hence, ω = {S 0 } ×ω 1 for someω 1 ⊂ R. Since Ψ t (ω) =ω, for all t ≥ 0, we have whereΨ t is the solution semi-flow associated with (17) . Thus,Ψ t (ω 1 ) =ω 1 . By the asymptotic global stability of (6) as t → ∞ when R E < 1 by using the method of fluctuations [56, 57] . Let Hence, when n ≥ 1 and n → ∞, it follows from (6) that From the first and second inequalities in (18), we have S ∞ ≤ S 0 and V ∞ ≤ V 0 . Otherwise, we have contradictions (0 < 0). Thus, it follows from the third inequality in (18) and Since R E < 1, we have R E − 1 < 0, and hence, I ∞ = 0 because I(t) ≥ 0. This proves Claim 1. It follows from Claim 1 that lim t→∞ I(t) = 0 when R E < 1, and hence, the system (6) is asymptotic to the limiting system (7) . Recall that (S 0 , V 0 ) is globally attractive in the limiting system (7), see Theorem 4.2. To lift the dynamics of the limiting system (7) to the main system (6), we use the theory of internally chain transitive sets to prove lim t→∞ (S(t), V (t), I(t)) = S 0 , V 0 , 0 . and ω = ω(φ 1 , φ 2 , φ 4 ) be the omega limit set for the solution semi-flow z t (φ 1 , φ 2 , φ 4 ) of (6). Hence, ω is an internally chain transitive set for z t , see e.g., [35, whereẑ t is the solution semi-flow associated with the equation Notice thatω becomes an internally chain transitive set forẑ t (ẑ t (ω) =ω) because ω is an internally chain transitive set for z t . When R E < 1, then lim t→∞ I(t) = 0 in (19) . Suppose the solutions of (19) take the form I(t) = ce λt where λ satisfies the characteristic equation Assume there exists a zero in (20) with Re(λ) then which is a contradiction because λ µ + γ + δ + 1 > 1 and R E e −λτ < 1 when R E < 1. Thus, all roots have negative real part. Therefore, lim t→∞ I(t) = 0. This proves Claim 2. Let W s (0) be the stable manifold of 0, then it follows from Claim 2 thatω ∩ W s (0) ≠ ∅. Hence by [35, Theorem 1.2.1] we haveω = {0}. Therefore, we have ω = (S 0 , V 0 , 0), and hence lim t→∞ (S(t), V (t), I(t)) = S 0 , V 0 , 0 , i.e., (S 0 , V 0 , 0) is globally attractive in system (6). (3) and the reverse Fatou lemma (see e.g., [58] ) that Thus, lim t→∞ E(t) = 0. Hence, we obtain the limiting system of (2): It is clear that R 0 is globally attractive in (22) . By using the theory of internally chain transitive sets, we lift the dynamics of the limiting system (22) to the main system (2), see the proofs of Since E 0 is the local stability when R E < 1, see Theorem 4.1, we obtain that E 0 is globally asymptotically stable in (2) when R E < 1. 8.5. Proof of Theorem 4.5. Fix a small 0 < σ ≪ 1. Since S(t)+ξV (t) f (I(t)) → S 0 + ξV 0 as (S, V, I) → (S 0 , V 0 , 0), in a neighborhood of (S 0 , V 0 , 0), we have S 0 + ξV 0 − σ < S(t) + ξV (t) f (I(t)) < S 0 + ξV 0 + σ. Claim 3. There exists an (σ) ∶= , such that for any φ ∈ X 1 lim sup By contradiction, suppose that u t (ψ) − E 0 < for some ψ ∈ X 1 . Thus, there exists t 0 > 0 such that S(t) − S 0 < , V (t) − V 0 < and I(t) < for t > t 0 + τ . Hence, (23) is satisfied. From the fourth equation of (2), we have dI dt > β(S 0 + ξV 0 − σ)e −µτ I(t − τ ) − (µ + γ + δ)I(t). For sufficiently small σ, the equation obtained from (24), by replacing > with =, is quasimonotone. Hence, it suffices to study the real roots of the characteristic equation ([54, Theorem 5.5.1]) ∆ 1 (λ) ∶= λ + (µ + γ + δ) − β(S 0 + ξV 0 − σ)e −µτ e −λτ = 0. Let σ = 0. Then, ∆ 1 (0) = (µ + γ + δ)(1 − R E ) < 0 when R E > 1. Notice that ∆ 1 is continuous, increases for λ > 0 and goes to ∞ when λ → ∞. Hence, there exists a positive rootλ > 0 satisfying (25) . Let λ 0 (σ) be the principle eigenvalue. Then, λ 0 (0) > 0, and hence, due the continuity of λ 0 , λ 0 (σ) > 0 for sufficiently small σ > 0. Thus, there exists a solution U (t) = ce λ 0 (σ)t with c > 0. Since I(t) ≥ 0 for t > 0, by the comparison theorem [54, Theorem 5.1.1], there exists a small K > 0 such that I(t) ≥ Kce λ 0 (σ)t for all t ≥ t 0 + m. Thus, I(t) → ∞ as t → ∞ due to the positivity of λ 0 (σ). which is a contradiction to the boundedness of (2). This proves Claim 3. Let ω 1 (φ) be the omega limit set of the orbit u t (φ) through φ ∈ X and define M ∂ = {φ ∈ X ∶ u t (φ) ∈ ∂X, t ≥ 0}. Let φ ∈ X and define a continuous function p ∶ X → R + by p(φ) = φ 4 (0). Therefore, p −1 (0, ∞) ⊂ X 1 and p(u t (φ)) > 0 for all t > 0 whenever p(φ) > 0. It follows from Claim 4 that any forward orbit of u t in M ∂ converges to E 0 . Let W s (E 0 ) is the stable manifold of E 0 . Then it follows form that Claim 3 that W s (E 0 ) ∩ X = ∅, and hence, there is no cycle in M ∂ from E 0 to E 0 . By [59] , there exists κ 1 > 0 such that lim inf t→∞ I(t) ≥ κ 1 for all φ ∈ X 1 , which implies the uniform persistence. 8.6. Proof of Theorem 4.6. From Theorems 2.1 and 4.5, we have κ 1 ≤ I(t) ≤ Λ µ. Consequently, from the first equation of (2) and (Q 1 ), we have When we replace ≥ by = in (26),κ is globally asymptotically stable in (26) . Hence, by applying the comparison arguments [55, Lemma 1.2], we have S(t) ≥κ 1 . Parallelly, from the equations of V and R in (2) we have V (t) ≥κ 2 and R(t) ≥κ 3 , respectively, whereκ 2 = ψκ 1 f (Λ µ) βξΛ µ + (η + µ) f (Λ µ) andκ 3 = ηκ 2 + γκ 1 µ . It follows from Theorem 4.5 and the integral form of E(t) in (3) Choose κ 2 = min{κ 1 ,κ 1 ,κ 2 ,κ 3 ,κ 4 }. This completes the proof. On vaccine efficacy and reproduction numbers Dynamics of vaccination in a time-delayed epidemic model with awareness Dynamics of cholera epidemics with impulsive vaccination and disinfection Complicated endemics of an SIRS model with a generalized incidence under preventive vaccination and treatment controls SVIR epidemic models with vaccination strategies An SVEIR model for assessing potential impact of an imperfect anti-sars vaccine An epidemiology model that includes a leaky vaccine with a general waning function The asymptotic behavior of a stochastic vaccination model with backward bifurcation A comprehensive probabilistic solution of random SIS-type epidemiological models using the random variable transformation technique Probabilistic solution of random SI-type epidemiological models using the random variable transformation technique The stochastic SI model with recruitment and deaths I. comparison with the closed SIS model A primer on stochastic epidemic models: Formulation, numerical simulation, and analysis The threshold of a stochastic SIS epidemic model with vaccination Influence of nonlinear incidence rates upon the behavior of SIRS epidemiological models Global analysis of an epidemic model with nonmonotone incidence rate Dynamical behavior of an epidemic model with a nonlinear incidence rate Bifurcation analysis of an SIRS epidemic model with a generalized nonmonotone and saturated incidence rate Stability analysis for a vector disease model A periodic disease transmission model with asymptomatic carriage and latency periods Modeling diseases with latency and relapse Analysis of an SEIRS epidemic model with two delays Threshold dynamics in an SEIRS model with latency and temporary immunity A time-delayed epidemic model for Ebola disease transmission Global stability of an SVEIR epidemic model with ages of vaccination and latency Global stability of a delayed epidemic model with latent period and vaccination strategy A delay SIR epidemic model with pulse vaccination and incubation times Pulse vaccination of an SEIR epidemic model with time delay Modeling the initial transmission dynamics of influenza a H1N1 in guangdong province, china Infectious AIDS: Have we been misled? A note on generation times in epidemic models Infectious disease epidemiology: theory and practice Backward bifurcation of an epidemic model with treatment Mathematics for life science and medicine Basic reproduction ratios for periodic compartmental models with time delay Estimation of effective reproduction numbers for infectious diseases using serological survey data The basic reproduction number (R 0 ) of measles: a systematic review Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission A methodology for performing global uncertainty and sensitivity analysis in systems biology Percentage of children aged 19 − 35 months who are vaccinated against measles, mumps and rubella in the U.S. from 1994 to 2017 Questions about measles Nursing: 2nd Year E-Book Measles virus, immune control, and persistence Life expectancy at birth Plug-and-play inference for disease dynamics: measles in large and small populations as a case study Epidemiology and Prevention of Vaccine-Preventable Diseases A dictionary of epidemiology A phase model with large time delayed coupling Determining important parameters in the spread of malaria through the sensitivity analysis of a mathematical model Reproduction numbers of infectious disease models Mathematical modeling in systems biology: an introduction Sensitivity analysis of a nonlinear lumped parameter model of HIV infection dynamics Introduction to functional differential equations Monotone Dynamical Systems: An Introduction to the Theory of Competitive and Cooperative Systems Ordinary differential equations and dynamical systems Global dynamics in sea lice model with stage structure Differential equation models of some parasitic infections: methods for the study of asymptotic behavior Mathematical Analysis: An Introduction to Functions of Several Variables Robust persistence for semidynamical systems The author would like to thank the referees for their careful reading and helpful suggestions.