# The Global Dynamics of an Age-Structured Hand-Foot-Mouth Disease Model with Saturation Incidence and Time Delay

Authors: Dongxue Yan; Mengqi Zhang; Jiashan Tang
Journal: Discrete Dynamics in Nature and Society
 (2023)

---

## Abstract

The global dynamic behavior of an age-structured hand-foot-mouth disease (HFMD) model with saturation incidence and time delay is investigated in the work. The time delay occurs during the transition from latent to infectious individuals. Firstly, the model is expressed as an abstract Cauchy problem. The presence of equilibria is then pointed; meanwhile, we define the model’s basic reproduction numberR0. The conclusions show that a threshold of R0=1 can be used to evaluate whether the disease is on the verge of extinction or is still present. When R0<1, the disease-free equilibrium is globally stable. While R0>1, there exists an endemic equilibrium, and the global stability of the endemic equilibrium is also demonstrated. Finally, some numerical examples are given to demonstrate the obtained conclusions.


---

## Body

## 1. Introduction

Hand-foot-mouth disease (HFMD) is a viral infection primarily caused by enteroviruses. Coxsackie Asckievirus (Cox A16) and enterovirus (EV71) are the most prevalent HFMD virus strains [1]. HFMD mostly affects children younger than the age of ten. But older children and adults may also be infected with HFMD [2]. One of the most substantial components of death in children is HFMD, which has become a main public health worry that may have unquantifiable societal and economic consequences. Thus, researches concerns about the spread of HFMD have risen in recent years.Because there is no anti-HFMD medical therapy or vaccination available, a huge number of scholars have built mathematical models to study the transmission dynamics and illness control of HFMD. Mathematical models of HFMD infection have advanced significantly over the last few decades. A simple SIR (susceptible-infected-recovered) model was developed to forecast the number of infectious individuals and the length of the outbreak in [3]. After that, more real factors have been included into this simple SIR model, and more complicated mathematical models were proposed in studying the dynamics of HFMD. For instance, Roy and Halder [4] established an SEIR model in 2010 by dividing the infectious group into asymptomatic (E (t)) and symptomatic (I (t)) infectious individuals. Furthermore, Liu [5] formulated an SEIQR model to investigate the effect of quarantine Q (t) on the disease development. Considering the asymptomatic individuals A and contaminated environments W, Wang et al. [6] discussed a SEIARW system and found that both factors are critical to delay and avoid epidemic outbreaks. Tan and Cao [7] developed a SEIVT epidemic model with vaccination to investigate HFMD’s transmission in 2018. Recently, Sun et al. in [8] proposed a diffusive HFMD model with a fixed latent period and nonlocal infections. Besides, these mathematical models are also used in the research of other infectious diseases, such as COVID-19 [9–11]; coffee berry disease [12]; and banana black Sigatoka disease [13].Very recently, Samanta [14] constructed and described an HFMD infection model with incidence rates β1SI/1+σ1S and β2SI/1+σ2S. Standard epidemiological models almost utilize a bilinear incidence rate βSI [15–17] on account of the law of mass action. There are three types of incidence forms used in epidemiological models when the community is saturated with infectious diseases: proportionate mixing incidence βSI/N [16, 18], nonlinear incidence βSqIq [19, 20], and saturation incidence βSI/1+σS [21] or βSIp/1+σIp [22]. In this paper, the saturation incidence rates βSI/1+σS will be considered.In the transmission of HFMD, individuals who are latently infected have no symptoms and cannot convey the virus to others until several days after exposure (usually 3–7 days [23]). As a result, time delay (latent period) is an important factor in the research of the dynamic behaviors of HFMD systems. Many scholars have investigated delayed HFMD models, for instance [14, 24, 25]. The active HFMD individuals will disseminate the bacterium in the population once the latently infectious persons become active. Thus, infected people’s infectivity is determined by their infection age, which is a significant factor for infectious diseases, particularly the ones with a prolonged latency and fluctuating infection rate like HFMD.Note that as time progresses, the disease develops within the individuals with the different infectivity or many times one may easily observe vary infectiousness of an infected individual at different stages of infection. Particularly in the transmission of HFMD, latent individuals have no symptoms and cannot convey the virus to others. When asymptomatic individuals can infect others (that is to say they have infectiousness), depending on how long they have been infected, which is usually called infection age [26]. Moreover, the infectivity of the host might continuously change with time and infection age. This means that infection age may be one of the informative factors to model for the diseases like HFMD. Hence, age-structured epidemic models should be extensively examined in order to provide better knowledge and further insights into transmission mechanisms. The authors of [26–34] have widely explored the dynamic behavior of age-structured epidemic systems.Inspired by the preceding discussion, in this work, we will investigate the HFMD model with age structure and time delay. LetEa be the infected individuals’ transmission rate at age a. To establish the fundamental form of Ea, we utilize the data (see Table 1) from the China Center for Disease Control (CDC) [35] in 2018 and the least square approach to match the transmission rate Ea at varied infection age a. The numerical simulation shows that the transmission rate Ea of HFMD follows the exponential function 6.275e−0.3873a with a≥τ (see Figure 1), after the shortest period τ necessary during the initial infection to become infectious class. Therefore, when we discuss HFMD, the transmission rate Ea is taken as the following exponential function:(1)Ea=α∗e−κa,a≥τ,0,a<τ,where τ represents the shortest time required from initial infection to become infectious class, Ea∈L+∞0,+∞,R, κ and α∗ are constant coefficients of (1), κ>0 and α∗>0. The equation (1) means that if the infection age a of the infected individuals is less that the time delay τ, the infected individuals do not have the infectivity; otherwise, the transmission rate is α∗e−κa.Table 1
The number of latent HFMD individuals in China.
AgeNumber×10−314.262722.544732.468241.571250.857060.480570.203680.084190.0807100.0395110.0038120.0037130.0036140.0035150.0034Figure 1
The fitting curve of the transmission rateEa of latent HFMD individuals who are 0–15 years old.The form of the age-structured HFMD model with saturation incidence and time delay is investigated as follows:(2)dStdt=Λ−μSt−βStIt1+σSt+rRt,∂Et,a∂t+∂Et,a∂a=−δaEt,a,t≥0,a≥0,dItdt=∫0+∞EaEt,ada−μ+ι1+γ1+qIt,t≥0,dQtdt=qIt−μ+ι2+γ2Qt,t≥0,dRtdt=γ1It+γ2Qt−μ+rRt,t≥0,Et,0=βStIt1+σSt,t≥0,where δa:=μ+Ea, and with the initial conditions(3)S0=S0≥0,E0,a=E0a∈L+10,+∞,I0=I0≥0,Q0=Q0≥0,R0=R0≥0,where St, It, Qt, and Rt represent the number of susceptible persons, infectious persons, quarantined persons, and recovered persons at time t, respectively. The number of latent persons with age a is expressed by t,a. As a result, system (2) is an SEIQR-type HFMD system, and the infection process flowchart is shown in Figure 2.Figure 2
Flow chart of the SEIQR model for HFMD.The parameters in model (2) are considered biologically relevant, and all of them are positive. The biological interpretations of the parameters of the system are described in Table 2.Table 2
The biological meanings of parameters in the system (2).
ParameterBiological meaningEaThe infected individuals’ transmission rate at ageaΛThe susceptible individuals’ continual influx rateβThe ratio at which the infectious individuals infect the susceptible onesμAll individuals’ natural mortalityι1, ι2The mortality due to disease of an infectious and quarantined individualγ1, γ2Infectious and quarantined individuals’s recovery ratesqThe quarantine raterThe loss of immunity rateIt is worth pointing out here that, in comparison to the model in [14], system (2), which considers the infection age and infectious delay, is obviously more general and realistic. The asymptotic behavior of solutions to system (2) is the focus of this work. Firstly, we will determine the basic reproduction number R0, and then study the infection-free equilibrium E0’s asymptotic stability for system (2) when R0<1 and R0>1, respectively, through discussing the position of the corresponding characteristic equations’ eigenvalues. Meanwhile, when R0>1, the global stability of a positive equilibrium E∗ is also discussed by constructing a suitable Lyapunov function. Finally, we are particularly interested in how age-structure, time delay, and saturation incidence affect the dynamics of the system under consideration.The following is the remainder of this work. In Section2, we give some preliminary conclusions and prove the well-posedness of the system (2). In Section 3, the existence of the equilibria, particularly the existence and uniqueness of the endemic equilibrium, as well as the linearized system of (2) surrounding the equilibria are discussed. The stability of the equilibrium is demonstrated, in Section 4. Ultimately, in Section 5, we summarize this paper and provide some numerical simulations to demonstrate our theoretical results.


