Assessing the Role of Vaccination in the Control of Hand, Foot, and Mouth Disease Transmission

1. Introduction

This involves incorporating vaccination into the study of Hand, Foot, and Mouth Disease transmission dynamics to understand the effects of different immunization strategies on outbreaks within populations. Vaccination has the potential to significantly lower the basic reproduction number ($R_0$), the number of secondary infections produced by the average infected person. The immunization of a large part of the population susceptible to infection breaks the chain of transmission, which will slow or reduce the spread of the virus. In communities where vaccination is highly prevalent, HFMD outbreaks may remain impressively low due to fewer people susceptible to the virus [1,2]. Besides, vaccination may assist in herd immunity, whereby the unvaccinated populations tend to benefit from the reduced virus circulation. In addition, vaccine efficacy and coverage rates, and even age distribution among the vaccinated population [3], might be other factors influencing the effectiveness of vaccination in HFMD control. This can be accomplished using mathematical modeling and epidemiological studies to determine what different vaccination scenarios would yield an effect on, thereby assisting program designers in developing specific immunization programs that would limit the burden of HFMD in at-risk populations.

The modes of transmission of infectious diseases represent a variety of mechanisms, generally classified into two broad categories: direct and indirect. Direct transmission may occur through actual physical contact, respiratory droplets, or several body fluids to permit the transfer of infectious agents from one host to another. Indirect transmission involves intermediate agents such as fomites (inanimate infected objects), vectors like mosquitoes, and air particles [4,5]. Mathematical modeling further elucidates the complex nature of disease spread and provides critical tools for describing and predicting transmission dynamics. The basic notions include the basic reproduction number, $R_0$, defined as the average number of secondary infections produced by one infected individual in a wholly susceptible population, indicating a pathogen’s contagiousness. SIR and SEIR compartment models have been one of the most used compartment models to simulate dynamics of diseases within populations by partitioning individuals into different compartments depending on their disease status [6,7,8,9,10,11,12,13,14,15].

Agent-based models explicitly simulate the individual-to-individual interactions within a population and hence provide yet further detail, and are thus much better suited to capturing heterogeneous mixing patterns and individual behavior. Stochastic models incorporate randomness in the transmission process; this is appropriate for modeling the natural variability and, in particular, the unpredictability of real-world disease transmission. The network model approaches that explicitly represent populations as interconnected networks of individuals offer insight into how the structure of social interactions drives disease dynamics. Besides these mathematical methods, a significant contribution is seen in modern epidemiology in data-driven approaches. Using data around surveillance analysis and genomic epidemiology, model-based approaches use real-time data and large datasets to track disease spread, identify risk factors, and optimize intervention strategies.

Another key aspect involves intervention modeling, such as different vaccination scenarios, simulations for quarantine and isolation policies, and social distancing measures. These models enable the modeling of the potential impacts of various interventions on disease containment and control. In reality, such models are used in planning pandemic responses, and governments and health agencies use them to forecast outbreaks, assess the effectiveness of interventions, and inform resource allocation. But once again, such models are hostage to the quality of data, the accuracy of assumptions behind them, and their ability to capture uncertainty and behavior changes in the population. Nonetheless, model techniques will remain invaluable for scientific insight and management of infectious disease spread.

2. Preliminaries

This section presents basic definitions, theorems, and results connected with the stability, existence, and uniqueness of the model developed based on the Caputo fractional derivative. Such a mathematical toolkit will provide a complete setting where the dynamic behavior of the model will be scrutinized. Moreover, we start by recalling some relevant definitions of the Caputo fractional derivatives, fixed point theorems, and stability criteria. These are critical elements to help validate the robustness and reliability of our model for capturing the dynamics of HFMD.

Definition 1 

([16]). The integral formula for the Caputo-fractional derivative of order $(\alpha \in (m-1,m\left]\right)$ for a function $\left(h\right(x\left)\right)$ that is $\left(m^{th}\right)$ -differentiable is given by:

$${}_{0}D_t^{\alpha }h\left(x\right)={\displaystyle \frac{1}{\Gamma (m-\alpha )}}{\int }_0^x{(x-t)}^{m-\alpha -1}{h}^{\left(m\right)}\left(t\right),dt,m=\lceil \alpha \rceil$$

Theorem 1 

([17]). Any contractive operator $(T:Y\to Y)$ that maps a complete metric space onto itself has a unique fixed point $(T\left(z^{*}\right)=z^{*})$. Furthermore, $\left(T\right)$ satisfies the following condition:

$$\mathrm{dist}(T\left(z^{*}\right),T\left(w\right))<K,\mathrm{dist}(z^{*},w),0<K<1.$$

Definition 2 

([18]). For integrable function $f:\mathbb{R}⟶\mathbb{R}$ and $0<\alpha \le 1$, the fractional-integral for the function f of order α is given by

$${\mathcal{I}}^{\alpha }f\left(y\right)={\displaystyle \frac{1}{\Gamma \left(\alpha \right)}}{\int }_0^y{(y-t)}^{\alpha -1}f\left(t\right)dt$$

Lemma 1 

([19]). Representing the Laplace transform associated with the Caputo fractional derivative:

$$\mathcal{L}{\{}_0{D}_t^{\alpha }f\left(t\right)\}=s^{\alpha }\mathcal{L}\left\{f\left(t\right)\right\}-s^{\alpha -1}f\left(0\right),0<\alpha <1$$

Theorem 2 

([20]). The equilibrium solutions $x^{*}$ of the Caputo-fractional differential equations system

$${}_{0}D_t^{\alpha }x\left(t\right)=f(t,x),x\left(t_0\right)=x_0$$

The Jacobian matrix ${\displaystyle \frac{\partial f}{\partial x_j}}$ evaluated at equilibrium points ensures local asymptotic stability (LAS) if its eigenvalues ${\lambda }_j$ satisfy

$$\left|arg\left({\lambda }_j\right)\right|=\pi >\alpha \pi /2,0<\alpha <1,\forall j=1,\cdots ,n.$$

Theorem 3 

([21]). Let $x\left(t\right)\in {\mathbb{R}}^{+}$ be a continuous and derivable function. If $x^{*}\in {\mathbb{R}}^{+}$ and $0<\alpha <1$, then for any time instant $t\ge t_0$, we have

$${}_{0}D_t^{\alpha }\left[x\left(t\right)-x^{*}-x^{*}ln{\displaystyle \frac{x\left(t\right)}{x^{*}}}\right]\le \left(1-{\displaystyle \frac{x\left(t\right)}{x^{*}}}\right){}_{0}D_t^{\alpha }x\left(t\right)$$

