# Modelling the Periodic Outbreak of Measles in Mainland China

Authors: Yuyi Xue; Xiaoe Ruan; Yanni Xiao
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

In mainland China, measles infection reached the lowest level in 2012 but resurged again after that with a seasonally fluctuating pattern. To investigate the phenomenon of periodic outbreak and identify the crucial parameters that play in the transmission dynamics of measles, we formulate a mathematical model incorporating periodic transmission rate and asymptomatic infection with waning immunity. We define the basic reproduction number as the threshold value to govern whether measles infection dies out or not. Fitting the reported measles cases from 2013 to 2016 to our proposed model, we estimate the basic reproduction numberR0 with immunization to be 1.0077. From numerical simulations, we conclude asymptomatic infection does not cause much new infections and the key parameters affecting the transmission of measles are vaccination rate, transmission rate, and recovery rate, which suggests the public to enhance vaccination and protection measures to reduce effective contacts between susceptible and infective individuals and treat infected individuals timely. To minimize the number of infected individuals at a minimal cost, we formulate an optimal control system to design optimal control strategies. Numerical simulations show the effectiveness of optimal control strategies and recommend us to implement the control strategies as soon as possible. In particular, enhancing vaccination is especially effective in lowering the initial outbreak and making disease recurrence less likely.


---

## Body

## 1. Introduction

Measles is an acute, viral infectious disease characterized by high fever, cough, and a maculopapular rash [1]. In 1980, 2.6 million people died from measles. In 2005, the World Health Organization Western Pacific Region set a goal to eliminate measles before 2012. For achieving it, China has taken some immunization measures including continuing a two-dose routine measles vaccine strategy (first and second dose at 8 and 18–23 months) while using supplementary immunization activities (SIAs) [2]. Because of these measures, the number of reported measles cases significantly decreased from approximately 140000 in 2008 to 6678 in 2012. However, in 2013, a resurgence of measles occurred and measles cases obviously increased to 30000 and retained at around 28000 in 2016 with obvious periodic feature. Consequently, measles still remains a serious public health problem in mainland China and how to eliminate it becomes challenging.Mathematical models play an essential role in investigating the transmission dynamics of the infectious diseases and helping to identifying the key parameters or process that significantly affect disease spread. A number of mathematical models have been formulated to study the transmission and control of measles infection. Hamer in 1906 [3] proposed a model to study the regular occurrence of measles in London. Then, a prediction model was developed which successfully predicted the 1997 measles epidemic of New Zealand. Roberts and Tobias [4] extended this model by including age structure to examine optimal timing of vaccine doses and compare the effects of various potential vaccine strategies. In 2005, Pang et al. [5] proposed a SEIR model with vaccination to investigate oscillations of measles and fit the case data of U.S. from 1951 to 1962 to design a control strategy. Considering some developing countries with high measles burden and limited resources, Verguet et al. [6] proposed a model to estimate optimal scheduling of SIAs without second routine vaccination. In contrast, for some countries where measles has been eliminated, Choi et al. [7] estimated susceptibility to measles by age and Mckee et al. [8] investigated the optimal age target for two routine dose without changing vaccination coverage in order to maintain elimination. Furthermore, Bai and Liu [9] fitted the model to the case data in China and Huang et al. [10] investigated the effect of various interventions on measles infection. Thorrington et al. [11] attempted to quantify the loss of QALY (quality-adjusted life years) due to measles at a population level in England and provided important parameters for future control interventions. In addition, some researchers found critical community size and herd immunity influence the extinction of measles in communities [12–16], and obtained that human mobility has an impact on the periodicity of measles outbreaks [17].There are studies showing the evidence of waning immunity in vaccinated individuals [18–20], indicating individuals with sufficiently low antibody levels are at risk of mild or subclinical measles infections [21, 22]. Although transmission of virus between such individuals with asymptomatic infection has not been demonstrated, serological studies have indicated that this may occur [23, 24]. In China, measles cases can still be observed in recent years even if we have initiated enhanced immunization and got high coverage rate of vaccination. Possible reasons may include mature women with low antibody can either be infected by their infected children without showing any symptom and in turn infect others or have babies who have low antibody and can be infected. Hence, the exploration of the influence of asymptomatic infection on the transmission of measles and identification of the key factors or processes that significantly influence the measles outbreak in mainland China remain unclear and fall within the scope of this study.The main purpose of this paper is to formulate a dynamic model incorporating periodic transmission and asymptomatic infection with waning immunity to identify the key parameters which influence seasonal fluctuation of measles and discuss optimal control measures to minimize the number of infected individuals and costs. In the next section, a periodic model with vaccination and asymptomatic infection due to waning immunity is introduced and then the existence and stability of periodic solution are discussed. By fitting monthly measles data in mainland China from 2013 to 2016, we estimate some unknown parameters and calculate the basic reproduction number. We then numerically evaluate effect of asymptomatic infection and other important parameters on the transmission dynamics of measles. In Section3, we apply optimal control theory to design optimal control strategies such that the number of infected individuals and costs are minimum. Finally, a discussion is presented in Section 4.


## 2. Periodic Model

### 2.1. Model Formulation

The underlying structure of the model comprises of eight compartments. We assume that susceptible individuals (S) receive the vaccine and enter into the compartment V1 at rate of p. Then, the vaccinated individuals (V1) progress to the class (V2) with waning of immunity at rate of ω. Susceptible individuals (S) become infected and enter into the latent undetectable class (E) at rate of β(t) (I + ξIv), where β(t) is the baseline periodic transmission rate and ξ (0 ≤ ξ ≤ 1) quantifies the relative transmissibility of asymptomatic infections. Then, the exposed individuals (E) progress to infectious class (I) at rate of σ and recover with rate of γ. Note that most individuals with waning immunity may not produce symptoms, once infected again; hence, it is negligible that latent individuals (E) move to asymptomatic infections class (Iv). The individuals in V1 class are those who are effectively vaccinated with high antibody titer, so we assume they cannot be infected as mentioned in [19, 20]. But because the antibody titer is low, individuals located in V2 may be infected and become exposed individuals Ev at rate of ηβ(t) (I + ξIv) with the modification factor η (0 ≤ η ≤ 1), denoting these individuals having low probability to be infected. Then, these latent individuals may move to symptomatic or asymptomatic infectious class (I) and (Iv) at rate of θσ or (1−θ)σ, where θ (0 ≤ θ ≤ 1) is the proportion of latent individuals becoming symptomatic infections. The recovery rate from Iv to R is γv. In general, asymptomatic infections recover faster than symptomatic infection, which means γv > γ. A flow diagram of the model is described in Figure 1 and the definitions of variables and parameters are given in Table 1. The specific dynamic model is as follows:(1)dSdt=Λ−pS−βtI+ξIvS−μS,dV1dt=pS−μ+ωV1,dV2dt=ωV1−ηβtI+ξIvV2−μV2,dEdt=βtI+ξIvS−σ+μE,dEvdt=ηβtI+ξIvV2−σ+μEv,dIdt=σE+θσEv−γ+μI,dIvdt=1−θσEv−γv+μIv,dRdt=γI+γvIv−μR,where  β(t) is nonnegative and periodic function with period  T (T > 0) and Λ and μ are recruitment rate and natural death rate per unit time. All parameters are assumed to positive constants.Figure 1
The flow diagram of the measles transmission.Table 1
Interpretation and values of the parameters used in the simulation.
Parameters
Description
Default value
Sources
μ
Natural death rate (month−1)
1.142 × 10−3
Calculated
σ
Rate of progression to infected stage from the exposed (month−1)
2.17
[11]
γ
Recovery rate for the infected stage (month−1)
3
[11]
Λ
Recruitment rate (month−1)
1589357
[9]
ω
Rate of waning immunity (month−1)
0.0033
[19]
γ
v
Recovery rate for the asymptomatic infected stage (month−1)
6
LS
η
Correction factor
0.01
LS
ξ
Correction factor
0.001
LS
θ
Proportion of progression to symptomatic infected stage
0.004
LS
p
Coverage rate of vaccine for entire population (month−1)
0.00014
LS
β(t)
Periodic transmission rate of symptomatic infection (month−1)
a
1.346
+
sin
2
π
t
/
12
+
2
Assumed
a
Baseline transmission rate (month−1)
1.9616 × 10−9
LS


### 2.2. Threshold Dynamics

It is clear that any solution of this system from nonnegative initial values is nonnegative. Then, we first discuss bound of solution of this system.Proposition 1.
The solution of system (1) is uniformly and ultimately bounded, i.e., there exists M>0,t⌣>0 such that (S, V1, V2, E, Ev, I, Iv, R) ≤ (M, M, M, M, M, M, M, M) for t>t⌣.Proof.
LetS(t) + V1(t) + V2(t) + E(t) + Ev(t) + I(t) + Iv(t) +R(t) = N; then,(2)N′=Λ−μN.
Hence,limt⟶+∞Nt=Λ/μ, which implies for any ε > 0, there exists t⌣>0 for t>t⌣, then Nt<Λ/μ+ε=˙M. Subsequently, We get (S, V1, V2, E, Ev, I, Iv, R) ≤ (M, M, M, M, M, M, M, M), i.e., the solution of system (1) is uniformly and ultimately bounded.


#### 2.2.1. Disease-Free Equilibrium

