A reaction-diffusion HFMD model with nonsmooth treatment function

Hand, foot, and mouth disease (HFMD) is a contagious viral illness that commonly affects infants and children. In some areas with high incidence of this disease, the relevant departments often use some strategies to strengthen treatment when the number of infected individuals exceeds a certain threshold. To assess the effectiveness of strengthening treatment strategies which depend on a certain threshold, we propose a new reaction-diffusion model with nonsmooth treatment function to investigate the spread of HFMD. In the case of the spatial domain being bounded, by defining the basic reproduction number R0, we use Lyapunov theory to prove that the disease-free equilibrium is globally asymptotically stable as R0 < 1, and the positive equilibrium is globally asymptotically stable as R0 > 1. In the case of the spatial domain being linear and unbounded, in order to study how the movement of children impacts the spatial spread of HFMD, we further consider the traveling waves. Finally, numerical simulations demonstrate the effectiveness of the theoretical analysis.


Introduction
HFMD, caused by multiple intestinal viruses of the Picornaviridae family, is a common infectious disease among infants and children [1,2]. Since HFMD was first reported in New Zealand in 1957 [3], it has spread rapidly all over the world [4][5][6][7][8]. Then, the transmission mechanism of HFMD has become a common concern in recent years.
Dynamic models are often used to solve practical problems [9][10][11] such as compartment dynamics models for analyzing the spread of infection [2,12]. Wang and Sung (2006) [12] and Tiing and Labadin (2008) [2] first applied the classical SIR model to the epidemic analysis of HFMD in a specific area. On the basis of the above, many factors affecting HFMD, such as isolation measures [13][14][15][16][17][18][19], enterovirus contaminated environment [17,20,21], seasonal variation [18,19,21,22], are considered in the dynamic modeling. Recently, we analyzed the effects of EV71 vaccination on HFMD in mainland China from 2016 to 2017 [23]. It should be noted that all of the above studies were based on time evolution. To the best of our knowledge, fewer mathematical models for HFMD were based on spatial evolution. The spread of infectious diseases is inseparable from the activities of the hosts in Figure 1 The treatment function T(I, I c ) in (1.1) spatially environments [24,25]. In addition, the research results in [26,27] showed that human spatial behavior has an important impact on the spread of HFMD. For instance, people with more active activities are more likely to contract diseases, such as boys than girls, scattered children than school children. Recently, having considered the impact of the spatial impact of the epidemic, we proposed a reaction-diffusion HFMD model [28].
The treatment is an important method to prevent and control the spread of diseases. In general, sudden changes in treatment strategies often lead to nonsmooth phenomena. Therefore, based on the previous studies in [29][30][31][32], we have established that the spatially nonsmooth (continuous but not differentiable) treatment function needs to meet two conditions (A1), (A2) [28]. Based on conditions (A1), (A2), we seek intensive treatments when the number of infected individuals exceeds a given threshold, denoted by I c . Then we give the following treatment (see where I := I(x, t) represents the number of infectious individuals at location x and time t. k 1 , k 2 are positive constants, k 2 > k 1 , and k 1 I, k 2 I represent general treatment and intensive treatment at x ∈ for t ≥ 0, respectively. In this paper, we propose a reaction-diffusion HFMD model with the nonsmooth treatment function (1.1) to investigate the effect of the movement of individuals on the spatial spread. Assuming that the spatial domain is bounded, we analyze the dynamic behaviors of the equilibria for our model. If it is assumed that the spatial domain is linear and unbounded, we prove the nonexistence of traveling waves under certain conditions. The organization of the paper is as follows: In Sect. 2, a novel reaction-diffusion model is built. In Sect. 3, the dynamic behaviors of the equilibria of the model are analyzed theoretically. In Sect. 4, a traveling wave solution of the model is discussed. In Sect. 5, simulations of the model and some prevention and control measures are performed. In Sect. 6, we discuss and summarize our conclusions.

Model formulation
In this section, a reaction-diffusion model with the strengthening treatment function (1.1) is proposed. The underlying structure of the model comprises the following classes of individuals: susceptible S(x, t), infectious I(x, t), infectious and hospitalized or quarantined or isolated Q(x, t), and recovered R(x, t). The model is represented by the following ordinary differential equations: where is the recruitment rate, β is the contact transmission rate, ζ is progression rate leaving the children group below six year old, m 1 , m 2 are the per capita disease induced death rates, and χ is the proportion of hospitalized individuals, T(I, I c ) is defined in Sect. 1. All the parameters of model (2.1) are positive. Due to the biology meaning of (2.1), we assume that the initial condition of (2.1) satisfies In the case of the spatial domain being bounded, we take into account the distribution of the individuals in spatial locations within a fixed bounded domain = [0, π] with smooth or Neumann boundary conditions at time t ∈ [0, ∞). The Neumann boundary conditions are where ∂ ∂n represents the differentiation along the outward normal n to ∂ . (b) In the case of the spatial domain being unbounded, we assume that = R for t ∈ [0, ∞).

Basic properties
In this section, we discuss the existence of solution for system (2.1). We note that Q, R do not appear in the first and second equations, which means that Q, T are decoupled from other equations of system (2.1). Here, we only need to discuss the dynamics for the following subsystem for variable S, I of system (2.1): with the Neumann boundary conditions and the initial conditions As a similar proof process of Theorem 3.1 in [28], we give the following result.
Hence, according to Theorem 3.1, we easily have the following result.

Equilibria
In this section, we first consider the distribution of equilibria for system (2.1). An equilibrium of system (2.1) satisfies  From (4.1), system (2.1) always has a disease-free equilibrium In order to find other equilibria, for system (2.1), we denote a basic reproduction number and a threshold value Next, we investigate the two cases: (a) 0 ≤ I ≤ I c for x ∈ , t ≥ 0; (b) 0 ≤ I c < I for x ∈ , t ≥ 0.
(a) The case 0 ≤ I ≤ I c for x ∈ , t ≥ 0. When 0 ≤ I ≤ I c , we have T(I, I c ) = k 1 I. Hence, it follows from (4.1) that βSI -(ζ + m 1 )Ik 1 I = 0, By (4.5), we have Therefore, if 1 < R 0 ≤ R c , then E * 1 is a unique positive equilibrium for system (2.1). (b) 0 ≤ I c < I for x ∈ , t ≥ 0. When 0 ≤ I ≤ I c , we have T(I, I c ) = k 2 I -(k 2k 1 )I c , and thus (4.1) becomes Next, it follows from the first equation of (4.6) that Substituting (4.7) into the second equation of (4.6), we have Hence, one has that Therefore, equation (4.8) has two real roots: It is clear that I 1 > 0 and I 2 < 0. For biological meaning, I 2 is meaningless. Hence, letting It follows from I * 2 > I c that After some complex computations, it follows from (4.9) that Hence, if R 0 > R c , then E * 2 is a unique positive equilibrium for system (2.1).

The stability of equilibrium
In this section, we analyze the stability of equilibria. We just need to discuss system (3.1). Based on our analyses on the equilibria of system (2.1), it is clear that system (3.1) has a disease-free equilibrium (S * 0 , I * 0 ) = ( /ζ , 0) with respect to E * 0 , a unique positive equilibrium (S * 1 , I * 1 ) with respect to E * 1 as 1 < R 0 ≤ R c , and a unique positive equilibrium (S * 2 , I * 2 ) with respect to E * 2 as 1 < R c < R 0 . To avoid confusion, we denote (S * 0 , I * 0 ) as E * 0 , (S * 1 , I * 1 ) as E * 1 , and (S * 2 , I * 2 ) as E * 2 for system (3.1). First, we discuss the local stability of equilibria. Assume that E * = (S * , I * ) is an arbitrary equilibrium of system (3.1). Let u 1 = S -S * , u 2 = I -I * . The linearization system of system where i = 1 or 2. Denoting u = (u 1 , u 2 ) T , the linearization system of system (5.1) is Therefore, we finally obtain the characteristic equation as follows: To analyze the local stability for each equilibrium E * , we just need to discuss the distribution of the roots for the characteristic equation (5.3). It follows from (5.3) that For system (3.1), we give the following local stability result.
It is obvious that if R 0 < 1, then (5.5) has two negative roots Hence, all the roots of the characteristic equation (5.6) at the equilibrium E * 0 are negative. That is, the disease-free equilibrium E * 0 of system (5.1) is locally asymptotically stable as R 0 < 1.
(c) At the equilibrium E * 2 , from (5.4), we have the following characteristic equation: It is obvious that b 1 > 0, b 2 > 0, and b 1 b 2 > 0. According to the Hurwitz criterion, one has that E * 2 is locally asymptotically stable as 1 < R c < R 0 . This completes the proof.
Theorem 5.2 (a) If R 0 < 1, then the disease-free equilibrium E * 0 is globally asymptotically stable; (b) If 1 < R 0 < R c , then the positive equilibrium E * 1 is globally asymptotically stable; (c) If 1 < R c < R 0 , then the positive equilibrium E * 2 is globally asymptotically stable.
Proof (a) To show the global stability of E * 0 , we consider the following Lyapunov function: Calculating the time derivative of V (x, t) along the trajectories of system (3.1), we have Therefore, if R 0 < 1, then E * 0 is globally asymptotically stable for system (3.1). (b) When 1 < R 0 < R c , system (3.1) has a unique positive equilibrium E * 1 . Moreover, since R 0 < R c is equivalent to I * 1 < I c , we only need to consider T(I, I c ) = k 1 I. Next, we consider the following Lyapunov function: Then we obtain that Moreover, we also have that  [33], for system (3.1) with the initial condition ϕ(x, 0) = (S 0 (x), I 0 (x)), we obtain that E * 1 is globally asymptotically stable. (c) Similarly, since 1 < R c < R 0 is equivalent to I * 2 > I c , then we consider T(I, I c ) = k 2 I -(k 2k 1 )I c for system (3.1).
Take the Lyapunov function as follows: Calculating the time derivative of V 3 (x, t) along the trajectories of system (3.1), we have 14) It follows that ∂V 3 (x,t) ∂t = 0 if and only if S = S * 2 , I = I * 2 . Therefore, the largest compact invariant set in 2 = {(S, I) | ∂V 2 (x,t) ∂t | (3.1) = 0} is just the singleton E * 2 . By the LaSalle invariance principle [33], if we take the initial condition ϕ(x, 0) = (S 0 (x), I 0 (x)) for system (3.1), then we have that E * 2 is globally asymptotically stable. This completes the proof.
Next, we consider this case R 0 = R c . If R 0 = R c , then system (3.1) has a unique equilibrium E * 1 = (S * 1 , I * 1 ) with I * 1 = I c . We note that system (3.1) is not smooth at I * 1 = I c . In this special case, we consider the following two systems (I < I c ) and (I < I c ): (5. 16) We note that when I * 1 = I c , the equilibrium of system (5.15) is the same as the equilibrium of system (5.16). From Theorem 5.1 and Theorem 5.2, if I * 1 = I c , then it is obvious that system (5.15) and system (5.16) are stable.