## 2. Preliminaries and Well-Posedness

For this second, model (2) will be converted as an abstract Cauchy problem (ACP); then we will prove the well-posedness of the system and discuss the nonnegativity and boundedness of solutions. For these intensives, we will first gather some background information about linear operators and C0-semigroup theory, and some notations will be utilized throughout the paper.We start by transforming the system (2) into an abstract evolution equation. Define that(4)X=R×L10,+∞,R×R4,and give the linear operator’s formulation A:DA⊆X⟶X(5)Axyzmn0=−μx−y′−δay−μ+ι1+γ1+qz−μ+ι2+γ2m−μ+rn−y0,with DA=R×W1,10,+∞,R×R×R×R×0. Then DA¯=R×L10,+∞,R×R×R×R×0, which implies that DA is not dense in X. A nonlinear operator L:DA¯⟶X given by the following:(6)Lxyzmn0=Λ−βxz1+σx+rn0∫0+∞Eayadaqzγ1z+γ2mβxz1+σx,and let Tt=St,Et,⋅,It,Qt,Rt,0T. Hence, based on the above notations, system (2) can be rewritten into an abstract Cauchy problem as follows:(7)ddtTt=ATt+LTt,t≥0,T0=T0,with T0=S0,E0a,I0,Q0,R0,0T.Generally speaking, solving an abstract differential problem with a strong solution such as (7) is challenging. As a result, we solve (7) in integrated form as follows:(8)Tt=T0+A∫0tTsds+∫0tLTsds.Set(9)X0=DA¯=R×L10,+∞,R×R×R×R×0,X0+=R+×L+10,+∞,R×R+×R+×R+×0.In Section3, the operator A,DA will be shown to be a Hille–Yosida operator, and thus, it produces a C0-semigroup on the closure of its domain, according to Theorem 3.8. in [36]. Next, the well-posedness result of the model (7) can be obtained.Theorem 1.
InX0+, model (2) has a unique continuous solution. Furthermore, the mapping Φ:0,+∞×X0+↦X0+ defined by Φt,T0=Tt,T0 for T0∈X0+ is a continuous semiflow, i.e., the mapping is continuous and meets the condition Φ0,⋅=I and Φt,Φs,⋅=Φt+s,⋅.From biological interpretation,St, It, Qt, and Rt represent the number of susceptible persons, infectious persons, quarantined persons, and recovered persons at time t, respectively. Et,a represents the number of latent persons with age a. Therefore, only nonnegative solutions St,Et,a,It,Qt,Rt of the system (2) are meaningful. The following result demonstrates that the nonnegative initial value solutions St,Et,a,It,Qt,Rt of (2) remain nonnegative and bounded eventually.Theorem 2.
For anyt≥0, all the solutions of system (2) with nonnegative initial values stay nonnegative and are eventually bounded.Proof.
To begin, by integrating the second equation of (2) along the characteristic line, we can get the following result:(10)Et,a=Et−a,0e−∫0aδθdθ,a≤t,E0a−te−∫a−ttδθdθ,a>t.
Similar to the proof of Lemma 2.2 in [27], it is easy to check that Et,a remains nonnegative for nonnegative initial data.
Next, before provingIt is nonnegative for ∀t≥0, assume that p=μ+ι1+γ1+q and Mt=∫0+∞EaEt,ada. Thus, the third equation of model (2) can be replaced as follows:(11)dItdt=−pIt+Mt.
Then we haveIt=e−pt∫0tMsept−sds+I0≥0, hence, It≥0 for ∀t≥0.
Through the fourth equation of system (2), we have(12)dQtdt=qIt−μ+ι2+γ2Qt≥−μ+ι2+γ2Qt,thus, Qt≥Q0e−μ+ι2+γ2t≥0.
Similarly, from the fifth equation, we can obtain(13)dRtdt=γ1It+γ2Qt−μ+rRt≥−μ+rRt,which indicates that Rt≥R0e−μ+rt≥0.
Then, we demonstrate for∀t≥0, St is nonnegative. If we assume there ∃t1>0 such that St1=0, and St>0 for ∀t∈0,t1. In fact, by the first equation of system (2), there is S′t1=Λ+rRt1>0, which obviously is inconsistent. Thus, St≥0 for any t≥0.
Finally, the boundness of the solutions of model (2) will be proved. The total number of latent persons is represented by E¯t=∫0+∞Et,ada. There is a biologically defined maximum age; hence, lima⟶+∞Et,a=0 is a plausible assumption. From system (2), we deduce(14)St+E¯t+It+Qt+Rt′=Λ−μSt−βStIt1+σSt+rRt+∫0+∞−∂Et,a∂a−δaEt,ada+∫0+∞EaEt,ada−μ+ι1+γ1+qIt+qIt−μ+ι2+γ2Qt+γ1It+γ2Qt−μ+rRt=Λ−μSt−μ∫0+∞Et,ada−μ+ι1It−μ+ι2Qt−μRt≤Λ−μSt−μE¯t−μIt−μQt−μRt=Λ−μSt+E¯t+It+Qt+Rt.
Therefore,(15)limsupt⟶+∞St+E¯t+It+Qt+Rt≤Λμ.
As a result, the omega limit set of model (2) is constrained in the below feasible region:(16)Ψ=S,E⋅,I,Q,R:S,I,Q,R,E⋅≥0,S+∫0+∞Et,ada+I+Q+R≤Λμ,which is obviously a positive invariant for model (2), implying that system (2) is ultimately bounded.