It is easy to obtain that system (1) has unique disease-free equilibrium E0=S¯,V¯1,V¯2,0,0,0,0,0, where S¯=Λ/p+μ,V¯1=Λp/p+μω+μ,V¯2=Λpω/p+μω+μμ. In the following, the stability of disease-free equilibrium will be discussed. First, we define the reproduction number of system (1) on the basis of [25].Denote(3)Ft=00βtS¯ξβtS¯00ηβtV¯2ξηβtV¯200000000,Vt=σ+μ0000σ+μ00−σ−θσγ+μ00−1−θσ0γv+μ.AssumeY(t, s), t ≥ s is the evolution operator of following linear system:(4)dydt=−Vty.That is to say,Y(t, s) satisfies(5)dYt,sdt=−VtYt,s,∀t≥s,Ys,s=I.Let Φ−V(t) be the fundamental solution matrix of linear system (4). Clearly, Φ−V(t) equals to Y(t, 0), t ≥ 0.Define a linear operatorL : CT ⟶ CT by(6)Lϕt=∫−∞tYt,sFsϕsds=∫0∞Yt,t−aFt−aϕt−ada,where CT  is the ordered Banach space of all T-periodic function from R to R4 with maximum norm ‖·‖ and positive cone CT+:=ϕ∈CT:ϕt≥0,∀t∈R. Suppose that ϕ ∈ CT is the initial distribution of infectious individuals in the periodic environment; then, F(s)ϕ(s) is the distribution of new infections produced by the infected individuals who were introduced at time  s. Given t (t ≥ s), then Y(t, s)F(s)ϕ(s) gives the distribution of those infected individuals who were newly infected at time s and remain in the infected compartments. Define the spectral radius of L as the basic reproduction number of system (1), i.e.,(7)R0=ρL.It can be verified that system (1) satisfies conditions (A1)–(A7) of periodic system [25]; then, the following theorem is established.Lemma 1.
Assume that (A1)–(A7) hold [25]. Then, the following statements are valid:(1)
R0 = 1 if and only if ρ(ΦF−V(T)) = 1(2)
R0 < 1 if and only if ρ(ΦF−V(T)) < 1(3)
R0 > 1 if and only if ρ(ΦF−V(T)) > 1Theorem 1.
IfR0 < 1, the disease-free equilibrium E0 is globally asymptotically stable; if R0 > 1, the disease-free equilibrium E0 is unstable.Proof.
According to Lemma1, it is easy to get E0 is local stable if R0 < 1. The following will show the global attractive.
By system (1), we obtain(8)dSdt≤Λ−pS−μS,dV1dt≤pS−ω+μV1,dV2dt≤ωV1−μV2.
By standard comparison theorem [26], we get for any ε > 0, there exists t1 > 0; if t > t1, then St<S¯+ε,V1t<V¯1+ε,V2t<V¯2+ε.
Consider auxiliary system(9)dE^dt=βtI^+ξI^vS¯+ε−σ+μE^,dE^vdt=ηβtI^+ξI^vV¯2+ε−σ+μEv^,dI^dt=σE^+θσE^v−γ+μI^,dI^vdt=1−θσE^v−γv+μI^v.
We can write (9) as(10)dXdt=Ft−Vt+εMtX,where X = (Ê, E^v, Î, I^v), and(11)Mt=00βtεξβtε00ηβtεξηβtε00000000.
SinceR0 < 1, that is ρ(ΦF(t)−V(t) (T)) < 1, then for any small ε, ρ(ΦF(t)−V(t)+εM(t) (T)) < 1 and consequently, μ1=1/TlnρΦFt−Vt+εMtT<0.
It follows from Lemma 2.1 of [27] that there exists a positive T-periodic function vt = (v1t, v2t, v3t, v4t) such that eμ1tvt is a solution of (9). Choose t¯>t1 and a small value α > 0 such that Xt¯≤αv0; then, we can get Xt≤αvt−t¯eμ1t−t¯,t≥t¯.
By standard comparison theorem [26], we get(12)Et,Evt,It,Ivt≤Xt≤αvt−t¯eμ1t−t¯.
Hence,limt⟶+∞Et=0,limt⟶+∞Evt=0,limt⟶+∞It=0,limt⟶+∞Ivt=0. By the last four equations of system (1), considering their limit system, i.e.,(13)dSdt=Λ−pS−μS,dV1dt=pS−ω+μV1,dV2dt=ωV1−μV2,dRdt=−μR.
We getlimt⟶+∞St=S¯,limt⟶+∞V1t=V¯1,limt⟶+∞V2t=V¯2,limt⟶+∞Rt=0. The global attractivity of E0 has been proved.


#### 2.2.2. Uniform Persistence of the System

In the following, we examine the uniform persistence and existence of positive periodic solutions of system (1).Theorem 2.
IfR0 > 1, system (1) is uniformly persistent, i.e., there is constant η¯>0, such that for all initial values S0,V10,V20,E0,Ev0,I0,Iv0,R0∈R+3×IntR+4×R+, the solution of system (1) satisfies limt⟶+∞infSt,V1t,V2t,Et,Evt,It,Ivt,Rt≥η¯,η¯,η¯,η¯,η¯,η¯,η¯,η¯. And there exists at least one positive T-periodic solution for system (1).Proof.
DefineX = {(S, V1, V2, E, Ev, I, Iv, R) : S, V1, V2, E, Ev, I, Iv, R ≥ 0}, X0=S,V1,V2,E,Ev,I,Iv,R:S,V1,V2,R≥0,E,Ev,I,Iv>0, ∂X0 = X\X0. Define a Poincaré map P:R+8⟶R+8 satisfying P(x0) = u(T, x0), for x0∈R+8, where u(t, x0) is the solution of system (1).
First, it is easy to see thatX, X0 are positively invariant and ∂X0 is relatively closed in X. It follows from the Proposition 1 that the solution of system (1) is uniformly and ultimately bounded. Thus, the semiflow P is point dissipative on R+8 and P:R+8⟶R+8 is compact. By Theorem 3.4.8 of [28], P admits a global attractor A.
DefineM∂ = {(S, V1, V2, E, Ev, I, Iv, R) ϵ ∂X0 : Pm(S, V1, V2, E, Ev, I, Iv, R) ∈ ∂X0, ∀m ≥ 0}. In the following, we claim M∂ = {(S, V1, V2, 0, 0, 0, 0, 0)}. First, it is obvious that {(S, V1, V2, 0, 0, 0, 0, 0)}⊂M∂. Next, certify M∂⊂{(S, V1, V2, 0, 0, 0, 0, 0)}. Assume there exists x0 ∈ M∂\{(S, V1, V2, 0, 0, 0, 0, 0)}. For example, let E0>0,Ev0=I0=Iv0=0; then, I′(0) = σE0 > 0, which implies there exists t2 > 0; when 0 < t < t2, we get I(t) > 0 and E(t) > 0. For 0 < t < t2, we obtain from equations of (1)(14)Evt=e−σ+μtEv0+∫0tηβsI+ξIvsV2seσ+μsds>0,Ivt=e−γv+μtIv0+∫0t1−θσEvseγv+μsds>0.
It indicates that (E(t), Ev(t), I(t), Iv(t), S(t), V1t, V2t, R(t))∉∂X0, for t ∈ (0, t2), which is a contradiction. Hence, we verify M∂ = {(S, V1, V2, 0, 0, 0, 0, 0)}.
Obviously,E0 is one fixed point of P in M∂. For any x0 ϵ M∂\E0, limm⟶+∞Pm(x0) = E0, which means E0 is isolated in M∂ and then isolated in ∂X0. Hence, we get E0 is acyclic covering in M∂.
The following showsWs(E0)∩X0 = ∅. First, by the continuity of solutions with respect to initial values, for any ε > 0, there exists δ0 > 0; if ‖x0 − E0‖ ≤ δ0, it follows ‖u(t, x0) − u(t, E0)‖ ≤ ϵ, for all t ∈ [0, T]. In the following, we first claim limm⟶+∞sup d(Pm(x0), E0) ≥ δ0.
If not, there existsx0 ∈ X0, such that limm⟶+∞sup d(Pm(x0), E0) < δ0. Without loss of generality, assuming d(Pm(x0), E0) < δ0, ∀m > 0. For any t>0,t=mT+t^, then(15)ut,x0−ut,E0=ut^,Pmx0−ut^,E0<ε,∀t≥0.
Let (S(t), V1t, V2t, E(t), Ev(t), I(t), Iv(t), R(t)) = u(t, x0); then, we get S¯−ε<St<S¯+ε,V¯1−ε<V1t<V¯1+ε,V¯2−ε<V2t<V¯2+ε,0<Et,Evt,It,Ivt<ε.
Hence, from system (1), we obtain(16)dEdt>βtI+ξIvS¯−ε−σ+μE,dEvdt>ηβtI+ξIvV¯2−ε−σ+μEv.
Consider auxiliary system(17)dE^dt=βtI^+ξI^vS¯−ε−σ+μE¯,dE^vdt=ηβtI^+ξI^vV¯2−ε−σ+μE^v,dI^dt=σE^+θσE^v−γ+μI^,dI^vdt=1−θσE^v−γv+μI^v.
For convenience, it can be rewritten as(18)dudt=Ft−Vtu−εMtu,where u = (Ê(t), E^v(t), Î(t), I^v(t)) and(19)Mt=00βtεξβtε00ηβtεξηβtε00000000.
By Lemma 2.1 of [27], there exists a positive T-periodic function vt such that eμ2tvt is the solution of (18), where μ2=1/TlnρΦF⋅−V⋅−εM⋅T. Since R0 > 1, that is, ρ(ΦF(·)−V(·) (T)) > 1, we can choose ε small enough such that ρ(ΦF(·)−V(·)−εM(·) (T)) > 1, and thus μ2 > 0. Choose α>0,t¯>0 such that ut¯>αv0; then, ut>αvt−t¯eμ2t−t¯.
By standard comparison theorem [26], we get(20)Et,Evt,It,Ivt≥ut≥αvt−t¯eμ2t−t¯.
Hence,(21)limm⟶+∞Et=∞,limm⟶+∞Evt=∞,limm⟶+∞It=∞,limm⟶+∞Ivt=∞.
This is contradiction toE(t), Ev(t), I(t), Iv(t) < ɛ. Thus, we get Ws(E0)∩X0 = ϕ and E0 is isolated in X.
Above all, by Theorem 3.11 [29], we get P is uniformly persistent with respect to (X0, ∂X0).
Next, prove that forR0 > 1, there is a positive periodic solution of system (1). It follows from Theorem 1.3.6 of [29] that the Poincaré map P has a fixed point X∗=S∗0,V1∗0,V2∗0,E∗0,Ev∗0,I∗0,Iv∗0,R∗0∈X0. The following shows that S∗0>0,V1∗0>0,V2∗0>0,R∗>0. Suppose not assuming S∗(0) = 0, then the solution S˜t through initial point S∗(0) satisfies equation(22)dS˜tdt=Λ−pS˜t−βtI˜t+ξI˜vtS˜t−μS˜t.
We getS˜t=e−∫0tμ+p+βsI˜s+ξI˜vsdsS∗0+∫0tΛe∫0sμ+p+βτI˜τ+ξI˜vτdτds. Since X∗ is a fixed point of P, we can obtain S˜T=Λe−∫0Tμ+p+βsI˜s+ξI˜vsds∫0Te∫0sμ+p+βτI˜τ+ξI˜vτdτds=S∗0>0, which is contradiction with S∗(0) = 0. It also denotes that the fixed point X∗ is positive.
LetX˜t=S˜t,V˜1t,V˜2t,E˜t,E˜vt,I˜t,I˜vt,R˜t be a periodic solution through X∗. the uniform persistence of the solution with respect to (X0, ∂X0), i.e., limm⟶+∞infSt,V1t,V2t,Et,Evt,It,Ivt,Rt>η¯,η¯,η¯,η¯,η¯,η¯,η¯,η¯, we get ΩX˜t=X˜t>0,∀t≥0.


