Modelling the Transmission Dynamics and Control of Mumps in Mainland China

Mumps is a common childhood viral disease and children have been vaccinated throughout the world since 1967. The incidence of mumps has increased with more than 300,000 young people infected with mumps annually in mainland China since 2005. Therefore, we designed and analyzed long-term mumps surveillance data in an SVEILR (susceptible–vaccinated–exposed–severely infectious–mildly infectious–recovered) dynamic transmission model with optimized parameter values to describe the dynamics of mumps infections in China. There were 18.02% of mumps infected young adults seeking medical advice. The vaccine coverage has been insufficient in China. Young adults with frequent contact and mild infection were identified as a major driver of mumps epidemics. The reproduction number of mumps was determined 4.28 in China. Sensitivity analysis of the basic reproduction number and the endemic equilibrium was conducted to evaluate the effectiveness of mumps control measures. We propose to increase the vaccine coverage and make two doses of MMR (Measles, mumps and rubella) vaccines freely available in China.


Introduction
Mumps is best known as a common childhood viral disease and is highly contagious to human beings. Initial signs and symptoms often include fever, muscle pain and headache, then usually followed by painful swelling of one or both parotid salivary glands [1]. The disease caused by the mumps virus, the causative agent of mumps infection, is an enveloped RNA virus [2]. Since the disease is generally benign and self-resolving, its mortality is rare, but aseptic meningitis can affect 10% of case-patients [3]. Mumps is a significant cause of pediatric deafness, and up to 37% of post-pubertal males develop orchitis, 13% of whom have impaired fertility [1].
Due to the lack of vaccines, most teenagers have been mainly infected by those patients aged 4-15 years. Transmission of the virus is by direct physical contact, droplet spread, or contaminated fomites [4][5][6]. The incubation period is about 15 to 24 days (median is 19 days) [7]. Infected patients become the most contagious in 1-2 days before onset of clinical symptoms and continue so for a few days afterwards. Generally speaking, the infectious period is about eight days [8], and the patients will recover between 10 to 14 days [9,10].
Mumps has periodic outbreaks and no specific antiviral therapy, and treatment is just mostly symptomatic and supportive [11]. In 1967, the attenuated mumps virus vaccine was licensed in the United States [12]. In addition, mumps is preventable by two doses of the mumps vaccine in many developed countries nowadays [13]. They includes it in their immunization programs, often in combination with measles and the rubella vaccine [14,15]. Nowadays, children aged 18-24 months just routinely receive one dose of the measles-mumps-rubella (MMR) vaccine free of charge in China, and National Health and Family Planning Commission of the People's Republic of China (NHFPC) doesn't emphasize inoculating the second dose of MMR. In recent years, there are more than 300,000 young people infected with mumps every year in China [16,17].
The number of mumps patients is just smaller than hand, foot and mouth disease in all pediatric infectious diseases [16,17]. Compared with other usual vaccine-preventable diseases, such as measles and pertussis, mumps is more common. Due to its less severe acute manifestation, mumps have been somewhat neglected. Nevertheless, the UK and USA have been inspired a new interest in mumps by some large outbreaks. In the UK, a large mumps epidemic began in 2004 and reached the peak in 2005 with about 56,000 reported cases [18,19]. Most of these cases were in young adults attending colleges or universities [20]. In 2006, the USA underwent a multi-state outbreak involving 6584 reported cases, with the highest attack rate among persons 18-24 years old, many of whom were college students [12]. In affected colleges, most case-patients had been inoculated with a second dose of the MMR vaccine more than 10 years previously [21,22]. This was the first large-scale US mumps outbreak among the population with two-dose vaccines. From then on, people realized that even the two doses of the vaccine could not completely control mumps either.
Mathematical models have become vital tools in understanding the spread and control of infectious diseases well. By setting up a suitable epidemic model, we can put forward a lot of practical prevention and control measures to curb the epidemic of disease. So far, there have been a few papers using dynamic models to study mumps [23,24]. Qu [23] proposed a non-autonomous SVEILHR (susceptible-vaccinated-exposed-mild infectious-severe infectious-hospitalized-recovered) model with a seasonal varying transmission rate to describe the epidemic of mumps, and suggested improving vaccine coverage and providing two doses of MMR (Measles, mumps and rubella) vaccine by the government in China. Liu [24] studied the effects of heterogeneity on mumps spread by constructing a multi-group SVEIAR (susceptible-vaccinated-exposed -symptomatic-asymptomatic-recovered) mumps model with asymptomatic infection, general vaccinated and exposed distributions and established the threshold dynamics of the model.
The basic structures of this paper are as follows. In the next section, an SVEILR mumps model is formulated, and we give the basic reproduction number R 0 and the existence of equilibria. Section 3 discusses the global stability of the model. In Section 4, we advance the optimal parameters, simulation of the real data from 2009 to 2014, and the prediction results of 2014 and 2015. Sensitivity analysis of the basic reproduction number R 0 and endemic equilibrium P * are carried out in Section 5. We conclude with some discussions in Section 6 about the role of vaccine and the preventive measures on mumps. We close with a conclusions section.

