A realistic model for the periodic dynamics of the hand-foot-and-mouth disease

1.   Introduction

Hand-foot-and-mouth disease (HFMD) is a disease that mainly affects children younger than five years of age. Symptoms start within three to six days of contact according to the Center for Disease Control and Prevention (CDC), an initial fever and generalized malaise followed by painful sores in the mouth that blister and ulcerate ^[1].

HFMD causes a major public health problem, especially in Asia. Although most cases of HFMD are mild, there are some severe cases that may develop heart and lung failure and probably leads to death ^[1,2,3]. In the last decade, many severe HFMD outbreaks have occurred in East Asia and Southeast Asia. A vaccine for HFMD was approved and applied in China in 2016. However, several HFMD outbreaks occurred recently ^[4,5]. This guides the decision makers to focus more on effective prevention and control strategies.

A distinct seasonal pattern in tropical and subtropical regions with more than two peaks may occur within a year ^[6,7,8]. Many research studies have confirmed that climatological changes in humidity and temperature have an important effect on the incidence of HFMD despite the inconsistent results from different studies ^[9,10,11]. Another factor that affects the variations in the number of HFMD incidences is the opening and closing of the daycare centers or schools ^[12,13].

Many previous efforts have been done to study the dynamics of HFMD. A number of these studies used models without any vaccination. For example, ^[14] and ^[15] used an SIR model to analyze the epidemic of HFMD in Taiwan and Malaysia respectively and they found that the disease spread quite rapidly and the number of susceptible persons is the only controllable parameter. Moreover, some studies for the dynamical behaviors of HFMD considered more complicated models including SEIQR models with and without vaccination ^{[16,17,18,19,20,21,22,23]}. For example, Zhu ^[23] used a SEIQR with sinusoidal periodic contact rate, to investigate the spread of seasonal HFMD in Wenzhou and he found that HFMD becomes an endemic disease in Wenzhou. Some other previous studies ^{[17,18,24,25,26,27,28]} used the sinusoidal function of period one year ($\beta(t) = \beta_0 + \beta sin(\omega t + \phi)$) for simulating the seasonal varying transmission rate. Unfortunately, they failed to simulate the one year multi-peaks pattern which appeared in the reported data of HFMD.

However, our study focuses on analysing and simulating the effect of the above factors and clarify their relations with the observed multi-peaks, (in a one year) pattern. We study an SEIQRS model and consider these factors in the dynamics of HFMD. We try to determine the effect of the periodic transmission and vaccination rates on the spread of HFMD.

In this paper, we formulate an SEIQRS model with seasonal contact rate and periodic vaccination strategy. We find that this model has no disease free equilibrium point. But a disease free periodic solution (DFPS) is possible for this model. The reproduction number $R_0$ is derived. An upper bound $R^{sup}_{0}$ and a lower bound $R^{inf}_{0}$ of the reproduction number $R_0$ are derived as well. This DFPS is globally asymptotically stable if $R^{sup}_{0} < 1$ and unstable when $R^{inf}_{0} > 1.$

In our simulation we use a novel form of periodic contact and vaccination rates with period one year as follows:

and

We designed these functions to fit the multi-peaks pattern of the real data. Our new form of periodic functions shows a nicer fitting with the real data, compared with the previous simulation results ^[23]. These more realistic simulations give a wider insight to the future of the spread of the HFMD. We plot and compare our obtained simulation results with both previous works and real data. On the other hand, our simulation results confirm our analytical results obtained in this paper.

2.   The model

In the paper, we shall study a more realistic SEIQRS model with periodic vaccination and loss of immunity. Periodic vaccination gives a more applicable realistic strategy to control the dynamics of the disease.