### 2.3. Numerical Simulation

In this section, we will numerically analyze model (1) and identify the effects of some parameters we are interested in on the outbreak of measles. According to reference [9], we obtain the recruitment rate Λ = 1589357 per month. Assuming the current average life expectancy of Chinese people is 73 years, then we calculate the natural death rate μ = 0.001142 per month. From [11], we derive that the exposed and infectious periods are about 13.8 and 10 days, which indicates σ = 2.17 month−1, γ = 3 month−1. We let the number of initial infectious individuals be the sum of the numbers of reported measles cases for last 10 days (one infectious period) in December of 2012, that is, I(0) = 300. By fitting model (1) to the monthly reported measles data from 2013 to 2016, we estimate some unknown parameters, which are listed in Table 1. Figure 2 shows a high goodness of fit (R2 = 0.8378) which demonstrates that the model can capture main trends of measles outbreak events between 2013 and 2016. Note that the recovery rate γv is estimated as 6 month−1, that is, 5 days, which means symptomatic patients will take twice the time to recover than asymptomatic infectious individuals. Furthermore, the estimated value of θ is small, demonstrating that most individuals in group Ev move into the class Iv. In other words, individuals in V2 class, once infected, mainly become asymptomatic. Using the estimated parameter values, we calculate the basic reproduction number R0 as 1.0077, which indicates measles will persist.Figure 2
Fitting to reported measles cases from 2013 to 2016. The initial values are (S(0), V1(0), V2(0), E(0), Ev(0), I(0), Iv(0), R(0)) = (1.1439 × 109, 9.9395 × 107, 3.23 × 107, 599, 1851, 300, 7316, 3.2652 × 107).It is worth noting that two routine immunizations are objected to children under 2 years old in mainland China; hence, the vaccinated individuals in our model actually belong to this age group. LetS1 be the number of susceptible children younger than 2 years old and pc be the effective coverage rate corresponding these children, while parameter p in the model (1) represents the coverage rate of vaccination for general population; then, we actually have pS = pcS1, and in the following simulations, we mainly use pc to investigate the impact of vaccination on disease outbreak and make discussion.To determine the influence of asymptomatic infection and other important parameters on the measles infection, we investigate how the basic reproduction number and the number of symptomatic infections vary with parameters. It follows from Figures3(a) and 4(a)–4(c) that reducing baseline transmission rate (a) or increasing coverage rate (pc) and recovery rate (γ) significantly decreases the basic reproduction number and consequently the number of symptomatic infections declines. This indicates limited contact, extensive vaccination, and timely treatment can effectively mitigate the outbreak and transmission of measles. Then, considering parameters θ, ξ, η relating to asymptomatic infection, we find that decreasing the values of θ, ξ, η has slight impact on reducing the new infections and controlling the transmission of measles except for minutely decreasing the peak size of symptomatic infections, as shown in Figures 3(b) and 4(d)–4(f). It indicates that asymptomatic measles infection in population does not cause much new infections.Figure 3
Change of the basic reproduction number with varying values of ratek. (a) We decrease baseline transmission rate (represented by (1 − k)a) and increase coverage rate of children under 2 years old and recovery rate (represented by (1 + k) pc, (1 + k)0, 0 ≤ k) ≤ 0.04). (b) We decrease parameters related to asymptomatic infection (represented by (1 − k)ξ, (1 − k)η, (1 − k)θ, 0 ≤ k) ≤ 0.9).
(a)
(b)Figure 4
Plots of the number of measles against time at varying parameter value: (a) the baseline transmission ratea; (b) recovery rate γ for the symptomatic infected; (c) coverage rate pc of children under 2 years old; (d) proportion θ of population from Ev to I; (e, f) correction factors ξ and η of asymptomatic infection.
(a)
(b)
(c)
(d)
(e)
(f)Further, in order to access the dependence ofR0 on parameter variations, we perform sensitivity analysis using the Latin hypercube sampling (LHS) method and calculate partial rank correlation coefficients (PRCCs) [30, 31] for various input parameters against the basic reproduction number. We choose uniform distribution in all parameters as listed in Table 2 since there is limited information on the distribution of each parameter. Figure 5 shows that, in contrast to parameters relative to asymptomatic infection, coverage rate pc, the baseline transmission rate a, and recovery rate γ are the first three parameters that mostly influence R0. These results indicate that reducing transmission rate via limiting contacts (e.g., personal protective and social distancing), treating infected individuals timely, and enhancing vaccination greatly lead to decline in new infections.Table 2
Basic reproduction number for PRCCs.
Parameter
Distribution
PRCC
p value
a
U(1.16 ∗ 10−9, 2 ∗ 10−9)
0.6142
0
γ
U(2.9, 4)
−0.4493
0
p
c
U(0.002, 0.9)
−0.6904
0
ξ
U(0.0005, 0.0015)
0.0195
0.7866
η
U(0.009, 0.015)
0.0528
0.4633
θ
U(0.004, 0.0055)
0.0604
0.4013Figure 5
PRCCs of the six parameters forR0 from 2013 to 2016. ∗ PRCCs are significant (p < 0.01).


## 2.1. Model Formulation

The underlying structure of the model comprises of eight compartments. We assume that susceptible individuals (S) receive the vaccine and enter into the compartment V1 at rate of p. Then, the vaccinated individuals (V1) progress to the class (V2) with waning of immunity at rate of ω. Susceptible individuals (S) become infected and enter into the latent undetectable class (E) at rate of β(t) (I + ξIv), where β(t) is the baseline periodic transmission rate and ξ (0 ≤ ξ ≤ 1) quantifies the relative transmissibility of asymptomatic infections. Then, the exposed individuals (E) progress to infectious class (I) at rate of σ and recover with rate of γ. Note that most individuals with waning immunity may not produce symptoms, once infected again; hence, it is negligible that latent individuals (E) move to asymptomatic infections class (Iv). The individuals in V1 class are those who are effectively vaccinated with high antibody titer, so we assume they cannot be infected as mentioned in [19, 20]. But because the antibody titer is low, individuals located in V2 may be infected and become exposed individuals Ev at rate of ηβ(t) (I + ξIv) with the modification factor η (0 ≤ η ≤ 1), denoting these individuals having low probability to be infected. Then, these latent individuals may move to symptomatic or asymptomatic infectious class (I) and (Iv) at rate of θσ or (1−θ)σ, where θ (0 ≤ θ ≤ 1) is the proportion of latent individuals becoming symptomatic infections. The recovery rate from Iv to R is γv. In general, asymptomatic infections recover faster than symptomatic infection, which means γv > γ. A flow diagram of the model is described in Figure 1 and the definitions of variables and parameters are given in Table 1. The specific dynamic model is as follows:(1)dSdt=Λ−pS−βtI+ξIvS−μS,dV1dt=pS−μ+ωV1,dV2dt=ωV1−ηβtI+ξIvV2−μV2,dEdt=βtI+ξIvS−σ+μE,dEvdt=ηβtI+ξIvV2−σ+μEv,dIdt=σE+θσEv−γ+μI,dIvdt=1−θσEv−γv+μIv,dRdt=γI+γvIv−μR,where  β(t) is nonnegative and periodic function with period  T (T > 0) and Λ and μ are recruitment rate and natural death rate per unit time. All parameters are assumed to positive constants.Figure 1
The flow diagram of the measles transmission.Table 1
Interpretation and values of the parameters used in the simulation.
Parameters
Description
Default value
Sources
μ
Natural death rate (month−1)
1.142 × 10−3
Calculated
σ
Rate of progression to infected stage from the exposed (month−1)
2.17
[11]
γ
Recovery rate for the infected stage (month−1)
3
[11]
Λ
Recruitment rate (month−1)
1589357
[9]
ω
Rate of waning immunity (month−1)
0.0033
[19]
γ
v
Recovery rate for the asymptomatic infected stage (month−1)
6
LS
η
Correction factor
0.01
LS
ξ
Correction factor
0.001
LS
θ
Proportion of progression to symptomatic infected stage
0.004
LS
p
Coverage rate of vaccine for entire population (month−1)
0.00014
LS
β(t)
Periodic transmission rate of symptomatic infection (month−1)
a
1.346
+
sin
2
π
t
/
12
+
2
Assumed
a
Baseline transmission rate (month−1)
1.9616 × 10−9
LS


