On mathematical modelling of measles disease via collocation approach

1.   Introduction

Mathematical epidemiology is the use of mathematical models and methodologies to gain insights into the transmission and behavior of contagious diseases among populations. This discipline holds a crucial place in epidemiological investigations and contributes to the formulation of public health strategies. These mathematical models simulate disease transmission, predict outbreaks, assess the impact of interventions (such as vaccination campaigns or social distancing measures), and evaluate the effectiveness of various control strategies. We refer to [1]–[4]. These models incorporate variables such as population size, disease transmission rates, incubation periods, and other epidemiological parameters, giving important information on how infectious illnesses behave and directing public health initiatives. Mathematical modeling has proven to be valuable in the investigation and control of diseases like COVID-19, TB virus, influenza and Ebola disease, waterborne disease, hepatitis B virus and rabies virus, etc. For the mentioned diseases which have been studied by using mathematical models, we refer to [5]–[7].

Fractional calculus (FC) distinguishes itself from classical differential and integral calculus in how it applies integral and differential operators. Fractional order differential equations (FODEs) have been crucial in efforts to improve our understanding of many diseases, resulting in the development and investigation of mathematical models for a variety of diseases [8]. The use of FC has increased its value as a useful research tool in a number of fields across various domains in the realms of basic sciences and engineering. One key feature of fractional derivatives is their ability to manage integrals and derivatives of any order, regardless of whether they are real or complex. This unique property gives them a nonlocal quality, meaning that the future condition depends not just on the present state but also on all past states. Within the framework of FC, three distinct types of fractional differential operators exist. The Riemann-Liouville (RL) and Caputo relies on the power-law kernel [9]. The Caputo-Fabrizio derivative is rooted in decay processes [10], [11]. Finally, the Mittag-Leffler law encompasses both the power law and exponential decay. These properties can be effectively described by using Atangana–Baleanu fractional-order derivatives. These remarkable characteristics have been utilized in various fields, encompassing mathematics, applied sciences, engineering, biology and physics [12]. Authors [13] derived and simulated the numerical solutions for several control techniques in different fractional orders by using the iterative fractional-order Adams-Bashforth methodology. Researchers have made significant attempts to investigate this specific topic, exploring it from many perspectives which include both theoretical and numerical investigations. They developed an extensive number of methodologies and processes, each specifically to provide theoretical, analytical, and numerical results such as unraveling pine wilt disease by using a spectral method [14], the transmission dynamics and sensitivity analysis of pine wilt disease through the use of a fractal fractional operator [15], the sensitivity analysis of COVID-19 [16], the development of a fractional order model with non-local kernels [17], the fractional mathematical modeling of malaria disease and typhoid fever disease [18], [19], analysis and optimal control of SEIRI epidmic model [20]. Also, the numerical solution of the nonlinear delay integrodifferential equation using a wavelet [21], the variable-order fractional differential equation using Haar wavelet [22], Ulam- type stability of the impulsive delay integrodifferential equation [23] the computational modeling of a measles epidemic in human population [24], the co-dynamics of measles and dysentery disease [25], and the modeling of childhood disease outbreak in a community with the inflow of susceptible and vaccinated new-born [26]. These approaches have proven to be quite useful in the study of a wide range of problems, including difficulties relating to existence, approximation, and stability theory, among others [27]. Solving fractional differential equations (FDEs) to obtain accurate or analytical solutions can be time-consuming. In response to this difficulty, numerous tools and methodologies have arisen in recent decades to handle this topic and achieve approximations or analytical results. Techniques such as the decomposition method [28], transform method, and perturbation method [29] are examples of analytical and semi-analytical procedures. Numerical approaches have been created to identify approximate solutions to various FDEs problems. Spectral methods [30] are examples of notable numerical techniques.