Traveling waves and minimum wave speed
In this section, we discuss the traveling waves for system (2.1) without spatial boundary constraints, i.e., = R. Similarly, we only need to discuss system (3.1) without spatial boundary constraints. The purpose is to determine the minimum wave speed c * connecting the disease-free equilibrium to the endemic equilibrium, which may predict the spatial spreading of the disease [34,35]. To discuss the existence of traveling wave solutions of system (2.1), we will use the general results based on Schauder's fixed point theorem in [36].
We have known that system (3.1) has a unique positive equilibrium E * 1 or E * 2 if and only if R 0 > 1. In this section, we assume that R 0 > 1 always holds. Then it follows from Theorem 5.2 that system (3.1) has an unstable disease-free equilibrium E * 0 = (S 0 , 0) and a stable positive equilibrium E * = (S * i , I * i ), i = 1 or 2. Next, we will study the nonexistence of traveling wave solution for system (3.1). The traveling wave solution of system (3.1) connecting E * 0 to E * with the wave speed c > 0 has the following form: Rewriting system (3.1) with the form ξ = xct, we obtain the following system: which admits an equilibriumÊ 0 = (S 0 , 0, 0, 0). Linearizing system (6.3) atÊ 0 , we obtain the following Jacobian matrix: (6.4) and the characteristic equation It follows from (6.5) that Next, we will discuss the distribution of the roots for f 1 (λ) = 0 and f 2 (λ) = 0, separately. Case (i): for f 1 (λ) = 0, i.e., we havē By using the Vieta theorem, equation (6.7) has a positive root and a negative real root as follows: Case (ii): for f 2 (λ) = 0, i.e., -4(ζ + m 1 + k 1 )(R 0 -1).