## 2.2. Threshold Dynamics

It is clear that any solution of this system from nonnegative initial values is nonnegative. Then, we first discuss bound of solution of this system.Proposition 1.
The solution of system (1) is uniformly and ultimately bounded, i.e., there exists M>0,t⌣>0 such that (S, V1, V2, E, Ev, I, Iv, R) ≤ (M, M, M, M, M, M, M, M) for t>t⌣.Proof.
LetS(t) + V1(t) + V2(t) + E(t) + Ev(t) + I(t) + Iv(t) +R(t) = N; then,(2)N′=Λ−μN.
Hence,limt⟶+∞Nt=Λ/μ, which implies for any ε > 0, there exists t⌣>0 for t>t⌣, then Nt<Λ/μ+ε=˙M. Subsequently, We get (S, V1, V2, E, Ev, I, Iv, R) ≤ (M, M, M, M, M, M, M, M), i.e., the solution of system (1) is uniformly and ultimately bounded.


### 2.2.1. Disease-Free Equilibrium

It is easy to obtain that system (1) has unique disease-free equilibrium E0=S¯,V¯1,V¯2,0,0,0,0,0, where S¯=Λ/p+μ,V¯1=Λp/p+μω+μ,V¯2=Λpω/p+μω+μμ. In the following, the stability of disease-free equilibrium will be discussed. First, we define the reproduction number of system (1) on the basis of [25].Denote(3)Ft=00βtS¯ξβtS¯00ηβtV¯2ξηβtV¯200000000,Vt=σ+μ0000σ+μ00−σ−θσγ+μ00−1−θσ0γv+μ.AssumeY(t, s), t ≥ s is the evolution operator of following linear system:(4)dydt=−Vty.That is to say,Y(t, s) satisfies(5)dYt,sdt=−VtYt,s,∀t≥s,Ys,s=I.Let Φ−V(t) be the fundamental solution matrix of linear system (4). Clearly, Φ−V(t) equals to Y(t, 0), t ≥ 0.Define a linear operatorL : CT ⟶ CT by(6)Lϕt=∫−∞tYt,sFsϕsds=∫0∞Yt,t−aFt−aϕt−ada,where CT  is the ordered Banach space of all T-periodic function from R to R4 with maximum norm ‖·‖ and positive cone CT+:=ϕ∈CT:ϕt≥0,∀t∈R. Suppose that ϕ ∈ CT is the initial distribution of infectious individuals in the periodic environment; then, F(s)ϕ(s) is the distribution of new infections produced by the infected individuals who were introduced at time  s. Given t (t ≥ s), then Y(t, s)F(s)ϕ(s) gives the distribution of those infected individuals who were newly infected at time s and remain in the infected compartments. Define the spectral radius of L as the basic reproduction number of system (1), i.e.,(7)R0=ρL.It can be verified that system (1) satisfies conditions (A1)–(A7) of periodic system [25]; then, the following theorem is established.Lemma 1.
Assume that (A1)–(A7) hold [25]. Then, the following statements are valid:(1)
R0 = 1 if and only if ρ(ΦF−V(T)) = 1(2)
R0 < 1 if and only if ρ(ΦF−V(T)) < 1(3)
R0 > 1 if and only if ρ(ΦF−V(T)) > 1Theorem 1.
IfR0 < 1, the disease-free equilibrium E0 is globally asymptotically stable; if R0 > 1, the disease-free equilibrium E0 is unstable.Proof.
According to Lemma1, it is easy to get E0 is local stable if R0 < 1. The following will show the global attractive.
By system (1), we obtain(8)dSdt≤Λ−pS−μS,dV1dt≤pS−ω+μV1,dV2dt≤ωV1−μV2.
By standard comparison theorem [26], we get for any ε > 0, there exists t1 > 0; if t > t1, then St<S¯+ε,V1t<V¯1+ε,V2t<V¯2+ε.
Consider auxiliary system(9)dE^dt=βtI^+ξI^vS¯+ε−σ+μE^,dE^vdt=ηβtI^+ξI^vV¯2+ε−σ+μEv^,dI^dt=σE^+θσE^v−γ+μI^,dI^vdt=1−θσE^v−γv+μI^v.
We can write (9) as(10)dXdt=Ft−Vt+εMtX,where X = (Ê, E^v, Î, I^v), and(11)Mt=00βtεξβtε00ηβtεξηβtε00000000.
SinceR0 < 1, that is ρ(ΦF(t)−V(t) (T)) < 1, then for any small ε, ρ(ΦF(t)−V(t)+εM(t) (T)) < 1 and consequently, μ1=1/TlnρΦFt−Vt+εMtT<0.
It follows from Lemma 2.1 of [27] that there exists a positive T-periodic function vt = (v1t, v2t, v3t, v4t) such that eμ1tvt is a solution of (9). Choose t¯>t1 and a small value α > 0 such that Xt¯≤αv0; then, we can get Xt≤αvt−t¯eμ1t−t¯,t≥t¯.
By standard comparison theorem [26], we get(12)Et,Evt,It,Ivt≤Xt≤αvt−t¯eμ1t−t¯.
Hence,limt⟶+∞Et=0,limt⟶+∞Evt=0,limt⟶+∞It=0,limt⟶+∞Ivt=0. By the last four equations of system (1), considering their limit system, i.e.,(13)dSdt=Λ−pS−μS,dV1dt=pS−ω+μV1,dV2dt=ωV1−μV2,dRdt=−μR.
We getlimt⟶+∞St=S¯,limt⟶+∞V1t=V¯1,limt⟶+∞V2t=V¯2,limt⟶+∞Rt=0. The global attractivity of E0 has been proved.


### 2.2.2. Uniform Persistence of the System

