key: cord-0934403-hk8ziliv
authors: Al-Darabsah, Isam
title: Threshold dynamics of a time-delayed epidemic model for continuous imperfect-vaccine with a generalized nonmonotone incidence rate
date: 2020-07-27
journal: Nonlinear Dyn
DOI: 10.1007/s11071-020-05825-x
sha: 7da15ad52b22cada80795e400d19086dd4e77657
doc_id: 934403
cord_uid: hk8ziliv

In this paper, we study the dynamics of an infectious disease in the presence of a continuous-imperfect vaccine and latent period. We consider a general incidence rate function with a non-monotonicity property to interpret the psychological effect in the susceptible population. After we propose the model, we provide the well-posedness property and calculate the effective reproduction number [Formula: see text] . Then, we obtain the threshold dynamics of the system with respect to [Formula: see text] by studying the global stability of the disease-free equilibrium when [Formula: see text] and the system persistence when [Formula: see text] . For the endemic equilibrium, we use the semi-discretization method to analyze its linear stability. Then, we discuss the critical vaccination coverage rate that is required to eliminate the disease. Numerical simulations are provided to implement a case study regarding data of influenza patients, study the local and global sensitivity of [Formula: see text] , construct approximate stability charts for the endemic equilibrium over different parameter spaces and explore the sensitivity of the proposed model solutions.

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 [1] . Since then, many authors have incorporated time delays in epidemic models in different scenarios, such as vaccination period [2] , asymptomatic carriage period [3] , immune period [4] and incubation period or latent period [3] [4] [5] [6] [7] . More precisely, in [3] , a disease transmission model with two delays in incubation and asymptomatic carriage periods is investigated. In [4] , the authors study an SEIRS epidemic model with constant latent and immune periods. In [5] , a latent period and relapse are considered in a general mathematical model for disease transmission. In [7] , the authors studied a time-delayed SIR model with nonlinear incidence rate and Holling functional type II treatment rate for epidemic transmission. Also, many authors studied time-delayed epidemic models with vaccination [8] [9] [10] [11] . For example, the authors in [9] study a vaccination model with a time delay to represent the time that an unaware susceptible individual takes to become aware of the infection. Due to the inherent complexity of epidemiological transmission, other works studied epidemic complex network models. For example, in [12] , the authors studied a semi-random epidemic network and discussed the relationship between its topological structure and the optimal control of the epidemic. In [13] , the authors used the concept of epidemiology to analyze data from real computer virus epidemics by using complex network models. They studied an SIS model on scale-free graphs by largescale simulations and analytical methods.

Vaccines are considered to be one of the most significant medical means of disease control and prevention [14] . They have played a major role in the spread and eradication of many infectious diseases, such as smallpox, or partially, like measles. Many authors in the literature have studied the dynamics of epidemics models with different types of vaccination schedules [15] [16] [17] [18] [19] . For instance, in [16] , an SIR model with a generalized incidence under preventive vaccination and treatment controls is proposed. In [20] , the authors developed an SIVS epidemic model with degree-related transmission rates and imperfect vaccination on scale-free networks.

In [17] , the authors establish two SVIR models: one with continuous vaccination strategy and another one with pulse vaccination strategy. An epidemic model to study the potential impact of a SARS vaccine when it is imperfect is proposed in [18] . The dynamics of cholera epidemics with impulsive vaccination is studied in [15] .