Wavelet analysis has become a popular topic of research in several scientific and technical domains. Wavelets are viewed as a fresh basis for functions by a number of scholars and as a tool for time-frequency analysis by others. Given that wavelets are a flexible tool with many mathematical components and several potential uses, it is obvious that all of these are accurate [31], [32]. By using wavelet techniques, we can break down a complicated function into several smaller ones and study each one separately at various scales. This feature, along with a fast wavelet approach, makes these methods very interesting for analysis and synthesis. Wavelet-based collocation techniques have become more popular in numerical analysis because of their fast convergence, low computational cost, and straightforward procedure [33]. Wavelet techniques are a relatively recent addition to the family of orthogonal functions, with notable and attractive properties such as orthogonality, compact support, unconstrained regularity and good localization. As a result, they are frequently used to derive the numerical solutions to several mathematical models that are being developed in the fields of biology, chemistry, and physical sciences. Various types of wavelets, such as Chebeshev [34], the Haar method [35], etc. Hermite technique were used by researchers [36] to compute numerical solutions to infectious disease model. Authors [37] applied Bernstein wavelets to study a biological model. These wavelet-based approaches not only possess a strong mathematical foundation that also exhibit the capability to tackle nonlinear problems effectively. Among the various types of wavelets, Haar wavelets hold a distinct place. They are characterized by a pairwise constant function that forms the Haar wavelet basis series, making them one of the simplest wavelet series in mathematics. The Haar wavelet's orthogonality, local support, and simplicity make it an effective choice for use in the derivation of the numerical solution of FODEs. Recently, the Haar wavelets approach has been used to solve the HIV infection fractional model [38], the SEIR epidemic model [39]. The ability of Haar wavelets to efficiently capture complex nonlinear dynamics while maintaining high accuracy makes them particularly well-suited for handling the intricacies of fractional-order models. Moreover, the robustness of Haar wavelet methods when dealing with nonlinearities, coupled with their computational speed, significantly contributes to efforts to model the dynamic nature of measles disease. Given the inherently nonlinear nature of epidemiological systems and the fractional-order dynamics involved in measles transmission, researchers may improve the efficacy of measles dynamics simulations and analyses by making use of Haar wavelet approaches, which not only minimize computer cost they also increase model correctness and dependability. Due to the advantages offered by Haar wavelets when applied to solve the nonlinear models of fractional order, we are motivated to use the Haar wavelet collocation method (HWCM) to address the nonlinear fractional model of measles disease. The aim was to analyze the dynamics of the system with the utmost precision and minimal error. This work's major goal was to present and explore an HWCM that uses the Haar wavelet basis to understand the numerical and geometrical behavior of a nonlinear mathematical model. The findings in this research are original and have not previously been reported in the literature. This technique eliminates complex numerical methods and yields valuable insights into the numerical behavior of the model. The proposed solution also uses numerical computing to tackle the problem. Researchers have also used other numerical and analytical techniques to study various biological models, we refer to [40]–[42]. Also researchers have deduced various results devoted to mathematical analysis including the existence theory, and Ulam-Hyers stability criteria which have now a days extended to mathematical problems of physical importance. For mentioned results, we refer to [43]–[46].

1.1.   Structure of the model

In several infectious diseases, there is an initial latent period following exposure corresponding to the period before individuals become infectious. This latent period is a crucial aspect of disease progression and cannot be overlooked when analyzing infectious stages. Consequently, it is sensible to incorporate an initial compartment into epidemiological models. We have created a clear and predictable mathematical model to explain the process of measles transmission. In order to construct the model, the entire population (N) is segmented into four distinct categories: susceptible $(\mathbb{S})$, exposed $(\mathbb{E})$, infected $(\mathbb{I})$, and recovered $(\mathbb{R})$. In Figure 1, we have detailed the changes that occurred between these groups. The susceptible class, denoted as $(\mathbb{S})$, experiences an increase due to births or immigration at a rate represented by $\mathcal{B}$. At a rate π, natural death also has an impact. At a rate η, infection results from contact with infected people. At a rate η, interaction with infected individuals results in the generation of a class $\mathbb{E}$, which is the group of exposed individuals. This class decreases as a results of testing and treatment for measles at a rate of ω, transitioning to the infected class at a rate of Ω. In addition to these factors, the population in this class is also influenced by the natural mortality rate denoted as π. The group of infected people, denoted as class $\mathbb{I}$, emerges as a result of the transition of exposed peoples at a rate of Ω. It is decreased at a rate of ρ through infection recovery and at a rate of π due to natural mortality. The model requires the assumption that both recovered exposed persons and recovered infected persons develop a lifelong immunity to the illness. As a result, a class $\mathbb{R}$ is created, which is made up of those who have total immunity to the illness. Natural mortality affects this class of recovered people at a rate of π [40], [41].