The Mumps Model and Basic Reproduction Number R 0
We propose a mathematical model to understand the transmission dynamics and prevalence of mumps in mainland China. Since we will simulate real data for half a year, we think the total population is constant in a short period of time. We assume that the birth rate equals the natural death rate, denoted by µ. The model is constructed based on the characteristics of mumps transmission; therefore, the population associated with mumps is divided into six epidemiological sub-classes: the proportion susceptible to total population (S), the proportion of vaccinated to total population (V), the proportion exposed to total population (E, infected but not infectious), the proportion severely infectious to total population (I, severely infectious requiring medical attention), the proportion mildly infectious to total population (L, mild infections, including both asymptomatic and those with mild symptoms and self-care), and the proportion recovered to total population (R), subject to the restriction S + V + E + I + L + R = 1. The transmission dynamics associated with these sub-classes are illustrated in Figure 1.
Considering the vertical transmission, we assume that a fraction ρ of the offspring from the exposed parents is born into the exposed class E. Consequently, the birth flux into the exposed class is given by ρµE and the birth flux into the susceptible class is given by µ − ρµE. The other defined parameters in the model Equation (2) are listed below: • β: transmission rate, • λ: loss of vaccination rate, • ε: vaccine coverage of the susceptible, • ε 1 : vaccine coverage of the exposed, • κ: invalid vaccination rate, • α: rate moving from exposed to severe or mild infectious, • γ: proportion of the severe infections seeking medical advice, • δ 1 : rate moving from severe infectious to recovered, • δ 2 : rate moving from mild infectious to recovered. To simplify the research, we don't consider the spatial stratified heterogeneity of the population and relevant determinants, and consider that the mumps model has homogeneous mixing, and an individual has an equal chance of contacting any individual among the population. The SVEILR mumps model is given by six ordinary differential equations: A person was infected with the virus and then fully recovered. After that, he is typically immune for life [20]. The first five equations are independent of R in Equation (1), and then it suffices to study the following sub-system: The biologically feasible region of Equation (2) is Ω = {(S, V, E, I, L) ∈ R 5 + : 0 ≤ S + V + E + I + L ≤ 1}, which can be verified as positively invariant (i.e., given non-negative initial values in Ω, all solutions to Equation (2) have non-negative components and stay in Ω for t ≥ 0) and globally attracting in R 5 + with respect to Equation (2). Therefore, we restrict our attention to the dynamics of Equation (2) in Ω.
It is easy to see that model Equation (2) always has a disease-free equilibrium P 0 = (S 0 , V 0 , 0, 0, 0), where According to the next generation matrix developed by van den Driessche and Watmough [25], we define the basic reproduction number of model Equation (2) as where , which represent the average numbers of the infected individuals by a single exposed parents, severely infectious, and mildly infectious individuals in a fully susceptible population, respectively.
The basic reproduction number R 0 represents for the average number of new infections brought out by one infectious during the initial patient's infectious (not sick) period [26]. If R 0 > 1, model Equation (2) has a unique endemic equilibrium P * = (S * , V * , E * , I * , L * ). Here, we do not give the exact expression of P * , the detailed analysis can be seen in the Appendix A. Then, we have the following proposition: Proposition 1. Model Equation (2) always has a disease-free equilibrium P 0 ; and, if R 0 > 1, model Equation (2) has a unique endemic equilibrium P * .