Theorem 4 

([22]). Let $\left(z^{*}\right)$ be the equilibrium point of ${(}_0{D}_t^{\alpha }z\left(t\right)=g(t,z))$, where$(\Omega \subset {\mathbb{R}}^m)$ is a domain containing $\left(z^{*}\right)$. By defining $(U(t,z):{\mathbb{R}}^{+}\cup 0\times \Omega \to \mathbb{R})$, and assuming U is a continuously differentiable (Lyapunov candidate) map such that:

• $W_1\left(z\right)\le U(t,z)\le W_2\left(z\right),$
• ${}_{0}D_t^{\alpha }U(t,z)\ge -W_3\left(z\right)$, for all nonnegative $\left(t\right)$ and for all $z\in \Omega ,\alpha \in (0,1)$,

for given continuous positive definite functions$(W_1\left(z\right),W_2\left(z\right)),and\left(W_3\left(z\right)\right)$ on $\left(\Omega \right)$. Then $\left(z^{*}\right)$ is globally asymptotically stable.

Our article is organized as follows: We begin with the proposed model, introducing the fractional model and the proofs of existence and uniqueness. This is followed by a section dedicated to Equilibrium Stability Analysis. Next, we include a section on Endemic Equilibrium. Finally, we conclude with a section on Discussion and Numerical Results.

3. Proposed Model for HFMD Dynamical Spreading

In this section, we introduce a system of differential equations incorporating fractional derivatives to model the dynamical spreading of Hand, Foot, and Mouth Disease (HFMD). Using fractional derivatives allows for a more accurate representation of memory and hereditary properties in the disease transmission process. This approach captures the complex dynamics of HFMD more effectively than traditional integer-order models. The equations presented here reflect these advanced modeling techniques, providing a robust framework for analyzing the spread of HFMD under various realistic scenarios, including vaccination efficiency. Figure 1 illustrates the interconnections between compartments, with all parameters clearly defined, while the system of equations 1 mathematically describes the dynamic behavior depicted in Figure 1.

$$\begin{array}{cc}\hfill & {}_{0}D_t^{\alpha }S={\Lambda }^{\alpha }-{\beta }^{\alpha }SI-\left({\tau }^{\alpha }+{\mu }^{\alpha }\right)S\hfill \\ & {}_{0}D_t^{\alpha }{S}_v=-{\beta }^{\alpha }\left(1-m^{\alpha }\right)S_{v}I+{\tau }^{\alpha }S-{\mu }^{\alpha }{S}_v\hfill \\ & {}_{0}D_t^{\alpha }I={\beta }^{\alpha }SI+{\beta }^{\alpha }\left(1-m^{\alpha }\right)S_{v}I-\left({\delta }^{\alpha }+{\omega }^{\alpha }+{\gamma }^{\alpha }+{\mu }^{\alpha }\right)I\hfill \\ \hfill & {}_{0}D_t^{\alpha }H={\omega }^{\alpha }I-\left({\pi }^{\alpha }+{\sigma }^{\alpha }+{\mu }^{\alpha }\right)H\hfill \\ \hfill & {}_{0}D_t^{\alpha }R={\gamma }^{\alpha }I+{\sigma }^{\alpha }H-{\mu }^{\alpha }R\hfill \\ \hfill & {}_{0}D_t^{\alpha }D={\delta }^{\alpha }I+{\pi }^{\alpha }H\hfill \end{array}$$

(1)

With the initial conditions:

$$S\left(0\right)\ge 0,S_v\left(0\right)\ge 0,I\left(0\right)\ge 0,H\left(0\right)\ge 0,R\left(0\right)\ge 0,D\left(0\right)\ge 0$$

Table 1 summarizes the detailed descriptions of the variables and parameters used in the proposed model. The model tracks different groups of the population based on their health state. S is the class of the susceptible population, and $S_v$ accounts for the vaccinated susceptible individuals. I denotes the proportion of those infected, H - those hospitalized, and R - those recovered. The parameters describe basic model dynamics: ${\Lambda }^{\alpha }$ describes new inflows into the population, ${\mu }^{\alpha }$ describes the mortality rate, and ${\beta }^{\alpha }$ describes the rate at which an infection spreads in the population. The vaccination efforts are described through ${\tau }^{\alpha }$, while through ${\gamma }^{\alpha }$ and ${\sigma }^{\alpha }$, the recovery rate models describe the infected and hospitalized populations. Hospitalization and death rates are described through ${\omega }^{\alpha }$, ${\delta }^{\alpha }$, and ${\pi }^{\alpha }$. Finally, $m^{\alpha }$ describes the efficiency of the vaccine to impede the infection. The above variables and parameters define the dynamics of infection, recovery, hospitalization, and population mortality.

3.1. Positivity and Boundedness

In a mathematical model for epidemiology, positivity and boundedness are the two major features that guarantee the model operates in a way that shows real population dynamics. Positivity guarantees that no model variable, representing the number of susceptible, infected, or recovered individuals at any given time, assumes a negative value, which does not make biological sense. Boundability means that the model variables do not grow boundless but rather stay within the confines of the total size of a finite population. These characteristics are crucial for ensuring the model’s validity and reliability in accurately portraying the dissemination and management of diseases.

Define the region $\Omega =\{(S\left(t\right),S_v\left(t\right),I\left(t\right),H\left(t\right),R\left(t\right))\in {\mathbb{R}}_{+}^5:S,S_v,I,H,R\le E_{\alpha ,1}\left(\varphi t^{\alpha }\right)\}$. Where $E_{\alpha ,1}\left(\varphi t^{\alpha }\right)$ is the Mittag-Leffler function. The following theorem shows that $\Omega$ is positively invariant.

Theorem 5. 

The non-negative region Ω that includes the solutions of the model’s (1) equations is positively invariant for $t\ge 0$.

Proof. 

Write the total population $T\left(t\right)$ to be the sum of the model’s components i.e. $T\left(t\right)=S\left(t\right)+S_v\left(t\right)+I\left(t\right)+H\left(T\right)+R\left(t\right)$. Then