We assume that the population is mixed uniformly and homogeneously with a constant population size $A.$ The population number $A$ is divided into five groups, $S(t)$, $E(t)$, $I(t)$, $Q(t)$, and $R(t)$ which represent the susceptible, exposed, infected, quarantined, and recovered individuals respectively. Figure 1 illustrates the flowchart of the dynamics of the spread of HFMD. This figure describes the flow of inputs/outputs into each of the five compartments of our model. A proportion of the susceptible individuals who have received the vaccine move to the recovered group, the rest get infected and join the latent period or died for any reason. on the other hand, there is a flux of newborns entering the susceptible class. The latent population in their latency period do not show any symptom. Some of the latent died naturally but after about $3 \sim 7$ days the latent individuals become infectious. A part of the infectious people will be hospitalized for treatment and isolated then they get recovered or died. Another proportion of those infectious moves to the recovery class, while some have died. Finally, some of the recovered lose their immunity and join the susceptibles again, another proportion died.

The model parameters are defined and listed as follows :

Note that, $\beta (t)$ and $r(t)$ are time dependent bounded seasonally varied functions with period one year and all parameter values are non negative. Using our assumptions and following ^[9,18,19] we can write our $SEIQRS$ model as a set of five coupled non-linear ordinary differential equations as follows:

Equations (2.1)–(2.5) represent a nonlinear first order system of differential equations that describes the spread of HFMD disease.

The disease transmission rate $\beta(t)$ and the vaccination rate $r(t)$ are nonconstant, positive, and bounded periodic functions. These periodic functions $\beta(t)$ and $r(t)$ are designed to imitate the seasonal variation in the reported cases for HFMD infectious disease due to climatological changes or opening and closing schools ^[12,13].

As we assume that, the vaccination rate $r(t)$ is a nonconstant function therefore, the system (2.1)–(2.5) could not have any equilibrium point. On the other hand, a disease free periodic solution (DFPS) is still possible when $E (t) = I (t) = Q (t) = 0$.

Define the region $D \subseteq R^5$ where,

The five dimensional set $D$, is obviously positively invariant, as the right-hand side of the system (2.1)–(2.5) is differentiable.

3.   Disease free periodic solution (DFPS)

In our model, the vaccination rate $r(t)$ is a periodic function, this force the system (2.1)–(2.5) to have no equilibrium point. On the other hand we expect that, there is a disease free periodic solution (DFPS) for the system (2.1)–(2.5) at $(\hat{S}(t), 0, 0, 0, \hat{R}(t))$. Similar to ^[29,30,31] this DFPS is given by

with the initial condition $S(t_0)\in \mathbb{R^{+}}$, integrating equation (3.1) we get,

So,

Equation (3.3) represents a recurrence relation between the susceptibles at times ($t_0 + n T; n = 1, 2, ....$). Now define the mapping ${G}$ such that

where $S_n = S(t_0+n T)$. If $S_{n1} \neq S_{n2}$ we have that,

Therefore, the mapping ${G}$ is a contraction mapping and it has a unique fixed point $S^*$ ^[32]. Moreover $S^*$; is depending on $t_0$; such that,

Hence, $S^*(t_0 +T) = S^*(t_0)$. So, $S^*$ is a periodic function of t. From (3.4), we have that, $S^*(t_0)$ is continuously differentiable with respect to $t_0$. On the other hand, we have

$\text { where } r_{min} = \min _{t \in[0, T]} r(t)$. Substituting (3.4) and (3.5) into (3.2), we obtain

Hence, $\hat{S}(t) = S^*(t), \hat{E}(t) = \hat{I}(t) = \hat{Q}(t) = 0$, and $\hat{R}(t) = R^*(t) = A -S^*(t)$ is a disease-free periodic solution of the system (2.1)–(2.5).

4.   The basic reproduction number

The basic reproduction number $R_0$ is defined as the average number of new cases produced by a single infected individual entered into a pure susceptible population. Now, we can drive a form of the basic reproduction number $R_0$ of the system (2.1)–(2.5) using a similar scenario to ^[33,34]. By linearizing the system (2.1)–(2.5) around $(\hat{S}(t), 0, 0, 0, \hat{R}(t))$, we find that:

which can be written in the matrix form as:

where,

Now, let $y(t) = \beta(t)\hat{S}(t)$. Note that $y(t)$ is a nonzero positive periodic function with period $T$, so $\overline{y(t)} = \frac{1}{T}\int_0^T\beta(\tau)\hat{S}(\tau) d\tau$ is average of $Y(t)$ over the time $T$. Now, arguing similar to ^[33], we can obtain the value of the basic reproduction number $R_0$ as follows:

where $\hat S (t)$ is a periodic solution of period $T$ and it is the only disease free solution of the system of differential equations (2.1)–(2.5).

Now, we can define $R_{0}^{sup}$ as an upper bound and $R_{0}^{inf}$ a lower one for the value $R_0$ respectively where,

and

5.   Stability of the (DFPS)

In this section, we extend the results obtained by ^[29,30,31]. We show first that if $R_{0}^{sup} < 1$ the DFPS $(\hat{S}(t), 0, 0, 0, \hat{R}(t))$ is globally asymptotically stable. In another section we shall prove that the DFPS is unstable when $R_{0}^{inf} > 1$.

5.1.   Stability of the DFPS when $R_{0}^{sup} < 1$

Now, we show that the DFPS is globally asymptotically stable when $R_{0}^{sup} < 1.$ We use an argument similar to that used in the proof of Theorem 1 ^[31] with some modifications to fit our SEIQRS model.

Lemma 1.

Proof

From Equation (2.1)

As $(\hat{S}(t), 0, 0, 0, \hat{R}(t))$ is a solution of the system (2.1)–(2.5) when $E(t) = I(t) = Q(t) = 0$ then arguing similar to ^[31] we deduce that,

and the proof of Lemma 1 is completed.

Now we go directly to the first main result which is the global stability of the disease free solution when $R_{0}^{sup} < 1$. Theorem 1 proves that DFPS is globally asymptotically stable (GAS) when $R_{0}^{sup} < 1$. Following Theorem 1 in ^[31], we give proof of the next Theorem.

Theorem 1. The DFPS, $(\hat{S}, 0, 0, 0, \hat{R})$ is GAS if $R_{0}^{sup} < 1$,

Proof

From Equations (2.3) and (2.4) we have

and

arguing similar to ^[31] and ^[35] we can easily find that

and

The idea behind this proof is simple as that, from Lemma 1 we have that, given $\epsilon > 0$ there exists $t_1$ such that $S(t)\le \hat {S} (t) + \epsilon$, $I(t) \le I^\infty + \epsilon$ and $Q (t)\le Q^\infty + \epsilon$ for all $t\ge t_1$.

Using Eq (2.2), $E(t)$ can be bounded above, then we use this upper bound to show that $E^\infty = 0$. Now similar to Theorem 1 ^[31], suppose that $E^\infty > 0.$ By integrating Eq (2.2) we have that for $t\ge t_0 > t_1,$

Now, recall that, $y(\eta) = \beta(\eta)\hat{S}(\eta)$. Therefore,

Suppose that $(k+1)T\ge t-t_0 \ge kT,$ then we have,

So, for $t \geq t_0$ we find that

Choose $t_2 > t_0$ large enough so that for $t\geq t_2$; $A e^{-(\xi +c_{2})(t-t_0)} < \epsilon$: Then for $t \geq t_2$,

where $\beta_{\max } = sup_{u\in[0, T ]} \beta(u)$. Now choose $\epsilon$ small enough so that,