## 3. Equilibria and Properties of Linearization at Equilibria

The nonlinear system (2) will be linearized around the equilibrium solutions in this section. In order to do this, we first consider whether or not the equilibrium exist. System (2) always maintains a disease-free equilibrium E¯0=S0,0,0,0,0 with S0=Λ/μ. For the purpose of finding the nontrivial equilibrium E∗=S∗,E∗a,I∗,Q∗,R∗ of system (2), the following equations are proposed:(17)Λ−μS∗−βS∗I∗1+σS∗+rR∗=0,dE∗ada=−δaE∗a,∫0+∞EaE∗ada−μ+ι1+γ1+qI∗=0,qI∗−μ+ι2+γ2Q∗=0,γ1I∗+γ2Q∗−μ+rR∗=0,E∗0=βS∗I∗1+σS∗.We solve the above equation in (17), assume p=μ+ι1+γ1+q, M=∫0+∞Eae−∫0aδθdθda, and then obtain(18)E∗a=E∗0e−∫0aδθdθ:=E∗0fa,τ,S∗=pβM−pσ,Q∗=qμ+ι2+γ2I∗,R∗=1μ+rγ1+γ2qμ+ι2+γ2I∗,with fa,τ=e−∫0aδθdθ=e−α∗e−κτ/κ+α∗e−κa/κ−μa.On the other hand, by the first equation (17) and conducting some computations, we obtain that(19)Λ−μpβM−pσ−pMI∗+rμ+rγ1+γ2qμ+ι2+γ2I∗=0.Thus,(20)I∗=Λ−μp/βM−pσp/M−r/μ+rγ1+γ2q/μ+ι2+γ2.Note that0<M≤1, assume that C=p/M−r/μ+rγ1+γ2q/μ+ι2+γ2. Because p/M=μ+ι1+γ1+q/M>γ1+q, and r/μ+r<1, γ1+γ2q/μ+ι2+γ2<γ1+q, thus r/μ+rγ1+γ2q/μ+ι2+γ2<γ1+q, hence C>0.We propose the definition of the basic reproduction number as follows:(21)R0≔ΛμβM−pσp=S0β∫0+∞Eae−∫0aδθdθdaμ+ι1+γ1+q−σ.R0 is defined to be the expected number of secondary cases produced in a completely susceptible population by a typical infected individual during its entire period of infection [23]. R0 can also be calculated by the recipe in van den Driessche and Watmough [37]. Moreover, Diekmann et al. [38] have proposed that R0 is the spectral radius of the next-generation matrix.In fact, each phrase inR0 has a distinct epidemiological meaning. e−∫0aδθdθ is the likelihood of a latent person living to the age of a. In addition, ∫0+∞Eae−∫0aδθdθda is the total number of infectious persons produced during a latent individual’s lifetime. 1/μ+ι1+γ1+q is the average infection period. β is the transmission rate of the infectious individual. S0 is the overall number of susceptible individuals. σ is the saturation incidence. As a result, R0 represents new cases generated by the average number of typical infected members during the course of the infection period. As determined by the preceding study, the basic reproduction number of system (2) is R0.Then, from (19), we have I∗=Λ1−1/R0/C. Thus, I∗>0 if R0>1, then system (2) has a single positive solution. As a result of the aforesaid study, we arrive at the following conclusion.Theorem 3.
The model (2) always has a disease-free equilibrium E0=S0,0,0,0,0. Furthermore, when R0>1, a unique endemic equilibrium E∗=S∗,E∗a,I∗,Q∗,R∗ exists.LetSt=xt+S¯, Et,a=yt,a+E¯a, It=zt+I¯, Qt=mt+Q¯, Rt=nt+R¯, where E¯=S¯,E¯a,I¯,Q¯,R¯ is any consistent condition of the model (2), and u˜t=xt,yt,a,zt,mt,nt,0, u¯=S¯,E¯a,I¯,Q¯,R¯,0. System (8) is identical to the following Cauchy problem(22)ddtu˜t=Au˜t+Lu˜t+u¯−Lu¯,t≥0,u˜0=u0−u¯.Through computations, the linearized system of (8) can be obtained as follows:(23)ddtu˜t=Au˜t+DLu¯u˜t,t≥0,u˜0=u0−u¯,in which(24)DLu¯xtyt,aztmtnt0=−βI¯1+σS¯2xt+S¯1+σS¯zt+rnt0∫0+∞Eayt,adaqztγ1zt+γ2mtβI¯1+σS¯2xt+S¯1+σS¯zt.Thus, on X the compactness of the bounded linear DLu¯ is obviously derived.Notice thatΩ:=λ∈C:Reλ>−μ, then the statement below can be proven.Theorem 4.
The operatorA,DA is a Hille–Yosida operator.Proof.
Forϕ,φ,ω,ν,η,ψ∈X,ϕ˜,φ˜,ω˜,ν˜,η˜,0∈DA,λ∈Ω, there is(25)λI−A−1ϕφωνηψ=g1g2g3g4g50⇔λ+μg1=ϕ,g2′=−λ+δag2+φ,λ+μ+ι1+γ1+qg3=ω,λ+μ+ι2+γ2g4=ν,λ+μ+rg5=η,g20=ψ.
Thus, we know(26)g1=ϕλ+μ,g2=e−∫0aλ+δθdθψ+∫0ae∫asλ+δθdθφsds,g3=ωλ+μ+ι1+γ1+q,g4=νλ+μ+ι2+γ2,g5=ηλ+μ+r.
Integrating the second equation (26) with regard to the age variable a and summating all the equations, we get(27)g1+g2L1+g3+g4+g5≤1λ+μϕ+φL1+ω+ν+η+ψ.
Thus, we have(28)λI−A−1≤1λ+μ,forallλ∈Ω,which implies that A,DA is a Hille ̶Yosida operator.
Based on Theorem 1.3 in [36] and Theorem 4, it shows that the operator A+DLu¯ is a Hille ̶Yosida operator. Furthermore, according to Theorem 3.8 in [36], we get two strongly continuous semigroups Ttt≥0 and Ttt≥0 generated by the part of A,DA and A+DLu¯,DA+DLu¯ on X0, respectively.
By the proof of Theorem4, the Hille   ̶Yosida is estimated; thus, we know St≤e−ξt. Moreover, DLu¯St:X0⟶X is compact for any t>0. Since(29)Tt=eDLu¯St,=St+∑k=1+∞DLu¯tk/k!St,which implies Ttt≥0 is a quasi-compact C0-semigroup. Thus, Theorem B.1 in [39] deduces that, for some η>0, eηtTt⟶0 as t⟶+∞ when any eigenvalue of A+DLu¯ has a negative real part.
Now, we may derive the result below, based on the previous arguments.Theorem 5.
The solution semiflowΦt,u0 of the system (2) as given in Theorem 1 satisfies the following characteristics.(i)
The steady-stateu¯ is locally asymptotically stable if any eigenvalue of A+DLu¯ has a strictly negative real part.(ii)
If at least one of the eigenvalues ofA+DLu¯ has a strictly positive part, then the steady state u¯ is not stable.