$$\begin{array}{cc}\hfill {}_{0}D_t^{\alpha }T& ={}_{0}D_t^{\alpha }S+{}_{0}D_t^{\alpha }{S}_v+{}_{0}D_t^{\alpha }I+{}_{0}D_t^{\alpha }H+{}_{0}D_t^{\alpha }R\hfill \\ \hfill & ={\Lambda }^{\alpha }-{\mu }^{\alpha }S-{\mu }^{\alpha }{S}_v-{\mu }^{\alpha }I-{\pi }^{\alpha }H-{\mu }^{\alpha }H-{\mu }^{\alpha }R\hfill \\ \hfill & ={\Lambda }^{\alpha }-{\mu }^{\alpha }T-{\pi }^{\alpha }H\hfill \\ \hfill & \le \varphi T\left(t\right),\mathrm{where},\varphi \mathrm{is}\mathrm{the}\mathrm{population}\mathrm{coefficient}.\hfill \end{array}$$

Therefore, we conclude that

$${}_{0}D_t^{\alpha }T\le \varphi T\left(t\right)$$

(2)

To proceed with the proof we, Let $\mathcal{L}\left\{T\left(t\right)\right\}=\widehat{T}\left(s\right)$ be the Laplace transform of $T\left(t\right)$ with the initial condition $T\left(0\right)=T_0$. By the application of the Laplace transform on Equation (2), we have

$$\begin{array}{cc}\hfill \mathcal{L}{\{}_0{D}_t^{\alpha }T\}& \le \mathcal{L}\left\{\varphi T\right(t\left)\right\}\hfill \\ \hfill s^{\alpha }\widehat{T}\left(s\right)-s^{\alpha -1}{T}_0& \le \varphi \widehat{T}\left(s\right)\hfill \\ \hfill & \widehat{T}\left(s\right)\le {\displaystyle \frac{s^{\alpha -1}{T}_0}{s^{\alpha }-\varphi }}\hfill \\ \hfill & ={\displaystyle \frac{T_0}{s-\frac{\varphi }{S^{\alpha -1}}}}\hfill \\ \hfill \widehat{T}\left(s\right)& \le {\displaystyle \frac{T_0}{s\left(1-\frac{\varphi }{s^{\alpha }}\right)}}\hfill \end{array}$$

(3)

To find $T\left(t\right)$, we take Laplace inverse of both sides of Equation (3), to get:

$$\begin{array}{cc}\hfill {\mathcal{L}}^{-1}\left\{\widehat{T\left(s\right)}\right\}& \le {\mathcal{L}}^{-1}\left\{T_0\sum _{n=0}^{\infty }{\displaystyle \frac{{\Phi }^n}{s^{n\alpha +1}}}\right\}\hfill \\ \hfill & =T_0\sum _{n=0}^{\infty }{\displaystyle \frac{{\left(\varphi t^{\alpha }\right)}^n}{\Gamma (n\alpha +1)}}\hfill \end{array}$$

Thus, $T\left(t\right)\le T_0{E}_{\alpha ,1}\left(\varphi t^{\alpha }\right)$, which completes the proof.

□

3.2. Existence and Uniqueness

We can prove the existence and uniqueness of the solution for the fractional epidemic model of foot, hand, and mouth disease by applying mathematical methodologies, mainly fixed-point theorems. The accent is on the Banach fixed-point theorem. This model, which combines fractional derivatives, better describes the dynamics of the disease due to the presence of memory in the population. The existence of a solution guarantees a well-defined trajectory of the disease spread for any given set of initial conditions. On the other hand, uniqueness guarantees that the behavior is deterministic; the course of the disease progression can thus be predicted and is not dependent on arbitrary choices. These properties are usually demonstrated by converting the system to an appropriate integral form and showing that the associated operator satisfies the fixed-point theorem conditions to guarantee the solution’s existence and uniqueness. The conclusion here is that the model produces a reliable prediction about the dynamics of the epidemic over time.

Theorem 6. 

Equations of System (1) are satisfying the Lipschitz continuous for $K_j\ge 0,j=1,2,3,4,5$

Proof. 

Rewrite the equations of System (1) in a steady-state mode as follows:

$$\begin{array}{cc}& g_1\left(S\right)={\Lambda }^{\alpha }-({\tau }^{\alpha }+{\mu }^{\alpha })S-{\beta }^{\alpha }SI\hfill \\ & g_2\left(S_v\right)={\tau }^{\alpha }S-{\mu }^{\alpha }{S}_v-\beta (1-m^{\alpha })S_{v}I\hfill \\ & g_3\left(I\right)={\beta }^{\alpha }SI+{\beta }^{\alpha }(1-m^{\alpha })S_{v}I-({\delta }^{\alpha }+{\omega }^{\alpha }+{\mu }^{\alpha }+{\gamma }^{\alpha })I\hfill \\ & g_4\left(H\right)={\omega }^{\alpha }I-({\pi }^{\alpha }+{\sigma }^{\alpha }+{\mu }^{\alpha })\hfill \\ & g_5\left(R\right)={\gamma }^{\alpha }I+{\sigma }^{\alpha }H-{\mu }^{\alpha }R\hfill \end{array}$$

Upon applying Theorem 1, we obtain:

$$\begin{array}{cc}\hfill ∥g_1\left(S\right)-g_1\left(S^{\prime }\right)∥& =∥({\tau }^{\alpha }+{\mu }^{\alpha })\left(S-S^{\prime }\right)+{\beta }^{\alpha }I\left(S-S^{\prime }\right)∥\hfill \\ & \le |{\tau }^{\alpha }+{\mu }^{\alpha }|∥S-S^{\prime }∥+|{\beta }^{\alpha }|\parallel I\parallel ∥S-S^{\prime }∥\hfill \\ & \le \left(\right|{\tau }^{\alpha }+{\mu }^{\alpha }|+|{\beta }^{\alpha }|{\max }_{t\in [0,T]}\parallel I\parallel )∥S-S^{\prime }∥\hfill \\ & =K_1∥S-S^{\prime }∥\hfill \end{array}$$

Where $K_1=|{\tau }^{\alpha }+{\mu }^{\alpha }|+|{\beta }^{\alpha }|{\max }_{t\in [0,T]}\parallel I\parallel <\infty$.

Similarly, the following results can be obtained:

$$\begin{array}{cc}\hfill ∥g_2\left(S_v\right)-g_2\left(S_v^{\prime }\right)∥& \le K_2\parallel S_v-S_v^{\prime }\parallel \hfill \\ \hfill ∥g_3\left(I\right)-g_3\left(I^{\prime }\right)∥& \le K_3\parallel I-I^{\prime }\parallel \hfill \\ \hfill ∥g_4\left(H\right)-g_4\left(H^{\prime }\right)∥& \le K_4\parallel H-H^{\prime }\parallel \hfill \\ \hfill ∥g_5\left(R\right)-g_5\left(R^{\prime }\right)∥& \le K_5\parallel R-R^{\prime }\parallel \hfill \end{array}$$

Where,

