New global dynamical results and application of several SVEIS epidemic models with temporary immunity.

This work applies a novel geometric criterion for global stability of nonlinear autonomous differential equations generalized by Lu and Lu (2017) to establish global threshold dynamics for several SVEIS epidemic models with temporary immunity, incorporating saturated incidence and nonmonotone incidence with psychological effect, and an SVEIS model with saturated incidence and partial temporary immunity. Incidentally, global stability for the SVEIS models with saturated incidence in Cai and Li (2009), Sahu and Dhar (2012) is completely solved. Furthermore, employing the DEDiscover simulation tool, the parameters in Sahu and Dhar'model are estimated with the 2009-2010 pandemic H1N1 case data in Hong Kong China, and it is validated that the vaccination programme indeed avoided subsequent potential outbreak waves of the pandemic. Finally, global sensitivity analysis reveals that multiple control measures should be utilized jointly to cut down the peak of the waves dramatically and delay the arrival of the second wave, thereinto timely vaccination is particularly effective.

Introduction

Immunization is believed to be one of the most successful and cost-effective public health interventions [1], for instance in worldwide eradication of small-pox and sharp reduction in the annual morbidity of most other vaccine-preventable diseases, such as polio, measles, hepatitis B, yellow fever [2], cholera [3], mumps [4] and influenza [5], [6], [7], [8]. Currently, immunization saves 2–3 million lives yearly and prevents debilitating illness, disability and death from the diseases. However, it is estimated that 19.4 million infants failed to be reached with routine immunization services in 2018 [1]. Due to low vaccination rate, the 2017–2018 seasonal influenza caused estimated 45 million illnesses, 21 million medical visits, 810,000 hospitalizations and 61,000 deaths in the United States [9], and now burden is not optimistic. Fortunately, timely vaccination programme played a core part in mitigating the pandemic (H1N1) 2009 [8] (pH1N1). Take Hong Kong China for instance: the subsequent potential waves of the pandemic [10] might be effectively mitigated with the launch of the pH1N1 vaccination programme for several priority groups [11], although the first wave failed to be timely contained due to the unavailability of the vaccine against the novel influenza strain [12] (see Fig. 1 ).

Fig. 1

Epidemic curve of the reported pH1N1 cases in Hong Kong China, 2009–2010.

Admittedly, immunization may not be once and for all because vaccine-induced immunity is generally temporary, and so are disease-acquired and natural immunity, which becomes one of major obstacles eliminating such these infectious diseases. Vaccines rarely provide the recipients with almost life-long immunity against re-infection. After being infected, susceptible individuals first become exposed but not infectious and then become infectious. The successfully recovered individuals acquire disease-induced immunity. Additionally, by virtue of natural immunity [13], [14], [15], a part of exposed individuals fail to develop disease but acquire temporary immunity. For example, the efficient innate immunity protects more than 90% of individuals infected with Mycobacterium tuberculosis [14]. A recent study [15] has showed that, similar to seasonal influenza, most infection (up to 75%) of the pandemic H1N1 strain was asymptomatic and gave the infected individuals temporary immunity.

The nonlinear epidemic dynamical models incorporating both temporary immunity and latency such as SEIRS, SVEIS models in [16], [17], [18], [19], have been developed to better understand the transmission dynamical behaviors of infectious diseases qualitatively and quantitatively. The exploitation of their global asymptotic stability has been of great interest and challenging to researchers in infectious disease modelling aimed at finding out the effective control interventions, seeing [16], [17], [18], [19]. While the Lyapunov function methods may become unsuitable to prove their global stability, the classical geometric approach for nonlinear autonomous differential equations based on additive compound matrix theory developed by Li and Muldowney [20], [21], [22] has been succeed in applying to these epidemic models [16], [17], [20], [21], [22]. For example, Cai and Li [16] proposed the following nonlinear SEIV epidemic model with temporary immunity:where the total population N consists of susceptible (S), latent (E), infectious (I) and vaccinated-recovered (V) classes. The nonlinear incidence βSI/φ(I), with and φ′(I) ≥ 0, generalizes saturated incidence and nonmonotone incidence capturing psychological effect [23], [24]. Along the work of [16], Sahu and Dhar [17] further developed a nonlinear SVEIS model with partial temporary immunity as follows:where susceptible class are vaccinated with certain vaccine at constant rate α, different from model (1.1) with a fraction of vaccinated newborns (denoted by p). We always assume that the same parameter represents the identical biological meaning throughout this paper, and the detailed biological descriptions of the parameters for model (1.2) are demonstrated in Table 1 . Note that [16], [17] applied the geometric approach based on the second additive compound matrix theory of [20] to the responding limiting systems and achieved global stability of the unique endemic equilibrium (EE) under the vaccination reproduction number and some additional restrictions. More recently, Lu and Lu [18], [19] improved the classical geometric approach of [20], [21], [22] and generalized the geometric criterion on global-stability problem and applied it to several nonlinear SEIRS models, successfully removing some restrict conditions on global stability of their EE.