Consequently, the deterministic model's diagram is as shown in Figure 1.

The following forms are used to present the fractional order SEIR model [27], [40], [42].

$N(\tau) = \mathbb{S}(\tau) + \mathbb{E}(\tau) + \mathbb{I}(\tau) + \mathbb{R}(\tau), \forall \tau$. According to (1.1), $(\mathbb{S} + \mathbb{E} + \mathbb{I} + \mathbb{R})^\prime = 0$; thus, $N(\tau)$ is constant and equal to N. System (1.1) is in a feasible region because $\Delta = \{(\mathbb{S} + \mathbb{E} + \mathbb{I} + \mathbb{R}) : 0 \leq \mathbb{S}, \mathbb{E}, \mathbb{I}, \mathbb{R} \leq N\}$. All associated parameters and state variables in the model stay non-negative while it follows the human population, where $\tau\geq 0$.

The remaining sections of this article are structured as follows. Section 2 discusses the preliminaries and fractional model formulation. Section 3 elaborates on the theoretical properties associated with the fractional model and qualitative analysis. In Section 4, we establish essential conditions for the Ulam-Hyers stability (UHS) of the solution within the context of the model under consideration. Section 5 presents the numerical scheme, along with graphical results and a corresponding discussion. Lastly, in Section 6, the conclusion is drawn.

2.   Preliminaries

Fractional-order models have garnered a significant amount of attention across various scientific disciplines and have been the focus of extensive research. Our exploration begins with the concepts of fractional-order integration and differentiation, as described in [9], [10]. We provide a brief overview of key lemmas and definitions from FC that are essential for studying the proposed model.

Definition 2.1. [1] For any function $\Theta \in L^1([0, \infty), \mathbb{R})$ the RL integral with order $\chi \in (0, 1)$ is given by

Definition 2.2. [1] For any function Θ the fractional order $\chi$ Caputo derivatives is defined as follows:

Definition 2.3. The subsequent equation holds true:

where r=0,1,…,n−1, with $b_r \in \mathbb{R}$, and $n = \lfloor \Theta \rfloor + 1$.

Lemma 2.4. Given g as a compact and continuous mapping from the Banach space $\mathcal{B}$ → $\mathcal{D}$, where $\mathcal{B}$ is characterized by elements $\mathcal{X} \in \mathcal{B}$ such that $\mathcal{X} = \Upsilon g \mathcal{X}$ with $\Upsilon$ in the interval [0, 1], when $\mathcal{D}$ is bounded, it implies that for the function g there exists at least one fixed point.

2.1.   Haar wavelet

Haar wavelet function $\mathcal{H}(\tau)$, along with its corresponding Haar scaling function denoted by $\widetilde{\mathcal{H}}_0(\tau)$ defined as:

On [0, 1), multi-resolution analysis produces a range of Haar wavelets, each of which can be represented as $\widetilde{\mathcal{H}}_m(\tau)$ [34]. As a result, this leads to the subsequent relationship:

where s=2i+q, q=0,1,…,2i−1, and i=0,1,…. Furthermore, we perform a translation of the Haar functions on the interval v−1≤τ<v, as follows

From [34] we conclude that the sequence $\left\langle \widetilde{\mathcal{H}}_s(\tau) \right\rangle_{s=0}^\infty$ constitutes a comprehensive orthonormal system within the space L2[0,1). Meanwhile, the sequence

is orthonormal in L2[0,ρ). This signifies that any function g(τ) within the domain L2[0,ρ) could represent the Haar orthonormal basis functions in series, as follows

Moreover, when this series is truncated, we obtain an approximate equivalent, denoted as yq(τ), for g(τ):

Here, the coefficients denoted by Cv,s can be calculated through the use of the inner product as for s=1,2,…,(q−1),v=1,2,…,ρ,

and the superscript T indicates the transpose of a matrix.

2.2.   Formulation of fractional model