$K_2=|{\mu }^{\alpha }|+|{\beta }^{\alpha }(1-m^{\alpha })|{max}_{t\in [0,T]}\parallel I\parallel <\infty$,

$K_3=|{\beta }^{\alpha }|{max}_{t\in [0,T]}\parallel S\parallel +|{\beta }^{\alpha }(1-m^{\alpha })|{max}_{t\in [0,T]}∥S_v∥+|{\delta }^{\alpha }+{\omega }^{\alpha }+{\gamma }^{\alpha }+{\mu }^{\alpha }|<\infty$,

$K_4=|{\pi }^{\alpha }+{\sigma }^{\alpha }+{\mu }^{\alpha }|<\infty ,\mathrm{and}{K}_5=\left|{\mu }^{\alpha }\right|<\infty$, which completes the proof. □

Lemma 2. 

A Volterra-integral equations form can represent the system of equations of the model (1).

Proof. 

Consider the system of equations of the Model (1)

$${}_{0}D_t^{\alpha }{x}_j\left(t\right)=g_j(t,x_j),j=1,\cdots ,5.$$

(4)

By the application of Definition (2) to both sides of the Equation (4), we get:

$${\mathcal{I}}^{\alpha }{(}_0{D}_t^{\alpha }{x}_j\left(t\right))={\mathcal{I}}^{\alpha }\left(g_j(t,x_j)\right),j=1,\cdots ,8.$$

(5)

The left-hand side reduced to $x_j\left(t\right)-x_j\left(0\right)$ and hence Equation (5) becomes

$$x_j\left(t\right)-x_j\left(0\right)={\displaystyle \frac{1}{\Gamma \left(\alpha \right)}}{\int }_0^t{(t-u)}^{\alpha -1}{g}_j(u,x_j)du,\mathrm{for}\mathrm{all}j.$$

(6)

Thus, the proof is completed, and the system of equations of the Model (1) is converted to the following equivalent system of Volterra-integral equations:

$$x_j\left(t\right)=x_j\left(0\right)+{\displaystyle \frac{1}{\Gamma \left(\alpha \right)}}{\int }_0^t{g}_j(u,x_j){(t-u)}^{\alpha -1}du,\mathrm{for}\mathrm{all}j.$$

(7)

□

Theorem 7 

([28]). For $0<\alpha <1$ and $G=[0,u]\subseteq \mathbb{R},J=[x_j\left(0\right)-k,x_j\left(0\right)+k]$. Let $g_j:G\times J⟶\mathbb{R}$ be Lipschitz condition and continuous bounded $i.e.$$\exists !M_j>0,\left|\right|g_j\left(x\right)\left|\right|<M_j$, where $M_j={max}_{t\in [0,u]}\left|\right|g_j\left(x\right)\left|\right|$. If ${\displaystyle \frac{K_j{\ell }_j}{M_j}}<1$, then ∃ only one solution for the initial value problem (1), namely, $x_j\in C[0,u^{*}]$, where $u^{*}=min\left\{u,{\left({\displaystyle \frac{{\ell }_j\Gamma (\alpha +1)}{M_j}}\right)}^{{\displaystyle \frac{1}{\alpha }}}\right\}$. $j=1,2,\cdots ,5$.

Proof. 