where $R_{0}^{sup}+\psi < 1$ and $\psi > 0.$ Hence for $t \geq t_{2},$ we find that $E(t) \leq\left(R_{0}^{sup}+\psi\right) E^{\infty}.$ Thus, $0 \leq E^{\infty} \leq\left(R_{0}^{sup}+\psi\right) E^{\infty}.$ then $, E^{\infty} = 0.$ Hence also $I^{\infty} = 0$ this leads to $Q^{\infty} = 0$ Moreover $E(t) \rightarrow 0$, $I(t) \rightarrow 0$ and $Q(t) \rightarrow 0$ as $t \rightarrow \infty$. Finally we go to show that $(S(t)-\hat{S}(t)) \rightarrow 0$ and $(R(t)-\hat{R}(t)) \rightarrow 0$ as $t \rightarrow \infty.$ As $\hat{R}(t)$ is a solution of Eq (2.5) when $E(t) = I(t) = Q(t) = 0,$ we find that,

Given $\epsilon_{1} > 0$ using Lemma $1,$ there exists $t_{3}$ such that $E(t)+I(t)+Q(t) \leq \epsilon_{1}$ and $S(t) \leq S(t)+\epsilon_{1}$ for all $t \geq t_{3}.$ So for $t \geq t_{3}$

where $r_{max} = \sup _{u \in [0, T]} r(u).$ By integrating this inequality, we have that

Now we can show that, given $\epsilon_{2} > 0,$ there exists $t_{4}$ such that $R(t)-\hat{R}(t) \leq \epsilon_{2}$ for $t \geq t_{4}.$ Therefore, $S(t) = A-R(t)-Q(t)-I(t)-E(t) \geq A-\hat{R}(t)-\epsilon_{1}-\epsilon_{2}$ for $t \geq t_{4}.$ From Lemma $1,$ we find that $S(t) \rightarrow \hat{S}(t)$ as $t \rightarrow \infty.$ Since $R(t) = A-S(t)-I(t)-Q(t)-E(t),$ then we must have, $R(t) \rightarrow A-\hat{S}(t) = \hat{R}(t)$ at $t \rightarrow \infty.$ This completes the proof of Theorem $1.$ Thus if $R_{0}^{sup} < 1,$ then DFPS is GAS.

5.2.   Instability of the DFPS when $R_{0}^{inf} > 1$

Here $R_{0}^{inf}$ is a lower bound for $R_0$. We use an argument consisting of a mixture of those used in ^[30], ^[36], and ^[37]. So, we can show that if $R_{0}^{inf} > 1$ the (DFPS) is not stable and the disease will fire up.

The first case when the infection is not present initially, $I(0) = E(0) = Q(0) = 0$ and $S+R = A$, then by Theorem 1, the system (2.1)–(2.5) tends to the DFPS as $t\rightarrow \infty$ whatever the value of $R_0.$ The second case is that, all of $E(0) > 0$, $I(0) > 0$, and $Q(0) > 0$ we shall prove that DFPS is unstable by using an argument similar to ^[30] and ^[31].

Theorem 2. If $R_{0}^{inf} > 1$, then the DFPS is unstable.

Proof

In this case, we suppose that the DFPS is stable when $R_{0}^{inf} > 1$ and we get a contradiction with our assumption. So if $(S, E, I, Q, R)$ starts near $(\hat S, 0, 0, 0, \hat R)$ we have that $(S, E, I, Q, R)$ must stay near $(\hat S, 0, 0, 0, \hat R)$ for the time being large. Therefore $(S, E, I, Q, R) \to (\hat S, 0, 0, 0, \hat R)$ as $t \to \infty.$ So, given $\varepsilon = (1/2)(1-(1/R_{0}^{inf})) > 0$ there exists $t_1 > 0$ such that $|S(t)-\hat S | \le \varepsilon \hat S,$ $|E(t)-0| \le \varepsilon,$ $|I(t)-0 | \le \varepsilon$, $|Q(t)-0| \le \varepsilon$ and $|R(t)-\hat R | \le \varepsilon$ for all $t \ge t_1$. Choose $\varepsilon_1 > 0$ such that $(R_{0}^{inf}+1)(1-\varepsilon_1) > 2.$

From Eq (4.6) recall $R_{0}^{inf}$,

Therefore, there exists an integer $k$ such that