Table 1

Notation description for model (1.2) and their values.

Notation	Description	Units	Range	Baseline	Source
Π	Recruitment rate	m ^−1	[0,6748]	130	Assumed
μ	Natural death rate	m ^−1	–	9.815×10^−4	[25]
1/γ	The mean infectious period	m	[0.1333,0.3333]	0.2333	[5], [26], [27], [28]
1/ω	Average time of immunity waning	m	[6,12.1655]	12.1655	[6], [7]
α	Vaccination rate	m^−1	[0,1]	0	[5], [7]
ξ	The recovery rate of exposed class
	due to natural immunity	m^−1	[3,30]	4.2857	fitted
β	The disease transmission coefficient	m^−1·p^−1	[0,1]	7.0219×10^−5	Fitted
1/σ	The latent period	m	[0.0333,0.1667]	0.1116	Fitted
κ	The inhibition effect	–	[0,1]	1.3458×10^−13	Fitted
q	Fraction of recovered individuals
	from disease developing immunity	–	[0,1]	0.9287	Fitted
S(0)	Initial value for susceptible class	p	[0, 7 × 106]	1.2959 × 105	Fitted
V(0)	Initial value for vaccinated class	p	[0, 7 × 106]	2.7970 × 105	Fitted
E(0)	Initial value for exposed class	p	[0, 7 × 106]	10	Assumed
I(0)	Initial value for infectious class	p	–	23	[12]

[Note: m,  p  represent month and person, respectively.]

Borrowing the ideas of [16], [17], [23], [24], we establish the following SVEIS epidemic model with general nonlinear incidence:in which, it is assumed that vaccine-induced, disease-acquired and natural immunity may last the nearly same time for some diseases like influenza, and the differential infectious force function g possesses the following properties reflecting some biological significances:

(P1) satisfies g(I) > 0 for I > 0.

(P2) is monotonously nonincreasing for I > 0, and .

(P3) I|g′(I)| ≤ g(I) for I > 0.

It is worth highlighting that saturated and nonmonotone incidences in [23], [24], [29] and [30], but not confined to them, fulfill (P1)-(P3), thus we lift restrictions on monotonicity of g(I) in spite of the introduction of (P3). With this geometric criterion in [18], we shall thoroughly address global threshold dynamics of models (1.3) and (1.2), characterized by their vaccination reproduction numbers. Incidentally, the unnecessary restrictions both in Theorem 4 in [16] and Theorem 5.5 in [17] are completely removed since model (1.3) reduces to model (1.1) if and . Of particular note is that we achieve global asymptotic stability for model (1.1) of [16] with nonmonotone incidence reflecting psychological effect, which also reserves threshold dynamics.

Furthermore, as an application of model (1.2), the reported pH1N1 case data of Hong Kong China [12] are utilized to estimate its parameters, aimed at accounting for the avoidance of the subsequent potential waves of the pandemic in 2010 (as predicted by WHO [10]) with the pH1N1 vaccination programme. Meanwhile, several disease-control measures are evaluated in terms of global sensitivity analysis for the vaccination reproduction number. In particular, this study arrives at a conclusion that joint usage of multiple control measures such as isolation, vaccination and treatment, can more effectively cut down the peak of the waves and dramatically delay the arrival of the second wave at the same time.