Lemma 6.1 The following statements are valid:
(a) If c < c * , then (6.5) has two distinct real roots and two complex roots; (b) If c = c * , then (6.5) has four real roots and two of them are equal. One is positive, and the others are negative; (c) If c > c * , then (6.5) has four real roots. One is positive, and the others are negative.
By Lemma 6.1, for the nonexistence of traveling wave solution of system (3.1) as c < c * , we have the following result. Proof According to Lemma 6.1, to ensure the existence of system (3.1) satisfying condition (6.1), we need to make sure that all the roots of (6.5) are real. Otherwise, a spiral solution near E 0 will destroy the positivity of S, I. Hence, from Lemma 6.1, we obtain that, for 0 < c < c * , system (3.1) does not have traveling wave solution. This completes the proof.
For c ≥ c * , since T(I, I c ) is discontinuous, we cannot prove the existence of traveling waves for system (3.1) theoretically. Hence, we will simulate the traveling wave solution for system (3.1) numerically in Sect. 7.2.

Simulations
In this section, we check the above main results with numerical simulations. First, we show the stability of the equilibria of system (2.1) and analyze the influence of threshold I c selection on HFMD control. Next, we simulate the traveling wave solutions of system (2.1) numerically.
In order to investigate the effect of threshold I c on HFMD control, we only consider the case R 0 > 1. Note that the positive equilibrium E * 1 or E * 2 depends on threshold level I c . Therefore, we fix β = 0.000080 (get R 0 = 1.4333) and vary threshold I c to analyze the disease control. In the case 1 < R 0 < R c , taking I c = 60 (R c = 1.4457), Fig. 3 shows the trends of HFMD with enhanced disease control and without enhanced disease control by simulating model (2.1). In the case 1 < R c < R 0 , taking I c = 40 (R c = 1.2971), Fig. 4 shows the trends of HFMD with enhanced disease control and without enhanced disease control by simulating system (2.1). In Fig. 3 and Fig. 4, comparing with without enhanced disease control, the solution I(x, t) of system (2.1) with respect to with enhanced disease control decreases rapidly and reaches steady state earlier. For biological meaning, introducing threshold therapy control can relieve HFMD outbreaks and earlier control of infectious diseases. In addition, in the case 1 < R c < R 0 , Fig. 4 shows that the state variable I(x, t) has the lower stable state in the long term behavior. Furthermore, Fig. 5(a) shows that the lower the threshold I c is, the lower the number of infected individuals is, and Fig. 5(a) shows that if I c < 58, i.e., in the case 1 < R c < R 0 , then the lower the threshold I c is, the lower the variable I * 1 of the positive equilibrium E * 1 is. The conclusion also indicates that improving threshold therapy control, i.e., setting a low threshold level I c , has a positive impact on the spread of disease.
To find some effective prevention and control measures for HFMD by using system (2.1), we further analyze the influence of the key parameter β on the solution of equation system (2.1), where β corresponds to the average transmission ability of an infectious individual. Taking I c = 50 and varying threshold β (the other parameter values and the initial conditions are the same as in Fig. 2 for system (2.1)), Fig. 6(a) shows the influence of β on the basic reproduction number R 0 and the key threshold value R c , and Fig. 6(b) shows the influence of β on the basic reproduction number R 0 the total number of the infectious individuals at time t, i.e., I(x, t) dx in terms of different β. In Fig. 6(a), R 0 , R c decrease with decreasing the contact transmission rate β, in which lower values of R 0 and R 0 means lower spread of the epidemics (see Fig. 6(b)). Moreover, when the contact transmission rate β < 0.00052, we obtain R 0 < 1, which means the disease is extinct (see Theorem 5.2). Therefore, we should improve health-care education, such as washing hands after using the toilet and before meals, making air fresh indoors, and so on, to reduce the viral transmission capacity.