Temporal memory effects are a common feature of biological processes, and particularly epidemiological dynamics, and they provide important new insights into nonlocal dynamics. Fractional derivatives provide a more efficient way to solve these difficult problems since time-varying kernels are intrinsic to non-integer order derivatives. Fractional derivatives appear in diverse forms in the literature, with the Caputo fractional derivative emerging as the most frequently encountered form. The Caputo operator offers a distinct benefit, as it does not require fractional initial values, unlike classical derivatives. Given these advantageous characteristics, we have opted to employ the Caputo operator in our computational model (1.1).

We add a time-varying kernel in the way described below in order to achieve the power correlation:

In integral form, the system (1.1) may be expressed as follows:

By inserting the χ−1 order Caputo derivative in (2.2), the following is obtained

Next, The operators $^CD^{\chi-1}_\tau$ and $I^{-(\chi-1)}$ exhibit an interesting property they mutually cancel each other out. Further, to keep the dimension balance on both sides, we re-write the model as follows:

Here, the basic reproductive number for the model (2.3) is given as follows:

The sensitivity index can be computed as follows [19]:

where p represents a parameter of expression (2.4). Using 2.5, we can compute the sensitivity index as follows:

We present the sensitivity index graphically in figure 2 as follows:

3.   Qualitative analysis

We evaluate the suggested model's well-posedness in this section [43], [44]. To achieve this, we employ methods from fixed-point theory to examine the solutions for the proposed system. The expression on the left-hand side of equation (1.1) assumes the following structure:

This allows us to systematically assess the model's solution stability and characteristics through well-established mathematical methods. Let the Banach space $\xi = C([0, T] \times \mathcal{R}^4,\mathcal{R})$, and $0 \leq \tau \leq \mathcal{T} < \infty$; then,

From (3.1), the proposed system (1.1) can takes the following form:

where

The Caputo initial value problem (3.1) along with definition (2.3) give

We rely on the following assumptions to demonstrate the existence of the (1.1):

(H1): $\exists$ positive constants $Q_{\mathcal{X}}$ and $\Psi_{\mathcal{X}}$ as $\forall$ $M \in \xi$:

(H2): $\exists$ a positive constant $\Psi_{\mathcal{X}}$ as $\forall$ $W, W' \in \xi$:

Theorem 3.1. [32] Assume that the conditions in (H1) hold and $\mathcal{X} : [0, \mathcal{T}] \times \mathcal{R}^4 \rightarrow \mathcal{R}$ is a continous mapping then there is at least one solution for equation (3.3). As a result, at least one solution exists for (1.1) with $vQ_{\mathcal{X}} < 1$, where $v = \frac{\mathcal{T}^\chi}{\Gamma(\chi+1)}$.

Proof. Assuming that (H1) is satisfied, for $\tau \in [0, \mathcal{T}]$, we define:

as a closed subset of ξ with convex properties, where $\zeta \geq \frac{v_0+v \Psi_\mathcal{X}}{1-vQ_\mathcal{X}}$. Further, define

as

Assume that

This implies that $\|T W(\tau)\|_\chi \leq \zeta$; hence, $T(L) \subset L$.

We examine τ1<τ2 within the interval $[0, \mathcal{T}]$, and we can conclude that T exhibits continuity. Our estimation becomes:

and so

Now, the right-hand side of (3.4) approaches 0, as τ₂ approaches τ₁. Therefore, $\|T W(\tau_2) - T W(\tau_1)\|_\xi \rightarrow 0$; clearly, T is uniformly continuous and bounded.

As a result T is completely continuous according to the Arzela–Ascoli theorem. Consequently, utilizing Schauder's fixed-point theorem, (1.1) possesses at least one solution.   □

Theorem 3.2. Assuming that (H2) is valid and $\mathcal{T}^\chi \Psi_\mathcal{X} < \Gamma(\chi + 1)$, the measles model (1.1) possesses a singular, unique solution.

Proof. Let W and W′ be two solutions in ξ, and consider T:ξ→ξ as the operator we have

The operator T exhibits continuity, and consequently, the Banach principle ensures the uniqueness of the solution to (1.1).   □

4.   Stability criteria

To conduct a comprehensive stability analysis of the proposed model, we revisit several definitions [45]. Consider T:ξ→ξ as a self-map defined by:

We will say that equation (4.1) exhibits Ulam-Hyers stability if, for every $\mathcal{E} > 0$ and W∈ξ, the solution satisfies:

Furthermore, there exists at most one solution W of (4.1) with gq>0, satisfying

Definition 4.1. [?] If, for all $W \in C(\mathbb{R})$ satisfying that W(0)=0 and for every solution W of eq. (4.2), with W being a solution of eq. (4.1), then the following inequality is satisfied:

then (4.1) exhibits generalized UHS.

Remark 4.2. [?] Let $Z(\tau) \in C([0, T];\mathbb{R})$, where W ¯∈ξ satisfies (4.3) under the following conditions

(i) $|Z(\tau)| \leq \mathcal{E}$,

(ii) TW ¯(τ)=W ¯+Z(τ).

For our analysis, we examine the perturbed initial value problem denoted in (3.2), i.e.,

with the initial condition W(0)=W0.

Lemma 4.3. [32]. The below inequality satisfies (4.4):

where

Theorem 4.4. [46] From Lemma 4.3, the solution to equation (3.2) exhibits UHS when TχLpΓ(χ+1)<1, which implies that the solution for the system (1.1) possesses a generalized UHS.

Proof. Consider an arbitrary solution W∈ξ, and let W ¯∈ξ be another solution (at most) for (3.2) we have:

we conclude that

Demonstrating the UHS of (3.2) also yields the generalized derivation of the UHS.   □

Definition 4.5. Consider (4.1) to demonstrate the stability in the sense of Ulam-Hyers-Rassias for a function $W \in C([0, \mathcal{T}], \mathbb{R})$ for any given $\mathcal{E} > 0$, with $W \in \xi$ as a solution of

there exists a solution W of (4.1) with gq>0 that satisfies the below condition

Definition 4.6. For $W \in C([0, \mathcal{T}], \mathbb{R})$, assume the existence of a constant Cq,W. Given that $\mathcal{E} > 0$, consider W as a solution of (4.5) and W as another solution of (4.1); then,

So (4.1) is generalized Ulam–Hyers–Rassias stable.

Lemma 4.7. The following inequality holds true for (4.3):

such that

Lemma 4.8. [?] According to Lemma (4.7), the solution of (4.3) has UHS and generalized UHS whenever $\frac{\mathcal{T}^\chi L_p}{\Gamma(\chi+1)} < 1$.

Proof. Consider an arbitrary solution, W∈ξ, and another solution, W ¯∈ξ, for (4.3). Here, we have:

This gives

As a result, (3.2) exhibits UHS, making it a case of generalized UHS    □

5.   Numerical scheme

Consider that $\mathbb{S}(\tau)$, $\mathbb{E}(\tau)$, $\mathbb{I}(\tau)$, and $\mathbb{R}(\tau)$ belong to L2[0,1), and can be represented by using Haar wavelets series, as follows

where' denotes the derivative, Xj,Yj,Zj,Wj are Haar series coefficients and Hj(τ) is the discretize Haar function [34]. We integrate the equations governing the transitions of individuals between these compartments to create a model that represents the progression of the epidemic over time. This integration yields a set of equations as follows:

Here, the operational matrix of integration is denoted by Pj,1(τ) for the Haar wavelet as explained in [34]. Now applying the Caputo derivative, we have

Assuming $\chi$ to be within the range (0, 1), it follows that n = 1, and hence

Next, from Haar approximations, the above becomes

After some calculation we get

Here, to approximate the integral in the prior system we have applied Haar's integration formula, as follows [35]:

Therefore, we obtain

and

and

and

Now, let

and

and

and

Combining the nodal points results in this nonlinear system yields

and

and

and

Utilizing Broyden's method, we can solve this system. The Jacobian matrix is expressed as follows:

The Jacobian matrix is determined by evaluating the following partial derivatives:

where

and

and

and

The unknown coefficients Xj,Yj,Zj,Wj are obtained from the solution of this above system; also, by substituting these coefficients into (5.1), we obtain the required solutions at nodal positions for $\mathbb{S}(\tau), \ \mathbb{E}(\tau),\ \mathbb{I}(\tau)$ and $\mathbb{R}(\tau)$. The rate of convergence in the experiments, as defined by the formula rρ(N) [34], [35], were computed by using the subsequent equation:

Next, we will present the graphical outcomes.

5.1.   Graphical results