To incorporate the effect of behavioral changes on the disease spreading dynamics, the authors in [21] introduce a nonlinear incidence rate of the form Sg(I ) = β S I q 1 1 + d I q 2 (1)

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+d I q 2 ) represents the inhibition effect from the behavioral change of the susceptible individuals when the number of infectious individuals increase [22] . It also can be used to describe the crowding effect of infectious individuals [23] . There are three types of incidence functions g(I ) = β I q 1 /(1 + d I q 2 ) based on the values of q 1 and q 2 [22, 24] : (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 [22] consider a nonmonotone incidence rate to represent the psychological effect with q 1 = 1 and q 2 = 2.

In [23] , a saturated incidence rate is considered with q 1 = q 2 = 2. For more details and examples, we refer the reader to [24] and references within. In this paper, we consider a general form of an nonmonotone 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 toward reducing the number of contacts among individuals per unit time [22] . In the real world, the latent period may vary from days, as in influenza A H1N1 [25] , to years, as in AIDS [26] . 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 [27] . Thus, the latent period has an influence on the epidemic growth [28] . The purpose of the work is to investigate the dynamical behavior of a time-delayed SEIR model with continuous imperfect-vaccine and discuss the effect of the latent period on the epidemic. The paper is organized as follows. In the next section, we present the mathematical model and the study its well-posedness. Then, in Sect. 3, we calculate the effective reproduction number R E . In Sect. 4, we use various methods, such as the Lyapunov functional technique and the method of fluctuations, to establish the global stability of the diseasefree equilibrium when R E < 1. Then, in Sect. 5, we study the system persistence when R E > 1. Moreover, we use the semi-discretization method of order one to study the local stability of endemic equilibrium. In Sect. 6, we discuss the critical vaccination coverage that is required to eliminate the disease. In Sect. 7, we consider an application to influenza transmission. Numerically, we also study the local and global sensitivities of R E , discuss the stability of the endemic equilibrium and examine the sensitivity of model solutions. Finally, we discuss our results in Sect. 8.

Motivated by [29] , 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),

is the infected population and R(t) is the recovered population. We assume there is a recruitment rate into the population and the natural death rate is the same for all the compartments. We make the following assumptions to describe the interaction among the five classes:

• When the susceptible individuals receive a vaccine, they move from class S to V . On the other side, 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., [30] ); • The vaccine is not perfect, in the sense that, the individuals in V are not on a fully protective level. Consequently, when the vaccinated individuals become infected, they move into E for τ period with a reduced transmission rate ξβ, where ξ ∈ [0, 1] is the reduction coefficient [30] . When individuals gain immunity, they move into the class R [17] ; • After the latent period τ , infected individuals become infectious and move from E to I . Then when they recover, they move to R. Individuals in R may lose immunity at rate α and become susceptible again, that is, they move to the class S (i.e., the vaccine is continuous).

These 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 e.g., [3] . The population flux among the five compartments is given in Fig. 1 and a description of the parameters is given in Table 1 . Denote the infection force function by β I / f (I ). Following [16, 31] , we assume that 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 policies. Consequently, we consider the following assumptions:

There is a ζ > 0 such that d dI

Here, ζ is the critical level of invectives, that is, the incidence rate is increasing when I ∈ (0, ζ ] and it is decreasing when I > ζ. Notice that the nonlinear incidence rate has the form β(S+ξ V )I f (I ) .

Consider the Banach space C := C([−τ, 0), R 5 ) with the maximum norm