Mathematical Analysis Results
Regarding the stability of the disease-free equilibrium P 0 and endemic equilibrium P * , we have the following Theorems. The detailed proof process can be seen in the Appendix B. Theorem 1. If R 0 < 1, then P 0 is stable, and, if R 0 > 1, then P 0 is unstable. Theorem 2. The disease-free equilibrium P 0 of model Equation (2) is globally asymptotically stable if R 0 ≤ 1.
The above results show that mumps will be eliminated from the community if the epidemiological threshold R 0 can be brought to a value less than unity. (2) has a unique endemic equilibrium P * = (S * , V * , E * , I * , L * ), which is globally asymptotically stable.
When ρ = 0, ε 1 = 0, we do not testify the global stability of the model Equation (2). Nevertheless, using parameter estimation in the next section, the values of ρ and ε 1 are very small, almost negligible. It shows that Theorem 3 still have great significance.
Mathematical analysis shows that the basic reproduction number R 0 is a sharp threshold completely determining the global dynamics of model Equation (2): if R 0 ≤ 1, the disease-free equilibrium P 0 is globally asymptotically stable and the mumps dies out; if R 0 > 1, the endemic equilibrium P * is globally asymptotically stable, namely, the mumps persists at the endemic equilibrium state if only it initially exists.

Data Analysis
It was reported by the Ministry of Health of the People's Republic of China that mumps included made the list as a Category C Infectious Disease (Monitoring and Managing of Infectious Disease) in 1989 [27]. On the basis of the Chinese Center for Disease Control and Prevention (China's CDC), using month as the period does statistical work on the patients who were infected by mumps (see Figure 2) [16,17]. China's CDC publishes the surveillance data of mumps from each province in mainland China except Hong Kong, Macao and Taiwan. From the network direct reporting surveillance data, we can acquire a lot of detailed information from cases, such as the area code, gender, occupation, date of birth, date of onset, date of diagnosis, address, and, especially, classification of disease, which is labeled as clinically diagnosed cases.
According to the data of mumps analysis, it is not difficult to find that the epidemic of mumps has certain regularity. Judging by larger number of mumps cases, we can find a peak between February and September, another peak appears between September of this year to February of the next year, and the trough appears in February and September each year [28]. In fact, it is associated with the holiday scheme of the elementary and secondary schools in China. Every year in February, which is just after the lunar New Year, all primary and secondary schools will end more than a month of winter vacation. Similarly, every year in July to September, all schools have a summer vacation to rest for two months. This time happens to be when the mumps go through the trough. On account of the patients with mumps being mainly 4 to 15 years old, we can conclude that the main cause of the mumps epidemic in China is close contact among the students in their schools. Considering Chinese summer vacation and winter vacation, in order to facilitate the fitting, we therefore divide the reported data of each year into two parts: one is the first half of the year (from February to September), denoted by FHY (the first half of the year), and the other is the latter half of the year (from September of this year to February of the next year), denoted by LHY (the latter half of the year). Therefore, the data of the February and September are simulated twice.
By applying a chi-square test to fit annual data (six data points in FHY, eight data points in LHY), at most six parameters (or initial values) every time can be fit into the first half of the year, and at most four parameters (or initial values) every time can be fit into the latter half of the year. We thus consider ε = ε 1 , The remaining parameters and initial values (β, λ, ε, γ, S(0) and V(0)) in model Equation (2) are estimated through calculating the minimum sum of chi-square value (MSCV) [30]. (5) Assume that the population about mumps is 140,000 every month. Before fitting, one put the real data reducing 140,000 times in the simulation as the proportion severely infectious seeking medical advice to total population. In addition, after fitting, we can multiply by 140,000 to the estimation valves.
The data presented in Figure 2 refers to the clinical data from China's CDC, denoted by I, and we can see the optimal parameter values and initial values of FHY and LHY from 2009 to 2014 in Tables A1 and A2. In addition, the chi-square test of goodness of fit is shown in Appendix C. Here, we also calculate the basic reproduction number, R 0 , the median and the arithmetic mean of the optimal parameters, which are estimated by MSCV as shown in Table 1. The difference between the same parameters is very small in FHY and LHY of different years. In addition, R 0 of each year is stable around 4, the range of R 0 is 2.6529-5.5563, the median of R 0 is 4.7497 and arithmetic mean of R 0 is 4.2806. This suggests that the endemic equilibrium of model Equation (2) is globally asymptotically stable. It is said that, in this current situation, mumps will continue to spread.
As we know, the vaccination rate of the susceptible is ε = η ϕ, and the vaccination rate of the exposed is ε 1 = η 1 ϕ 1 . Here, η and η 1 are defined to be the percentage of the number of vaccinations, and ϕ are ϕ 1 to defined to be the rate of progression to the vaccinated. Unfortunately, we can't get the respective parameter values of η, η 1 , ϕ and ϕ 1 . We just obtain their products ε and ε 1 in model Equation (2) by MSCV.
Using model Equation (2) and the parameters from Tables A1 and A2, one carries on the data fitting to the number of severely infectious seeking medical advice individuals, as shown in Figure 2, and the numerical results are found to be a good match with the data of mumps from 2009 to 2014. From 1960 to 1980, R 0 in the Netherlands, England and Wales was 11-14, and in various states of USA, R 0 was 4-7 [9]. In recent years, the literature estimates the basic reproduction number of mumps was about 6.5428 in China [23]. Edmunds [31] used the default matrix configuration to estimate R 0 for mumps was about 3.6-4.5 in many European countries (England and Wales, the Netherlands, Finland, Denmark, East Germany and Italy) from 1970-1990. Kanaan [32] used the matrix models with individual heterogeneity to estimate R 0 for mumps was 19.3 (95% credible region (CR) 4.0-31.5) in UK in 1986, if heterogeneity is not considered, R 0 was much smaller (about 4.44). In addition, by our model, we estimate that the arithmetic mean of R 0 of mumps is about 4.2806. In the current literature, R 0 for mumps varies quite substantially between 3.6 and 19.3. This may be partly due to different source populations and timing of the studies. As our R 0 was at the lower end of this variation, we think that this was mostly due to a large proportion of once vaccinated (with partial immunity; i.e., asymptomatic or mild infection; hence, lower subsequent infections due to lower virus transmission).
Nowadays, mumps becomes a endemic disease in China and the number of cases tend to be stable each year. There are some connections between the epidemic of mumps and the vacation in China when it is summer and winter. Therefore, the fitting parameters and initial condition of the latter half of year 2013 and the first half of year 2014 may be used to predict the number of clinical cases of the latter half of year 2014 and the first half of year 2015. In fact, the prevalence status of mumps from the last year can be used to predict the epidemic situation of next year.