The outline of this paper is summarized as follows. In Section 2, we offer insight into global threshold dynamics for model (1.3), including the existence, local and global asymptotic stability of its equilibria. Section 3 completely addresses the global dynamics of model (1.2). Section 4 performs parameter estimation and global sensitivity analysis for the vaccination reproduction number of model (1.2) with the purpose of seeking for effective control measures. Finally, we close the paper with a conclusion and discussion section.

Global threshold dynamics for model (1.3)

The existence of the equilibria

For model (1.3), one can easily obtain that the biologically feasible regionis the positively invariant set by similar arguments in [16]. Apparently, the disease-free equilibrium (DFE) of model (1.3) always exists, where . Thus, by application of the next generation matrix approach in [31], the vaccination reproduction number (e.g., seeing [32], [33]) is calculated asclearly remaining the same with the model in [16] when .

By some direct but tedious algebra operations, it can be deduced that the I* component in the EE is determined by the following equationwhereIn what follows, we are going to focus mainly on analyzing the positive real solution of Eq. (2.2). A simple induction then showsIt deduces from (P2) that .

In the case of together with G′(0) > 0, and it can be revealed that G(I) > 0 as I is sufficiently small, guaranteeing the existence of positive real root for Eq. (2.2) from Fig. 2 , denoted by I*. And its uniqueness is verified by reduction to absurdity as follows. Provided that another positive solution I * of (2.2) nearest to I*, if it exists, must satisfy G′(I *) ≥ 0 owing to the continuity of G(I). Actually, together with g′(I *) ≤ g(I *)/I * deduced from (P3), we arrive atwhere one utilizes the equality derived by the equations that the EE satisfies. An obvious contradiction exists as shown in Fig. 2. Thus, the positive solution I* is unique, which can lead to the uniqueness of S*, V*, E* from the analysis above.

Fig. 2

The existence and uniqueness of positive real root for Eq. (2.2).

In the case of Eq. (2.2) must admit no positive solution. Otherwise, let I ⋆ be its smallest one. Combining and G′(0) ≤ 0 yields that G(I) ≤ 0 for sufficiently small I. Since the function G(I) continuously increases to 0 from the non-positive value, it is clear to see that G′(I ⋆) ≥ 0, which contradicts with G′(I ⋆) < 0 in (2.4).

To sum up, model (1.3) has a unique EE P* if and only if (iff) .

For model (1.3) , a DFE P 0 always exists and the EE P* exists uniquely iff .

Local stability

(i) The DFE P 0 is local asymptotically stable (LAS) if but becomes unstable if ; (ii) The EE P* is LAS iff .

The Jacobian matrix of model (1.3) takes the following form of(i) The characteristic equation at P 0 isObviously, its all eigenvalues possess negative real parts when that is, P 0 is LAS. If there exists a positive root, so the DFE becomes unstable.

(ii) Calculating the characteristic equation at P*, one reacheswhere . Clearly, .

Case I. Let g′(I*) > 0. One asserts that all eigenvalues of the following equationsatisfy Reλ < 0. Suppose, for contradiction, that there exists one eigenvalue with . From (2.8) and (P3), the following contradiction is attained Case II. Let g′(I*) ≤ 0. Equality (2.8) is recast as . For the Routh-Hurwitz conditions can be ensured byWe thus infer that all eigenvalues obey Reλ < 0. Combining Cases I and II leads to local stability of P* for . □

Global stability

The DFE P 0 of model (1.3) is GAS in Ω if .

By the first equation of (1.3) and it is easy to ascertain thatwhich asserts that S ≤ S 0 (similar to [4]). Otherwise, let us suppose that S > S 0, thus dS/dt < 0. It follows that S ≤ S 0 when S(0) ≤ S 0, which is absurd as our assumption. Hence, our claim S ≤ S 0 is valid. Observe that g(I)/I ≤ β for I > 0 can be ensured by (P3) (seeing, e.g., [34]). Construct Lyapunov function and its time derivative of W(t) along the solutions of model (1.3) is estimated asprovided that . From the LaSalle’s Invariance Principle [35] and local stability of P 0 in Theorem 2.2, we can derive its global asymptotic stability for . □