In the following, we examine the uniform persistence and existence of positive periodic solutions of system (1).Theorem 2.
IfR0 > 1, system (1) is uniformly persistent, i.e., there is constant η¯>0, such that for all initial values S0,V10,V20,E0,Ev0,I0,Iv0,R0∈R+3×IntR+4×R+, the solution of system (1) satisfies limt⟶+∞infSt,V1t,V2t,Et,Evt,It,Ivt,Rt≥η¯,η¯,η¯,η¯,η¯,η¯,η¯,η¯. And there exists at least one positive T-periodic solution for system (1).Proof.
DefineX = {(S, V1, V2, E, Ev, I, Iv, R) : S, V1, V2, E, Ev, I, Iv, R ≥ 0}, X0=S,V1,V2,E,Ev,I,Iv,R:S,V1,V2,R≥0,E,Ev,I,Iv>0, ∂X0 = X\X0. Define a Poincaré map P:R+8⟶R+8 satisfying P(x0) = u(T, x0), for x0∈R+8, where u(t, x0) is the solution of system (1).
First, it is easy to see thatX, X0 are positively invariant and ∂X0 is relatively closed in X. It follows from the Proposition 1 that the solution of system (1) is uniformly and ultimately bounded. Thus, the semiflow P is point dissipative on R+8 and P:R+8⟶R+8 is compact. By Theorem 3.4.8 of [28], P admits a global attractor A.
DefineM∂ = {(S, V1, V2, E, Ev, I, Iv, R) ϵ ∂X0 : Pm(S, V1, V2, E, Ev, I, Iv, R) ∈ ∂X0, ∀m ≥ 0}. In the following, we claim M∂ = {(S, V1, V2, 0, 0, 0, 0, 0)}. First, it is obvious that {(S, V1, V2, 0, 0, 0, 0, 0)}⊂M∂. Next, certify M∂⊂{(S, V1, V2, 0, 0, 0, 0, 0)}. Assume there exists x0 ∈ M∂\{(S, V1, V2, 0, 0, 0, 0, 0)}. For example, let E0>0,Ev0=I0=Iv0=0; then, I′(0) = σE0 > 0, which implies there exists t2 > 0; when 0 < t < t2, we get I(t) > 0 and E(t) > 0. For 0 < t < t2, we obtain from equations of (1)(14)Evt=e−σ+μtEv0+∫0tηβsI+ξIvsV2seσ+μsds>0,Ivt=e−γv+μtIv0+∫0t1−θσEvseγv+μsds>0.
It indicates that (E(t), Ev(t), I(t), Iv(t), S(t), V1t, V2t, R(t))∉∂X0, for t ∈ (0, t2), which is a contradiction. Hence, we verify M∂ = {(S, V1, V2, 0, 0, 0, 0, 0)}.
Obviously,E0 is one fixed point of P in M∂. For any x0 ϵ M∂\E0, limm⟶+∞Pm(x0) = E0, which means E0 is isolated in M∂ and then isolated in ∂X0. Hence, we get E0 is acyclic covering in M∂.
The following showsWs(E0)∩X0 = ∅. First, by the continuity of solutions with respect to initial values, for any ε > 0, there exists δ0 > 0; if ‖x0 − E0‖ ≤ δ0, it follows ‖u(t, x0) − u(t, E0)‖ ≤ ϵ, for all t ∈ [0, T]. In the following, we first claim limm⟶+∞sup d(Pm(x0), E0) ≥ δ0.
If not, there existsx0 ∈ X0, such that limm⟶+∞sup d(Pm(x0), E0) < δ0. Without loss of generality, assuming d(Pm(x0), E0) < δ0, ∀m > 0. For any t>0,t=mT+t^, then(15)ut,x0−ut,E0=ut^,Pmx0−ut^,E0<ε,∀t≥0.
Let (S(t), V1t, V2t, E(t), Ev(t), I(t), Iv(t), R(t)) = u(t, x0); then, we get S¯−ε<St<S¯+ε,V¯1−ε<V1t<V¯1+ε,V¯2−ε<V2t<V¯2+ε,0<Et,Evt,It,Ivt<ε.
Hence, from system (1), we obtain(16)dEdt>βtI+ξIvS¯−ε−σ+μE,dEvdt>ηβtI+ξIvV¯2−ε−σ+μEv.
Consider auxiliary system(17)dE^dt=βtI^+ξI^vS¯−ε−σ+μE¯,dE^vdt=ηβtI^+ξI^vV¯2−ε−σ+μE^v,dI^dt=σE^+θσE^v−γ+μI^,dI^vdt=1−θσE^v−γv+μI^v.
For convenience, it can be rewritten as(18)dudt=Ft−Vtu−εMtu,where u = (Ê(t), E^v(t), Î(t), I^v(t)) and(19)Mt=00βtεξβtε00ηβtεξηβtε00000000.
By Lemma 2.1 of [27], there exists a positive T-periodic function vt such that eμ2tvt is the solution of (18), where μ2=1/TlnρΦF⋅−V⋅−εM⋅T. Since R0 > 1, that is, ρ(ΦF(·)−V(·) (T)) > 1, we can choose ε small enough such that ρ(ΦF(·)−V(·)−εM(·) (T)) > 1, and thus μ2 > 0. Choose α>0,t¯>0 such that ut¯>αv0; then, ut>αvt−t¯eμ2t−t¯.
By standard comparison theorem [26], we get(20)Et,Evt,It,Ivt≥ut≥αvt−t¯eμ2t−t¯.
Hence,(21)limm⟶+∞Et=∞,limm⟶+∞Evt=∞,limm⟶+∞It=∞,limm⟶+∞Ivt=∞.
This is contradiction toE(t), Ev(t), I(t), Iv(t) < ɛ. Thus, we get Ws(E0)∩X0 = ϕ and E0 is isolated in X.
Above all, by Theorem 3.11 [29], we get P is uniformly persistent with respect to (X0, ∂X0).
Next, prove that forR0 > 1, there is a positive periodic solution of system (1). It follows from Theorem 1.3.6 of [29] that the Poincaré map P has a fixed point X∗=S∗0,V1∗0,V2∗0,E∗0,Ev∗0,I∗0,Iv∗0,R∗0∈X0. The following shows that S∗0>0,V1∗0>0,V2∗0>0,R∗>0. Suppose not assuming S∗(0) = 0, then the solution S˜t through initial point S∗(0) satisfies equation(22)dS˜tdt=Λ−pS˜t−βtI˜t+ξI˜vtS˜t−μS˜t.
We getS˜t=e−∫0tμ+p+βsI˜s+ξI˜vsdsS∗0+∫0tΛe∫0sμ+p+βτI˜τ+ξI˜vτdτds. Since X∗ is a fixed point of P, we can obtain S˜T=Λe−∫0Tμ+p+βsI˜s+ξI˜vsds∫0Te∫0sμ+p+βτI˜τ+ξI˜vτdτds=S∗0>0, which is contradiction with S∗(0) = 0. It also denotes that the fixed point X∗ is positive.
LetX˜t=S˜t,V˜1t,V˜2t,E˜t,E˜vt,I˜t,I˜vt,R˜t be a periodic solution through X∗. the uniform persistence of the solution with respect to (X0, ∂X0), i.e., limm⟶+∞infSt,V1t,V2t,Et,Evt,It,Ivt,Rt>η¯,η¯,η¯,η¯,η¯,η¯,η¯,η¯, we get ΩX˜t=X˜t>0,∀t≥0.


## 2.2.1. Disease-Free Equilibrium

It is easy to obtain that system (1) has unique disease-free equilibrium E0=S¯,V¯1,V¯2,0,0,0,0,0, where S¯=Λ/p+μ,V¯1=Λp/p+μω+μ,V¯2=Λpω/p+μω+μμ. In the following, the stability of disease-free equilibrium will be discussed. First, we define the reproduction number of system (1) on the basis of [25].Denote(3)Ft=00βtS¯ξβtS¯00ηβtV¯2ξηβtV¯200000000,Vt=σ+μ0000σ+μ00−σ−θσγ+μ00−1−θσ0γv+μ.AssumeY(t, s), t ≥ s is the evolution operator of following linear system:(4)dydt=−Vty.That is to say,Y(t, s) satisfies(5)dYt,sdt=−VtYt,s,∀t≥s,Ys,s=I.Let Φ−V(t) be the fundamental solution matrix of linear system (4). Clearly, Φ−V(t) equals to Y(t, 0), t ≥ 0.Define a linear operatorL : CT ⟶ CT by(6)Lϕt=∫−∞tYt,sFsϕsds=∫0∞Yt,t−aFt−aϕt−ada,where CT  is the ordered Banach space of all T-periodic function from R to R4 with maximum norm ‖·‖ and positive cone CT+:=ϕ∈CT:ϕt≥0,∀t∈R. Suppose that ϕ ∈ CT is the initial distribution of infectious individuals in the periodic environment; then, F(s)ϕ(s) is the distribution of new infections produced by the infected individuals who were introduced at time  s. Given t (t ≥ s), then Y(t, s)F(s)ϕ(s) gives the distribution of those infected individuals who were newly infected at time s and remain in the infected compartments. Define the spectral radius of L as the basic reproduction number of system (1), i.e.,(7)R0=ρL.It can be verified that system (1) satisfies conditions (A1)–(A7) of periodic system [25]; then, the following theorem is established.Lemma 1.
Assume that (A1)–(A7) hold [25]. Then, the following statements are valid:(1)
R0 = 1 if and only if ρ(ΦF−V(T)) = 1(2)
R0 < 1 if and only if ρ(ΦF−V(T)) < 1(3)
R0 > 1 if and only if ρ(ΦF−V(T)) > 1Theorem 1.
IfR0 < 1, the disease-free equilibrium E0 is globally asymptotically stable; if R0 > 1, the disease-free equilibrium E0 is unstable.Proof.
According to Lemma1, it is easy to get E0 is local stable if R0 < 1. The following will show the global attractive.
By system (1), we obtain(8)dSdt≤Λ−pS−μS,dV1dt≤pS−ω+μV1,dV2dt≤ωV1−μV2.
By standard comparison theorem [26], we get for any ε > 0, there exists t1 > 0; if t > t1, then St<S¯+ε,V1t<V¯1+ε,V2t<V¯2+ε.
Consider auxiliary system(9)dE^dt=βtI^+ξI^vS¯+ε−σ+μE^,dE^vdt=ηβtI^+ξI^vV¯2+ε−σ+μEv^,dI^dt=σE^+θσE^v−γ+μI^,dI^vdt=1−θσE^v−γv+μI^v.
We can write (9) as(10)dXdt=Ft−Vt+εMtX,where X = (Ê, E^v, Î, I^v), and(11)Mt=00βtεξβtε00ηβtεξηβtε00000000.
SinceR0 < 1, that is ρ(ΦF(t)−V(t) (T)) < 1, then for any small ε, ρ(ΦF(t)−V(t)+εM(t) (T)) < 1 and consequently, μ1=1/TlnρΦFt−Vt+εMtT<0.
It follows from Lemma 2.1 of [27] that there exists a positive T-periodic function vt = (v1t, v2t, v3t, v4t) such that eμ1tvt is a solution of (9). Choose t¯>t1 and a small value α > 0 such that Xt¯≤αv0; then, we can get Xt≤αvt−t¯eμ1t−t¯,t≥t¯.
By standard comparison theorem [26], we get(12)Et,Evt,It,Ivt≤Xt≤αvt−t¯eμ1t−t¯.
Hence,limt⟶+∞Et=0,limt⟶+∞Evt=0,limt⟶+∞It=0,limt⟶+∞Ivt=0. By the last four equations of system (1), considering their limit system, i.e.,(13)dSdt=Λ−pS−μS,dV1dt=pS−ω+μV1,dV2dt=ωV1−μV2,dRdt=−μR.
We getlimt⟶+∞St=S¯,limt⟶+∞V1t=V¯1,limt⟶+∞V2t=V¯2,limt⟶+∞Rt=0. The global attractivity of E0 has been proved.


## 2.2.2. Uniform Persistence of the System