Sensitivity Analysis of the Basic Reproduction Number and Endemic Equilibrium
The beginning of a disease transmission is directly related to the basic reproduction number R 0 , while disease prevalence is closely linked with the endemic equilibrium P * [34]. Sensitivity analysis is conducted to identify how closely input parameters that are related to predictor parameter, and determine level of change necessary for an input parameter to acquire the ideal value of a predictor parameter. In order to study the effectiveness of mumps control strategies, we compute the sensitivity indices [34][35][36], which intrinsically measure the relative changes in R 0 or the state variables at P * with changes of model parameters. Definition 1. The normalized forward sensitivity index of a variable, u, depending differentiably on a parameter, p, is defined as [34]: ∆ u p = p u × ∂u ∂p .
In view of that, it is impossible to control the values of µ, we only compute sensitivity indices of R 0 and P * with respect to the eight parameters p i (i = 1, 2, ..., 8): β, δ 1 , δ 2 , ε, ε 1 , λ, α, γ. Here, we consider the parameters from the first half of 2014 (as shown in Table A1). By the way, we consider δ 1 = δ 2 = 1/12, according to the Equation (4), it is obvious that γ has no influence on R 0 . Furthermore, the sensitivity indices of R 0 to the seven parameters, ∆ R 0 p i are shown in Table 2.
In addition, one computes the sensitivity indices of the P * , which determines the relative impacts of different parameters on disease prevalence. Using the parameters (as shown in Table A1) from the first half of the year 2014 yields P * = (0.146893606, 0.192406471, 0.000461376, 0.000054189, 0.000238137).