Choose $t_2\ge t_1$ such that, $E(t_2) \ge E(0) e^{-(\xi+c_{2})t_2} > 0$, $I(t_2) \ge I(0) e^{-(\xi+c_{3}+c_4)t_2} > 0$ and $Q(t_2) \ge Q(0) e^{-(\xi+c_{5})t_2} > 0$ as $E(0) > 0$, $I(0) > 0$ and $Q(0) > 0$. Now define $\varepsilon_2$ and $\eta$ as follows:

By continuity $\eta > 0$ and if $\eta < \infty$ then $I(t_2 + \eta) = \varepsilon_2$, $Q(t_2 + \eta) = (c_3/(\xi+c_{5}))\varepsilon_2$ or $E(t_2 +\eta) = { ((\xi+c_{3}+c_4)/ c_2) \varepsilon_2 }$. This leads to a contradiction.

Now from Eq (2.3) and by the definition of $\varepsilon_2$ we find that,

similarly, from Eq (2.4) we find that,

but on the other hand from Eq (2.2) we find that

Now

If $\eta\le kT$ then

If $\eta\ge kT$ then

As $R_{0}^{inf} > 1$ and $\varepsilon = \frac{(R_{0}^{inf} - 1)}{2R_{0}^{inf}} > 0$ then using Eq (5.9)

Hence $\eta < \infty$ leads to a contradiction so $\eta = \infty$ and $I(t_2 +\delta)\ge \varepsilon_2$, $Q(t_2 + \delta)\ge (c_3/(\xi+c_5))\varepsilon_2$ and $E(t_2 + \delta)\ge ((\xi+c_3+c_4)/c_2)\varepsilon_2$ for all $\delta\ge 0.$ This contradicts the fact that the trajectory tends to the DFPS. Hence the DFPS cannot be stable for $R_{0}^{inf} > 1,$ and the proof is completed.

These results prove that the disease free solution for our SEIQRS model represented by the system (2.1)–(2.5) with a nonconstant periodic transmission and vaccination rates $r(t)$ in $[0, T]$, is globally asymptotically stable if $R_{0}^{sup} < 1$ and not stable if $R_{0}^{inf} > 1$. from these results we deduce that,

is the sufficient condition to keep the DFPS globally asymptotically stable.

The obtained results, for our SEIQRS model with nonconstant periodic vaccination rate $r(t)$, extend the corresponding results obtained by ^[30], ^[31], ^[37], and ^[35]. These results are original for an SEIQRS model with periodicity in the transmission and vaccination rates.

6.   Simulation results

Now we look numerically at the behavior of the system (2.1)–(2.5). The software package Mathematica12 is used to solve our system of nonlinear ordinary differential equations (2.1)–(2.5). Real data about HFMD in Egypt dose not available. So we use data taken from the district of Wenzhou China (Table 1 ^[23]). We choose Wenzhou China as a subtropical area that is similar to the climate of Egypt. Simulations are performed using parameters from the literature ^[23]. We compare our simulation results with the reported data of HFMD in Wenzhou from March 2010, to December 2013. Also, we compare our results with the simulation results of ^[23].

The parameters estimated in ^[23] are used in our simulation to give a fair comparison with their results. The parameters are set to the following values:

As stated before our novel periodic transmission rate $\beta(t)$ is designed to fit the yearly multi-peaks pattern in the reported data. Recall $\beta(t)$ as follows:

By using the least-square fitting of the numbers of $I(t)$, we can estimate the values of $\beta_{ij}: i, j = 0, 1, 2.$

Figure 2 shows that our simulation results are often closer to the real data than the results obtained by ^[23]. Clear multiple peaks, within a year, are appeared in our simulation which gives a good fitting to the real data. These results show that our model is more realistic than ^[23]. In Figure 2 we use the parameters' values of the transmission periodic function which are estimated from the least square fitting of $I(t)$ during the period from 2010 to 2013, as follows: $\beta_{10}$ = 0.64, $\beta_{11}$ = 0.665, $\beta_{12}$ = 0.0035, $\beta_{20}$ = 0.63, $\beta_{21}$ = 0.445, $\beta_{22}$ = 0.0017.