Simulations for traveling wave solutions
In this section, by using similar numerical methods as those in [35,36], some numerical simulations of the traveling wave solutions of system (2.1) are presented to support the analytic results obtained above. We note that traveling wave solution is a global concept, so there is no way to give the simulated images at infinity. However, it is feasible to give the variation trend of the traveling wave solution of system (2.1) in the bounded region. Therefore, we only give numerical simulations on the local area. Taking = 160, ζ = 0.01077, m 1 = 0.0001731, m 2 = 0.0001731, χ = 0.05, β = 0.00007, we get R 0 = 1.2541 and the minimum wave speed c * = 0.0092. In the case 1 < R 0 < R c , taking I c = 60 (R c = 1.39) and the following initial data (2.1): If x ∈ [0, π 4 ], then we set S 0 (x) = S * 1 , I 0 (x) = I * , π], then we set S 0 (x) = S * 0 , I 0 (x) = I * 0 = 0, Q 0 (x) = Q * 0 = 0, R 0 (x) = R * 0 . Then the numerical simulations support the existence of traveling wave solutions (see Fig. 7).
Similarly, in the case 1 < R c < R 0 , take I c = 60 (R c = 1.39), fix the other parameters as above and choose the following initial data (2.1): Fig. 8 shows the existence of traveling wave solutions. We investigate the effect of the threshold treatment on the spatial spread of HFMD, i.e., the effect of selection of threshold I c on traveling wave solutions of system (2.1) (see Fig. 9). In Fig. 9, we obtain that the varieties of I c only affect the magnitude of waves but do not change the speed of waves. In addition, Fig. 9 shows that the lower the threshold I c is, the lower the peak value of the wave is. In a word, reducing the value of the threshold I c can reduce the severity of the spread of the disease in space, but it cannot change the speed of the disease transmission.
We investigate the effect of the contact transmission rate β on the spatial spread of HFMD. Figure 10 shows that the varieties of β not only affect the magnitude of waves but also change the speed of waves. Moreover, we also obtain that the lower the threshold β is, the lower the peak value of the wave solution is. From this perspective, to reduce the severity of the spread of disease in space, we should reduce the value of the contact transmission rate β.