is a normal cone of C and the interior of C is not empty. Let r > 0 and u = (

The initial data set for system (2) is in form:

where the form for φ 3 (0) follows from the implicit solution of E(t) in system (2) which has form:

The following theorem shows the nonnegativity and boundedness of system (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. Furthermore, the region = (S(t), V (t), E(t), I (t), R(t)) ∈ R 5 + :

is a positive invariant set and globally attractive set for (2) .

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 

Thus, N (t) = μ is globally asymptotically stable on (4 The general solution of (4) can be written as

Therefore, when

Consequently, the set is the globally attractive set for (2).

In system (2), when

,

,

.

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

and I (t) at t = 0, then from (6) the distribution of the remaining population of classes E(t) and

The total number of newly infected individuals is

due to the nonsingularity of the matrix M 2 . Then it follows that, the next infection operator is

In the literature (see e.g., [35] ), the reproduction number R E for system (2) is the spectral radius of the matrix M 0 , which is

Since we introduce a vaccination program in system (2) , R E is called the effective reproduction number which gives the actual number of secondary infections per infectious person at any time [16, 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 susceptibleaccinated individuals that one infected individual can produce in a disease-free population S 0 \V 0 .

When R E is less than unity, the epidemiological interpretation is that an 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 establish the global stability of E 0 when R E < 1 .

As we mentioned in Sect. 3, E 0 exists for all parameters values. The following result indicates the instability and local stability of E 0 in (2) . See [35, Theorem 2.1 and Corollary 2.1].

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 , I and R are decoupled in (2) , it suffices to study the following system:

. To obtain the global stability of E 0 in (2), first, we prove the following results, Theorems 4.2 and 4.3, by using the Lyapunov functional technique, the method of fluctuations and the theory of limiting systems and chain transitive sets.

It follows from (9) that the largest invariant set in the set of dU dt = 0 is (S 0 , V 0 , R 0 ). By the LaSalle's invariance principle [37] ,

Now, we show that the limit supremum of I (t) is zero in (8) as t → ∞ when R E < 1 by using the method of fluctuations [38, 39] . Let

Claim 1 When R E < 1, then I ∞ = 0 in (8).

Hence, when n ≥ 1 and n → ∞, it follows from (8) that

Since

, and hence, 

when R E < 1, and hence, the system (8) is asymptotic to the limiting system (9) .

Recall that (S 0 , V 0 , R 0 ) is globally attractive in the limiting system (9), see Theorem 4.2. To lift the dynamics of the limiting system (9) to the main system (8), we use the theory of internally chain transitive sets to prove

and ω = ω(φ 1 , φ 2 , φ 4 , φ 5 ) be the omega limit set for the solution semi-flow z t (φ 1 , φ 2 , φ 4 , φ 5 ) of (8). Hence, ω is an internally chain transitive set for z t , see e.g., [35, 

Notice thatω becomes an internally chain transitive set forẑ t (ẑ t (ω) =ω) because ω is an internally chain transitive set for z t .

Claim 2 When R E < 1, then lim t→∞ I (t) = 0 in (13) .

Suppose the solutions of (13) take the form I (t) = ce λt where λ satisfies the characteristic equation

Assume there exists a zero in (14) with Re(λ) then

which is a contradiction because λ μ + γ + δ + 1 > 1 and |R E | e −λτ < 1

when 

i.e., S 0 , V 0 , 0, R 0 is globally attractive in system (8) .

The following result shows the global stability of E 0 for system (2) .

When R E < 1, the disease-free equilibrium E 0 is globally asymptotically stable for system (2) in X .

Proof It follows from Theorem 4.3, the integral form of E(t) in (3) and the reverse Fatou lemma (see e.g., [40] 

Thus, lim

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.

In this section, we prove the persistence of system (2) when R E > 1. Define

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 . 

Proof Fix a small 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 (16) is satisfied. From the fourth equation of (2), we have

For sufficiently small σ , the equation obtained from (17) , by replacing > with =, is quasimonotone. Hence, it suffices to study the real roots of the characteristic equation ([33, Theorem 5.5.1])

Let σ = 0. Then, 1 

Notice that 1 is continuous, increases for λ > 0 and goes to ∞ when λ → ∞. Hence, there exists a positive rootλ > 0 satisfying (18) . 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 [33, Theorem 5.1.1], there exists a small K > 0 such that I (t) ≥ K ce λ 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 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 from that Claim 3 that W s (E 0 ) ∩ X = ∅, and hence, there is no cycle in M ∂ from E 0 to E 0 . By [41] , there exists κ 1 > 0 such that lim inf t→∞ I (t) ≥ κ 1 for all φ ∈ X 1 , which implies the uniform persistence.

Furthermore, we can prove the uniform persistence of system (2) with respect to X 1 . (2), i.e., there is a positive number κ 2 > 0 such that every solution in system (2) 

Proof From Theorems 2.1 and 5.1, we have κ 1 ≤ I (t) ≤ /μ. Consequently, from the first equation of (2) and (Q 1 ), we have

When we replace ≥ by = in (19),

is globally asymptotically stable in (19) . Hence, by applying the comparison arguments [34, Lemma 1.2], we have S(t) ≥κ 1 . Parallely, from the equations of V and R in (2) we have V (t) ≥κ 2 and R(t) ≥κ 3 , respectively, wherê

and

It follows from Theorem 5.1 and the integral form of

Choose κ 2 = min{κ 1 ,κ 1 ,κ 2 ,κ 3 ,κ 4 }. This completes the proof.

Since C([−τ, 0], R 4 + ) is a convex set and the system (8) is ultimately bounded and uniformly persistent with respect to X 1 when R E > 1, it follows from [42, Theorem 3.1] that (8) has at least a positive equilibrium point (S * , V * , I * , R * ) when R E > 1. Consequently, the system (2) has the positive equilibrium point E 1 = (S * , V * , E * , I * , R * ). The value of E * can be found from the integral form of E(t) in (3). It is not possible to derive an explicit formula for the components of E 1 or guarantee its uniqueness due to the presence of the exponential terms in the model.

Regarding the stability of E 1 , a self-contained proof seems to have a tedious calculation due to the fourthorder transcendental characteristic equation. However, we use the semi-discretization method of order one to study the linear stability of the endemic equilibrium [43, 44] .

Let R E > 1 and assume that E 1 exists. By setting x = (S, V, E, I, R) − E 1 , the linearized system of (2) about E 1 is

where

Now, define the solution operator U : C → C of (21) by

When all of the nonzero elements of the spectrum of the monodromy operator U (the Floquet multipliers of system (21) ) are within the unit circle of the complex plane, the zero solution of (21) is stable. While when one or more of the Floquet multipliers are on the unit circle and the rest of them are inside the unit circle, the zero solution may undergo a bifurcation [45] .

To study the location of the Floquet multipliers, we use the semi-discretization method which is an efficient numerical method based on a special kind of discretization technique with respect to the past effects only [46] . By employing this method, we define a Floquet transition matrix, which is an approximation to the infinite-dimensional monodromy operator U corresponding to the linear delayed system (21) . Usually, the semi-discretization method is used to study the linear stability when the system is non-autonomous (contains time-dependent periodic delays or coefficients functions). However, since system (21) is autonomous, we can choose an arbitrary period [43] . Consequently, we assume a period T for the system (21) , and hence, the length of the discretization interval is h = T /K where K is the number of subintervals of [0, T ]. Let t i = ih. Then, in each discretization interval [t i , t i+1 ], the firstorder semi-discretization approximate the delayed term x(t − τ ) by the Lagrange first-order polynomial

where

and r = int(τ/ h + 1/2) with int denoting the integerpart function, see [43, 44] . The scheme of the approximation in (23) is shown in Fig. 2 , and more details are provided in [44, Chapter 3 ].

Consequently, system (21) can be approximated by a system of ordinary differential equations

where i = 0, . . . , K − 1. By using the variation of constants formula, the general solution of (24) can be written as

Using the notation x(t i ) = x i , when t = t i+1 . Then, the solution over one discrete step can be formulated as where P = e Ah and

If A −1 exists, then (25) can be written as

Now, we define the augmented state vector as

Combining z i and (24) leads to the discrete map

Utilizing that T = K h and applying (26) K times with initial state z 0 gives the monodromy mapping

which represents a finite-dimensional approximation of the monodromy operator U associated with E 1 of (2). The rate of convergence for the first-order semi-

When all the eigenvalues of are inside the unit circle of the complex plane, then E 1 is asymptotically stable. While when one or more of the eigenvalues are on the unit circle and the rest of them are inside the unit circle, E 1 may undergo a bifurcation [45] . In the numerical simulations, Sect. 7.4, we implement the above algorithm to construct an approximate stability region for E 1 . We notice that when E 1 is unstable, the model exhibits periodic oscillations as it is expected from SIRS models with delay [47] .

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 becomes

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 [48, 49] . In the case of R 0 , everyone is susceptible while in R E not all contacts will become infected due immunity, hence, R E is less than R 0 from epidemiological point of view. Notice that, R E (ψ) can be written as

Thus,

we have that ξ(α+μ) α+η+μ 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 when 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 Fig. 3 . 

Biologically, Lemma 1 indicates that ψ * is the vaccination coverage rate to eradicate of the disease.

From V 0 and Theorem 2.1, we have

Hence,

Hence, in a disease-free population, the proportion of vaccinated individuals is

Therefore, when R 0 > 1, the critical proportion of the population that should be vaccinated when the vaccination is imperfect and ψ = ψ * is given by Figure 4a and b shows the contour plot of ρ ε in R 0and τ -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 fit the model with data of influenza patients as a case study. Secondly, we study the local and global sensitivity of R E with respect to the parameters of system (2) . Thirdly, we discuss the stability of endemic equilibrium. Finally, we investigate the sensitivity system of the system (2) with respect to main parameters. Through this section, we take f (I ) = 1 + 0.001I 5 .

We use the system (2) to simulate the data of influenza patients (weekly percentage) in North Carolina from January to April, 2011 [50] . A range for (2) parameters is given in Table 2 . Figure 5 shows that the numerical solution of (2) provides a good agreement with the real data.

The local sensitivity analysis of R E provides insight into the proportional change in R E responding to a small variation 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 [56, 57] :

From the explicit formula of R E in (7), we derive an analytical expression for ϒ R E p to each parameter described in Table 1 . From (7) and (33), we have

,

First, we notice for the parameters and β, the sensitivity indices ϒ R E and ϒ R E β are independent of any other parameters; hence, they are locally and globally valid. Also, both parameters are equally important for

Consequently, when or β increases by 100%, then R E increases by 100%. For the other parameters p ∈ {ψ, μ, ξ, η, τ, γ , δ}, they have different impacts on R E due to the different absolute value of the forward sensitivity indices ϒ R E p in (34) . We use the values in Table 2 to calculate numeric values for ϒ R E p , see Fig. 6 . For example, when ψ increases by 100%, R E decreases by 86% while when α increases by 100%, R E increases by 5%; hence, due to absolute value of the sensitivity index we have φ is more important for R E . From Fig. 6 , we notice that the order of parameters from the highest importance to the lowest is μ, (and β), ψ, γ , τ , α, ξ and η.

When there are large perturbations in all parameters, global sensitivity analysis is typically used which includes sampling a given range of parameter values [58] . 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 [59] . The PRCC values vary in the interval [−1, 1] such that there is a perfect negative\positive correlation when the value is −1\1, also, the PRCC value is statistically significant when |PRCC value| > 0.5. Figure 7 shows PRCC values of R E where the parameters sampling are produced from LHS with uniform distribution over the parameter values in Table 2 with 1,000,000 samples. We notice from Fig. 7 that the most influential parameters on R E , ordered from highest to lowest, are γ , ξ , β, η, τ , ψ and α.

In this section, we fix = 300, η = 1, μ = 0.13, ξ = α = 0.01.

We use the package SemiDiscretizationMethod.jl on Julia to implement the algorithm in Sect. 5.1 and con- Fig. 7 PRCC results, the orange boxes represent the partial rank correlation coefficients of R E . The small figures in the right side and bottom show the PRCC scatter plots struct an approximate stability region for E 1 in two parameters space. Figure 8 shows the stability charts in βγ -plane, βψ-plane and γ ψ-plane. In Fig. 8b , when ψ is fixed in the interval [0.4, 0.6] and β increases, a stability switches occur where endemic equilibrium losses its stability and then becomes stable for larger β, see Fig. 9a . While when β is fixed and ψ increases, a Hopf bifurcation occurs and the endemic equilibrium losses its stability and becomes unstable, see Fig. 9b . The latter case also occurs when β or ψ is fixed and γ increases, see Fig. 8a and c. Figure 9c shows a one-parameter bifurcation diagram, by varying the value of τ . For small τ , the endemic equilibrium E 1 is stable (Fig. 10a) . By increasing the value of τ , E 1 loses the stability and a Hopf Bifurcation occurs around τ ≈ 2.94. For τ > 2.94, a unique stable periodic solution of system (2) exists (Fig. 10c) . In Fig. 10 , we plot the phase portrait of the system (2) with different initial condition and various value of τ , we notice that system (2) exhibits global asymptotic stability behavior when R E > 1.

The sensitivity system of (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 is represented by p X p X . The detailed method is given in [60] . Figures 11 and 12 show the semi-relative and logarithmic sensitivity curves for S, V, I and R with respect to the parameter β, ψ, α and τ . From Figs. 11 and 12, we can interpret that the perturbation of τ has a big influence over S, V , I and R. For t ∈ (0, 10), a remarkable positive affect of the parameter β occurs on I and R, while an opposite affect appears on the variables S and a (912.618, 161. Fig. 11 The semi-relative sensitivity curves of S, V, I and R with respect to the parameter β, ψ, α and τ Fig. 12 The logarithmic sensitivity curves of S, V, I and R with respect to the parameter β, ψ, α and τ Fig. 13 The semi-relative (row 1) and logarithmic (row 2) sensitivity curves for S, V, I and R with respect to the parameter d (red) and q (green) in f (I ) = 1 + d I q V . Moreover, at t ≈ 4(7), we notice that both parameters τ and β have the largest effects on S and V (I and R). The perturbation of ψ and α has a noticeable positive affect on V and R.

Finally, since the function f (I ) is independent of the value of R E , we consider f (I ) = 1 + d I q and study the influence of the parameters d and q on the solution of the system (2). In Fig. 13 , the 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. More precisely, q has an oscillatory influence on {S, V, I, R}. There are small affects (around zero) of d on the four variables.

Nowadays, epidemiological modeling plays a key role in providing strategies for the prevention and control of many communicable diseases. Vaccination, meanwhile, is considered to be one of the most favored and effective methods of mitigation and elimination of epidemics. In the paper, we have studied infectious disease transmission dynamics in the presence of an imperfect vaccine by a stage-structured mathematical model. To interpret the "psychological effect" when an infectious disease being spread in a population, we have considered a general nonmonotone nonlinear incidence rate function. We have shown that the solutions of the proposed model exist (uniquely determined) and they are nonnegative and bounded, that is, the model is biologically well-posed.

Due to the existence of time delay, 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 [16, 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 used the semi-discretization method to analyze the linear stability of the endemic equilibrium. Also, 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 noticed any influence of the latent time τ on ρ when the vaccine efficacy is fixed. However, we have observed that the value of ρ increases as R 0 increases when is fixed. The quantity 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. In other words, R E cannot become below the unity even when ψ becomes high. Through the numerical simulations:

• We have fitted the model with data of influenza patients as a case study. We have noticed that at the peak level of infection, nearly 6% of the population is infected; • We have carried out global and local sensitivities analysis for R E . We have found that the latent time τ has a noticeable effect on R E . Regarding 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; • We have constructed an approximate stability region for the endemic equilibrium E 1 and noticed that when E 1 loses its stability, a unique stable periodic solution exists via Hopf Bifurcation. For example, for small τ , the endemic equilibrium is stable and as the value of τ increases, E 1 loses the stability and a Hopf Bifurcation occurs and a unique stable periodic solution exists. Moreover, we have observed that the model exhibits global asymptotic stability behavior when R E > 1. 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 [3] . For certain infectious diseases, asymptomatic carriers are a potential source for transmission such as Typhoid Fever, HIV and, most recently, the COVID-19.

Stability analysis for a vector disease model. Rocky Mt

Global stability of an SVEIR epidemic model with ages of vaccination and latency

A periodic disease transmission model with asymptomatic carriage and latency periods

Analysis of an SEIRS epidemic model with two delays

Modeling diseases with latency and relapse

A time-delayed epidemic model for ebola disease transmission

Stability behavior of a nonlinear mathematical epidemic transmission model with time delay

Global stability of a delayed epidemic model with latent period and vaccination strategy

Dynamics of vaccination in a time-delayed epidemic model with awareness

A delay SIR epidemic model with pulse vaccination and incubation times

Pulse vaccination of an SEIR epidemic model with time delay

Suboptimal control and targeted constant control for semi-random epidemic networks

Epidemic spreading in scale-free networks

On vaccine efficacy and reproduction numbers

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 antisars vaccine

An epidemiology model that includes a leaky vaccine with a general waning function

Dynamical analysis and control strategies of an SIVS epidemic model with imperfect vaccination on scale-free networks

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

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

Global results for an epidemic model with vaccination that exhibits backward bifurcation

An Introduction to

Backward bifurcation of an epidemic model with treatment

Introduction to Functional Differential Equations

Monotone Dynamical Systems: An Introduction to the Theory of Competitive and Cooperative Systems

Ordinary Differential Equations and Dynamical Systems

Basic reproduction ratios for periodic compartmental models with time delay

Estimation of effective reproduction numbers for infectious diseases using serological survey data

The Stability of 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

Permanence implies the existence of interior periodic solutions for FDEs

Delay Differential Equations

Semi-discretization for Time-delay Systems: Stability and Engineering Applications

Applied Nonlinear Dynamics: Analytical, Computational, and Experimental Methods

Semi-discretization method for delayed systems

Nonlinear oscillations in epidemic models

The basic reproduction number (r 0 ) of measles: a systematic review

Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission

Optimal two-phase vaccine allocation to geographically different regions under uncertainty

CDC: Estimates of influenza vaccination coverage among adults in the United States

IAC: Ask the experts-influenza

CDC: Past seasons vaccine effectiveness estimates

WHO: How flu spreads

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

A methodology for performing global uncertainty and sensitivity analysis in systems biology

Sensitivity analysis of a nonlinear lumped parameter model of HIV infection dynamics

Acknowledgements The author would like to thank the referees for their careful reading and helpful suggestions.

Conflict of interest The author declares that he has no conflict of interest.