In the sequel, we shall employ the general criterion for global stability for autonomous differential equations developed by [18] to establish global stability of the EE P* of model (1.3). A brief outline on this geometrical approach [18], [20], [21], [22] is presented as follows.

Let us consider the nonlinear autonomous dynamical system:where the function and Q is an open set. For (2.9), the solution with is defined as x(t, x 0) and its equilibrium as x*. Moreover, let us assign satisfying when . We assume that system (2.9) admits a dimensional invariant manifold defined by .

The following three hypotheses are satisfied:

(H1) Γ is simply connected.

(H2) There is a compact absorbing set D ⊂ Q ⊂ Γ.

(H3) System (2.9) admits a unique equilibrium x* in Γ.

The general geometric criterion of Lu and Lu is recapped as follows.

(see Theorem 2.6 in ). The unique equilibrium x* of (2.9) is globally asymptotically stable (GAS) in Γ provided that

(x(0, x 0)) of system (2.9) , there are a matrix C(t), a sufficiently large τ 1 > 0 and constants such that and where b(t) and c(x(0, x 0)) and C(t), respectively.

Denote the interior, the boundary of Ω by and ∂Ω, respectively. Uniform persistence in of model (1.3) for can be deduced from the instability of P 0 and P 0 ∈ ∂Ω.

Model (1.3) is uniform persistent in if .

The EE P* of model (1.3) is GAS in if .

The third additive compound matrix of [22] for model (1.3) acquires the formwhere .

Assign with . The invariant manifold for (1.3) is . Following [22], it turns out to be and . In the sequel, let and be the 4 × 4 identity matrix. Then the coefficient matrix readswhere .

Meanwhile, model (1.3) can be recast into

Note that Theorem 2.4 implies that there is a constant π 0 > 0 such that π 0 ≤ S, V, E, I ≤ Π/μ. It follows from (P1) that there are constants l, L > 0 such that l ≤ g(I) ≤ L. Assign π ≔ μπ 0/Π. By I|g′(I)| ≤ g(I) in (P3) and (2.12), c(t) are respectively estimated asChoose the matrix C(t) in Lemma 2.1 as . It is easy to check that where . By Lemma 2.1, the EE is GAS in . □

Let and then model (1.3) reduces to the model with saturated incidence of [16], which retains global threshold dynamics from Theorem 2.5, improving Theorem 4 in [16]. More importantly, the sharp threshold dynamics result is extended to the model with nonmonotone incidence capturing psychological effect of [16].

Global threshold dynamics for model (1.2)

In this section, for simplicity, we take satisfying (P1), (P2) and

(P3)′ Ig′(I) ≤ g(I) for I > 0.

Following the same reasoning as the proof of Theorems 2.1–2.2 in Subsections 2.1–2.2, one easily draws the following conclusions on the existence, local stability of the DFE and the EE for model (1.2), where .

Model (1.2) always has a DFE and the EE is unique if the vaccination reproduction number

The DFE is LAS when and it is unstable but the EE is LAS when .

In what follows, we make a thorough inquiry into global stability of model (1.2). Using the similar arguments as the analysis of Theorems 2.3-2.4 in Subsection 2.3 can lead to global stability of the DFE and persistence of model (1.2) as follows.

If the DFE is GAS in Ω.

If model (1.2) is uniform persistent in .

In order to achieve global stability of the EE, we focus mainly on the significant differences and skip the repeated parts with the proof of Theorem 2.3 in Subsection 2.3. The coefficient matrix B(t) for model (1.2) is calculated aswhere . And model (1.2) can be transformed into

Clearly, g(I) meets (P1), (P2) and (P3)′, and g′(I) > 0 for I > 0. Uniform persistence ensures that there exists positive constants π 0, l, L such that π 0 ≤ S, V, E, I ≤ Π/μ, and l ≤ g(I) ≤ L. Let π ≔ μπ 0/Π. Two cases will be considered to estimate c 1(t).

Case 1. . Employing (3.2), g′(I) > 0 and (P3)′ results in Case 2. . Similar proof in Theorem 2.3 gives S ≤ S 0, being equivalent to . By g′(I) > 0, we can arrive at

We can similarly infer thatBy applying Lemma 2.1, the above is concisely stated into Theorem 3.5.

