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.


Introduction
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 [2]. However, adults are immune to the disease due to the antibodies in their bodies, although most of them are exposed to the virus [3].
The virus of HFMD spreads easily through coughing, sneezing, and infected stool. It usually takes 3 − 7 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 [4]. 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 [5].
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][7][8][9][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 [11]. 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 [14]. Roy and Halder proposed a deterministic SEIR model of HFMD and did numerical simulations [15]. 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 [2].
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.

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 [1]. 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 1 infected with EV71 and infectious individuals 2 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 1 ( ) and 2 ( ), vaccination individuals ( ), and recovery individuals ( ). It is clear that ( ) = ( )+ ( )+ 1 ( )+ 2 ( )+ ( )+ ( ). The dynamical model for HFMD transmission in children below the age of 10 years is in the following: where > 0 is the birth rate of the population; ≥ 0 is the vaccine rate of the population; 1 > 0 is the transmission coefficient of the infectious individuals infected with EV71; 2 > 0 is the transmission coefficient of the infectious individuals infected with CVA16; > 0 is the natural death rate; > 0 is the per-capita rate of the progression from latent individuals to infectious individuals; ≥ 0 is the percentage of individuals infected with EV71 from latent individuals to infectious individuals; accordingly, 0 < 1 − < 1 is the percentage of individuals infected with CVA16 or other human enteroviruses serotypes from latent individuals to infectious individuals; 1 , 2 > 0 is the disease induced death rate of infectious individuals 1 , 2 , respectively; 1 , 2 ≥ 0 is the treatment rate of the infectious individuals 1 , 2 , respectively; ≥ 0 is the removal rate of population; 1 ≥ 0 and 2 ≥ 0 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 = 1 = 2 , = 1 = 2 , = 1 = 2 , and ( ) = 1 ( ) + 2 ( ); we simplify model (1) to the following model: In the next section, we will discuss dynamics of system (2). It is obvious that any solution of system (2) with nonnegative initial values is nonnegative.
Following van den Diessche and Watmough [17,18], we can obtain the basic reproduction number: Each term in 0 has clear epidemiological interpretation. /( + ) is the proportion that latent individuals progress to infectious class. 1/( + + ) 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 0 , for model (2) in the absence of controls, i.e., in the case = = 0, which means that model (2) does not have vaccination individuals ( ) and recovery individuals ( ), is proportional to the transmission coefficient and is given by It is clear that which implies that the vaccination and treatment have contributed to decrease of the 0 . That is, the vaccination and treatment help to slow down the HFMD spread. Three parameters have a high impact on 0 : and decrease 0 , respectively, and increases 0 .
In the following, we will discuss the stability of equilibria of system (5). The stability of disease-free equilibrium of system (5) firstly was proved. Theorem 3. If 0 ≤ 1, the disease-free equilibrium 0 of system (5) is globally asymptotically stable, while if 0 > 1, the disease-free equilibrium 0 of system (5) is unstable.
Proof. The Jacobian matrix of system (5) at the disease-free equilibrium 0 is ) .
In the following, we study the global stability when 0 ≤ 1.
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. ) .

The HFMD Model with Optimal Controls
In this section, we present the optimal control problem by adding to the model (2) two control functions 1 ( ) and 2 ( ).
The aim is to find the optimal values * 1 and * 2 of the controls 1 and 2 , such that the associated state trajectories 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 , = 1, 2. The controls are bounded between 0 and 1. We consider state system (25) of ordinary differential equations in 5 with the set of admissible control functions given by The objective functional is given by where the constants 1 and 2 are a measure of the relative cost of the interventions associated with the controls 1 and 2 , respectively. We consider the optimal control problem of determining ( ( ), ( ), ( ), ( ), ( )), associated with an admissible control pair ( * 1 , * 2 ) ∈ on the time interval [0, ], satisfying (25) and the initial conditions (0), (0), (0), (0), and (0) and minimizing cost functional (27); that is, Theorem 5. There exists an optimal control pair ( * 1 , * 1 ) such that 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 [22].
According to the Pontryagin Maximum Principle [23], 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 , = 1, 2, 3, 4, 5, are the adjoint variables.
with zero transversality. To get the characterization of the optimal control given by [23], we solve the equations on the interior of the control set: Using bounds on the controls, we obtain the desired characterization.

Numerical Results and Discussion
In this section, the numerical simulation results of the opti-  Figure 2 shows the control variables 1 and 2 when 0 > 1. By Figure 2, we see that, to minimize the number of infectious and latent individuals, the control 1 keeps the increasing trend during 5 days; during the remaining 15 days, it decreases to the lower bound. The control 2 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 1 on the basic regeneration number 0 in the following. Because   Figure 3. We obtain Figure 4. The higher the immune loss rate 1 , the greater the number of infections and latent individuals. At the same time, higher immunization loss rate 1 indicates that the immunization control measures 1 are weakened. Accordingly, the treatment control measures 2 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 1 and 2 are implemented simultaneously, the blue dots indicate that only control measure 1 is implemented, and the black dotted line indicates that only control measure 2 is implemented. The numerical simulation results show that the number of infections is minimal when the control measures 1 and 2 are implemented simultaneously. If only one control was implemented, the treatment control 2 would be more effective than vaccination control 1 in controlling the number of infectious and latent individuals. The trend of the controls 1 and 2 is displayed in Figure 6 under different control measures.

Data Availability
No data were used to support this study.

Conflicts of Interest
All authors declare that they have no conflicts of interest.