Discussion
It has been observed that the strengthening treatment is an important method to prevent and control the spread of HFMD. Moreover, the spread of infectious diseases depends on the movements of the hosts in spatial environment [24,25]. Therefore, we propose a reaction-diffusion HFMD model with a continuous and nonsmooth treatment function to study the effect of threshold-dependent interventions on the spread of HFMD in spatial environment. It is worth emphasizing that, as introduced in the descriptions in Sect. 1, the previous studies for HFMD modeling rarely considered spatial factors, and most used the linear treatment functions which do not contain the threshold-dependent interventions for HFMD. In addition, comparing with other models for other diseases considering the discontinuous treatment functions [29], we do not need to study the dynamics in a transformed Filippov system.
The theoretical analysis of our model is divided into two cases. In the case of the spatial domain being bounded, we have investigated the dynamic behaviors of our model. By defining two critical parameters as the basic reproduction number R 0 and the threshold parameter R c , the main dynamic behaviors of system (2.1) are (see Theorem 5.2): (i) If R 0 < 1, then system (2.1) has a unique disease-free equilibrium E 0 and E 0 is globally asymptotically stable; If R 0 > 1, then the disease-free equilibrium E 0 is unstable; (ii) If 1 < R 0 < R c , then the unique positive equilibrium E * 1 of system (2.1) is globally asymptotically stable; (iii) If 1 < R c < R 0 , then the unique positive equilibrium E * 2 of system (2.1) is globally asymptotically stable.
Moreover, in the case of the spatial domain being linear and unbounded, we further determine the nonexistence and existence of traveling wave solutions. Finally, numerical simulations have supported our theories.