The EE of model (1.2) is GAS in if .

An immediate consequence of Theorem 3.5 yields global threshold dynamics of model (1.2), getting rid of the unnecessary restrictions in Theorem 5.5 from [17]. Additionally, model (1.2) with the incidence satisfying (P1),(P2) and (P3)′, e.g., [29], [30], also reserves global threshold stability by the same proof.

From the analysis in Sections 2 and 3, it can be similarly verified that the following SVEIS model with temporary immunity and nonlinear incidence satisfying (P1)-(P3) is a sharp threshold system characterized by its vaccination reproduction number,

An application of model (1.2)

Vaccination was the most cost-effective intervention for mitigating the 2010 influenza A(H1N1) pandemic. On 28 August 2009, WHO advised that the countries in the northern hemisphere should prepare for a second wave of pandemic spread [10]. Fortunately, the pH1N1 vaccination programme for five priority groups was launched, such as medical workers, pregnant women, people over 65 or with chronic illness, children aged between 6 months to 6 years [11]. Because the susceptible individuals aged over 6 months were vaccinated with the pH1N1 vaccine instead of newborns and up to 75% of H1N1 infection was asymptomatic due to nature immunity [15], model (1.2) is applied to illustrate that vaccination effectively contained subsequent potential waves of the pandemic (H1N1) 2009 in Hong Kong China in this section.

Data

At the end of every month from May 2009 to October 2010, the pH1N1 case data of Hong Kong were released by official website of Center for Health Protection, Hong Kong China (available at https://www.chp.gov.hk/sc/statistics/data/10/26/43/416.html [12]), and the data from May 2009 to June 2010 are chosen to fit the parameter values of model (1.2) owing to its high smooth degree (see Fig. 1). Indeed, the prevalence level of from July to October 2010 showed the small fluctuations and kept low (also seeing [8]). The first wave of the pandemic failed to be avoided (see Fig. 1) since there was no available vaccine against the novel influenza strain before 21 December 2009. It was on that day, the pH1N1 vaccination programme for five priority groups was launched and started [11] to minimize any potential second wave and 4182 doses of pH1N1 vaccine were administered [36]. Notice that the vaccine recipients will develop immunity in about 15 days [7] (delayed vaccination, e.g.,[2]), so the start time of generating vaccine-induced immunity can be approximated as 1 January 2010 as shown in Fig. 3 (a).

Fig. 3

(a) Comparison of the reported pH1N1 case data in Hong Kong China and the simulated solution I(t) of model (1.2); (b) The second wave of the H1N1 pandemic is observed through simulation using the estimated parameter values if the pH1N1 vaccination programme had not been carried out.

Parameter estimation

The intervals or values of parameters and initial condition of model (1.2) are estimated (as shown in Table 1) and explained as follows.

According to Subsection 4.1, we set vaccination rate during the 2009 pandemic, but α in (0,1] during the 2010 pandemic from [5]. The vaccine effectiveness is up to 99% [37], thus the vaccine is considered to be perfect.

Since life expectancy is about 83.74 years in Hong Kong in 2010 [25], natural death rate per month (m).

Following [5], [26], [27], [28] and [6], [7], the infectious duration varies from 4 to 10 days and the immunity period changes in the scope of 180 days to 2 years, respectively, so and 1/ω ∈ [6, 24.3333]. Let us take the infectious duration and the immunity period as 7 days [27], [28] and 1 years [6], respectively, then m and m.

The latent period (1/σ) ranges from 1 day to 5 days according to Refs[5], [26], [27], [28]., then 1/σ ∈ [0.0333, 0.1667]. From [5], [26], [28], it may be realistic for the influenza A(H1N1) to consider that exposed individuals recover after 1–10 days due to natural immunity, namely, ξ in [3,30]. It is not hard to obtain that the values of parameters q, β and κ belong to [0,1] based on some existing works (e.g., [17], [23]).

The number of births of Hong Kong in 2009 [38] was 8.21 × 104 per year, namely, 6748 . Considering that vast majority of newborns were taken protective measures, about 2% of the number of births is chosen as recruitment rate of S class, so m. The number of total population of Hong Kong during 2009–2010 [38]) was about 7.0 × 106, thus . Together with the case data [12], the initial value is fixed. And we assume that .