In the following, we examine the uniform persistence and existence of positive periodic solutions of system (1).Theorem 2.
IfR0 > 1, system (1) is uniformly persistent, i.e., there is constant η¯>0, such that for all initial values S0,V10,V20,E0,Ev0,I0,Iv0,R0∈R+3×IntR+4×R+, the solution of system (1) satisfies limt⟶+∞infSt,V1t,V2t,Et,Evt,It,Ivt,Rt≥η¯,η¯,η¯,η¯,η¯,η¯,η¯,η¯. And there exists at least one positive T-periodic solution for system (1).Proof.
DefineX = {(S, V1, V2, E, Ev, I, Iv, R) : S, V1, V2, E, Ev, I, Iv, R ≥ 0}, X0=S,V1,V2,E,Ev,I,Iv,R:S,V1,V2,R≥0,E,Ev,I,Iv>0, ∂X0 = X\X0. Define a Poincaré map P:R+8⟶R+8 satisfying P(x0) = u(T, x0), for x0∈R+8, where u(t, x0) is the solution of system (1).
First, it is easy to see thatX, X0 are positively invariant and ∂X0 is relatively closed in X. It follows from the Proposition 1 that the solution of system (1) is uniformly and ultimately bounded. Thus, the semiflow P is point dissipative on R+8 and P:R+8⟶R+8 is compact. By Theorem 3.4.8 of [28], P admits a global attractor A.
DefineM∂ = {(S, V1, V2, E, Ev, I, Iv, R) ϵ ∂X0 : Pm(S, V1, V2, E, Ev, I, Iv, R) ∈ ∂X0, ∀m ≥ 0}. In the following, we claim M∂ = {(S, V1, V2, 0, 0, 0, 0, 0)}. First, it is obvious that {(S, V1, V2, 0, 0, 0, 0, 0)}⊂M∂. Next, certify M∂⊂{(S, V1, V2, 0, 0, 0, 0, 0)}. Assume there exists x0 ∈ M∂\{(S, V1, V2, 0, 0, 0, 0, 0)}. For example, let E0>0,Ev0=I0=Iv0=0; then, I′(0) = σE0 > 0, which implies there exists t2 > 0; when 0 < t < t2, we get I(t) > 0 and E(t) > 0. For 0 < t < t2, we obtain from equations of (1)(14)Evt=e−σ+μtEv0+∫0tηβsI+ξIvsV2seσ+μsds>0,Ivt=e−γv+μtIv0+∫0t1−θσEvseγv+μsds>0.
It indicates that (E(t), Ev(t), I(t), Iv(t), S(t), V1t, V2t, R(t))∉∂X0, for t ∈ (0, t2), which is a contradiction. Hence, we verify M∂ = {(S, V1, V2, 0, 0, 0, 0, 0)}.
Obviously,E0 is one fixed point of P in M∂. For any x0 ϵ M∂\E0, limm⟶+∞Pm(x0) = E0, which means E0 is isolated in M∂ and then isolated in ∂X0. Hence, we get E0 is acyclic covering in M∂.
The following showsWs(E0)∩X0 = ∅. First, by the continuity of solutions with respect to initial values, for any ε > 0, there exists δ0 > 0; if ‖x0 − E0‖ ≤ δ0, it follows ‖u(t, x0) − u(t, E0)‖ ≤ ϵ, for all t ∈ [0, T]. In the following, we first claim limm⟶+∞sup d(Pm(x0), E0) ≥ δ0.
If not, there existsx0 ∈ X0, such that limm⟶+∞sup d(Pm(x0), E0) < δ0. Without loss of generality, assuming d(Pm(x0), E0) < δ0, ∀m > 0. For any t>0,t=mT+t^, then(15)ut,x0−ut,E0=ut^,Pmx0−ut^,E0<ε,∀t≥0.
Let (S(t), V1t, V2t, E(t), Ev(t), I(t), Iv(t), R(t)) = u(t, x0); then, we get S¯−ε<St<S¯+ε,V¯1−ε<V1t<V¯1+ε,V¯2−ε<V2t<V¯2+ε,0<Et,Evt,It,Ivt<ε.
Hence, from system (1), we obtain(16)dEdt>βtI+ξIvS¯−ε−σ+μE,dEvdt>ηβtI+ξIvV¯2−ε−σ+μEv.
Consider auxiliary system(17)dE^dt=βtI^+ξI^vS¯−ε−σ+μE¯,dE^vdt=ηβtI^+ξI^vV¯2−ε−σ+μE^v,dI^dt=σE^+θσE^v−γ+μI^,dI^vdt=1−θσE^v−γv+μI^v.
For convenience, it can be rewritten as(18)dudt=Ft−Vtu−εMtu,where u = (Ê(t), E^v(t), Î(t), I^v(t)) and(19)Mt=00βtεξβtε00ηβtεξηβtε00000000.
By Lemma 2.1 of [27], there exists a positive T-periodic function vt such that eμ2tvt is the solution of (18), where μ2=1/TlnρΦF⋅−V⋅−εM⋅T. Since R0 > 1, that is, ρ(ΦF(·)−V(·) (T)) > 1, we can choose ε small enough such that ρ(ΦF(·)−V(·)−εM(·) (T)) > 1, and thus μ2 > 0. Choose α>0,t¯>0 such that ut¯>αv0; then, ut>αvt−t¯eμ2t−t¯.
By standard comparison theorem [26], we get(20)Et,Evt,It,Ivt≥ut≥αvt−t¯eμ2t−t¯.
Hence,(21)limm⟶+∞Et=∞,limm⟶+∞Evt=∞,limm⟶+∞It=∞,limm⟶+∞Ivt=∞.
This is contradiction toE(t), Ev(t), I(t), Iv(t) < ɛ. Thus, we get Ws(E0)∩X0 = ϕ and E0 is isolated in X.
Above all, by Theorem 3.11 [29], we get P is uniformly persistent with respect to (X0, ∂X0).
Next, prove that forR0 > 1, there is a positive periodic solution of system (1). It follows from Theorem 1.3.6 of [29] that the Poincaré map P has a fixed point X∗=S∗0,V1∗0,V2∗0,E∗0,Ev∗0,I∗0,Iv∗0,R∗0∈X0. The following shows that S∗0>0,V1∗0>0,V2∗0>0,R∗>0. Suppose not assuming S∗(0) = 0, then the solution S˜t through initial point S∗(0) satisfies equation(22)dS˜tdt=Λ−pS˜t−βtI˜t+ξI˜vtS˜t−μS˜t.
We getS˜t=e−∫0tμ+p+βsI˜s+ξI˜vsdsS∗0+∫0tΛe∫0sμ+p+βτI˜τ+ξI˜vτdτds. Since X∗ is a fixed point of P, we can obtain S˜T=Λe−∫0Tμ+p+βsI˜s+ξI˜vsds∫0Te∫0sμ+p+βτI˜τ+ξI˜vτdτds=S∗0>0, which is contradiction with S∗(0) = 0. It also denotes that the fixed point X∗ is positive.
LetX˜t=S˜t,V˜1t,V˜2t,E˜t,E˜vt,I˜t,I˜vt,R˜t be a periodic solution through X∗. the uniform persistence of the solution with respect to (X0, ∂X0), i.e., limm⟶+∞infSt,V1t,V2t,Et,Evt,It,Ivt,Rt>η¯,η¯,η¯,η¯,η¯,η¯,η¯,η¯, we get ΩX˜t=X˜t>0,∀t≥0.


## 2.3. Numerical Simulation

In this section, we will numerically analyze model (1) and identify the effects of some parameters we are interested in on the outbreak of measles. According to reference [9], we obtain the recruitment rate Λ = 1589357 per month. Assuming the current average life expectancy of Chinese people is 73 years, then we calculate the natural death rate μ = 0.001142 per month. From [11], we derive that the exposed and infectious periods are about 13.8 and 10 days, which indicates σ = 2.17 month−1, γ = 3 month−1. We let the number of initial infectious individuals be the sum of the numbers of reported measles cases for last 10 days (one infectious period) in December of 2012, that is, I(0) = 300. By fitting model (1) to the monthly reported measles data from 2013 to 2016, we estimate some unknown parameters, which are listed in Table 1. Figure 2 shows a high goodness of fit (R2 = 0.8378) which demonstrates that the model can capture main trends of measles outbreak events between 2013 and 2016. Note that the recovery rate γv is estimated as 6 month−1, that is, 5 days, which means symptomatic patients will take twice the time to recover than asymptomatic infectious individuals. Furthermore, the estimated value of θ is small, demonstrating that most individuals in group Ev move into the class Iv. In other words, individuals in V2 class, once infected, mainly become asymptomatic. Using the estimated parameter values, we calculate the basic reproduction number R0 as 1.0077, which indicates measles will persist.Figure 2
Fitting to reported measles cases from 2013 to 2016. The initial values are (S(0), V1(0), V2(0), E(0), Ev(0), I(0), Iv(0), R(0)) = (1.1439 × 109, 9.9395 × 107, 3.23 × 107, 599, 1851, 300, 7316, 3.2652 × 107).It is worth noting that two routine immunizations are objected to children under 2 years old in mainland China; hence, the vaccinated individuals in our model actually belong to this age group. LetS1 be the number of susceptible children younger than 2 years old and pc be the effective coverage rate corresponding these children, while parameter p in the model (1) represents the coverage rate of vaccination for general population; then, we actually have pS = pcS1, and in the following simulations, we mainly use pc to investigate the impact of vaccination on disease outbreak and make discussion.To determine the influence of asymptomatic infection and other important parameters on the measles infection, we investigate how the basic reproduction number and the number of symptomatic infections vary with parameters. It follows from Figures3(a) and 4(a)–4(c) that reducing baseline transmission rate (a) or increasing coverage rate (pc) and recovery rate (γ) significantly decreases the basic reproduction number and consequently the number of symptomatic infections declines. This indicates limited contact, extensive vaccination, and timely treatment can effectively mitigate the outbreak and transmission of measles. Then, considering parameters θ, ξ, η relating to asymptomatic infection, we find that decreasing the values of θ, ξ, η has slight impact on reducing the new infections and controlling the transmission of measles except for minutely decreasing the peak size of symptomatic infections, as shown in Figures 3(b) and 4(d)–4(f). It indicates that asymptomatic measles infection in population does not cause much new infections.Figure 3
Change of the basic reproduction number with varying values of ratek. (a) We decrease baseline transmission rate (represented by (1 − k)a) and increase coverage rate of children under 2 years old and recovery rate (represented by (1 + k) pc, (1 + k)0, 0 ≤ k) ≤ 0.04). (b) We decrease parameters related to asymptomatic infection (represented by (1 − k)ξ, (1 − k)η, (1 − k)θ, 0 ≤ k) ≤ 0.9).
(a)
(b)Figure 4
Plots of the number of measles against time at varying parameter value: (a) the baseline transmission ratea; (b) recovery rate γ for the symptomatic infected; (c) coverage rate pc of children under 2 years old; (d) proportion θ of population from Ev to I; (e, f) correction factors ξ and η of asymptomatic infection.
(a)
(b)
(c)
(d)
(e)
(f)Further, in order to access the dependence ofR0 on parameter variations, we perform sensitivity analysis using the Latin hypercube sampling (LHS) method and calculate partial rank correlation coefficients (PRCCs) [30, 31] for various input parameters against the basic reproduction number. We choose uniform distribution in all parameters as listed in Table 2 since there is limited information on the distribution of each parameter. Figure 5 shows that, in contrast to parameters relative to asymptomatic infection, coverage rate pc, the baseline transmission rate a, and recovery rate γ are the first three parameters that mostly influence R0. These results indicate that reducing transmission rate via limiting contacts (e.g., personal protective and social distancing), treating infected individuals timely, and enhancing vaccination greatly lead to decline in new infections.Table 2
Basic reproduction number for PRCCs.
Parameter
Distribution
PRCC
p value
a
U(1.16 ∗ 10−9, 2 ∗ 10−9)
0.6142
0
γ
U(2.9, 4)
−0.4493
0
p
c
U(0.002, 0.9)
−0.6904
0
ξ
U(0.0005, 0.0015)
0.0195
0.7866
η
U(0.009, 0.015)
0.0528
0.4633
θ
U(0.004, 0.0055)
0.0604
0.4013Figure 5
PRCCs of the six parameters forR0 from 2013 to 2016. ∗ PRCCs are significant (p < 0.01).