Discussion
According to sensitivity analysis of R 0 and P * with regard to parameters, several effective measures for mumps' control and prevention can be put forward and put into practice.
(1) Cut off transmission routes as soon as possible in order to restrain the disease spread among the crowd (reduction β). β is the most sensitive parameter to R 0 . In addition, β is the most sensitive parameter for S * , the second sensitive parameter for V * , E * , L * , and the third sensitive parameter for I * , respectively. Consequently, parents and teachers should frequently urge students to popularize health knowledge and preserve good personal hygiene habits. Students should wash hands before meals and, after using the toilet, minimize the possibility of getting in touch with other students.
(2) Increasing vaccine coverage (increasing ε = η ϕ). ε is the third sensitive to R 0 and L * , and the most sensitive parameter for V * . Nevertheless, it is difficult to separately obtain the mumps vaccination rate η and the average vaccination time 1/ϕ. For most cities and provinces in China, according to the National Immunization Schedule by National Health and Family Planning Commission of the People's Republic of China (NHFPC) [37], MMR is one of the vaccinations supplied for free by the government, and children just get vaccinated by one dose.
However, as we know, the most common preventative measure against mumps is a vaccination with two doses of the mumps vaccine, applying in many developed countries [22,23]. In the United Kingdom, it is conventional to give children at age 13 months with another dose of vaccine at 3-5 years (preschool). This confers immunity for life. The efficacy of the vaccine depends on the strain of the vaccine, but it is usually around 80% [38]. The American Academy of Pediatrics proposes the routine administration of MMR vaccine at ages 12-15 months and at 4-6 years old [39]. In Canada, provincial governments and the Public Health Agency of Canada have all took part in awareness campaigns to encourage students ranging from the kindergarten to college to get vaccinated. Thus, the Chinese Government, health departments across the country and hospitals should promote teenagers of the right age to continue inoculating with the mumps vaccination.
Of course, if we want to explicate two doses of vaccine accurately, we will need to use pulse model and so on, and this will be our future work.
(3) Recover as soon as possible (decreasing incubation 1/δ 1 and 1/δ 2 ). δ 2 is the most sensitive parameter to L * , and the second sensitive to R 0 , S * . δ 1 is the most sensitive parameter for I * , the third sensitive parameter for S * , respectively. Try to cut down the source of infection. Patients should be treated actively so that they can recover early. Meanwhile, hospitals are supposed to enhance infection control practices to avoid nosocomial cross infection.
(4) Since δ 1 = δ 2 , γ has no influence on R 0 , but it is a sensitive parameter to I * and L * . Thus, we suggest that symptomatic patients should seek medical attention and take necessary precautions to prevent further spread. Actually, 18.02% of persons infected with the mumps virus seek medical attention (seeing Table 1, γ = 0.1802), therefore forming the group of severely infectious (I). The vast majority, 81.98% of infected mumps patients, do not seek medical attention and form the group (L). This group consisted of both mild infections and those completely asymptomatic. From [40], 20-40% of mumps infections may be asymptomatic. That is to say, about 60-80% infections have mild (part of group L) or more severe symptoms (mostly in group I). From our results, the proportion of patients seeking medical attention (γ) was not high. Therefore, a large proportion of mumps infections are self-treated, thus promoting further spread of the disease.
Owing to the huge and high heterogeneity of the country, in fact, we should consider the spatial stratified heterogeneity of the population and relevant determinants (e.g., climate, temperature, humidity, longitude and latitude, even people's behavior habits, and so on).
Taking different spatial stratifications, we have reconsidered a new multi-group model Equation (5). The total population N is divided into n groups, and each subgroup represents different areas. Some parameters in the model Equation (5) are listed below: • β ij : rate of disease transmission between susceptible individuals in group i and infectious individuals in group j, • σ ik : transfer rate move from the k-th susceptible group to i-th susceptible group, • ζ il : transfer rate move out the i-th susceptible group into l-th susceptible group, • ϕ ik : transfer rate move from the k-th vaccinated group to i-th vaccinated group, • υ il : transfer rate move out the i-th vaccinated group into l-th vaccinated group, • ω ik : transfer rate move from the k-th recovered group to i-th recovered group, • φ il : transfer rate move out the i-th recovered group into l-th recovered group.
Other parameters of different subgroups are the same as model Equation (1). Due to the incubation period and the duration of the mumps being shorter, we only consider that three categories' compartments (S, V and R) have transfer of each other, and the multi-group model as follows: According to the clinical data from different provinces of China's CDC [16,17], we can also simulate the parameters of the model. From the analyses of the parameters and the basic reproduction number, we can obtain a lot more meaningful conclusions by studying the data of different spatial and natural environment in different regions. The results of this study will appear in the later papers.

Conclusions
We have described here an improved autonomous model with a bilinear incidence rate. The model has been developed both in practical and theoretical terms compared to our previous model [23]. The model was concise and we performed extensive sensitivity analysis. The model proved the following theoretical results: (i) by Routh-Herwitz criteria, we proved that, if R 0 < 1, the disease-free equilibrium is stable, and R 0 > 1, the disease-free equilibrium is unstable; (ii) by LaSalle's Invariance Principle, we proved the disease-free equilibrium is globally asymptotically stable if R 0 ≤ 1; and, (iii) with the help of the Lyapunov function, we proved that, if R 0 > 1, for ρ = ε 1 = 0, the unique endemic equilibrium is globally asymptotically stable. Therefore, we would propose the use of this model for further studies.
Using the arithmetic mean of the parameters of model Equation (2), we can calculate R E = 0, R I = 0.8326, and R L = 3.7879. Combining the parameter analysis, one can find that vertical transmission is not the important factor that causes mumps' sustained outbreaks (ρ was relatively small and almost equal to 0 by our simulation), and young adults with frequent contacts and mild infection were identified as a major driver of mumps epidemics. It can be concluded that a mild infection that carries the virus but doesn't receive medical attention plays an important part in leading to the epidemic of mumps.
In this work, an SVEILR mumps model incorporating imperfect vaccination, vertical transmission, and mild cases are formulated to describe the dynamics of mumps transmission. As far as we can, we conduct statistical assessments on the sensitivity analysis of R 0 and the endemic equilibrium P * to parameters, and put forward several corresponding measures to control the spread of mumps. Because of the huge and high heterogeneity of the country, it is important to consider the spatial stratified heterogeneity of the population with relevant determinants (e.g., climate, temperature, humidity, longitude and latitude, even people's behavior habits, and so on) and choose certain strata of the population. However, in such a huge population it is very difficult to make the strata compatible. We will try to study this problem in the future. Author Contributions: Yong Li, Xianning Liu and Lianwen Wang participated in study conception and design, and data collection; Yong Li and Lianwen Wang performed the mathematical analyses and statistical analyses; Yong Li wrote the manuscript. All authors contributed to the interpretation of the results, revised the manuscript critically and approved the final version of the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.

Proof of Theorem 2. Consider the following Lyapunov function:
The derivative of H along the solutions of model Equation (2), together with S ≤ S 0 and V ≤ V 0 , yieldsḢ Since all of the parameters are positive, it follows thatḢ ≤ 0 for R 0 ≤ 1. There exists a singleton P 0 , as the maximal compact invariant set in the set {(S, V, E, I, L) ∈ Ω :Ḣ = 0}. Therefore, by LaSalle's Invariance Principle [41], every solution of model Equation (2), with initial conditions in Ω approaches P 0 as t → ∞, whenever R 0 ≤ 1. This completes the proof.

Proof of Theorem 3. Let
When ρ = ε 1 = 0, the model Equation (2) is transformed into the following form: It is easy to find that model Equation (A10) has a unique endemic equilibrium P 1 * (1, 1, 1, 1, 1), and the global stability of P 1 * is the same as that of P * . Thus, we will investigate the global stability Defining the Volterra-type Lyapunov function For the equilibrium state P * , we have the following equations: Then, differentiating H with respect to t along solution curves of model Equation (A10) and considering Equation (A11) gives     Table A3. Chi-square values and degrees of freedom for the first half of the year (FHY, e.g., 09(2-9) stands for February to September of 2009, and υ denotes degrees of freedom and AV denotes the accepting value at a 5% significant level with degrees of freedom 1).  Table A4. Chi-square values and degrees of freedom for the latter half of the year (LHY, e.g., 09(9)-10(2) stands for September of 2009 to February of 2010, and υ denotes degrees of freedom and AV denotes the accepting value at a 5% significant level with degrees of freedom 1). 09(9)-10(2) 10(9)-11(2) 11(9)-12(2) 12(9)-13(2) 13(9)-14 (