Above all, the values of the remaining parameters β, ξ, σ, q, κ and the initial values S(0), V(0) are estimated (seeing Table 1) with the 8 cases data from May to December 2019 by the DEDiscover simulation tool [39], where we choose the method of Hybrid DESQP Optimization Algorithm, combining global differential evolution and local sequential quadratic programming. From the parameter estimation results above, the values of tend to 0 and 1, respectively. This entails that several standard model selection criteria are employed to evaluate the superiority of models fitting the data [40], including Akaike information criterion (AIC) and Bayesian information criterion (BIC), and their variations such as AICc, with their smaller values corresponding to a better model. It can be observed from Table 2 that model (1.2) with and is selected as the best model by the criteria above, and its simulation results are presented in Fig. 3(a). This suggests that the simple mass action incidence βSI may appropriately reflect the short-term transmission process of the emerging influenza A(H1N1) virus and partial temporary immunity should be incorporated into the influenza models. Furthermore, we analyze the error of fitting to evaluate the performances and reliability of model (1.2) with and and MAPE (the mean absolute percentage error) and RMSPE (the root mean square percentage error) are computed as MAPE =37.36%, RMSPE =6.76%, respectively. Based on the criteria of MAPE and RMSPE in [41], [42], our model can yield reasonable forecasting results. Lastly, from Table 1 it can be checked that the latent period and are in agreement with the reality.

Table 2

Model selection for the pH1N1 case data.

Model	AIC	BIC	AICc
κ=1.3458×10^−13,q=0.9287	113.0446	113.1240	113.7112
κ=0,q=0.9287	113.0446	113.1240	113.7112
κ=0,q=1	113.8155	113.8950	114.4822

