1. Introduction
Hepatitis A Virus (HAV) is one type of hepatitis viruses, it generally causes self-limiting infection and does not result in chronic infection. However, several serious complications can occur, especially in older persons or in combination with risk factors and a small portion of infected individuals may die. It spreads through the fecal-oral route, typically by consuming contaminated food or water or through close personal contact with an infected person. HAV spread more rapidly in the areas with poor sanitation and hygienic measures. Around 90% of children are infected before the age of ten. The risk of infection is also high among those persons who inject drugs (PWID) and those of men who have sex with men (MSM) [35], [36], [47], [57].
The HAV has a worldwide distribution and causes about 1.5 million clinical cases each year. In 2016, WHO reported about 7134 deaths globally [56], [57] . About 3864 cases were reported in European countries in 2021 [16] and almost 44926 cases were reported in United States from 2016-2023 by US center for disease control and prevention [7].
Mathematical modelling serves as a valuable tool and plays a crucial role in understanding and predicting the spread of infectious diseases. Most infectious diseases are modelled as systems of ordinary differential equations (ODEs). These models assume instantaneous transitions and lack memory effects. However, involvement of fractional-order derivatives in epidemiology addressed this limitation by incorporating memory and hereditary properties through non-integer derivatives, offering a more accurate representation of complex biological processes. The fractional derivatives extend the concept of ordinary differentiation to non-integer order and have been increasingly used in various fields due to providing flexible and accurate models for
*Corresponding Author
Received July \(9^{th}\), 2025, Revised September \(10^{th}\), 2025, Accepted for publication October \(25^{th}\), 2025. Copyright ©2025 Published by Indonesian Biomathematical Society, e-ISSN: 2549-2896, DOI:10.5614/cbms.2025.8.2.4
complex phenomena. The integer-order derivatives are limited to local characteristics, whereas fractionalorder derivatives have broad scope and are influenced by the history and memory of the phenomenon [41], [43].
Various types of fractional derivatives are employed for investigating such real-world phenomena, like Caputo, Caputo-Fabrizio, Atangana-Baleanu and other types of fractional derivatives [2], [21], [26]. Each type of the fractional derivatives has its specifications and applications. In this work Caputo fractional derivatives are employed due to the nature of the initial and boundary conditions. As the initial and boundary conditions for differential equations with the Caputo derivatives are analogous to the case of integer-order differential equations, so they can be interpreted in the same way [43].
The Caputo fractional derivatives have been extensively used in epidemiology, there exist a rich and growing body of literature highlighting its significance, such as [2], [28], [29], [38], [52]. Many other authors also have shown strong interest in the integer-order of the models to be reduced to fractional-orders, such as the work by Paul et al. [40], where the authors reduced classic SIR model into fractional-order SIR model in consideration of Caputo derivatives. Especially, researchers attempted to formulate epidemiological behaviour of HAV for forecasting its future behaviour, guiding public health policies, and optimizing strategies for disease prevention and control.
Recently, Mwaijande and Mpogolo [24] developed a mathematical model for dynamics of HAV in consideration with vaccination and sanitation as prevention measures. They used Routh's stability criteria for local stability of disease free equilibrium point and obtained a Lyapunov function for the global stability of endemic equilibrium. Additionally, they studied sensitivity of the model and revealed that the model exhibits a forward bifurcation. Another recent literature on modelling dynamics of HAV is authored by Ben Aribi et al. [6], where their core interest is global stability analysis of their developed model, using Lyapunov method and numerical analysis based on the available data of Tunisia. They also visualized numerical results of the model in consideration with and without vaccines.
In a mathematical model developed by Wameko et al. [55], co-infection of HAV and Typhoid fever is investigated. The authors firstly studied sub models and then full model, as well as usual qualitative analysis and numerical analysis of the model were carried out. They also considered optimal control in their study, they revealed that prevention strategy has significant impact in reducing transmission of the co-infection and eventually they concluded that it can be successfully reduced by applying control measures.
In addition to earlier studies, several recent works highlight the growing importance of fractional calculus in epidemic modeling and numerical simulations. For example, novel fractional-order epidemic frameworks and stability analyses have been presented in [33], [34], while computational approaches for fractional epidemic systems were discussed in [32]. Similarly, [31] demonstrated the predictive potential of fractional epidemic models applied to real data. Further methodological advances, such as new fractional differential operators and numerical schemes, have been proposed in [15], [42]. Together, these contributions reinforce the relevance of adopting fractional-order approaches in epidemiology and motivate our work on Hepatitis A transmission dynamics.
Although a few researchers have investigated mathematical models for HAV transmission, most existing works are based on integer-order derivatives. Since the disease dynamics are influenced by history and memory, it is important to employ fractional-order derivatives that can capture these effects more accurately. The novelty of this work lies in three aspects, (i) the formulation of fractional-order HAV model in the sense of Caputo derivatives that incorporates memory and hereditary properties, (ii) the incorporation of awareness and a Holling type-II vaccination function to reflect realistic control strategy [50] and (iii) the estimation of possible model parameters using HAV data from U.S. To the best of our knowledge, such a comprehensive fractional-order framework for HAV dynamics has not been studied previously.
Furthermore, stability analysis has been performed, values of possible parameters have been estimated with the help of mean absolute error (MAE) estimation, while sensitivity analysis and simulations were also conducted for validating and supporting the qualitative findings. Finally, it has been concluded that reliable and precise results can be obtained using Caputo fractional-order derivatives. Also, awareness against HAV has a significant impact on disease dynamics, it does not ensure recovery of the individuals but slows down spread of the disease. However intervention of precaution vaccines indicated inevitable impact to reducing new infections. Applying both awareness and vaccination reduce new infections remarkably.
1.1. Preliminaries
In this section, basic definitions, theorems and properties are presented that will be useful in this paper.
Definition 1.1. The Caputo fractional derivative of order \(\alpha\) of a continuous function f(t) is defined as
\[{}^{C}D_{t}^{\alpha}f(t) = \frac{1}{\Gamma(n-\alpha)} \int_{0}^{t} (t-s)^{n-\alpha-1} f^{(n)}(s) ds,\] where \(\Gamma(*)\) is Gamma function, \(\alpha \in (n-1,n)\) and \(n \in \mathbb{N}\).
Particularly when \(\alpha \in (0,1)\), we have
\[{}^{C}D_{t}^{\alpha}f(t) = \frac{1}{\Gamma(1-\alpha)} \int_{0}^{t} (t-s)^{-\alpha}f'(s)ds.\]
Theorem 1.1. Let f(t) be n times continuously differentiable and \({}^CD_t^{\alpha}f(t)\) be piecewise continuous on \([0,\infty)\), where \(\alpha>0, \alpha\in(n-1,n)\) and \(n\in\mathbb{N}\) then Laplace transform of the Caputo derivative is defined as
\[\mathscr{L}\{^{C}D_{t}^{\alpha}f(t)\} = s^{\alpha}F(s) - \sum_{k=0}^{n-1} s^{\alpha-k-1}f^{(k)}(0),\] where \(\mathcal{L}\{f(t)\}=F(s)\).
Theorem 1.2. For any \(B \in \mathbb{C}^{n*n}\) and c, d > 0, let
\[\mathscr{L}\lbrace t^{d-1}E_{c,d}(Bt^c)\rbrace = \frac{s^{c-d}}{s^c - B},\] for \(\Re e(s) > \|s\|^{\frac{1}{c}}\), where \(\Re e(s)\) is real part of s and \(E_{c,d}(*)\) is Mittag-Liffler function.
Proposition 1.1. Let \(\alpha, \beta > 0\) and \(z \in \mathbb{C}\) then the Mittag-Liffler function satisfies
\[E_{c,d}(z) = zE_{c,c+d}(z) + \frac{1}{\Gamma(d)}.\]
Similar results can also be observed in [3], [11], [12], [25], [38].
2. MODEL FORMULATION
Since, it is pretty significant to understand dynamics of any phenomenon and consider assumptions for formulating it mathematically, hence in this section dynamics of HAV are brought into consideration and few assumptions are made for deriving a realistic model. Assumptions:
- 1) Deaths can only happen among infected individuals due to fulminant hepatitis.
- 2) Awareness is only considered among the susceptible individuals.
- Vaccination is only considered among aware susceptible individuals and it is not recommended during infection or after infection.
- Passing over treatment against HAV, as availability of no any specific and effective treatment against it is assured.
- 5) As hepatitis A (HA) is self-limited illness, therefore recovery of infected individual only occurs through self-reactivity of immune receptors neither through vaccination nor through treatment.
- 6) Each parameter in Table 1 is assumed for a specific dynamic and all have non-negative values.
For better understanding the dynamics of any disease among population, it is usually required to separate the populations into classes with same characteristics and then formulate the interactions between them as mathematical equations. Here, the attentive disease among population is HA. At the initial stage the whole population is considered susceptible, then a portion of susceptible population has forewarned about the disease, equipped them with the knowledge about HA and all the factors that affecting the disease. Later on, both the aware and unaware susceptible individuals acquire HA infection and they only recover through self limiting
Table 1: Involved parameter of the model with description.
| Parameter | Description |
|---|---|
| b | Recruitment rate |
| d | Natural death rate |
| \(d_1\) | Disease caused death rate |
| \(\rho\) | Infection rate of susceptible individuals |
| \(\omega\) | Precaution vaccination rate of susceptible individuals |
| a | Awareness rate of susceptible individuals |
| \(\gamma\) | Self-limiting recovery rate of infected individuals |
immunity as the vaccination is not effective during infection, as well as there is no any specific and effective treatment available for HA infection [46], [57]. Therefore three classes of population are considered, that are compartments of susceptible individuals S, HA infected individuals I and recovered individuals I, where the dynamics of the disease and interactions between individuals are formulated mathematically as the system below.
\[\begin{split} \frac{dS}{dt} &= b - \rho(1-a)SI - \rho aSI - af(S) - dS, \\ \frac{dI}{dt} &= \rho(1-a)SI + \rho aSI - \gamma I - (d+d_1)I, \\ \frac{dR}{dt} &= af(S) + \gamma I - dR, \end{split} \tag{1}\] with initial condition as \(S(0) = S_0\), \(I(0) = I_0\) and \(R(0) = R_0\), also
\[f(S) = \frac{\omega S}{1 + rS}.\]
Here f(S) is Holling type-II vaccination function and r is saturation constant of vaccines availability and supply, also \(\omega\) denotes vaccination rate. Since vaccination serves as a control mechanism, incorporating it as Holling type-II functional response can be regarded as density-dependent non-linear control. Several other control strategies for disease mitigation are discussed in [5], [27], [30], [48]. In this study, density-dependent non-linear control is taken into account to assess the impact of vaccination on susceptible individuals under resource constraints. The density-dependent non-linear control strategy captures both saturation and limitation effects, reflecting how vaccination efficacy varies with population density. In this strategy, the effectiveness of vaccination depends on the density of susceptible individuals. At low supply and availability recovery is slow due to resources constraints. Conversely, when vaccination resources are plentiful, recovery accelerates, which demonstrates how intervention efficiency is influenced by population density. [13], [39], [45].
The awareness analysis is performed through the parameter a which is awareness parameter and as well as ensuring that the condition \(\rho(1-a)>\rho a\) holds initially. Model (1) is generalized to a system of fractional differential equations (FDEs) as follows:
\[{}^{C}D_{t}^{\alpha}S = b - \rho(1 - a)SI - \rho aSI - af(S) - dS,\] \[{}^{C}D_{t}^{\alpha}I = \rho(1 - a)SI + \rho aSI - \gamma I - (d + d_{1})I,\] \[{}^{C}D_{t}^{\alpha}R = af(S) + \gamma I - dR,\] (2)
where \({}^CD_t^{\alpha}\) denotes Caputo derivative of order \(\alpha\) w.r.t. time.
In Model (2) the parameters have been retained in their classical form to preserve their original biological and epidemiological interpretations. Applying Caputo order on the involved parameters can increase model complexity and reduce identifiability. For theoretical analysis maintaining classical parameters simplifies the mathematical structure and enhances clarity [28], [51].
3. QUALITATIVE ANALYSIS
3.1. Positivity and Boundedness
In this section, positivity and boundedness are investigated and proven as the theorems below.
Theorem 3.1. If S(0) ≥ 0, I(0) ≥ 0 and R(0) ≥ 0 then the solutions S(t), I(t) and R(t) of the system are non-negative for all t ≥ 0.
Proof: Investigating positivity of each state variable individually. Positivity of S(t) is followed from the first equation of the model, which simplifies to
\[^{C}D_{t}^{\alpha}S = b - \rho SI - af(S) - dS.\]
For the purpose of positivity, only negative terms of the equation are considered
\[^{C}D_{t}^{\alpha}S \ge -\rho SI - af(S) - dS.\] (3)
Assuming that, intervention of vaccination is not applied to susceptible individuals i.e. ω = 0, as well as it is assumed that for time t > 0, there is fixed positive number of infected individuals m, then Inequality (3) reduces to
\[^{C}D_{t}^{\alpha}S\geq -pS,\] where p = ρm + d and is constant. Now using Laplace transform and recalling Theorem 1.1, we get
\[\mathcal{L}\lbrace^C D_t^{\alpha} S(t)\rbrace \ge -p \,\mathcal{L}\lbrace S(t)\rbrace,\]
\[\Rightarrow s^{\alpha} \mathcal{S}(s) - \sum_{k=0}^{n-1} s^{\alpha-k-1} S^{(k)}(0) \ge -p \mathcal{S}(s),\] where L {S(t)} = S (s). Since 0 < α < 1 then the previous inequality is obtained as
\[s^{\alpha} \mathcal{S}(s) - s^{\alpha - 1} S(0) \ge -p \mathcal{S}(s),\]
\[\Rightarrow \mathcal{L}\{S(t)\} = \mathcal{S}(s) \ge \frac{S(0)s^{\alpha - 1}}{s^{\alpha} + p}.\]
Taking inverse Laplace transform, we have
\[\mathscr{L}^{-1}\{\mathscr{S}(s)\} = S(t) \ge S(0)\mathscr{L}^{-1}\left\{\frac{s^{\alpha-1}}{s^{\alpha}+p}\right\}.\]
Using Laplace transform of Mittag Liffler function, we get
\[S(t) \ge S(0)E_{\alpha,1}(-pt^{\alpha}).\]
Since the Mittag Liffler function Eα,1(−ptα) ≥ 0 for 0 < α < 1, so S(t) ≥ 0 for all t ≥ 0. Similarly positivity of second and third equations of the model can be easily proven that I(t) ≥ 0, R(t) ≥ 0 ∀t ≥ 0. This completes the proof [3], [38].
Theorem 3.2. The feasible region of Model (2), defined as
\[\Omega = \left\{ (S, I, R) \in \mathbb{R}^3_+, N(0) \le N(t) \le \frac{b}{d} \right\},\tag{4}\] is positively invariant, where N(t) = S(t) + I(t) + R(t) is total population size.
Proof: Obtaining fractional derivative of total population by adding all equations of Model (2), we have CDα t N(t) ≤ b − dN(t).
Solving the Caputo fractional differential equation of order 0 < α < 1, using Laplace transforms as follows
\[\Rightarrow s^{\alpha} \mathcal{N}(s) - \sum_{k=0}^{n-1} s^{\alpha-k-1} f^{(k)}(0) \le \frac{b}{s} - d\mathcal{N}(s),\] where \(\mathcal{L}\{N(t)\} = \mathcal{N}(s)\). Since \(0 < \alpha < 1\), then we have
\[s^{\alpha} \mathcal{N}(s) - s^{\alpha - 1} N(0) \leq \frac{b}{s} - d \mathcal{N}(s),\] \[\Rightarrow (s^{\alpha + 1} + sd) \mathcal{N}(s) \leq b + N(0) s^{\alpha},\] \[\Rightarrow \mathcal{N}(s) \leq \frac{b}{s^{\alpha + 1} + sd} + \frac{N(0) s^{\alpha}}{s^{\alpha + 1} + sd}.\] (5)
By taking inverse Laplace transform the previous inequality is obtained as
\[N(t) \le b\mathcal{L}^{-1} \left\{ \frac{1}{s^{\alpha+1} + sd} \right\} + N(0)\mathcal{L}^{-1} \left\{ \frac{s^{\alpha}}{s^{\alpha+1} + sd} \right\}. \tag{6}\]
Using the Laplace transform of Mittag Liffler function given as Theorem 1.2, then Inequality (6), is obtained as
\[N(t) \le bt^{\alpha} E_{\alpha,\alpha+1}(-dt^{\alpha}) + N(0)E_{\alpha,1}(-dt^{\alpha}). \tag{7}\]
Applying Proposition 1.1 on Inequality (7) then we have
\[N(t) \leq bt^{\alpha} \left( \frac{1}{dt^{\alpha}} - \frac{1}{dt^{\alpha}} E_{\alpha,1}(-dt^{\alpha}) \right) + N(0)E_{\alpha,1}(-dt^{\alpha}),\] \[= \frac{b}{d} \left( 1 - E_{\alpha,1}(-dt^{\alpha}) \right) + N(0)E_{\alpha,1}(-dt^{\alpha}). \tag{8}\]
The Inequality (8) can be equivalently expressed as
\[N(t) \le \frac{b}{d} + \left(N(0) - \frac{b}{d}\right) E_{\alpha,1}(-dt^{\alpha}).\]
The Mittag Liffler function satisfies the condition \(0 \le E_{\alpha,1}(-dt^{\alpha}) \le 1\), so by multiplying it with the negative quantity \((N(0) - \frac{b}{d})\), we obtain
\[\left(N(0) - \frac{b}{d}\right) \le \left(N(0) - \frac{b}{d}\right) E_{\alpha,1}(-dt^{\alpha}) \le 0.\]
Since, \(\frac{b}{d}\) is positive, this gives
\[N(0) \le \frac{b}{d} + \left(N(0) - \frac{b}{d}\right) E_{\alpha,1}(-dt^{\alpha}) \le \frac{b}{d},\] which equivalently, be expressed as
\[N(0) \le N(t) \le \frac{b}{d}.\tag{9}\]
This implies that the feasible region \(\Omega\) is an invariant set of the system and all solutions of the system lie within the feasible region [3], [37], [38]. \(\blacksquare\) Since the model exhibits positivity and boundedness, it is biologically meaningful.
3.2. Equilibria and Basic Reproduction Number
The model exhibits two disease free equilibrium (DFE) points. First, when the disease does not exist at all in the considered population and precaution vaccines are not applied to susceptible individuals, then \(E_1\) is a DFE point. Second is DFE \(E_2\), when the disease dies out but precaution vaccines are applied to the susceptible population. The system also has an endemic equilibrium point \(E^*\), they are given by
\[E_1 = (\frac{b}{d}, 0, 0), \tag{10}\]
\[E_2 = \left(\frac{b(\omega r + 1)}{(dr + a)\omega + d}, 0, \frac{a\omega b}{d((dr + a)\omega + d)}\right),\tag{11}\]
\[E^* = (S^*, I^*, R^*), \tag{12}\] where
\[\begin{split} S^* &= \frac{d+d_1+\gamma}{\rho}, \\ I^* &= \frac{(d+d_1+\gamma)\Big(\omega(rd+a)+d\Big)-\rho b(\omega r+1)}{\rho(d+d_1+\gamma)(\omega r+1)}, \\ R^* &= \frac{b\rho\gamma(\omega r+1)+a\omega(d+d_1+\gamma)(d+d_1)-\Big(d\gamma(d+d_1+\gamma)(\omega r+1)\Big)}{d\rho(\omega r+1)(d+d_1+\gamma)} \end{split}\]
Basic reproduction number is usually denoted by \(\mathcal{R}_0\) and is determined for \(E_1\) and \(E_2\) respectively, given as follows
\[\mathcal{R}_0 = \max\{\mathcal{R}_1, \mathcal{R}_2\},\tag{13}\]
where
\[\mathscr{R}_1 = \frac{\rho b}{d(d+d_1+\gamma)}, \qquad \mathscr{R}_2 = \frac{\rho b(\omega r+1)}{((dr+a)\omega+d)(\gamma+d+d_1)},\] are basic reproduction numbers for \(E_1\) and \(E_2\) respectively.
3.3. Local Stability Analysis
Local stability of equilibria have been investigated and proved using a matrix corresponding to the system, given by
\[J = \begin{bmatrix} -\left(\rho(1-a)I + \rho aI + \frac{a\omega}{\omega r + 1} + d\right) & -\left(\rho(1-a)S + \rho aS\right) & 0\\ \rho(1-a)I + \rho aI & \rho(1-a)S + \rho aS - (d+d_1+\gamma) & 0\\ \frac{a\omega}{\omega r + 1} & \gamma & -d \end{bmatrix}. \tag{14}\]
Theorem 3.3. The DFE point \(E_1\) is locally asymptotically stable if \(\mathcal{R}_1 < 1\), otherwise it is unstable.
Proof: As mentioned, the DFE point \(E_1\) satisfies whenever vaccines are not applied to susceptible population i.e. \(\omega = 0\), then in this case the Jacobian matrix of the system at \(E_1\) is obtained as
\[J_{E_1} = \begin{bmatrix} -d & -\frac{b\rho}{d} & 0\\ 0 & \frac{b\rho - d(d+d_1+\gamma)}{d} & 0\\ 0 & \gamma & -d \end{bmatrix}.\] (15)
Third column of the matrix (15) is zero except diagonal entry, hence the non-zero diagonal entry in the column is an eigenvalue of the matrix i.e. \(\lambda_1 = -d\).
In complex plane, the eigenvalue \(\lambda_1\) can be written as \(\lambda_1 = -d + \iota \cdot 0\), then
\[|arg(\lambda_1)| = \tan^{-1}\left(\frac{0}{-d}\right) = \pi.\]
This satisfies that \(|arg(\lambda_1)| = \pi > \frac{\alpha\pi}{2}\) for d > 0. Other eigenvalues can be obtained from the reduced matrix given below
\[J_{E_{11}} = \begin{bmatrix} -d & \frac{b\rho}{d} \\ 0 & \frac{b\rho - d(d+d_1 + \gamma)}{d} \end{bmatrix}. \tag{16}\]
Since the matrix (16) is upper triangular, eigenvalues of the matrix are diagonal entries, given as
\[\lambda_2 = -d,\] \[\lambda_3 = \frac{b\rho - d(d + d_1 + \gamma)}{d}.\]
The eigenvalue \(\lambda_2 = \lambda_1\), which clearly satisfies that \(|arg(\lambda_2)| = \pi > \frac{\alpha\pi}{2}\) for d > 0.
Elaborating the eigenvalue \(\lambda_3\), it can be written as
\[\lambda_{3} = \frac{\mathcal{R}_{1}d(d+d_{1}+\gamma) - d(d+d_{1}+\gamma)}{d},\] \[\lambda_{3} = \frac{d(d+d_{1}+\gamma)(\mathcal{R}_{1}-1)}{d}.\] (17)
From (17), it is deduced that, the eigenvalue \(\lambda_3 < 0\) if \(\mathscr{R}_1 < 1\). Therefore \(|arg(\lambda_3)| = \pi > \frac{\alpha\pi}{2}\) holds. Since the condition \(|arg(\lambda_i)| > \frac{\alpha\pi}{2}\) is valid for \(\lambda_i\), i=1,2,3. Hence, based on the Routh Hurwitz stability criterion [1], [10], [14], the DFE point \(E_1\) is locally asymptotically stable and unstable otherwise.
Theorem 3.4. The DFE point \(E_2\) is locally asymptotically stable if \(\mathcal{R}_2 < 1\) and unstable otherwise.
Proof: Here, Jacobian matrix of the system computed at \(E_2\) is given as
\[J_{E_2} = \begin{bmatrix} -\left(\frac{a\omega}{\omega r + 1} + d\right) & -\frac{b\rho(\omega r + 1)}{(dr + a)\omega + d} & 0\\ 0 & \frac{\left(br\rho - rd^2 - \left(r(d_1 + \gamma) + a\right)d - a(d_1 + \gamma)\right)\omega - d^2 - (d_1 + \gamma)d + b\rho}{(dr + a)\omega + d} & 0\\ \frac{a\omega}{\omega r + 1} & \gamma & -d \end{bmatrix}.\](18)
Easily eigenvalues of matrix (18) can be obtained as \(\lambda_1=-d\) , \(\lambda_2=-\left(\frac{a\omega}{\omega r+1}+d\right)\) and
\[\lambda_3 = \frac{\left(br\rho - rd^2 - \left(r(d_1 + \gamma) + a\right)d - a(d_1 + \gamma)\right)\omega - d^2 - (d_1 + \gamma)d + b\rho}{(dr + a)\omega + d}.\]
Rearranging \(\lambda_3\) as
\[\lambda_{3} = \frac{\left(br\rho - rd^{2} - rd(d_{1} + \gamma) + ad - a(d_{1} + \gamma)\right)\omega - d(d + d_{1} + \gamma) + b\rho}{(dr + a)\omega + d},\] \[= \frac{br\rho\omega - \omega(dr + a)(d + d_{1} + \gamma) - d(d + d_{1} + \gamma) + b\rho}{(dr + a)\omega + d},\] \[= \frac{b\rho(\omega r + 1) - (d + d_{1} + \gamma)\left((dr + a)\omega + d\right)}{(dr + a)\omega + d}.\] (19)
Now replacing \(\mathcal{R}_2\) in Equation (19), then
\[\lambda_{3} = \frac{\mathscr{R}_{2}(d+d_{1}+\gamma)\Big((dr+a)\omega+d\Big) - (d+d_{1}+\gamma)\Big((dr+a)\omega+d\Big)}{(dr+a)\omega+d},\] \[= \frac{(d+d_{1}+\gamma)\Big((dr+a)\omega+d\Big)(\mathscr{R}_{2}-1)}{(dr+a)\omega+d}.\] (20)
Here, \(\mathcal{R}_2 < 1\), guarantees that \(\lambda_3 < 0\).
Since, all the eigenvalues \(\lambda_i < 0\), i=1,2,3, the condition \(|arg(\lambda_i)| = \pi > \frac{\alpha\pi}{2}\) holds. Therefore, the DFE point \(E_2\) is locally asymptotically stable but is unstable otherwise.
Theorem 3.5. The endemic equilibrium \(E^*\) is locally asymptotically stable if \(\mathcal{R}_2 > 1\) but is unstable otherwise.
_
Proof: Jacobian matrix of the system at \(E^*\) is obtained as
\[J_{E^*} = \begin{bmatrix} -\frac{b\rho}{d+d1+\gamma} & -(d+d_1+\gamma) & 0\\ \frac{\left(br\rho - rd^2 - d\left(r(d_1+\gamma) + a\right) - a(d_1+\gamma)\right)\omega - d(d+d_1+\gamma) + b\rho}{(d+d_1+\gamma)(\omega r + 1)} & 0 & 0\\ \frac{(d+d_1+\gamma)(\omega r + 1)}{\omega r + 1} & \gamma & -d \end{bmatrix}.\] (21)
Clearly an eigenvalue of matrix (21) is \(\lambda_1 = -d < 0\), which satisfies that \(|arg(\lambda_1)| = \pi > \frac{\alpha\pi}{2}\) for d > 0. Investigating stability of the system by observing either roots or coefficients of the characteristic equation of the reduced matrix, given as
\[J_E = \begin{bmatrix} -\frac{b\rho}{d+d1+\gamma} & -(d+d_1+\gamma) \\ \frac{(br\rho - rd^2 - d(r(d_1+\gamma) + a) - a(d_1+\gamma))\omega - d(d+d_1+\gamma) + b\rho}{(d+d_1+\gamma)(\omega r + 1)} & 0 \end{bmatrix}.\] (22)
Characteristic equation of the matrix \(J_E\) in (22) is written as
\[\lambda^2 + A\lambda + B = 0, (23)\]
where
\[A = \frac{b\rho}{d+d_1+\gamma},\] \[B = \frac{\left(br\rho - rd^2 - d\left(r(d_1+\gamma) + a\right) - a(d_1+\gamma)\right)\omega - d(d+d_1+\gamma) + b\rho}{\omega r + 1}.\]
The stability conditions for the quadratic polynomial in (23) are either Routh Hurwitz conditions or the conditions provided in [1]. According to Routh Hurwitz stability criterion, whenever the characteristic equation in (23) has positive coefficients, then it has negative roots, so the corresponding system is locally asymptotically stable [10], [14]. Here the coefficient of \(\lambda^2\) is positive and also the coefficient A > 0 for \(b, \rho, d, d_1, \gamma > 0\). However, for observing sign of coefficient B, rewriting it as
\[B = \frac{br\rho\omega + b\rho - \omega(dr+a)(d+d_1+\gamma) - d(d+d_1+\gamma)}{\omega r + 1},\]\[= \frac{b\rho(\omega r + 1) - (d+d_1+\gamma)\Big((dr+a)\omega + d\Big)}{\omega r + 1}.\]
Substituting \(\mathcal{R}_2\) from (13), then
\[B = \frac{\mathscr{R}_2(d+d_1+\gamma)\Big((dr+a)\omega+d\Big) - (d+d_1+\gamma)\Big((dr+a)\omega+d\Big)}{\omega r+1},\] \[= \frac{(d+d_1+\gamma)\Big((dr+a)\omega+d\Big)(\mathscr{R}_2-1)}{\omega r+1}.\] (24)
Here, \(\mathscr{R}_2 > 1\) implies that the coefficient B > 0 for at least d > 0, or \(d_1 > 0\) or \(\gamma > 0\). Since, the coefficients A > 0 and B > 0 with \(\mathscr{R}_2 > 1\), thus the polynomial has negative roots. It is summarized that, matrix (21) has negative eigenvalues, hence it clearly satisfies the condition \(|arg(\lambda_i)| > \frac{\alpha\pi}{2}\), for \(\lambda_i\), i = 1, 2, 3. Eventually, it is concluded that the endemic equilibrium \(E^*\) is locally asymptotically stable if \(\mathscr{R}_2 > 1\) but unstable otherwise.
3.4. Global stability analysis
Global stability analysis is a vital tool for determining the behaviour of complex systems and ensuring their resilience to deviations. In the context of epidemic models, global stability is performed by Lyapunov method and LaSalle's invariance principle, where the LaSalle's invariance principle extends the Lyapunov method and often requires identifying invariant sets, which may not be straightforward for complex systems. Recently, [4] proposed a simplified method of stability of non-linear systems. In contrast to these, the Lyapunov direct method is particularly well-suited for establishing global stability and the LaSalle's invariance principle is particularly effective when applied to systems with limited complexities, as they ensure asymptotic convergence toward one of the equilibria of the system. When a suitable Lyapunov function is constructed, it provides a clear and rigorous framework for demonstrating global stability, making it the most appropriate and effective choice for the objectives of this study. No any suitable method is available for constructing a Lyapunov function, however, some general forms of Lyapunov function are available [23], [54].
Theorem 3.6. The DFE point \(E_1\) is globally asymptotically stable if there exists a continuously differentiable function V, such that V is positive definite and its time Caputo derivative is negative definite at the equilibrium, further it is radially unbounded.
Proof: Let, the Lyapunov function \(V:\Omega\to\mathbb{R}\) defined as
\[V(S, I, R) = \left(S - \frac{b}{d}\right)^2 + I^2 + R^2.\] (25)
For showing that V is positive definite, one can easily verify \(V(E_1)=0\) by substituting \(E_1\in\Omega\) in the function V given in (25). Since, the function V involves square terms, hence V is always non-negative but it is strictly positive i.e. V(x)>0 for \(x\in\Omega\backslash E_1\), where \(x=(\bar S,\bar I,\bar R)\neq 0\). This validates positive definiteness of V.
Investigating negative definiteness of \({}^CD_t^{\alpha}V(S,I,R)\) by taking derivative of V, using the Definition 1.1, as
\[{}^{C}D_{t}^{\alpha}V(S,I,R) = {}^{C}D_{t}^{\alpha} \left[ \left( S - \frac{b}{d} \right)^{2} + I^{2} + R^{2} \right],\] \[= \frac{1}{\Gamma(1-\alpha)} \int_{0}^{t} (t-s)^{-\alpha} \frac{d}{ds} \left[ \left( S - \frac{b}{d} \right)^{2} + I^{2} + R^{2} \right] ds,\] \[= \frac{1}{\Gamma(1-\alpha)} \int_{0}^{t} (t-s)^{-\alpha} \frac{d}{ds} \left( S - \frac{b}{d} \right)^{2} ds + \frac{1}{\Gamma(1-\alpha)} \int_{0}^{t} (t-s)^{-\alpha} \frac{d}{ds} \left( I^{2} \right) ds\] \[+ \frac{1}{\Gamma(1-\alpha)} \int_{0}^{t} (t-s)^{-\alpha} \frac{d}{ds} \left( R^{2} \right) ds,\] \[= 2 \left( S - \frac{b}{d} \right) \frac{1}{\Gamma(1-\alpha)} \int_{0}^{t} (t-s)^{-\alpha} \frac{dS}{ds} ds + 2I \frac{1}{\Gamma(1-\alpha)} \int_{0}^{t} (t-s)^{-\alpha} \frac{dI}{ds} ds\] \[+ 2R \frac{1}{\Gamma(1-\alpha)} \int_{0}^{t} (t-s)^{-\alpha} \frac{dR}{ds} ds. \tag{26}\]
Applying Definition 1.1 on Equation (26), then
\[{}^{C}D_{t}^{\alpha}V(S,I,R) = 2\left(S - \frac{b}{d}\right){}^{C}D_{t}^{\alpha}S + 2I{}^{C}D_{t}^{\alpha}I + 2R{}^{C}D_{t}^{\alpha}R. \tag{27}\]
Simplifying Model (2) and substituting the model equations into (27), we get
\[{}^{C}D_{t}^{\alpha}V(S,I,R) = g + h + j, \tag{28}\]
where
\[g = 2\left(S - \frac{b}{d}\right)\left(b - dS - \rho SI - af(S)\right),\]
\[h = 2I\left(\rho SI - (\gamma + d + d_1)I\right),\]
\[j = 2R\left(af(S) + \gamma I - dR\right).\]
Analysing each term g, j and k individually, at \(x = (\bar{S}, \bar{I}, \bar{R})\). For analysing g, let \(K = \bar{S} - \frac{b}{d}\), then \(-dK = b - d\bar{S}\), which implies that
\[g = 2K \Big( -dK - \rho SI - af(\bar{S}) \Big).\]
Here, g<0 for every \(x\in\Omega\backslash E_1\). The term j<0 for \(x\in\Omega\backslash E_1\), if \(\rho\bar{S}< d+d_1+\gamma\). Also, the term k<0 for \(x\in\Omega\backslash E_1\), whenever \(-d\bar{R}\) dominates the term i.e. \(d\bar{R}>af(\bar{S})+\gamma\bar{I}\). Combining the results of g,j and k, then it is summarized that \({}^CD^\alpha_tV(x)<0\) for \(x\in\Omega\backslash E_1\) under the conditions \(\bar{S}<\frac{d+d_1+\gamma}{\rho}\) and \(\bar{R}>\frac{af(\bar{S})+\gamma\bar{I}}{d}\). Since, the considered function has positive terms with squares, then clearly V(x) increases as the point x goes farther from \(E_1\). Finally, for x with maximum length in the region \(\Omega\), the function V will reach to the upper bound of the region. Hence, it clearly hold \(V(x) \to \infty\) as \(||x|| \to \infty\) for \(x \in \Omega\), i.e., therefore V is radially unbounded in the feasible region.
Concluding that \(E_1\) is globally asymptotically stable due to the existence of V, which is positive definite for every \(x \in \Omega \backslash E_1\) and radially unbounded for every \(x \in \Omega\), also its time derivative in the sense of Caputo is negative definite for \(x \in \Omega \backslash E_1\), under the conditions \(\bar{S} < \frac{d+d_1+\gamma}{\rho}\) and \(\bar{R} > \frac{af(\bar{S})+\gamma\bar{I}}{d}\) [44].
Theorem 3.7. The DFE point \(E_2\) is globally asymptotically stable if there exists a continuously differentiable function V, such that V is positive definite and its time Caputo derivative \({}^CD_t^{\alpha}V\) is negative definite at the equilibrium, further it is radially unbounded.
Proof: Let, the Lyapunov function \(V:\Omega\to\mathbb{R}\) be defined as
\[V(S, I, R) = \left(S - \frac{b(\omega r + 1)}{\left((dr + a)\omega + d\right)}\right)^2 + I^2 + \left(R - \frac{a\omega b}{d\left((dr + a)\omega + d\right)}\right)^2.\] (29)
Time Caputo derivative of function V in Equation (29) is obtained as
\[^{C}D_{t}^{\alpha}V(S,I,R) = ^{C}D_{t}^{\alpha}\left[\left(S - \frac{b(\omega r + 1)}{\left((dr + a)\omega + d\right)}\right)^{2} + I^{2} + \left(R - \frac{a\omega b}{d\left((dr + a)\omega + d\right)}\right)^{2}\right],\] \[= 2\left(S - \frac{b(\omega r + 1)}{\left((dr + a)\omega + d\right)}\right)^{C}D_{t}^{\alpha}S + 2I^{C}D_{t}^{\alpha}I + 2\left(R - \frac{a\omega b}{d\left((dr + a)\omega + d\right)}\right)^{2}CD_{t}^{\alpha}R.\]
Substituting \({}^CD_t^{\alpha}S, {}^CD_t^{\alpha}I\) and \({}^CD_t^{\alpha}R\) from Model (2), then
\[{}^{C}D_{t}^{\alpha}V(S,I,R) = 2\left(S - \frac{b(\omega r + 1)}{\left((dr + a)\omega + d\right)}\right)\left(b - \rho SI - af(S) - dS\right),\] \[+ 2I^{2}\left(\rho S - (\gamma + d + d_{1})\right)\] \[+ 2\left(R - \frac{a\omega b}{d\left((dr + a)\omega + d\right)}\right)^{2}\left(af(S) + \gamma I - dR\right).\]
Now, it can be clearly seen, that V is positive definite as \(V(E_2) = 0\) and V(x) > 0 for every \(x \in \Omega \setminus \{E_2\}\). Furthermore, \({}^CD_t^{\alpha}V\) is negative definite as \({}^CD_t^{\alpha}V(E_2)=0\) and \({}^CD_t^{\alpha}V(x)<0\) for every \(x\in\Omega\backslash\{E_2\}\), while \(b< d(\bar{S}+\bar{R})\). Since \(V(x)\to\infty\) as \(||t||\to\infty\) for for all \(x\in\Omega\), hence V is radially unbounded.
Here, the positive definiteness of V and negative definiteness of its time Caputo derivative ensure asymptotic stability, however radially unboundedness show global behaviour of the function. Hence based on these it is concluded that the DFE point \(E_2\) is globally asymptotically stable.
Theorem 3.8. The endemic equilibrium point \(E^*\) is globally asymptotically stable if a continuously differentiable function V can be determined such that V is positive definite at \(E^*\) and is radially unbounded, additionally time Caputo derivative \({}^CD_t^{\alpha}V\) is negative definite at the equilibrium.
Proof: Considering the Lyapunov function \(V:\to \Omega \to \mathbb{R}\), defined as
\[V(S, I, R) = (S - S^*)^2 + (I - I^*)^2 + (R - R^*)^2.\] (30)
The function V is clearly continuously differentiable and is positive definite as as \(V(E^*) = 0\) and V(x) > 0 for every \(x = (S, I, R) \in \Omega \backslash E^*\) due to the square terms.
Applying Definition 1.1, the Caputo derivative of V w.r.t. time is obtained as
\[{}^{C}D_{t}^{\alpha}V(S,I,R) = 2(S-S^{*}){}^{C}D_{t}^{\alpha}S + 2(I-I^{*}){}^{C}D_{t}^{\alpha}I + 2(R-R^{*}){}^{C}D_{t}^{\alpha}R\] \[= 2(S-S^{*})[b-\rho SI - af(S) - dS]\] \[+ 2(I-I^{*})[\rho SI - (\gamma + d + d_{1})I]\] \[+ 2(R-R^{*})[af(S) + \gamma I - dR].\] (31)
Investigating each term of \({}^CD_t^{\alpha}V(S,I,R)\) individually. Let the first term be denoted by \(T_1\) as
\[T_1 = 2(S - S^*)[b - \rho SI - af(S) - dS]. \tag{32}\]
Substituting the equilibrium relation \(b = \rho S^*I^* + af(S^*) + dS^*\) in (32), we get
\[T_1 = 2(S - S^*)[\rho(S^*I^* - SI) + a(f(S^*) - f(S)) + d(S^* - S)].\] (33)
Case 1: If \(S^* < S\), then the factor \(2(S - S^*) > 0\), but the terms in the other factor will be
\[\rho(S^*I^* - SI) < 0, (34)\]
\[a(f(S^*) - f(S)) < 0, (35)\]
\[d(S^* - S) < 0. (36)\]
Multiplying the positive quantity \(2(S-S^*)>0\) with the negative quantities in (34-36), we obtain
\[2\rho(S - S^*)(S^*I^* - SI) < 0, (37)\]
\[2a(S - S^*) \Big( f(S^*) - f(S) \Big) < 0, \tag{38}\]
\[2d(S - S^*)(S^* - S) < 0. (39)\]
Combine (37-39), we get
\[2(S - S^*)[\rho(S^*I^* - SI) + a(f(S^*) - f(S)) + d(S^* - S)] < 0,\] \[\Rightarrow T_1 < 0.\] (40)
Case 2: If \(S^* > S\), then \(2(S - S^*) < 0\) but we get the terms as
\[\rho(S^*I^* - SI) > 0, (41)\]
\[a(f(S^*) - f(S)) > 0, (42)\]
\[d(S^* - S) > 0. (43)\]
Similarly, multiplying \(2(S - S^*) < 0\) with the positive quantities in (41-43) and combine them , we will get (40), hence \(T_1\) is also negative in this case.
Case 3: If \(S^* = S\), then \(2(S - S^*) = 0\), which implies that \(T_1 = 0\).
Now, investigating the second term. Let it be denoted as
\[T_2 = 2(I - I^*)[\rho SI - (\gamma + d + d_1)I]. \tag{44}\]
Case 1: If \(I^* < I\), then the factor \(2(I - I^*) > 0\), implies that \(T_2 < 0\) is valid, when the condition \(S < \frac{\gamma + d + d_1}{\rho}\) holds.
Case 2: If \(I^* > I\), then \(2(I - I^*) < 0\), which implies that \(T_2 < 0\), when \(S > \frac{\gamma + d + d_1}{\rho}\).
Case 3: If \(I^* = I\), then \(2(I - I^*) = 0\), hence \(T_2 = 0\).
Similarly, deliberating the third term. Let it be denoted by \(T_3\) as
\[T_3 = 2(R - R^*)[af(S) + \gamma I - dR]. \tag{45}\]
Case 1: If \(R^* < R\), then \(2(R - R^*) > 0\), so \(T_3 < 0\) while \(af(S) + \gamma I < dR\). Case 2: If \(R^* > R\), then \(2(R - R^*) < 0\), so \(T_3 < 0\) while \(af(S) + \gamma I > dR\). Case 3: If \(R^* = R\), then \(2(R - R^*) = 0\), so \(T_3 = 0\).
From above cases, it is summarized that \({}^CD_t^{\alpha}V(x) < 0\) for every \(x = (S, I, R) \in \Omega \backslash E^*\), holding the conditions mentioned and it vanishes at the equilibrium point \(E^*\), i.e. \({}^CD_t^{\alpha}V(E^*)=0\), so the singleton \(\{E^*\}\) is the only invariance set in the feasible region \(\Omega\), hence by LaSalle's invariance principle [23] the endemic equilibrium \(E^*\) is asymptotically stable. The function V is radially unbounded throughout \(\Omega\), thus it is concluded that \(E^*\) is globally asymptotically stable.
4. SIMULATION
4.1. Real-World Data
It is challenging to acquire comprehensive data of HAV. Despite this obstacle, it has been managed to obtain pertaining statistics of the infection and disease deaths from Center for Disease Control and Prevention-United States [8] available for period 2013-2022, shown in Table 2. Vaccination coverage statistics available for years 2004-2015, has been acquired from the recent work by Stroffolini and his co-author [49], shown in Table 3.
Table 2: Available statistics of infected individuals and disease deaths HAV pertaining to HAV [8].
| 10 | |||
|---|---|---|---|
| Year | Infected individuals | Disease deaths | 20 |
| 2013 | 1781 | 80 | 20 |
| 2014 | 1239 | 76 | 20 |
| 2015 | 1390 | 67 | 20 |
| 2016 | 2007 | 70 | 20 |
| 2017 | 3366 | 91 | 20 |
| 2018 | 12474 | 171 | 20 |
| 2019 | 18846 | 225 | 20 |
| 2020 | 9952 | 179 | 20 |
Table 3: Vaccination coverage statistics of HAV [49].
| Years | Coverage of vaccine (%) | Vaccinated individuals |
|---|---|---|
| 2004 | 7.7 | 7700 |
| 2005 | 16.86 | 16850 |
| 2006 | 16.86 | 16850 |
| 2007 | 20.034 | 20034 |
| 2008 | 20.034 | 20304 |
| 2009 | 17.05 | 17050 |
| 2010 | 18.267 | 18267 |
| 2011 | 17.05 | 17050 |
| 2012 | 17.05 | 17050 |
| 2013 | 17.05 | 1750 |
| 2014 | 26.4 | 26400 |
| 2015 | 26.4 | 26400 |
For enhancement of the data, Autoregressive Integrated Moving Average (ARIMA) and Exponential Smoothing (Holt's Model) have been used. ARIMA is a widely recognized statistical tool for time series analysis, traditionally used for forecasting future values based on historical data. However, its application to predict past values, often referred to as now-casting or hind-casting. ARIMA can estimate future and past values using recent data, ensuring continuity and completeness in the dataset, which is critical for trend analysis, [19], [20], [22]. Das and Muralidharan [9] have used a hybrid version of Holt's model to achieve accurate forecasts. As this technique is unrestricted by assumptions that usually bind ARIMA models, it has been used as a more robust alternative for hind-casting in specific cases, where assumptions pertaining to ARIMA models were observed to be violated.
In our study, we have used this technique to overcome hindrances faced due to lack of data pertaining to infected individuals, disease death and vaccinated individuals. The available data of infected individuals and deaths due to HAV have been improved by hind-casting the time points from 2004 to 2012, shown in Figure 1 and Figure 2. As well as, the data pertaining to vaccination has been enhanced by forecasting the time points for years 2016-2022, shown in the Figure 3.
Figure 1: Reported cases of HAV with hind-casted data points.
Figure 2: Reported cases of deaths due to HAV with hind-casted data points.
Figure 3: Vaccination coverage against HAV with forecasted data points.
4.2. Estimation
As mentioned earlier, that getting hand on comprehensive data for the purpose of fitting it into the model is a challenging task, therefore, the involved parameters b, d and a are assumed, while using the real world data given in Section 4.1, the parameter d1 is calculated (see Appendix A in [18]), the other parameters ρ, γ and ω are estimated by minimizing Mean Absolute Error (see Appendix B in [18]). Estimated values of the parameters for Model (1) are shown in Table 4.
| Parameter | Value | Source |
|---|---|---|
| b | 0.03408 | Assumed |
| d | 0.00616 | Assumed |
| d1 | 0.001579438 | Calculated |
| ρ | 0.007918408 | Estimated |
| ω | 0.857638 | Estimated |
| a | 0.428854899 | Assumed |
| γ | 0.783613446 | Estimated |
Table 4: Numerical values of the parameters.
4.3. Visualization
Since, it is very challenging to analytically solve a system of fractional-order differential equations due to the non-integer nature of the derivative, hence for validating the results, numerical simulations were conducted, where the findings are visualized as time series plots for all the state variables and different values of the fractional-order in this subsection. The results are visualized for the values of parameters given in Table 4, using the fde12 function [17] in MATLAB, which is based on the fractional Adams–Bashforth–Moulton predictor-corrector method for solving initial value problems involving Caputo fractional derivatives.
In Figure 4a and Figure 4b transmission of the population between compartments is shown in consideration with/without awareness and vaccination. From these figures, it can be observed that awareness and vaccination have significant impact over population. In presence of awareness, less number of susceptible individuals acquire the infection but the number of these individuals is comparatively high in absence of awareness. Additionally, in the first year the recovery of individuals through vaccination is faster in presence of vaccination but it is slower in absence of vaccination.
- (a) Transmission of individuals between compartments with awareness and vaccination.
- (b) Transmission of individuals between compartments without awareness and vaccination.
Figure 4: Transmission of individuals between compartments
In Figure 5, variations of γ can be observed. It shows that, higher natural immunity leads to high recovery and vice versa. This highlights the need of maintaining healthy life style and good hygiene. In Figure 6, effect of vaccination can be noticed, which shows that intensive HAV vaccination efforts targeting susceptible individuals leads to a significantly higher recovery rate during the early stages of infection. Since, no effective treatment exists against HAV, hence vaccination is a very crucial prevention strategy against HAV infection.
In Figure 7, transmission of individuals between compartments in consideration with Caputo order α is illustrated, it can be perceived that changes in the Caputo fractional-order cannot be neglected, where effect of varying fractional-order of the Caputo derivatives on susceptible, infected and recovered individuals can be seen in Figures 7a to 7c respectively.
Figure 5: Effect of natural immunity over HAV infection with different values of γ.
Figure 6: Effect of vaccination over HAV infection with different values of ω.
Figure 8: Infected individuals with different values of α, in presence of awareness.
Figures 8 and 9 show aware and unaware infected individuals in presence and absence of awareness respectively. The effect of Caputo order α can also be seen in the figures. From the figures it can be observed that, HAV infection is high among unaware individuals but is low among aware individuals. However, if awareness is not applied, then HAV infection will be at its peak, as it is illustrated Figure 9.
Figure 7: Dynamics of the infection in consideration with different values of α.
Figure 9: Infected individuals with different values of α, in absence of awareness.
In Figures 7 to 9, the corresponding changes in the disease dynamics are also illustrated for the case when α = 1, representing the classical integer-order model, allowing for direct comparison with the fractional-order cases. The comparison of the integer-order (α = 1) with fractional-orders (α < 1) reveals distinct dynamical features. In particular, the fractional-order model exhibits slower decay of the infected class and prolonged persistence of infection, indicating that memory effects delay disease eradication. By contrast, the integerorder model shows comparatively faster decay, which underscores the advantage of fractional derivatives in capturing hereditary and long-memory effects in HAV transmission.
The numerical simulations presented in this section provide important biological insights into HAV transmission. Figures 4- 6 demonstrate that increasing awareness and vaccination reduces the number of susceptible individuals entering the infected class, thereby lowering infection prevalence. In particular, vaccination shortens the infectious period and reduces outbreak peaks, consistent with public health evidence that HAV vaccines are highly effective in preventing community-level spread. Figures 7-9 further illustrate how varying the fractional-order α influences long-term epidemic outcomes. Smaller values of α correspond to stronger memory effects, which can sustain infection levels for a longer period even when transmission is reduced. This observation highlights the importance of incorporating population-level memory in HAV modeling, since it may explain why outbreaks sometimes persist despite interventions. Overall, the simulations confirm that the combined effects of awareness, vaccination, and fractional dynamics provide a more realistic framework for HAV control strategies.
So far, general insights on HAV dynamics are obtained from the visualized results of the model. Now performing sensitivity analysis to understand importance of parameters.
4.4. Sensitivity analysis
Sensitivity analysis helps to understand, how crucially a parameter is boosting transmission. Importance of sensitivity analysis is that it tells researchers, which parameter should be paid most numerical attention. By this it can be said, that a most sensitive parameter must be carefully estimated as it causes drastic changes on the dynamics. Here each involved parameter of R0 has been analysed for determining most sensitive parameter which increases newly infected individuals. This has been carried out with the help of normalized forward sensitivity of R0 defined as
\[\phi_{p_i}^{\mathcal{R}_0} = \frac{\partial \mathcal{R}_0}{\partial p_i} \times \frac{p_i}{\mathcal{R}_0},\tag{46}\] where pi is i th parameter of R0. The sensitivity indices are obtained using the normalized forward sensitivity formula, shown in the Table 5 and illustrated in Figures 10 and 11. The parameters b and ρ have positive indices and they are sensitive parameters of R1, whereas d, d1 and γ have no impact causing secondary infections. Similarly, b, ρ and r have positive indices and among them b and ρ are the most sensitive parameter of R2, whereas the other parameters d, d1, γ, a and ω are not sensitive parameters and have no impact to cause secondary infections. Any small variation of the sensitive parameters change numerical values of R1 and R2 drastically.
Table 5: Sensitivity indices of each parameter in R0 = max{R1, R2}.
| Parameter | Index(R1) | Index(R2) |
|---|---|---|
| b | +1 | +1 |
| d | -1.0078 | -0.0243 |
| d1 | -0.0020 | -0.0020 |
| ρ | +1 | +1 |
| γ | -0.9902 | -0.9902 |
| a | -0.9835 | |
| r | +0.00042 | |
| ω | -0.9831 |
For determining nature of the sensitive parameters, whether they increase or decrease new infections, the influence of the parameters over R1 and R2 have been visualized. Influence of ρ and b over R1 can be observed from Figure 12, where they can be described as promoters or risk factors as increasing or decreasing them can increase or decrease secondary infections drastically.
Figure 10: Sensitivity indices of R1. Figure 11: Sensitivity indices of R2.
Figure 12: R1 with respect to ρ and b. Figure 13: R2 with respect to ρ and b.
Figure 14: R2 with respect to ρ and r. Figure 15: R2 with respect to b and r.
Additionally, influence of ρ and b over R2 can be observed from Figure 13, where both the parameters are promoters, increasing/decreasing them will increase/decrease new infections. Figures 14 and 15 show the influence of r with respect to ρ and b over R2, it can be observed that effect of r is insignificant. In brief, it can be said that all the parameters with negative indices are always worthy in eliminating HAV, but the sensitive parameters b and ρ are risk factors, small deviation in b and ρ cause high number of new infections, remarkably in absence of vaccination.
CONCLUSION
This work presents a fractional-order SIR model for HAV dynamics, employing Caputo derivatives to effectively incorporate memory characteristics of the disease. The model uniquely integrates awareness campaigns and a non-linear vaccination strategy through a Holling type-II function. Through rigorous mathematical analysis, we established the model's well-posedness, local and global stability results. Further, possible parameters are estimated using real-world data from United States. The model is numerically solved using the predictor–corrector Adams–Bashforth–Moulton scheme, which is well-suited for fractional systems due to its accuracy and stability in handling memory-dependent dynamics and the simulation results highlighted the significant impact of fractional-order dynamics.
The findings show that fractional models offer superior flexibility and predictive capability compared to classical models. Both awareness and vaccination are shown to reduce infection levels with their combined application being most effective. Sensitivity analysis identified the infection rate (ρ) and recruitment rate (b) as the most influential parameters, suggesting that controlling these could significantly curb HAV transmission. Overall, this study underscores the utility of fractional calculus in epidemiology and provides a robust modelling framework for HAV. The model can guide public health interventions, particularly in areas with limited resources and high susceptibility due to inadequate sanitation or vaccine coverage.
ACKNOWLEDGEMENT
The authors extend their heartfelt thanks to all those who contributed to the improvment of this manuscript, especially grateful to Prof. Purnachandra Rao Koya for his insightful suggestions and generous support, which significantly enhanced the quality of our work. We also sincerely appreciate the anonymous reviewers for their constructive feedback and valuable comments that helped refine the manuscript. Furthermore, B.S. acknowledges the Indian Council for Cultural Relations (ICCR), Department of External Affairs, Government of India, New Delhi, for their generous sponsorship of his Ph.D. program.