## 4. Stability of Equilibria

Depending on the preceding analysis, the local stability of disease-free equilibriumE¯0 will be firstly analyzed in this section.Theorem 6.
WhenR0<1, the disease-free equilibrium E¯0=S0,0,0,0,0 of system (2) is locally asymptotically stable. Nevertheless, when R0>1,E¯0 is unstable.Proof.
Letxt=St−S0, yt,a=Et,a, zt=It, mt=Qt, nt=Rt. Linearizing system (23) at E¯0, we get the system as follows:(30)x′t=−μxt−βS01+σS0zt+rnt,∂yt,a∂t+∂yt,a∂a=−δayt,a,z′t=∫0+∞Eayt,ada−μ+ι1+γ1+qzt,m′t=qzt−μ+ι2+γ2mt,n′t=γ1zt+γ2mt−μ+rnt,yt,0=βS01+σS0zt.
Substitutext=x0eλt, yt,a=y0aeλt, zt=z0eλt, mt=m0eλt, nt=n0eλt into the (30), giving the equations(31)λ+μx0=−βS01+σS0z0+rn0,dy0ada=−λ+δay0a,λ+μ+ι1+γ1+qz0=∫0+∞Eay0ada,λ+μ+ι2+γ2m0=qz0,λ+μ+rn0=γ1z0+γ2m0,y00=βS01+σS0z0.
Solving equation (31), we obtain(32)y0a=y00e−∫0aλ+δθdθ,z0=y00λ+μ+ι1+γ1+q∫0+∞Eae−∫0aλ+δθdθda.
Substitutingz0 into the last equation (31) and noting y00 is arbitrary, then(33)βS0/1+σS0λ+μ+ι1+γ1+q∫0+∞Eae−∫0aλ+δθdθda−1=0,or∆0λ,τ≔βS0/1+σS0λ+μ+ι1+γ1+q∫0+∞Eae−∫0aλ+δθdθda−1=0,which is the characteristic equation of model (30) at E¯0.
Obviously, ifλ∈R, ∆0λ,τ is a real function of λ that is decreasing, and ∆0λ,τ satisfies that(34)limλ⟶+∞∆0λ,τ=−1,∆00,τ=R0−1.
Therefore, ifR0>1, the characteristic equation ∆0λ,τ=0 has at least one real positive solution, and the equilibrium E¯0 is not stable. However, when R0<1, then ∆0λ,τ=0 does not have complex roots with non-negative real parts. We assume λ=b+ωi is an arbitrary complex root with b>0; thus,(35)1=∆0λ,τ+1=βS0/1+σS0b+ωi+μ+ι1+γ1+q∫τ+∞Eae−∫0ab+ωi+δθdθda≤βS0/1+σS0b+ωi+μ+ι1+γ1+q∫τ+∞Eae−∫0ab+δθdθda≤βS0/1+σS0b+μ+ι1+γ1+q∫τ+∞Eae−∫0ab+δθdθda=∆0b,τ+1≤∆00,τ+1=R0<1,which is an inconsistency. As a result, the complex root of ∆0λ,τ=0 must have a negative real portion. Thus, when R0<1, E¯0 is locally asymptotically stable.Theorem 7.
WhenR0<1, the disease-free equilibrium E¯0=S0,0,0,0,0 of system (2) is globally asymptotically stable.Proof.
Based on the result of Theorem6, here we only need to demonstrate that for every nonnegative solution St,Et,a,It,Qt,Rt of Model (2),(36)limt⟶+∞St,Et,a,It,Qt,Rt=S0,0,0,0,0.
Indeed, (15) implies that(37)limsupt⟶+∞St≤Λμ=S0.
As a result, for∀ε>0, there is t0 such that St≤S0+ε, for ∀t≥t0, and thus,(38)Et,0≤βS0+ε1+σS0+εIt,fort≥t0.
Then, consider the following linear system:(39)∂E˘t,a∂t+∂E˘t,a∂a=−δaE˘t,a,t≥0,a≥0,dI˘tdt=∫0+∞EaE˘t,ada−μ+ι1+γ1+qI˘t,t≥0,dQ˘tdt=qI˘t−μ+ι2+γ2Q˘t,t≥0,dR˘tdt=γ1I˘t+γ2Q˘t−μ+rR˘t,t≥0,E˘t,0=βS0+ε1+σS0+εIt,t≥0.
We may conclude that system (39) permits a solution of the type(40)E˘t,a,I˘t,Q˘t,R˘t=E˘0aeλ0t,I˘0eλ0t,Q˘0eλ0t,R˘0eλ0t,where E˘0a, I˘0, Q˘0 and R˘0 are positive and λ0 is a root of the corresponding characteristic equation of this system by utilizing the similar way of proving Theorem 6. The formulation (10) of solution Et,a for the second equation in the system (2) yields that Et,a≤E˘t,a, for t≥t0; thus, It≤I˘t, which implies Qt≤Q˘t; therefore, Rt≤R˘t.
That is,(41)Et,a,It,Qt,Rt≤E˘0aeλ0t,I˘0eλ0t,Q˘0eλ0t,R˘0eλ0t,t≥t0.
SinceR0<1, all the eigenvalues of the system (39) have a negative real part, due to the continuous reliance of λ0 on ε, there exists ε>0 s.t. λ0<0, indicating that(42)limt⟶+∞Et,a,It,Qt,Rt=0,0,0,0.
Thus, the first equation in (2) is asymptotic to equation(43)dS˘tdt=Λ−μS˘.
Obviously, there islimt⟶+∞S˘t=S0. Using the asymptotic autonomous semiflow theory (see Corollary 4.3 in [40]), we can find out that(44)limsupt⟶+∞St=S0.
Thus, ifR0<1, then(45)limt⟶+∞St,Et,a,It,Qt,Rt=S0,0,0,0,0.
The demonstration is now complete.
In the following, we primarily concern about the stability of endemic equilibriumE∗=S∗,E∗a,I∗,Q∗,R∗. We linearize system (2) at endemic equilibrium E∗ to investigate its local stability. We do this by introducing the perturbation variables xt=St−S∗, yt,a=Et,a−E∗a, zt=It−I∗, mt=Qt−Q∗, nt=Rt−R∗, which lead to(46)x′t=−μxt−βI∗1+σS∗2xt+S∗1+σS∗zt+rnt,∂yt,a∂t+∂yt,a∂a=−δayt,a,z′t=∫0+∞Eayt,ada−μ+ι1+γ1+qzt,m′t=qzt−μ+ι2+γ2mt,n′t=γ1zt+γ2mt−μ+rnt,yt,0=βI∗1+σS∗2xt+S∗1+σS∗zt.
In order to analyze the stability ofE∗, we seek for solutions of the type xt=x0eλt, yt,a=y0aeλt, zt=z0eλt, mt=m0eλt, nt=n0eλt. Then, the following eigenvalue issue is obtained:(47)λ+μx0=−βI∗1+σS∗2x0+S∗1+σS∗z0+rn0,dy0ada=−λ+δay0a,λ+μ+ι1+γ1+qz0=∫0+∞Eay0ada,λ+μ+ι2+γ2m0=qz0,λ+μ+rn0=γ1z0+γ2m0,y00=βI∗1+σS∗2x0+S∗1+σS∗z0.
We obtain by solving the second and the third equations in (47):(48)y0a=y00e−∫0aλ+δadθ,z0=y00λ+μ+ι1+γ1+q∫0+∞Eae−∫0aλ+δθdθda.
By the fourth and the fifth equations in (47), there are(49)m0=qλ+μ+ι2+γ2z0,n0=γ1+γ2qλ+μ+ι2+γ2z0.
The first equation in (47) yields that(50)x0=−y00λ+μ+rγ1+γ2q/λ+μ+ι2+γ2λ+μz0.
By the last equation in (47), and then combining with (48)–(50), we obtain(51)λ+μ=βI∗1+σS∗2−1+rγ1+γ2q/λ+μ+ι2+γ2λ+μ+ι1+γ1+q∫0+∞Eae−∫0aλ+δθdθda+βS∗/1+σS∗λ+μλ+μ+ι1+γ1+q∫0+∞Eae−∫0aλ+δθdθda.
Thus, the characteristic equation (47) is derived as follows:(52)λ+μ+1−ϒβI∗/1+σS∗2λ+μ=βS∗/1+σS∗λ+μ+ι1+γ1+q∫0+∞Eae−∫0aλ+δθdθda,whereϒ:=rγ1+γ2q/λ+μ+ι2+γ2/λ+μ+ι1+γ1+q∫0+∞Eae−∫0aλ+δθdθda.
Forλ with Reλ≥0, we have ∫0+∞Eae−∫0aλ+δθdθda<M≤1, and rγ1+γ2q/λ+μ+ι2+γ2/λ+μ+ι1+γ1+q<rγ1+q/λ+μ+ι1+γ1+q<r; thus,(53)ϒ≤rγ1+γ2q/λ+μ+ι2+γ2λ+μ+ι1+γ1+q∫0+∞Eae−∫0aλ+δθdθda<r≤1.
Thus, it is easy to demonstrate that the left side of the characteristic (52) holds for λ with Reλ≥0:(54)λ+μ+1−ϒβI∗/1+σS∗2λ+μ>1.
For the right side of (52), when Reλ≥0, we obtain that(55)βS∗/1+σS∗λ+μ+ι1+γ1+q∫0+∞Eae−∫0aλ+δθdθda≤βS∗/1+σS∗λ+μ+ι1+γ1+q∫0+∞Eae−∫0aλ+δθdθda≤βS∗/1+σS∗μ+ι1+γ1+q∫0+∞Eae−∫0aδθdθda=βS∗1+σS∗⋅Mp=S∗S0R0=1,where the notations of p,M,S∗,R0 defined at the beginning of Section 3. Therefore, for λ with Reλ≥0 the left-hand side of (52) is strictly greater than 1, but the right-hand side of (52) is not greater than 1. That is an inconsistency, i.e., the solutions of the characteristic (52) must have a negative real part. As a result, the endemic equilibrium E∗ is locally asymptotically stable, which allows us to derive the following conclusion.Theorem 8.
WhenR0>1, the endemic equilibrium E∗ of system (2) is locally asymptotically stable.As the system (2) is a high-dimensional system, it is hard to estimate the global stability of the endemic equilibrium E∗. However, we can verify the global stability of E∗ for the situation r=0. In this case, the function gx,y and the Lyapunov functional are defined as follows:(56)gx,y=x−y−ylnxy,F=F1+F2+F3,where(57)F1=gS,S∗,F2=∫0+∞BagEt,a,E∗ada,F3=βS∗/1+σS∗μ+ι1+γ1+qgI,I∗,Ba=∫a+∞βS∗/1+σS∗Ea/μ+ι1+γ1+qe−∫0sδθdθds.Using (17), we take the derivative of each portion of the Lyapunov functional F described in (57) along with the solutions of system (2) separately. Firstly, we obtain(58)F1′=1−S∗S0Λ−μS−Et,0=1−S∗S0μS∗+E∗0−μS−Et,0=−μSS−S∗2+E∗0−Et,0−S∗S0E∗0+S∗S0Et,0.By using (10), we can rewrite(59)F2=∫0tBagEt−a,0e−∫0aδθdθ,E∗ada+∫t+∞BagE0a−te−∫a−taδθdθ,E∗ada=∫0tBt−rgEr,0e−∫0t−rδθdθ,E∗t−rdr+∫0+∞Bt+rgE0re−∫rt+rδθdθ,E∗t+rdr.Secondly, using (18) and the fact xgxx,y+ygyx,y=gx,y yield(60)F2′=B0gEt,0,E∗0+∫0tBt−rgEr,0e−∫0t−rδθdθ,E∗t−rdr−∫0tBt−rδt−rEr,0e−∫0t−rδθdθgxEr,0e−∫0t−rδθdθ,E∗t−rdr−∫0tBt−rδt−rE∗t−rgyEr,0e−∫0t−rδθdθ,E∗t−rdr+∫0+∞Bt+rgE0re−∫rt+rδθdθ,E∗t+rdr−∫0+∞Bt+rδt+rE0re−∫rt+rδθdθgxEr,0e−∫rt+rδθdθ,E∗t+rdr−∫0+∞Bt+rδt+rE∗t+rgyE0re−∫rt+rδθdθ,E∗t+rdr=B0gEt,0,E∗0+∫0+∞Ba−δaBagEt,a,E∗ada.Denote thatB0=S∗/S0R0=1, and(61)Ba=δaBa−βS∗/1+σS∗Eaμ+ι1+γ1+q.Then, we have(62)F2′=Et,0−E∗0−E∗0lnEt,0E∗0+∫0+∞βS∗/1+σS∗Ea/μ+ι1+γ1+qE∗a−Et,a−E∗alnEt,aE∗ada.DifferentiatingF3′ along the solutions of system (2) yields(63)F3′=βS∗/1+σS∗μ+ι1+γ1+q1−I∗I∫0+∞EaEt,ada−μ+ι1+γ1+qI=βS∗/1+σS∗μ+ι1+γ1+q∫0+∞EaEt,ada−I∗I∫0+∞EaEt,ada−μ+ι1+γ1+qI+μ+ι1+γ1+qI∗.Adding (58), (62), and (63) together yields(64)F′=F1′+F2′+F3′=−μSS−S∗2−S∗S0E∗0−E∗0lnEt,0E∗0+∫0+∞βS∗/1+σS∗Ea/μ+ι1+γ1+qE∗a2+lnEt,aE∗a−I∗IEt,aE∗ada=−μSS−S∗2+∫0+∞βS∗/1+σS∗Ea/μ+ι1+γ1+qE∗a2−S∗S0−lnEt,0E∗0+lnEt,aE∗a−I∗IEt,aE∗ada.Notice that(65)∫0+∞βS∗/1+σS∗Ea/μ+ι1+γ1+qE∗a1−E∗0SIEt,0S∗I∗da=0.Therefore, we have(66)F′=−μSS−S∗2+∫0+∞βS∗/1+σS∗Ea/μ+ι1+γ1+qE∗aHt,ada,where(67)Ht,a≔lS∗S0+lI∗Et,aIE∗a,lx≔1+lnx−x.For anyx>0, lx≤0 with equality holding if and only if x=1. Thus, Ht,a≤0. Hence, F′≤0. Then, it can be confirmed that the singleton E∗ is the biggest invariant set of F′=0. According to [41], we know that E∗ is globally asymptotically stable, as shown by the compact global attractor A0=E∗.The preceding considerations lead to the following conclusion.Theorem 9.
SupposeR0>1, and r=0, the endemic equilibrium E∗ of system (2) is globally asymptotically stable.Remark 1.
ForR0>1, and r>0, it is difficult to construct a suitable Lyapunov function to verify the global stability of the endemic equilibrium E∗. Here, we take r=0.025,0.05 and use numerical simulation to display the endemic equilibrium E∗. Compared with the case r=0, we find that the endemic equilibrium E∗ is still globally asymptotically stable (see Figure 3). Moreover, it can be seen that with the increasing of, the number of latent and infectious persons increase. Hence, the rate of loss of immunity r has a negative impact on the control of the disease.Figure 3
The solutions of system (2) with different r when R0>1. St represents the susceptible individuals, Et=∫0+∞Et,ada represents the total number of latent individuals, It represents the infectious individuals, Qt represents the quarantined individuals, and Rt represents the recovered individuals.