In this study, we employed the fractional SEIR model to visualize the dynamics of individuals in different states: susceptible, exposed, infected, and recovered. To ensure the reliability of our investigations, we conducted numerical simulations and graphical analyses for a range of $\chi$ values. To solve this model computationally, we utilized the HWCM. Specifically, we have considered the fractional SEIR model (2.3), where we have set the initial values as follows: $\mathbb{S}(0) = 600$, $\mathbb{E}(0) = 250$, $\mathbb{I}(0) = 100$, and $\mathbb{R}(0) = 50$. It is important to note that we have assumed a constant total population size (N). Additionally, we have defined the following parameter values for the numerical simulations and results: $\mathcal{B} = 0.32, \pi = 0.2, \Omega = 0.01, \rho = 0.2, \eta = 0.01,$ and $\omega = 0.25$ per day. The resulting figures, encompassing Figures 2 through 5, offer valuable insights into the behaviors of different groups, including susceptible, exposed, infected, and recovered people. Figure 6 displays the fractional-order derivatives of susceptible people for different values of the fractional parameter $\chi$ ranging from 0.75 to 1. It demonstrates a decreasing number of vulnerable people over time as a result of viral exposure, which is consistent with other epidemiological models. In Figure 7, we can observe a clear trend in the population of individuals who have been exposed, showing a consistent and rapid increase as the fractional derivative converges towards its classical solution. The rise can be attributed to an increased proportion of susceptible people being infected and moving into the exposed category. In Figure 8, we can observe an increase in the count of infected individuals as the fractional order approaches 1. This phenomenon is attributed to the heightened sensitivity of the fractional order to different values of the fractional parameter denoted as $\chi$. Figure 9 depicts that as the fractional-order derivative approaches the classical value, we observe a consistent increase in the count of individuals who have successfully recovered, owing mostly to the recovery of infected individuals and thereby playing a significant role in the containment of the transmission of the disease. The population growth in the recovered category will accelerate with an increase in the fractional order. In Figure 10, we visualize the dynamics of the entire SEIR model (2.3) at $\chi=1$. A more accurate representation of the data was made possible by the application of Haar wavelet collocation techniques, which has led to a more exact and trustworthy model. It is also to be emphasized that we have observed a high rate of infected population at the initial stage of infection, but after a certain time, the rate of increase of infected individuals density sloweed down. Further, we have observed that the recovery rate was slow at the initial stage but became high after a certain time. The comparative analysis between the susceptible, exposed, infected, and recovered individuals in the fractional SEIR model has been shown graphically for various values of fractional order. Graphical representations demonstrate how parameters and fractional orders affect the behaviors of susceptible, exposed, infected, and recovered individuals, providing valuable insights into population dynamics during disease outbreaks. This study emphasizes the necessity of employing modern techniques to gain a deeper understanding of the dynamics of the measles disease.

6.   Conclusion

We have successfully applied the HWCM to solve fractional measles models efficiently. We investigated the application of fractional derivatives within the Caputo framework to analyze the dynamics of the measles epidemic model. Our method accounts for a number of variables that influence the spread of viruses, and the incorporation of fractional derivatives represents a major advancement in terms of accuracy and effectiveness. To efficiently handle these fractional derivatives, we made use of the Broyden methodology and the HWCM. The Ulam-Hyers method was helpful in determining the stability of the system and offered a vital perspective on the dependability of the model. The major contributions of this study are the application of the HWCM, our investigation of fractional derivatives in the Caputo framework, advancements in accuracy and effectiveness, and our impact on understanding population dynamics. We have shown how the parameters and fractional orders affect the behaviours of susceptible, exposed, infected, and recovered individuals through graphical representations of the dynamics within these compartments across a range of fractional order values. This study emphasizes the need to utilize modern techniques to acquire a greater understanding of the complexities of the measles outbreak. Consequently, the suggested approach is very efficient and applicable to many mathematical models, such as models of cancer treatment, drug targeting systems, and biotherapy. This methodology can be easily implemented in computer programs, and with a small modification to the existing approach, it may be extended to higher degrees. The study's finding suggests that the methodology is not only applicable to measles models but also to various other mathematical models, such as influenza, tuberculosis, or HIV/AIDS, for the analysis and prediction of disease spread dynamics.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.