The results of parameter estimation above yield that the vaccination reproduction number of 2009 is computed into which is consistent with the conclusion in [28], [43] (ranging from 1.2 to 2.3). From Theorems 3.4 and 3.5, the disease may be persistent and become endemic. Without vaccination, as forecasted by WHO [10], the second wave is indeed observed through simulation using the estimated parameter values (see Fig. 3(b)), thus vaccination is imperative if the vaccine is available. Furthermore, vaccination rate is estimated with the case data from January to June 2010 (other parameter values remain the same with Table 1, and initial condition (93492,312020,627,1287) is the simulation result in December 2009, corresponding to such that the pandemic was contained quickly, as proved in Theorem 3.3 and shown in Fig. 3(a).

Sensitivity analysis

The vaccination reproduction number of model (1.2), measuring the average number of secondary cases that are caused when one index case is introduced into a disease-free population [32], [33] in which a vaccination programme is carried out, may determine the transmissibility, severity and outcome of the pandemic. In order to seek for effective disease-control measures, we therefore shall be concerned with the effects of input parameters ω, β, α, γ, ξ on . Based on Latin Hypercube Sampling (LHS) and partial rank correlation coefficients (PRCCs) [44], global uncertainty and sensitivity analysis for is conducted to reveal the influence degree on model outcomes. These interesting parameters are considered to obey normal distributions with means coming from baseline values given in Table 1. And their PRCC values are computed through 5000 simulations per run and demonstrated in Fig. 4 (a) and Table 3 .

Fig. 4

Sensitivity analysis: (a) PRCC values for the vaccination reproduction number ; (b)β is decreased by 10%, (c)α is increased by 10% and (d)γ is increased by 10% of their baseline values in Table 1, respectively.

Table 3

PRCC values for which are ranked from the most sensitive to the least.

Parameter	Mean	Standard deviation	PRCC	p-value
β	7.0219×10^−5	2.3406×10^−5	0.9176	0
α	0.3527	0.1175	−0.8824	0
ω	0.0822	0.0137	0.6496	0
γ	4.2857	0.4286	−0.5418	0
ξ	4.2857	0.4286	−0.3103	0.1979×10^−110

Finally, numerical simulations are carried out to evaluate the effectiveness of disease-control measures. In Table 3, input parameters β, α, ω, γ, ξ are ranked in descending order according to their influences on new infections. In fact, it seems difficult to prolong immunity duration related to the parameter ω. For this reason, we only consider the impacts of parameters β, α and γ. In detail, β has positive impact on and α, γ have negative impacts on it. Thus, we decrease the value of β by 10% and increase the value of γ by 10%, respectively. As discussed above, vaccination was such an effective health intervention, that the H1N1 pandemic was successfully curbed in 2010. In consideration of frequent outbreaks of current seasonal flu (including influenza A(H1N1), B and C) epidemics in many countries, such as the United States [9] and China with low vaccination rate, it may be interesting and significant to assume that the vaccine is available and vaccination is carried out at the begin of the pandemic. 10% and 20% of vaccination rate are used to study the effect of vaccination on the pandemic. And the other parameter values and initial values of Table 1 are fixed. Simulation results are presented in Fig. 4 (b)-(d).

Undoubtedly, reducing the disease transmission coefficient β, such as epidemic propaganda, isolation, sterilization and wearing a mask, cuts down the peak of the first wave and delays the arrival of the second wave, but its two peak values fail to decrease obviously even though parameter β is the first sensitive, seeing Fig. 4 (b). On the other hand, increasing vaccination rate α and shortening the disease course of disease γ (e.g., antiviral therapy) lower more dramatically the peak values of the first and second waves than reducing β, but the peak of the second wave arrives much earlier than reducing β (as shown in Fig. 4 (c) and (d)). Therefore, it is possible for policymaker to use multiple control measures jointly during the influenza pandemic. It is also acknowledged that timely vaccination is particularly effective at reducing the outbreak peaks than the other two measures.

Conclusion and discussion

Immunization has been bringing mankind great success to prevent the disease transmission every year [2], [3], [4], [5], [6], [7], [8], [1], and a long latent period of infectious disease may generate dramatically different model prediction and thus allows of no to neglect [26]. What’s more, nonlinear incidence can reproduce the inhibition effect from behavioral changes of individuals and the impact of other factors like severity and stage of the infection [16], [17], [45]. The current work formulates an SVEIS model with vaccination, latency, nonlinear incidence and temporary immunity and establishes its global threshold stability by a novel geometric criterion in [18]. Most pointedly, the open questions on global threshold stability of their EE for two nonlinear SVEIS models with saturated incidence in [16], [17] are also well addressed. Inspired by [18], the introduction of the property (P3) on the infectious force function g(I) leads us to successfully achieve global threshold dynamics for the SVEIS models with nonmonotone incidence reflecting psychological effect. Furthermore, let then an application of Theorem 2.5 yields that model (1.1) is a sharp threshold system provided that φ(I) meets and 0 ≤ Iφ′(I) ≤ 2φ(I), such as for 0 < r ≤ 2.

In 2009, the novel influenza A(H1N1) virus caused the first pandemic of 21st century. We apply model (1.2) to illuminate the avoidance of the potential second wave of the pandemic (H1N1) 2009 in Hong Kong, China (as predicted by [10]) with the pH1N1 vaccination programme, and it is revealed that timely vaccination is more effective at lowering the outbreak peaks than other measures. This offers solid support for implementation of immunization strategy to cope with current global seasonal influenza burden, measles cases surge and COVID-19 pandemic if the vaccines are available.

This research is also subject to several limitations as follows. In details, observe that HBV vaccine is administered to both newborns and susceptible individuals, so both two vaccination ways can be incorporated into these SVEIS models, which, together with [4] we guess, can still preserve the threshold dynamics since insights provided by several SVEIS models studied above, can inform us that vaccination for either newborns or susceptible individuals and temporary immunity fail to change their threshold stability (see Theorems 2.5, 3.5 and Remark 3.2). Additionally, we just consider the nonlinearity of incidence rate on I, perfect vaccines, constant total population and postulate that vaccine-induced and disease-acquired immunity last the same time. It would be interesting to introduce more general incidence S ϱ f(I) (ϱ > 0), distinct vaccinated class (V) and recovered class (R), incomplete vaccination and varying total population size (e.g., [4], [18], [19], [21], [45]) to improve the accuracy of model prediction. Certainly, more analytical techniques are needed, and these issues are left as future investigations.