## 3. Optimal Control Strategies

In this section, we extend the basic model (1) by including some feasible controls and formulate the optimal control system, which is defined by the following model:(23)dSdt=Λ−1+u3tpS−1−u1tβtI+ξIvS−μS,dV1dt=1+u3tpS−μ+ωV1,dV2dt=ωV1−η1−u1tβtI+ξIvV2−μV2,dEdt=1−u1tβtI+ξIvS−σ+μE,dEvdt=η1−u1tβtI+ξIvV2−σ+μEv,dIdt=σE+θσEv−1+u2tγI−μI,dIvdt=1−θσEv−γv+μIv,dRdt=1+u2tγI+γvIv−μR.The control functionu1(t) represents reducing contact between susceptible and infective individuals on the basis of personal protection (e.g., wearing masks and enhancing awareness about measles) or social distancing (e.g., Fengxiao) measures. The control functions u2(t) and u3(t) represent the enhanced treatment and vaccination, respectively.Our goal is to minimize the number of infectious individuals over a time interval [0,T] at a minimal cost. Then, the objective function J is defined as(24)Ju1,u2,u3=∫0TIt+12W1u12t+W2u22t+W3u32tdt,where W1, W2, and W3 are the weight constants for the control functions. The optimal control strategies can be obtained by finding an optimal pair (U∗, Y∗) such that(25)JU∗=minΩJUt,where Ω = {U(t) = (u1, u2, u3)ε L3(0, T)3|0 ≤ u1(t) ≤ b1, 0 ≤ u2(t) ≤ b2, 0 ≤ u3(t) ≤ b3} with positive constants bi (i = 1, 2, 3) and Y∗=S∗,V1∗,V2∗,E∗,Ev∗,I∗,Iv∗,R∗.


### 3.1. Optimality System

The existence of optimal controls follows from Corollary 4.1 of [32] since the state system satisfies the Lipschitz property with respect to the state variables and is a linear function of U, while the integrand of J is a convex function for u1(t), u2(t) and u3(t).Theorem 3.
There exist optimal controlsU∗(t) and corresponding solutions Y∗(t) that minimize J(U(t)) over Ω. In order for the above statement to be true, it is necessary that there exist continuous functions λi(t) such that(26)λ1′t=λ1−λ41−u1tβI+ξIv+λ1−λ21+u3tp1+μλ1,λ2′t=λ2−λ3ω+μλ2,λ3′t=λ3−λ51−u1tηβI+ξIv+μλ3,λ4′t=λ4−λ6σ+μλ4,λ5′t=λ7−λ6θσ+λ5−λ7σ+μλ5,λ6′t=−1+λ1−λ41−u1tβS+λ3−λ5η1−u1tβV2+λ6−λ81+u2tγ+μλ6,λ7′t=λ1−λ41−u1tξβS+λ3−λ51−u1tξηβV2+λ7−λ8γv+μλ7,λ8′t=μλ8,with the transversality conditions(27)λ1T=λ2T=λ3T=λ4T=λ5T=λ6T=λ7T=λ8T=0.Furthermore, the optimal control is given as follows:(28)u1∗t=maxminλ4−λ1βI∗+ξIv∗S∗+λ5−λ3ηβI∗+ξIv∗V2∗W1,b1,0,u2∗t=maxminλ6−λ8γI∗W2,b2,0,u3∗t=maxminλ1−λ2pS∗W3,b3,0.Proof.
The necessary conditions satisfied by optimal controls can be derived from Pontryagin’s maximum principles. Let Hamilton function be(29)H=I+12W1u12+12W2u22+W3u32+λ1dSdt+λ2dV1dt+λ3dV2dt+λ4dEdt+λ5dEvdt+λ6dIdt+λ7dIvdt+λ8dRdt.
Then, the adjoint equations with transversality satisfyλ1′=−∂H/∂S,λ2′=−∂H/∂Vk,λ3′=−∂H/∂V2,λ4′=−∂H/∂E,λ5′=−∂H/∂Ev,λ6′=−∂H/∂I,λ7′=−∂H/∂Iv,λ8′=−∂H/∂R, with  λ1(T) = λ2(T) = λ3(T) = λ4(T) = λ5(T) = λ6(T) = λ7(T) = λ8(T) = 0.
The Hamilton functionH is minimized with respect to the controls by differentiating H with respect to u1, u2, u3 on the set Ω, respectively, and then gives the following optimal conditions:(30)∂H∂u1=0,∂H∂u2=0,∂H∂u3=0.
Solving foru1∗,u2∗,u3∗ yields(31)u1∗=λ4−λ1βI∗+ξIv∗S∗+λ5−λ3ηβI∗+ξIv∗V2∗W1,u2∗=λ6−λ8γI∗W2,u3∗=λ1−λ2p1S∗W3.
By using the bounds 0 ≤u1(t) ≤ b1, 0 ≤ u2 ≤ b2, 0 ≤ u3 ≤ b3, we have the properties as Theorem 3.


### 3.2. Numerical Simulations

As mentioned in Section2.3, based on relationship between pc and p, we initially define the control function u3c(t) as the enhanced vaccination rate for children under 2 years old to replace u3(t) for general population and numerically explore the effects of optimal control measures on measles outbreak. The optimality system is solved by using the forward-backward sweep method. With an initial guess for the control variables, the state system and adjoint system are solved forward and backward, separately. The control variables are updated by using optimality condition, and this process is repeated until the system is converged. In the numerical simulations, we use parameter values given in Table 1. To balance the population and control functions in the objective functions, we choose weight constants values W1 = 20, W2 = 200, and W3 = 0.2 and the upper bound of controls b1 = 0.2, b2 = 0.2, and b3 = 5 for illustration.To study the effect of optimal control strategies on the transmission of measles infection, we initially plot the time series of the number of symptomatic infected individuals and corresponding optimal functions, as shown in Figure6. It follows from Figure 6(a) that the optimal control strategies significantly mitigate the outbreak of disease and decrease the number of symptomatic infections. Figures 6(b)–6(d) give the optimal control measures u1∗t,u2∗t and u3c∗t, which monotonically decrease after a period time of implementation with high intensity. This indicates we should enhance vaccination as soon as possible when measles begin outbreak, while measures such as personal protection and treatment should be implemented as much as possible throughout the whole epidemic.Figure 6
(a) The number of infected individuals with symptoms under the two situations: no control is used and three types of controls are implemented simultaneously; (b)–(d) The time series of corresponding optimal control functionsu1∗t,u2∗t,u3∗t.
(a)
(b)
(c)
(d)To further specify the extent of each control measure in reducing the outbreak of measles, we compare the number of symptomatic individuals with only one control measure and without any control measure, as shown in Figure7. The results illustrate each optimal control measure has great impact on decreasing the measles infections. In particular, the optimal enhanced vaccination for children u3c∗t, which has limited impact on lowering the first outbreak, results in a low probability of occurrence of repeated outbreaks via decreasing the susceptibility to measles.Figure 7
The number of infected individuals with symptoms under the two scenarios: without control and with only one type of control (top row). The time series of each optimal control function.u1∗t,u2∗t,u3c∗t.The timing of implementation of control measures is always a critical issue. Therefore, we study the optimal control functions and their influences on the reduction of infected individual under three different starting times: (1) the time when the number of infective individuals is at the lowest level, (2) the time when the number of infections increases but is relatively small, and (3) the time when measles spreads fast with severe level. Here, we can choose the 10th, 13th, and 15th month from January of 2013. Figure8(a) demonstrates that the relatively late implementation of control measures leads to a great increase in the number of infectious individuals. Comparing optimal control functions u1∗t,u2∗t, and u3c∗t with different starting times of control measure as shown in Figures 8(b)–8(d), we obtained that as the starting time is delayed, the optimal control strategy u1∗t has no obvious changes, while the optimal control function u2∗t must be implemented with maximum level for a long period. Different from u1∗t and u2∗t, Figure 8(d) shows that the optimal enhanced vaccination for children u3c∗t should be implemented with more maximum level (shown in green curves) and then be implemented with relatively low maximum level (shown in black curves) as the initiation time is delayed. These results indicate the earlier the control measures are implemented, the lower the outbreak is and the less (period or intensity) the optimal controls are needed.Figure 8
The effect of controls’ starting time on the number of infected individuals and the optimal control functionsu1∗t,u2∗t,u3c∗t under three cases: starting at 10th, 13th, and 15th month from January of 2013.
(a)
(b)
(c)
(d)