## 5. Conclusion and Discussions

In this research, an age-structured SEIQR model of HFMD has been developed and studied, with the infection age and time lag. The well-posedness of the system (2) is explored first, and following that the solutions’ positivity and boundedness for system (2) is discussed. Next, the basic reproduction number R0 is derived, which has been shown to be the threshold for determining disease extinction or survival. That is to say, when R0<1, the disease-free equilibrium E0 is globally asymptotically stable; otherwise, the system is unstable. Moreover, we explore the local stability of the infected equilibrium E∗ by analyzing the distribution of roots to the related characteristic equation and the global stability by constructing suitable Lyapunov functions for the epidemic model (2) when r=0, and demonstrate the E∗ is still globally stable when r≠0 by numerical simulation. Thus, by managing the basic reproduction number R0, we may control the spread of HFMD in the population.In the following, we utilize numerical simulations to demonstrate the theoretical results of the system (2) by Matlab software and reveal the influences of age-structure, saturation incidence σ and time delay τ on the disease dynamics. Firstly, according to Chinese real data about HFMD, we choose the parameters in the system (2). The entire population of China in 2020 is 1412120000 and the birth rate is 8.52 per thousand [42]; thus, Λ=1002605.2. And the natural mortality rate of China is 7.07 per thousand, thus μ=0.00707. From [35], the number of infected cases in 2018 is about 2353310. We assume that q=0.45, thus I0=1294321, Q0=1058989. The report shows that the exposure period for HFMD is about 4–7 days; thus, the initial latent population is E0=392218. We just supply an a priori estimate for the initial recovered population R0=4532432. The ratio of the entire youthful population aged from 1 to 14 is 17.9; therefore, S0=252770000. We may derive the death rate of infection ι1=0.000025. The decubation of infection was about 2 weeks [43, 44], that is, the recovery rate of infectious individuals γ1=0.0714. Besides, in all the following figures, we set Et:=∫0+∞Et,ada represents the total number of latent individuals.(1)
The maximum infection age we chose is 100,β=0.0005, ι2=0.000015, γ2=0.075, r=0.025, τ=0.5, θ=0, α∗=0.35 reference from [23]. Accordingly, we obtain R0<1. The numbers of the susceptible persons, the latent persons, the infectious persons, the quarantined persons, and the recovered persons with three different initial values are displayed in Figure 4, implies that the disease-free equilibrium E0 of system (2) is globally asymptotically stable when R0<1. It is compatible with Theorem 6, which indicates that if we control the fundamental regeneration number R0 to be kept below 1, the illness will be eradicated.(2)
WhenR0<1, the impact of delay “τ” on the solutions of system (2) are displayed in Figure 5. We take τ=0,0.5,0.8, respectively, and keep the other parameters the same with Figure 4. Figure 5 states that as the delay τ becomes longer, the time for achieving disease-free equilibrium E0 increases.(3)
UnderR0>1, the stability of positive equilibrium E∗ will be shown next. When we take the parameters β=0.22, and others keep the same with Figure 4, then we get R0>1. Figure 6 illustrates that the solutions of system (2) with three distinct initial values will trend to E∗ as t approaches infinity, revealing that the positive solution E∗ probably be globally asymptotically stable when R0>1.(4)
Fourthly, Figure7 shows the influence of delay “τ” on the solutions of system (2) when R0>1. The delay τ does not influence the stability of E∗, but from Figure 7 we can see that as the delay τ (in the incidence of latent infections) increases, the number of HFMD infections, quarantines, and recoveries decreases.(5)
Accordingly, the distribution of latent individuals with respect to infection age at the endemic equilibriumE∗, E∗a is shown in Figure 8(a), which corresponding to the second solution line in Figure 6, and the distributions with both infection age and time, Et,a is shown in Figure 8(b).(6)
Finally, we examine that the effect of the saturation incidence rate parameters of the system (2) on the dynamics of HFMD when R0>1. The impact of the saturation incidence rate on the solutions of the system (2) are displayed in Figure 9. When σ=0.195,0.2,0.205, respectively, and the other parameters are the same as in Figure 6, then R0>1. From Figure 9, we know that as the saturation incidence rate’s parameters σ increases, the number of the susceptible persons St raises, as well as the number of the latent persons Et, infectious persons It, quarantined persons Qt, and recovered persons Rt all decrease. In summary, the higher the saturation incidence rates’s parameter σ is, the faster the HFMD can be controlled.Figure 4
The graph of disease-free equilibriumE0 of system (2), where St represents the susceptible individuals, Et=∫0+∞Et,ada represents the latent individuals, It represents the infectious individuals, Qt represents the quarantined individuals, and Rt represents the recovered individuals. The three different color curves correspond to three sets of different initial values.Figure 5
The solutions of the system (2) with different τ when R0<1.Figure 6
The stability of endemic equilibriumE∗ of system (2), where St represents the susceptible individuals, Et=∫0+∞Et,ada represents the total latent individuals, It represents the infectious individuals, Qt represents the quarantined individuals, and Rt represents the recovered individuals.Figure 7
The solutions of system (2) with different τ when R0>1.Figure 8
The distributions of latent individuals atE∗. (a) The distributions E∗a at E∗. (b) The distributions Et,a.
(a)(b)Figure 9
The solutions of system (2) with different σ.Comprehensive to is it can know the growth of time delayτ adverse to the formation of disease-free equilibrium E0 when R0<1, and causes the increasing of the number of HFMD infections, quarantines and recoveries decrease when R0>1. Therefore, the increase in the saturation incidence rate σ is benefit for the control of HFMD.



---

*Source: 7090661-2023-01-24.xml*