The Dynamics and Optimal Control of a Hand-Foot-Mouth Disease Model
We build and study the transmission dynamics of a hand-foot-mouth disease model with vaccination. The reproduction number is given, the existence of equilibria is obtained, and the global stability of disease-free equilibrium is proved by constructing the Lyapunov function. We also apply optimal control theory to the hand-foot-mouth disease model. The treatment and vaccination interventions are considered in the hand-foot-mouth disease model, and the optimal control strategies based on minimizing the cost of intervention and minimizing the number of the infected people are given. Numerical results show the usefulness of the optimization strategies.
Hand, foot, and mouth disease (HFMD) is a common infectious disease caused by a group of viruses known as enteroviruses (EVs) [1, 2]. HFMD usually affects children, typically affecting children who are less than 10 years, but it can also affect adults . However, adults are immune to the disease due to the antibodies in their bodies, although most of them are exposed to the virus .
The virus of HFMD spreads easily through coughing, sneezing, and infected stool. It usually takes days for a person to get symptoms of HFMD disease after being exposed to the virus of HFMD. This is called the incubation period of HFMD. Although many HFMD infected people remain asymptomatic, the symptoms of HFMD include sores in or on the mouth and on the hands, feet, and sometimes the buttocks and legs. The sores may be painful, and these sores may be eased with the use of medication . In fact, there is no specific treatment for HFMD, and many doctors do not issue medicine for this illness since HFMD is a viral disease that has to run its course .
Numerous large outbreaks of HFMD have occurred in many areas of the world, such as the United States of America, Australia, Malaysia, Japan, and China since 1997, which have caused death and neurological sequelae, and have become a growing public health threat [6–10]. Fortunately, Chinese scientists have developed the first vaccine to protect children against enterovirus 71, or EV71, that causes the common and sometimes deadly HFMD in 2013 . However, the literature on the mathematical modeling of the transmission of HFMD is rather scant. In particular, there are fewer literatures on mathematical models of HFMD with vaccination. Urashima et al. and Wang and Sung tried to find the relationship between the outbreak of HFMD with the weather patterns in Taiwan and Tokyo, respectively [12, 13]. Tiing and Labadin used a deterministic SIR model to predict the number of infected cases and the duration of an outbreak when it occurs in Sarawak . Roy and Halder proposed a deterministic SEIR model of HFMD and did numerical simulations . Liu and Yang et al. used the SEIQRS model to take into account the quarantine measure [5, 16]. Recently, Samanta discussed a delay HFMD model with pulse vaccination strategy .
In this paper, we only consider the children below the age of 10 years since the children above the age of 10 years are immune to the disease because their immune systems are relatively perfect. The aim of our study is to use mathematical modeling to gain some insights into the transmission dynamics of HFMD when the population is vaccinated. The paper is organized as follows. In Section 2, we formulate the HFMD model with vaccination and define the basic reproduction number. In Section 3, we obtain the existence of equilibria of model, prove the global stability of disease-free equilibrium, and analyze the global stability of endemic equilibrium of model by constructing the Lyapunov function. In Section 4, we discuss the optimal control problem by adding two control functions. At last, we display the numerical simulation and give the conclusion.
2. Model Formulation
Enteroviruses (EVs) that are most frequently reported as causing HFMD outbreaks include enterovirus 71 (EV71) and coxsackievirus A16 (CVA16). Other human enteroviruses serotypes, such as CVA4, CVA5, CVA6, and CVA10, have also been reported in cases of HFMD . Because only EV71 vaccine was on market which could prevent the HFMD induced by EV71 infection, we will consider dividing the infectious individuals into two classes, which are infectious individuals infected with EV71 and infectious individuals infected with CVA16 or other human enteroviruses serotypes.
Let be total number of children below the age of 10 years at time . We divide children below the age of 10 years into five compartments, including susceptible individuals , latent individuals , infectious individuals and , vaccination individuals , and recovery individuals . It is clear that . The dynamical model for HFMD transmission in children below the age of 10 years is in the following:where is the birth rate of the population; is the vaccine rate of the population; is the transmission coefficient of the infectious individuals infected with EV71; is the transmission coefficient of the infectious individuals infected with CVA16; is the natural death rate; is the per-capita rate of the progression from latent individuals to infectious individuals; is the percentage of individuals infected with EV71 from latent individuals to infectious individuals; accordingly, is the percentage of individuals infected with CVA16 or other human enteroviruses serotypes from latent individuals to infectious individuals; is the disease induced death rate of infectious individuals , , respectively; is the treatment rate of the infectious individuals , , respectively; is the removal rate of population; and are the loss of immunity rate of vaccination individuals and recovery individuals, respectively.
In our paper, in order to make the qualitative mathematical analysis, let , , , and ; we simplify model (1) to the following model:
Proof. By using , from system (2), we haveIt is obvious that , which implies that . That is, every solution of system (2) eventually enters , and is positively invariant with respect to system (2). This proves the lemma.
The dynamics of system (2) will be investigated in the following bounded feasible region:
Using the relation , we may reduce system (2) to the following equivalent system:with , on the positively invariant setIn the following, since system (5) has the same dynamic as (2), we will discuss the dynamic of system (5) on .
Following van den Diessche and Watmough [17, 18], we can obtain the basic reproduction number:Each term in has clear epidemiological interpretation. is the proportion that latent individuals progress to infectious class. is the average infectious period. is the total amount of population in the case that the infected individuals in population do not exist. Thus, are average new cases generated by a typical infectious member in the entire infection period.
The basic reproduction number , for model (2) in the absence of controls, i.e., in the case , which means that model (2) does not have vaccination individuals and recovery individuals , is proportional to the transmission coefficient and is given byIt is clear thatwhich implies that the vaccination and treatment have contributed to decrease of the . That is, the vaccination and treatment help to slow down the HFMD spread.
Three parameters have a high impact on : and decrease , respectively, and increases .
3. The Existence and Stability of Equilibria
We first discuss the existence of equilibria of system (5). Directly calculating system (5), we obtain the disease-free equilibrium , where , and . In addition, there exists a endemic equilibrium when , whereSummarizing the above discussion, we can obtain the following result.
Proof. The Jacobian matrix of system (5) at the disease-free equilibrium isIt is clear that and are the eigenvalues of matrix . The rest of the eigenvalues of matrix satisfy the following equation:Obviously,It implies that (12) has two real roots, and , which satisfyIf , we have , which implies that the real parts of and are both negative. That is, the disease-free equilibrium is locally asymptotically stable. Meanwhile if , we obtain . It implies that the real part of or the real part of is positive. Therefore, is unstable.
For the critical case , the Jacobian matrix has three negative real eigenvalues , , and , and one zero eigenvalue.
We introduce the matrix of eigenvectorswith , such that , whereWe make the linear transformation , wherewith + / , .
The Jacobian matrix for the differential equations of about the zero equilibrium is exactly . To analyze the local asymptotic stability of this zero equilibrium, we need to calculate the restricted dynamical system on the center manifold for sufficiently small and , , . Note that ; from the second, third, and forth equations of system (5), we obtainNext, we make use of , , , and to obtainSince restricted system (19) is stable about , original system (5) is locally stable about the disease-free equilibrium when .
In the following, we study the global stability when . Let ; we have Furthermore, if and only if or . Therefore, the largest compact invariant set in is the singleton . LaSalle’s invariance principle  then implies that is globally stable in .
Next, we discuss the global asymptotical stability of the endemic equilibrium of system (5). The local stability of the endemic equilibrium firstly was discussed, and the global stability of the endemic equilibrium also was discussed by constructing the Lyapunov function.
Theorem 4. If , the endemic equilibrium of system (5) is locally asymptotically stable.
Proof. The Jacobian matrix of system (5) at the endemic equilibrium isIt is clear that the eigenvalues of matrix satisfy the following equation:whereBy directly calculating, we haveThe Routh-Hurwitz criterion  implies that all eigenvalues of characteristic equation (22) have negative real part; that is, is locally asymptotically stable when .
It is difficult to show the global stability of endemic equilibrium by the theoretical methods. We will use the numerical simulation to display the global stability of endemic equilibrium ; see Figure 1. The parameters are taken to be , , , , , , , , , and , respectively. Accordingly, the basic reproduction number . The simulation demonstrates that endemic equilibrium may be globally stable when .
4. The HFMD Model with Optimal Controls
In this section, we present the optimal control problem by adding to the model (2) two control functions and . The HFMD model with controls is given by the following equations:The aim is to find the optimal values and of the controls and , such that the associated state trajectories , , , , and are solution of system (25) in the time interval with initial conditions , , , , and and minimize the objective functional. Here the objective functional considers the number of latent individuals , the number of the infectious individuals , and the implementation cost of the strategies associated with the controls , . The controls are bounded between 0 and 1.
We consider state system (25) of ordinary differential equations in with the set of admissible control functions given byThe objective functional is given bywhere the constants and are a measure of the relative cost of the interventions associated with the controls and , respectively.
We consider the optimal control problem of determining , associated with an admissible control pair on the time interval , satisfying (25) and the initial conditions , , , , and and minimizing cost functional (27); that is,
Theorem 5. There exists an optimal control pair such thatsubject to state system (25) with initial conditions , , , , and .
Proof. The integrand of the objective functional given by (27) is convex on the closed, convex control set . The conditions for the existence of optimal control are satisfied as the model is linear in the control variables and is bounded by a linear system in the state variables .
According to the Pontryagin Maximum Principle , we now derive the necessary conditions that a pair of optimal controls and corresponding states must satisfy. To this purpose, we define the Hamiltonian function for the system:where , , are the adjoint variables.
Theorem 6. Given an optimal control on and corresponding state solution of corresponding state system (25) with initial conditions , , , , and , there exist adjoint variables , , satisfyingwith transversality conditions (or boundary conditions) beingFurthermore, the optimal controls and are given by
Proof. The adjoint system results from Pontryagin’s Principle :with zero transversality. To get the characterization of the optimal control given by , we solve the equations on the interior of the control set:Using bounds on the controls, we obtain the desired characterization.
5. Numerical Results and Discussion
In this section, the numerical simulation results of the optimized control measures for HFMD model (25) with vaccination are presented. First, we solve system (25) over the time interval using a forward fourth-order Runge-Kutta scheme and transversality conditions , . Then, system (31) is solved by a backward fourth-order Runge-Kutta scheme using the current iteration solution of (25). The controls are updated by using a convex combination of the previous controls and the values from (33). The iteration is stopped when the values of unknowns at the previous iteration are very close to the ones at the present iteration.
We first compare the number of the susceptible individuals, latent individuals, infectious individuals, vaccination individuals, and the recovery individuals with and without controls, respectively. Take , , , , , , , , , , , , and the initial conditions , we have the basic reproduction number , and the numerical results are depicted in Figures 2 and 3. Figure 2 shows the control variables and when . By Figure 2, we see that, to minimize the number of infectious and latent individuals, the control keeps the increasing trend during 5 days; during the remaining 15 days, it decreases to the lower bound. The control is at the upper bound during 17 days; during the remaining 3 days, it decreases to the lower bound. Figure 3 shows that the number of the susceptible individuals is higher, the numbers of the latent individuals, infectious individuals, and recovery individuals are lower, and the number of the vaccination individuals is with almost no change when controls are considered.
We discussed the influence of immune loss rate on the spread of disease. That is, we discuss the influence of on the basic regeneration number in the following. Because ofwe know the basic reproduction number increases as increases. In order to control the basic regeneration number less than 1, we need the immune loss rate to satisfy , under the condition , where . Let , respectively; other parameter values are the same as those in Figure 3. We obtain Figure 4. The higher the immune loss rate , the greater the number of infections and latent individuals. At the same time, higher immunization loss rate indicates that the immunization control measures are weakened. Accordingly, the treatment control measures need to be strengthened.
We also compared the number of infections and latent individuals under different control measures (see Figure 5). In Figure 5, the red dots indicate that control measures and are implemented simultaneously, the blue dots indicate that only control measure is implemented, and the black dotted line indicates that only control measure is implemented. The numerical simulation results show that the number of infections is minimal when the control measures and are implemented simultaneously. If only one control was implemented, the treatment control would be more effective than vaccination control in controlling the number of infectious and latent individuals. The trend of the controls and is displayed in Figure 6 under different control measures.
No data were used to support this study.
Conflicts of Interest
All authors declare that they have no conflicts of interest.
The authors also would like to thank Professor Jianquan Li for his suggestions, which helped them to deal with problems they met. This study was supported by NSFC Grants 11301314, 11501443, and 11671243, by Science and Technology Bureau of Weiyang District, Xi’an, Grant 201708, and by Shaanxi Provincial Education Department Grant 15JK1081.
G. P. Samanta, “Analysis of a delayed hand-foot-mounth disease epidemic model with pulse vaccination,” Systems Science Control Engineering, vol. 2, pp. 61–73, 2014.View at: Google Scholar
G. L. Gilbert, K. E. Dickson, M.-J. Waters, M. L. Kennett, S. A. Land, and M. Sneddon, “Outbreak of enterovirus 71 infection in victoria, australia, with a high incidence of neurologic involvement,” The Pediatric Infectious Disease Journal, vol. 7, no. 7, pp. 484–488, 1988.View at: Publisher Site | Google Scholar
L. G. Chan, U. D. Parashar, M. S. Lye et al., “Deaths of children during an outbreak of hand, foot, and mouth disease in Sarawak, Malaysia: clinical and pathological characteristics of the disease,” Clinical Infectious Diseases, vol. 31, no. 3, pp. 678–683, 2000.View at: Publisher Site | Google Scholar
T. Fujimoto, M. Chikahira, S. Yoshida et al., “Outbreak of central nervous system disease associated with hand, foot, and mouth disease in Japan during the summer of 2000: detection and molecular epidemiology of enterovirus 71,” Microbiology and Immunology, vol. 46, no. 9, pp. 621–627, 2002.View at: Publisher Site | Google Scholar
F.-C. Zhu, F.-Y. Meng, J.-X. Li et al., “Efficacy, safety, and immunology of an inactivated alum-adjuvant enterovirus 71 vaccine in children in China: A multicentre, randomised, double-blind, placebo-controlled, phase 3 trial,” The Lancet, vol. 381, no. 9882, pp. 2024–2032, 2013.View at: Publisher Site | Google Scholar
M. Urashima, N. Shindo, and N. Okabe, “Seasonal models of herpangina and hand-foot-mouth disease to simulate annual fluctuations in urban warming in Tokyo,” Japanese Journal of Infectious Diseases, vol. 56, no. 2, pp. 48–53, 2003.View at: Google Scholar
O. Diekmann, J. A. P. Heesterbeek, and J. A. J. Metz, “On the definition and the computation of the basic reproduction ratio R0 in models for infectious disease in heterogeneous populations,” Journal of Mathematical Biology, vol. 28, no. 4, pp. 365–382, 1990.View at: Publisher Site | Google Scholar | MathSciNet
J. P. LaSalle, “The stability of dynamical systems,” in Proceedings of the CBMS-NSF regional conf.ser in Appl. Mat, vol. 25, Philadelphia, 1976.View at: Google Scholar
M. Kot, Elements of Mathematical Biology, Wiley Interscience, 1962.View at: MathSciNet
W. H. Fleming and R. W. Rishel, Deterministic and Stochastic Optimal Control, vol. 1, Springer, New York, NY, USA, 1975.View at: MathSciNet
L. S. Pontryagin, V. G. Boltyanskii, R. V. Gamkrelidze, and E. F. Mishchenko, The Mathematical Theory of Optimal Processes, Wiley Interscience, 1962.View at: MathSciNet