## 3.1. Optimality System

The existence of optimal controls follows from Corollary 4.1 of [32] since the state system satisfies the Lipschitz property with respect to the state variables and is a linear function of U, while the integrand of J is a convex function for u1(t), u2(t) and u3(t).Theorem 3.
There exist optimal controlsU∗(t) and corresponding solutions Y∗(t) that minimize J(U(t)) over Ω. In order for the above statement to be true, it is necessary that there exist continuous functions λi(t) such that(26)λ1′t=λ1−λ41−u1tβI+ξIv+λ1−λ21+u3tp1+μλ1,λ2′t=λ2−λ3ω+μλ2,λ3′t=λ3−λ51−u1tηβI+ξIv+μλ3,λ4′t=λ4−λ6σ+μλ4,λ5′t=λ7−λ6θσ+λ5−λ7σ+μλ5,λ6′t=−1+λ1−λ41−u1tβS+λ3−λ5η1−u1tβV2+λ6−λ81+u2tγ+μλ6,λ7′t=λ1−λ41−u1tξβS+λ3−λ51−u1tξηβV2+λ7−λ8γv+μλ7,λ8′t=μλ8,with the transversality conditions(27)λ1T=λ2T=λ3T=λ4T=λ5T=λ6T=λ7T=λ8T=0.Furthermore, the optimal control is given as follows:(28)u1∗t=maxminλ4−λ1βI∗+ξIv∗S∗+λ5−λ3ηβI∗+ξIv∗V2∗W1,b1,0,u2∗t=maxminλ6−λ8γI∗W2,b2,0,u3∗t=maxminλ1−λ2pS∗W3,b3,0.Proof.
The necessary conditions satisfied by optimal controls can be derived from Pontryagin’s maximum principles. Let Hamilton function be(29)H=I+12W1u12+12W2u22+W3u32+λ1dSdt+λ2dV1dt+λ3dV2dt+λ4dEdt+λ5dEvdt+λ6dIdt+λ7dIvdt+λ8dRdt.
Then, the adjoint equations with transversality satisfyλ1′=−∂H/∂S,λ2′=−∂H/∂Vk,λ3′=−∂H/∂V2,λ4′=−∂H/∂E,λ5′=−∂H/∂Ev,λ6′=−∂H/∂I,λ7′=−∂H/∂Iv,λ8′=−∂H/∂R, with  λ1(T) = λ2(T) = λ3(T) = λ4(T) = λ5(T) = λ6(T) = λ7(T) = λ8(T) = 0.
The Hamilton functionH is minimized with respect to the controls by differentiating H with respect to u1, u2, u3 on the set Ω, respectively, and then gives the following optimal conditions:(30)∂H∂u1=0,∂H∂u2=0,∂H∂u3=0.
Solving foru1∗,u2∗,u3∗ yields(31)u1∗=λ4−λ1βI∗+ξIv∗S∗+λ5−λ3ηβI∗+ξIv∗V2∗W1,u2∗=λ6−λ8γI∗W2,u3∗=λ1−λ2p1S∗W3.
By using the bounds 0 ≤u1(t) ≤ b1, 0 ≤ u2 ≤ b2, 0 ≤ u3 ≤ b3, we have the properties as Theorem 3.


## 3.2. Numerical Simulations

As mentioned in Section2.3, based on relationship between pc and p, we initially define the control function u3c(t) as the enhanced vaccination rate for children under 2 years old to replace u3(t) for general population and numerically explore the effects of optimal control measures on measles outbreak. The optimality system is solved by using the forward-backward sweep method. With an initial guess for the control variables, the state system and adjoint system are solved forward and backward, separately. The control variables are updated by using optimality condition, and this process is repeated until the system is converged. In the numerical simulations, we use parameter values given in Table 1. To balance the population and control functions in the objective functions, we choose weight constants values W1 = 20, W2 = 200, and W3 = 0.2 and the upper bound of controls b1 = 0.2, b2 = 0.2, and b3 = 5 for illustration.To study the effect of optimal control strategies on the transmission of measles infection, we initially plot the time series of the number of symptomatic infected individuals and corresponding optimal functions, as shown in Figure6. It follows from Figure 6(a) that the optimal control strategies significantly mitigate the outbreak of disease and decrease the number of symptomatic infections. Figures 6(b)–6(d) give the optimal control measures u1∗t,u2∗t and u3c∗t, which monotonically decrease after a period time of implementation with high intensity. This indicates we should enhance vaccination as soon as possible when measles begin outbreak, while measures such as personal protection and treatment should be implemented as much as possible throughout the whole epidemic.Figure 6
(a) The number of infected individuals with symptoms under the two situations: no control is used and three types of controls are implemented simultaneously; (b)–(d) The time series of corresponding optimal control functionsu1∗t,u2∗t,u3∗t.
(a)
(b)
(c)
(d)To further specify the extent of each control measure in reducing the outbreak of measles, we compare the number of symptomatic individuals with only one control measure and without any control measure, as shown in Figure7. The results illustrate each optimal control measure has great impact on decreasing the measles infections. In particular, the optimal enhanced vaccination for children u3c∗t, which has limited impact on lowering the first outbreak, results in a low probability of occurrence of repeated outbreaks via decreasing the susceptibility to measles.Figure 7
The number of infected individuals with symptoms under the two scenarios: without control and with only one type of control (top row). The time series of each optimal control function.u1∗t,u2∗t,u3c∗t.The timing of implementation of control measures is always a critical issue. Therefore, we study the optimal control functions and their influences on the reduction of infected individual under three different starting times: (1) the time when the number of infective individuals is at the lowest level, (2) the time when the number of infections increases but is relatively small, and (3) the time when measles spreads fast with severe level. Here, we can choose the 10th, 13th, and 15th month from January of 2013. Figure8(a) demonstrates that the relatively late implementation of control measures leads to a great increase in the number of infectious individuals. Comparing optimal control functions u1∗t,u2∗t, and u3c∗t with different starting times of control measure as shown in Figures 8(b)–8(d), we obtained that as the starting time is delayed, the optimal control strategy u1∗t has no obvious changes, while the optimal control function u2∗t must be implemented with maximum level for a long period. Different from u1∗t and u2∗t, Figure 8(d) shows that the optimal enhanced vaccination for children u3c∗t should be implemented with more maximum level (shown in green curves) and then be implemented with relatively low maximum level (shown in black curves) as the initiation time is delayed. These results indicate the earlier the control measures are implemented, the lower the outbreak is and the less (period or intensity) the optimal controls are needed.Figure 8
The effect of controls’ starting time on the number of infected individuals and the optimal control functionsu1∗t,u2∗t,u3c∗t under three cases: starting at 10th, 13th, and 15th month from January of 2013.
(a)
(b)
(c)
(d)


## 4. Discussion

In mainland China, measles reached its lowest level in 2012 and then resurged again with seasonal transmission since 2013. These outbreaks have led to an interest in understanding the underlying biological mechanisms and public concern of effective control. To investigate the phenomenon of seasonal outbreak in China and identify the crucial parameters played in the transmission dynamics of measles, we propose a mathematical model with periodic transmission rate. In this system, we extend the existing dynamic model [9, 10] by considering asymptomatic infection due to vaccinated populations with low antibodies level being infected. We investigate the threshold dynamics by defining the basic reproduction number. Further, based on the baseline model (1), we propose three types of feasible control measures and formulate a optimal control system to discuss how to design the optimal measures to control epidemic outbreak.We theoretically analyze dynamic behavior of the periodic system by defining the thresholdR0. We prove the disease-free equilibrium E0 is globally asymptotically stable for R0 < 1, meaning the measles can be eliminated. If R0 > 1, the disease-free equilibrium E0 is unstable and system (1) is uniformly persistent, which indicates the disease cannot be effectively eliminated. Then, by fitting our proposed model to the reported cases in mainland China from year 2013 to 2016, we estimate some unknown parameters and calculate the basic reproduction number R0 with immunization as 1.0077. Compared with original basic reproduction number of measles (R0 = 12.8) [4], this result implies that current vaccination policies have dramatically reduced the spread of disease. Numerical simulations indicate that asymptomatic infection has slight impact on the repeated outbreaks of measles and the crucial factors that significantly affect new infections are vaccination rate, transmission rate, and recovery rate. Further, we propose three types of control measures (limiting contacts between susceptible and infective individuals and enhancing treatment and vaccination) and determine the optimal control strategies to minimize the number of infected individuals at minimum costs. It has been illustrated numerically that the optimal control measures have a very attractive effect on reducing the number of infected individuals and the interventions should be implemented as soon as possible. We further observed that enhancing vaccination is especially effective in lowering the initial outbreak and making disease recurrence less likely.Note that due to lack of data on measles cases as for different ages, we formulate a measles transmission model without considering age structure. However, the vaccination of measles mostly is implemented among children, and we shall formulate the dynamic model with age structure, conduct cost-effectiveness assessment, and explore the optimal vaccination age in near future.



---

*Source: 3631923-2020-03-20.xml*