Figure 3 shows that there is a decrease in the height of the peaks of the number of infected persons $I(t)$ over time. However, this solution has a one year multi-peaks similar to those of the reported real data. Figure 3 predicts that the number of the new HFMD cases will fluctuate around a decreasing level.

Figure 4 plots the number of infected against time in weeks when $R_{0} > 1$. This figure shows that our simulation results are always very close to the reported data. Our simulation results give a more precise fitting than the results of ^[23]. Our model produces a multi-peaks, in one period, solution while the model in ^[23] has a mono peak carve. In Figure 4 the parameters of the transmission periodic function are given estimated values from the least square fitting of $I(t)$ from 2011 to 2012, as follows: $\beta_{10}$ = 0.49, $\beta_{11}$ = 0.94, $\beta_{12}$ = 0.0031, $\beta_{20}$ = 0.36, $\beta_{21}$ = 0.70, $\beta_{22}$ = 0.00065.

Figure 5 shows that there is a stable endemic periodic solution of our model, with a period of one year, when the basic reproductive number ($R_{0} > 1$). This means that the HFMD will take off and becomes endemic when $R_{0}$ is larger enough than one in value.

Now, when we apply the vaccination strategy of the form:

Where the parameters of the vaccination function $r(t)$ are given the following values: $\rho_{10}$ = 0.00098, $\rho_{11}$ = 0.00188, $\rho_{12}$ = 0.0031, $\rho_{20}$ = 0.00072, $\rho_{21}$ = 0.0014, $\rho_{22}$ = 0.00065. In this case, we find that the susceptibles and the infected populations tend to the DFPS as shown in Figure 5.

Figure 6 (a) shows that in the absence of vaccination the HFMD fire up when $R_{0} > 1$. While applying our vaccination strategy forces the disease to die out and (b) show that without vaccination the susceptible population stays high enough to provide a flow of new case which give the disease the chance to persist in the population. On the other hand, applying our vaccination strategy forces the susceptibles to be small enough to fluctuate close to the DFPS and keep $R_{0}$ to be less than one.

7.   Conclusions

We study an SEIQRS model for the dynamics of HFMD. Our model introducing a seasonality in both the contact and vaccination rates. In section 1 we give a brief introduction to the dynamics of HFMD as a public health problem. In section 2 we formulate our model and define the model parameters. The exitance of the DFPS is proved in section 3 and the formula of the basic reproduction number $R_0$ is derived in section 4. We define an upper bound of $R_{0}^{sup}$ and a lower bound $R_{0}^{inf}$ in this section as well. We use these lower and upper bounds to prove that, the disease free periodic solution DFPS, is globally asymptotically stable when $R_{0}^{sup} < 1$ and unstable if $R_{0}^{inf} > 1$. Section 5 contains the stability results. We use our novel form of periodic contact and vaccination rates during all of our simulations in section 6. Previous works used a sinusoidal form of periodic transmission rate which is failed to simulate the yearly multi-peaks pattern in the reported data. We claim that our suggested transmission periodic function forces our model to have a multi-peaks pattern. This result gives our model more realistic behaviours and produces more accurate simulation results. Our simulation results are original for this kind of periodic functions.

Future work

In our study, the population is assumed to be homogeneously mixed and has a constant size. This assumption is unrealistic, so in future work, one can use network epidemic models to simulate the fact that, peoples have different chances to catch the disease ^[38,39,40]. Another future work is to study the behaviour of the endemic solution when the reproduction number $R_0 > 1$.

Acknowledgments

This project is funded by the Academy of Scientific Research and Technology (ASRT), Egypt, Grant No(6739), (ASRT) is the 2$^{nd}$ affiliation of this research.

Conflict of interest

The authors declare that they have no competing interests.