Let $X=\{x_j\left(t\right)\in C([0,u],{\mathbb{R}}^5):||x_j\left(t\right)-x_j\left(0\right)||\le {\ell }_j$, $j=1,2,\cdots ,5$. Since every sequence $x_j^n$ in X converges to $x_j\in C([0,u],{\mathbb{R}}^5)$ with respect to infinity norm, $\left|\right|·{\left|\right|}_{\infty }$, and x is continuous and $\left|\right|x_j\left(t\right)-x_j\left(0\right)\left|\right|\le {\ell }_j$, then $\forall n\left|\right|x_j^n\left(t\right)-x\left(0\right)\left|\right|\le {\ell }_j$. Thus, X is closed and hence a complete metric space.

Define $\mathcal{F}:X⟶X$ such that

${\displaystyle {\mathcal{F}}_j\left(x\left(t\right)\right)=x_j\left(0\right)+\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^t{k}_{\alpha }(u,t){\mathcal{F}}_j(u,x_j\left(t\right))du}$, where $k_{\alpha }(u,t)={(t-u)}^{\alpha -1}$ then

$$\begin{array}{cc}\hfill \left|\right|{\mathcal{F}}_j\left(x_j\left(t\right)\right)-x_j\left(0\right)\left|\right|& =∥{\displaystyle \frac{1}{\Gamma \left(\alpha \right)}}{\int }_0^t{k}_{\alpha }(u,t)g_j(u,x_j\left(t\right))du∥,j=1,2,\cdots ,5.\hfill \\ \hfill & \le {\displaystyle \frac{1}{\Gamma \left(\alpha \right)}}{\int }_0^t{k}_{\alpha }(u,t)\left|\right|g_j(u,x_j\left(t\right))\left|\right|du\hfill \\ \hfill & \le {\displaystyle \frac{M_j}{\Gamma \left(\alpha \right)}}{\int }_0^t{k}_{\alpha }(u,t)du\hfill \\ \hfill & ={\displaystyle \frac{M_j{t}^{\alpha }}{\Gamma (\alpha +1)}}\le {\displaystyle \frac{M_j{\left(u^{*}\right)}^{\alpha }}{\Gamma (\alpha +1)}}\le {\ell }_j\hfill \end{array}$$

Therefore, $\left|\right|{\mathcal{F}}_j\left(x_j\left(t\right)\right)-x_j\left(0\right)\left|\right|\le {\ell }_j$, and hence, $\mathcal{F}$ maps X onto itself.

To show that $\mathcal{F}$ is a contraction operator, we assume that $y,z\in X$, such that

$$\begin{array}{cc}\hfill \parallel {\mathcal{F}}_j\left(y_j\left(t\right)\right)-{\mathcal{F}}_j\left(z_j\left(t\right)\right)\parallel & =∥{\displaystyle \frac{1}{\Gamma \left(\alpha \right)}}{\int }_0^t{k}_{\alpha }(u,t)(g_j(u,y_j\left(t\right))-g_j(u,z_j\left(t\right))du∥,j=1,2,\cdots ,5.\hfill \\ \hfill & \le {\displaystyle \frac{1}{\Gamma \left(\alpha \right)}}{\int }_0^t{k}_{\alpha }(u,t)\left|\right|(g_j(u,y_j\left(t\right))-g_j(u,z_j\left(t\right))\left|\right|du\hfill \\ \hfill & \le {\displaystyle \frac{K_j}{\Gamma \left(\alpha \right)}}\left|\right|y_j\left(t\right)-z_j\left(t\right)\left|\right|{\int }_0^t{k}_{\alpha }(u,t)du\hfill \\ \hfill & ={\displaystyle \frac{K_j{t}^{\alpha }}{\Gamma (\alpha +1)}}\left|\right|y_j\left(t\right)-z_j\left(t\right)\left|\right|\hfill \\ \hfill & \le {\displaystyle \frac{K_j{\left(u^{*}\right)}^{\alpha }}{\Gamma (\alpha +1)}}\le {\displaystyle \frac{K_j{\ell }_j}{M_j}}\left|\right|y_j\left(t\right)-z_j\left(t\right)\left|\right|.\hfill \end{array}$$

Therefore, $\parallel {\mathcal{F}}_j\left(y_j\left(t\right)\right)-{\mathcal{F}}_j\left(z_j\left(t\right)\right)\parallel \le {\displaystyle \frac{K_j{\ell }_j}{M_j}}\left|\right|y_j\left(t\right)-z_j\left(t\right)\left|\right|$. Thus, by the assumption ${\displaystyle \frac{K_j{\ell }_j}{M_j}}<1$,$\mathcal{F}$ is a contraction mapping. Therefore, it possesses a unique fixed point. The system (1) has a unique solution. □

4. Equilibrium Stability Analysis

4.1. Basic Reproduction Number

The reproduction number, referred to as ${\mathcal{R}}_0$, is an essential threshold parameter that defines the stability of the foot, hand, and mouth disease epidemic model. This parameter quantifies the mean number of secondary infections that result from one infected individual within a completely susceptible population. If ${\mathcal{R}}_0<1$, then the disease-free equilibrium is stable, and the disease will die out eventually; otherwise, if ${\mathcal{R}}_0>1$, then the disease will spread in the population, and the endemic equilibrium becomes stable. The reproduction number can be obtained by the next-generation matrix method, which involves linearizing the differential equations around the disease-free equilibrium system and then analyzing the transmission and recovery rates. The principal eigenvalue of the resultant matrix yields ${\mathcal{R}}_0$, offering valuable information regarding the likelihood of disease outbreaks and guiding control measures [29]. Following [30], the infection states are divided into transmission and transition matrices as follows:

$$F=\left[\begin{array}{cc}{\beta }^{\alpha }S+{\beta }^{\alpha }(1-m^{\alpha })S_v& 0\\ 0& 0\end{array}\right]\mathrm{and}V=\left[\begin{array}{cc}{\delta }^{\alpha }+{\omega }^{\alpha }+{\gamma }^{\alpha }+{\mu }^{\alpha }& 0\\ -{\omega }^{\alpha }& {\pi }^{\alpha }+{\sigma }^{\alpha }+{\mu }^{\alpha }\end{array}\right]$$

By the application of theorem (2) in [30], the basic reproduction number is defined to be the spectral radius of the matrix $F{V}^{-1}$, i.e., ${\mathcal{R}}_0=\rho \left(F{V}^{-1}\right)$. Therefore, the basic reproduction number can be given by the following expression:

$${\mathcal{R}}_0={\displaystyle \frac{{\Lambda }^{\alpha }{\beta }^{\alpha }{\mu }^{\alpha }+{\Lambda }^{\alpha }{\beta }^{\alpha }{\tau }^{\alpha }(1-m^{\alpha })}{{\mu }^{\alpha }({\mu }^{\alpha }+{\tau }^{\alpha })({\delta }^{\alpha }+{\omega }^{\alpha }+{\gamma }^{\alpha }+{\mu }^{\alpha })}}$$

(8)

4.2. Existence of Equilibria

To find the Hand, foot, and mouth-free equilibrium point (HFMEEP), we set the right-hand side of the model’s (1) equations to Zero. to get

$$\begin{array}{cc}\hfill {\Lambda }^{\alpha }-{\beta }^{\alpha }SI-\left({\tau }^{\alpha }+{\mu }^{\alpha }\right)S& =0\hfill \\ \hfill -{\beta }^{\alpha }\left(1-m^{\alpha }\right)S_{v}I+{\tau }^{\alpha }S-{\mu }^{\alpha }{S}_v& =0\hfill \\ \hfill {\beta }^{\alpha }SI+{\beta }^{\alpha }\left(1-m^{\alpha }\right)S_{v}I-\left({\delta }^{\alpha }+{\omega }^{\alpha }+{\gamma }^{\alpha }+{\mu }^{\alpha }\right)I& =0\hfill \\ \hfill {\omega }^{\alpha }I-\left({\pi }^{\alpha }+{\sigma }^{\alpha }+{\mu }^{\alpha }\right)H& =0\hfill \\ \hfill {\gamma }^{\alpha }I+{\sigma }^{\alpha }H-{\mu }^{\alpha }R& =0\hfill \end{array}$$

(9)

The solution of the system (9) is given by:

$${\mathcal{E}}_0=\left({\displaystyle \frac{{\Lambda }^{\alpha }}{{\mu }^{\alpha }+{\tau }^{\alpha }}},{\displaystyle \frac{{\tau }^{\alpha }{\Lambda }^{\alpha }}{{\mu }^{\alpha }\left({\mu }^{\alpha }+{\tau }^{\alpha }\right)}},0,0,0\right)$$

4.3. Local Stability Analysis

This section investigates local stability analysis of the fractional mathematical model of hand, foot, and mouth disease. The system’s behavior around its equilibrium points, particularly the disease-free equilibrium, is dealt with. Considering the nonlinear system of fractional differential equations, a linear system is considered around the approximate equilibrium solution. Then, the stability of such an equilibrium can be determined by the eigenvalue analysis of the Jacobian matrix of the system. In most cases, a fractional-order model’s Jacobian matrix characteristic equation involves the fractional powers of the eigenvalues and the integer-order classical models. When viewed in the context of fractional derivatives, if all the eigenvalues have negative real parts, then the equilibrium point is locally stable; small perturbations will decay in time and lead the system back to the equilibrium point. This, in turn, gives such an analysis a high degree of importance because understanding the dynamics of initial outbreaks and under what conditions this disease can be controlled or even eradicated is a must. Thus, the Jacobian matrix for the Model (1) is given by

$$\mathcal{J}=\left[\begin{array}{cccc}-{\beta }^{\alpha }I-\left({\tau }^{\alpha }+{\mu }^{\alpha }\right)& 0& 0& 0\\ {\tau }^{\alpha }& -{\beta }^{\alpha }\left(1-m^{\alpha }\right)I-{\mu }^{\alpha }& 0& 0\\ {\beta }^{\alpha }I& {\beta }^{\alpha }\left(1-m^{\alpha }\right)I& -\left({\delta }^{\alpha }+{\omega }^{\alpha }+{\gamma }^{\alpha }+\mu \alpha \right)& 0\\ 0& 0& {\omega }^{\alpha }& -\left({\pi }^{\alpha }+{\sigma }^{\alpha }+{\mu }^{\alpha }\right)\end{array}\right]$$

Theorem 8. 

The hand, foot, and mouth disease-free equilibrium point, ${\mathcal{E}}_0$, is locally asymptotically stable.

Proof. 

At the HFMD free equilibrium point, ${\mathcal{E}}_0=\left({\displaystyle \frac{{\Lambda }^{\alpha }}{{\mu }^{\alpha }+{\tau }^{\alpha }}},{\displaystyle \frac{{\tau }^{\alpha }{\Lambda }^{\alpha }}{{\mu }^{\alpha }\left({\mu }^{\alpha }+{\tau }^{\alpha }\right)}},0,0,0\right)$, the Jacobian matrix is given as follows:

$$\begin{array}{c}\hfill \mathcal{J}\left({\mathcal{E}}_0\right)=\left[\begin{array}{cccc}-\left({\tau }^{\alpha }+{\mu }^{\alpha }\right)& 0& 0& 0\\ {\tau }^{\alpha }& -{\mu }^{\alpha }& 0& 0\\ 0& 0& -\left({\delta }^{\alpha }+{\omega }^{\alpha }+{\gamma }^{\alpha }+{\beta }^{\alpha }\right)& 0\\ 0& 0& {\omega }^{\alpha }& -\left({\pi }^{\alpha }+{\sigma }^{\alpha }+{\mu }^{\alpha }\right)\end{array}\right]\end{array}$$

The eigenvalues of the Jacobian matrix at the HFMD free equilibrium point, $\mathcal{J}\left({\mathcal{E}}_0\right)$, are the solutions of the characteristic equation of the $\mathcal{J}\left({\mathcal{E}}_0\right)$, i.e., $|\lambda Identity-\mathcal{J}({\mathcal{E}}_0\left)\right|=0$. Namely, ${\lambda }_1=-\left({\delta }^{\alpha }+{\beta }^{\alpha }+{\gamma }^{\alpha }+{\omega }^{\alpha }\right),{\lambda }_2=-{\mu }^{\alpha },{\lambda }_3=-\left({\mu }^{\alpha }+{\pi }^{\alpha }+{\sigma }^{\alpha }\right)\mathrm{and}{\lambda }_4=-\left({\mu }^{\alpha }+{\tau }^{\alpha }\right)$. Since all the eigenvalues of $\mathcal{J}\left({\mathcal{E}}_0\right)$ are negative real numbers, then $\left|arg\left({\lambda }_j\right)\right|=\pi >{\displaystyle \frac{\alpha \pi }{2}},0<\alpha <1$ for $j=1,2,3,4$.

Thus, by theorem (2), the HFMD free equilibrium point, ${\mathcal{E}}_0$ is locally asymptotically stable. □

5. Endemic Equilibrium

Endemic equilibrium of hand, foot, and mouth disease can be defined as a steady-state condition whereby the disease persists within a population at a constant level over time. This, in other words, means that the number of infected individuals at this equilibrium remains the same, neither increasing nor decreasing significantly, since the rates of new infections and recoveries balance out one another. Mathematical models of the disease use parameters on transmission rates and recovery rates for population dynamics in analyzing the endemic equilibrium. The endemic equilibrium is helpful for a long-term understanding of the disease and for assessing control strategies. This equilibrium is hugely important because the stability analysis will show whether small perturbations will cause the disease prevalence to go back to its steady state or give rise to outbreaks, and it may thus be central to informing public health interventions.

5.1. Existence of the HFMD Endemic Equilibrium

Let ${\mathcal{E}}^{*}=\left(S^{*},S_v^{*},I^{*},H^{*},R^{*}\right)$ be an equilibrium point of the Model (1). Thus, the following theorem follows.

Theorem 9. 

The proposed Model (1) has a unique endemic equilibrium point ${\mathcal{E}}^{*}$, whenever ${\mathcal{R}}_0>1$.

Proof. 

By solving the equations of the steady-state model (9) at the endemic equilibrium point ${\mathcal{E}}^{*}$, we get

$$S^{*}={\displaystyle \frac{{\Lambda }^{\alpha }}{{\beta }^{\alpha }{I}^{*}+{\tau }^{\alpha }+{\mu }^{\alpha }}},S_v^{*}={\displaystyle \frac{{\tau }^{\alpha }{\Lambda }^{\alpha }}{{\beta }^{\alpha }{I}^{*}+{\tau }^{\alpha }+{\mu }^{\alpha }}}·{\displaystyle \frac{1}{{\beta }^{\alpha }\left(1-m^{\alpha }\right)I^{*}+{\mu }^{\alpha }}}$$

$$I^{*}={\displaystyle \frac{-b\pm \sqrt{b^2-4ac}}{2a}},\mathrm{where}$$

$$\begin{array}{cc}\hfill a& =\left({\delta }^{\alpha }+{\omega }^{\alpha }+{\gamma }^{\alpha }+{\mu }^{\alpha }\right)\left(1-m^{\alpha }\right){\left({\beta }^{\alpha }\right)}^2\hfill \\ \hfill b& =\left[\left({\delta }^{\alpha }+{\omega }^{\alpha }+{\gamma }^{\alpha }+{\mu }^{\alpha }\right)\left({\tau }^{\alpha }+{\mu }^{\alpha }\right){\beta }^{\alpha }\left(1-{\mu }^{\alpha }\right)+{\mu }^{\alpha }{\beta }^2\right)-{\left({\beta }^{\alpha }\right)}^2{\Lambda }^{\alpha }\left(1-m^{\alpha }\right)\hfill \\ \hfill c& =\left({\delta }^{\alpha }+{\omega }^{\alpha }+{\gamma }^{\alpha }+{\mu }^{\alpha }\right)\left({\tau }^{\alpha }{\mu }^{\alpha }+{\left({\mu }^{\alpha }\right)}^2\right)-{\tau }^{\alpha }{\beta }^{\alpha }{\Lambda }^{\alpha }(1-m^{\alpha })-{\mu }^{\alpha }{\beta }^{\alpha }{\Lambda }^{\alpha }\hfill \end{array}$$

$$H^{*}={\displaystyle \frac{{\omega }^{\alpha }{I}^{*}}{\left({\pi }^{\alpha }+{\sigma }^{\alpha }+{\mu }^{\alpha }\right)}},R^{*}={\displaystyle \frac{{\delta }^{\alpha }{I}^{*}}{{\mu }^{\alpha }}}+{\displaystyle \frac{{\sigma }^{\alpha }{\omega }^{\alpha }{I}^{*}}{\left({\pi }^{\alpha }+{\sigma }^{\alpha }+{\mu }^{\alpha }\right){\mu }^{\alpha }}}$$

It is readily seen that all the compartments are positive whenever $I^{*}>0$. Moreover, $I^{*}>0$ iff $-b+\sqrt{b^2-4ac}>0$. $\mathrm{i}.\mathrm{e}{b}^2-4ac>b^2$ which implies that $ac<0$. That is $c<0$ since a is always positive. Thus,

$$\left({\delta }^{\alpha }+{\omega }^{\alpha }+{\gamma }^{\alpha }+{\mu }^{\alpha }\right)\left({\tau }^{\alpha }{\mu }^{\alpha }+{\left({\mu }^{\alpha }\right)}^2\right)-{\tau }^{\alpha }{\beta }^{\alpha }{\Lambda }^{\alpha }(1-m^{\alpha })-{\mu }^{\alpha }{\beta }^{\alpha }{\Lambda }^{\alpha }<0$$

which implies that

$${\displaystyle \frac{{\tau }^{\alpha }{\beta }^{\alpha }{\Lambda }^{\alpha }(1-m^{\alpha })+{\mu }^{\alpha }{\beta }^{\alpha }{\Lambda }^{\alpha }}{\left({\delta }^{\alpha }+{\omega }^{\alpha }+{\gamma }^{\alpha }+{\mu }^{\alpha }\right)\left({\tau }^{\alpha }{\mu }^{\alpha }+{\left({\mu }^{\alpha }\right)}^2\right)}}>1$$

Therefore, ${\mathcal{R}}_0>1$. Which completes the proof. □

5.2. Global Stability Analysis

Global stability analysis of a fractional mathematical model of hand, foot, and mouth disease deals with the system’s long-term behavior over the entire state space, so it goes beyond the local neighborhood of equilibrium points. This is important because it explains whether the system will always end in a disease-free or endemic equilibrium, regardless of initial conditions. Specifically, the global stability analysis in fractional models is often carried out with a specially constructed Lyapunov function, where the analysis demonstrates that the system’s total energy or "distance" from equilibrium decays monotonically. This means that a model is globally stable with a Lyapunov function being positive definite and its fractional derivative along the trajectories of the system being negative definite. The existence of such a function would ensure global convergence to a stable equilibrium of the system; otherwise, it would spread uniformly in the population or die out. Thus, this analysis plays a vital role in the design of control strategies, giving information on how the system behaves regarding its robustness and interventions. Thus, the following theorem was obtained:

Theorem 10. 

The hand, foot, and mouth endemic equilibrium point, ${\mathcal{E}}^{*}$, is global asymptotically stable if ${\mathcal{R}}_0>1$.

Proof. 

Let ${\mathcal{R}}_0>1$, then by Theorem (9) the Model (1) has a unique disease-present equilibrium point, ${\mathcal{E}}^{*}$. To proceed with the proof, we define the following Lyapunov function:

$$\begin{array}{cc}\hfill {\mathcal{F}}_2\left(S,S_v,I,H\right)& =\left[S-S^{*}-S^{*}ln{\displaystyle \frac{S}{S^{*}}}\right]+\left[S_v-S_v^{*}-S_v^{*}ln{\displaystyle \frac{S_v}{S-v^{*}}}\right]\hfill \\ \hfill & +\left[I-I^{*}-I^{*}ln{\displaystyle \frac{I}{I^{*}}}\right]+\left[H-H^{*}-H^{*}ln{\displaystyle \frac{H}{H^{*}}}\right]\hfill \end{array}$$

(10)

By taking the Caputo derivative for both sides of Equation (10), we get:

$$\begin{array}{cc}\hfill {}_{0}D_t^{\alpha }{\mathcal{F}}_2& {=}_0{D}_t^{\alpha }\left[S-S^{*}-S^{*}ln{\displaystyle \frac{S}{S^{*}}}\right]{+}_0{D}_t^{\alpha }\left[S_v-S_v^{*}-S_v^{*}ln{\displaystyle \frac{S_v}{S_v^{*}}}\right]\hfill \\ \hfill & {+}_0{D}_t^{\alpha }\left[I-I^{*}-I^{*}ln{\displaystyle \frac{I}{I^{*}}}\right]{+}_0{D}_t^{\alpha }\left[H-H^{*}-H^{*}ln{\displaystyle \frac{H}{H^{*}}}\right]\hfill \end{array}$$

(11)

Applying theorem (3) on Equation (11), yields

$${}_{0}D_t^{\alpha }{\mathcal{F}}_2\le \left[1-{\displaystyle \frac{S^{*}}{S}}\right]{}_{0}D_t^{\alpha }S+\left[1-{\displaystyle \frac{S_v^{*}}{S_v}}\right]{}_{0}D_t^{\alpha }{S}_v+\left[1-{\displaystyle \frac{I^{*}}{I}}\right]{}_{0}D_t^{\alpha }I+\left[1-{\displaystyle \frac{H^{*}}{H}}\right]{}_{0}D_t^{\alpha }H$$

(12)

Substituting the Caputo derivatives from Model (1) into Equation (12), we get:

$$\left[1-{\displaystyle \frac{S^{*}}{S}}\right]\times \left[{\Lambda }^{\alpha }-{\beta }^{\alpha }SI-\left({\tau }^{\alpha }+{\mu }^{\alpha }\right)S\right],\mathrm{but}\Lambda ={\beta }^{\alpha }{I}^{*}{S}^{*}+({\mu }^{\alpha }+{\tau }^{\alpha })S^{*}$$

Therefore,

$$\left[1-{\displaystyle \frac{S^{*}}{S}}\right]{}_{0}D_t^{\alpha }S\le {\beta }^{\alpha }{I}^{*}{S}^{*}\left[1-{\displaystyle \frac{S^{*}}{S}}-{\displaystyle \frac{SI}{S^{*}{I}^{*}}}+{\displaystyle \frac{I}{I^{*}}}\right]$$

(13)

Similarly, we will have the following:

$$\begin{array}{cc}\hfill \left[1-{\displaystyle \frac{S_v^{*}}{S_v}}\right]{}_{0}D_t^{\alpha }{S}_v& \le {\beta }^{\alpha }(1-m^{\alpha })S_v^{*}{I}^{*}\left[{\displaystyle \frac{S}{S^{*}}}\left(1-{\displaystyle \frac{S_v^{*}}{S_v}}\right)+{\displaystyle \frac{I}{I^{*}}}\left(1-{\displaystyle \frac{S_v}{S_v^{*}}}\right)\right]\hfill \end{array}$$

$$\begin{array}{cc}\hfill & +{\mu }^{\alpha }{S}_v^{*}\left[\left(1-{\displaystyle \frac{S_v}{S_v^{*}}}\right)+{\displaystyle \frac{S}{S^{*}}}\left(1-{\displaystyle \frac{S_v^{*}}{S_v}}\right)\right]\hfill \\ \hfill \left[1-{\displaystyle \frac{I^{*}}{I}}\right]{}_{0}D_t^{\alpha }I& \le -{\beta }^{\alpha }(1-m^{\alpha })S_v^{*}{I}^{*}\left[{\displaystyle \frac{S_v}{S_v^{*}}}\left(1-{\displaystyle \frac{I}{I^{*}}}\right)\right]-{\beta }^{\alpha }{S}^{*}{I}^{*}\left[{\displaystyle \frac{S}{S^{*}}}\left(1-{\displaystyle \frac{I}{I^{*}}}\right)\right]\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill & +\left({\delta }^{\alpha }+{\omega }^{\alpha }+{\gamma }^{\alpha }+{\mu }^{\alpha }\right)I^{*}\left(1-{\displaystyle \frac{I}{I^{*}}}\right)\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill \left[1-{\displaystyle \frac{H^{*}}{H}}\right]{}_{0}D_t^{\alpha }H& \le \left({\pi }^{\alpha }+{\sigma }^{\alpha }+{\mu }^{\alpha }\right)H^{*}\left[\left(1-{\displaystyle \frac{H}{H^{*}}}\right)+{\displaystyle \frac{I}{I^{*}}}\left(1-{\displaystyle \frac{H^{*}}{H}}\right)\right]\hfill \end{array}$$

(16)

Upon adding Equations (13) to (16) and by the application of the Arithmetic-Geometric mean inequality, which is the geometric mean is less than or equal to the arithmetic mean, we conclude that

$${}_{0}D_t^{\alpha }{\mathcal{F}}_2\le 0.$$

(17)

Therefore, by theorem (10), the hand, foot, and mouth endemic equilibrium point, ${\mathcal{E}}^{*}$, is globally asymptotically stable, and the proof is completed. □

6. Discussion and Numerical results

In this section, we validate the proposed model and the theorems presented in the previous section. Our results demonstrate the alignment between the theoretical analysis and the observed dynamics of both pandemic and endemic scenarios, particularly focusing on the reproduction number threshold. Using MATLAB, we simulate the model to visualize the results and analyze the sensitivity of the reproduction number. This allows us to observe the effects of various parameters on the spread of Hand, Foot, and Mouth Disease (HFMD). The simulations provide insights into how changes in specific parameters can influence the disease dynamics, offering valuable information for effective disease management and control strategies.

Figure 2a,b showcases two scenarios of disease progression with varying levels of fractional derivatives—alpha at $0.3$ and $0.8$, illustrating how changes in this parameter can impact the long-term behavior of an endemic illness. Each sub-figure represents a fractional derivative value, and an inset focuses on the long-term disease life cycle. Figure 2c illustrates the dynamics of Foot, Hand, and Mouth Disease using an alpha fractional derivative of $0.5$. Notably, the basic reproduction number ($R_0$) is calculated to be $1.98$, which exceeds the threshold of 1, indicating pandemic potential. Additionally, the inset figure focuses on the long-term behavior of disease progression.

Figure 3 illustrates the dynamics of Foot, Hand, and Mouth Disease using varying alpha fractional derivatives. The left figure represents short-term trends, while the right figure focuses on long-term behavior. In the short term, infection cases peak around $96-98$, declining sharply with higher alpha values. In the long term, all infection cases are gradually decreased over time.

The three-dimensional Figure 4 depicts the dynamics of Foot, Hand, and Mouth Disease using varying alpha fractional derivatives (specifically, $\alpha =0.5$). The axes are labeled as x-axis: $S+S_v$ (representing susceptible vaccinated individuals). y-axis: $I\left(t\right)$ (indicating the number of infected individuals over time). z-axis: $H+R$ (likely referring to hospitalized plus recovered individuals). Notably, all trajectories converge to a common point during the pandemic. This convergence occurs regardless of initial conditions, emphasizing that any initial disease discovery eventually leads to a consistent outcome.

Figure 5 consists of two graphs related to Foot, Hand, and Mouth Disease dynamics. The right graph illustrates the impact of varying recovery rates ($\gamma$) on infection cases, while the left graph focuses on vaccination efficiency (m). Increasing recovery rates and vaccination efforts contribute to reducing infection cases and controlling disease spread.”

Figure 6 represents Partial Rank Correlation Coefficients (PRCC) for the basic reproduction number, denoted as $R_0$. The graph displays vertical bars above and below a horizontal axis (which represents zero). The bars above the axis indicate positive PRCC values, while those below represent negative PRCC values. Each bar is labeled with a different Greek letter or symbol (such as $\beta$, $\lambda$, $\mu$, $\tau$, m, $\Lambda$, $\omega$, and $\gamma$). This visualization shows how different variables correlate with $R_0$ when other variables are held constant. PRCC is relevant in fields like epidemiology, where $R_0$ represents the basic reproduction number of an infection. It helps us understand how changes in specific factors impact disease spread. Our simulations validate the theoretical model, showing the importance of the reproduction number in HFMD dynamics. Sensitivity analysis reveals how parameter changes affect disease spread, aiding effective intervention strategies.

7. Conclusions

This study investigated the impact of vaccination on the transmission dynamics of Hand, Foot, and Mouth Disease (HFMD) using a fractional-order derivative model. The mathematical analysis section established the proposed model’s existence, uniqueness, and stability. Numerical simulations, implemented using MATLAB, were conducted to validate the model’s behavior under various scenarios. The results demonstrated the effectiveness of vaccination in reducing disease prevalence. Sensitivity analysis using the Partial Rank Correlation Coefficient (PRCC) identified key parameters influencing disease transmission. Overall, this research provides valuable insights into the dynamics of HFMD and underscores the importance of vaccination strategies for disease control.