Biomechanical responses of the cornea after small incision lenticule extraction (SMILE) refractive surgery based on a finite element model of the human eye

Yinyu Song; Lihua Fang; Qinyue Zhu; Ruirui Du; Binhui Guo; Jiahui Gong; Jixia Huang; Yinyu Song; Lihua Fang; Qinyue Zhu; Ruirui Du; Binhui Guo; Jiahui Gong; Jixia Huang

doi:10.3934/mbe.2021211

Mathematical Biosciences and Engineering

2021, Volume 18, Issue 4: 4212-4225. doi: 10.3934/mbe.2021211

Previous Article Next Article

Research article

Biomechanical responses of the cornea after small incision lenticule extraction (SMILE) refractive surgery based on a finite element model of the human eye

Key Laboratory of Nondestructive Test (Ministry of Education), Nanchang Hangkong University, Nanchang 330063, China

Academic Editor: Yang Kuang

Received: 20 March 2021 Accepted: 12 May 2021 Published: 17 May 2021

Purpose

To investigate the biomechanical responses of the human cornea after small incision lenticule extraction (SMILE) procedures, especially their effects of SMILE surgery on stress and strain.

Methods

Based on finite element analysis, a three-dimensional (3D) model of the human eye was established to simulate SMILE refractive surgery procedures. Stress and strain values were calculated by inputting the intraocular pressure (IOP).

Results

After SMILE refractive surgery procedures, the stress and strain of the anterior and posterior corneal surfaces were significantly increased. The equivalent stress and strain on the anterior and posterior corneal surfaces increased with increasing diopter and were concentrated in the central area, whereas the values of stress and strain at the incision site on the anterior surface of the cornea were approximately 0. Compared with the anterior corneal surface, the stress and strain of the posterior surface were larger. Increasing IOP caused an approximately linear change in stress and a nonlinear increase in corneal strain. In addition, we found that the incision sizes and direction had less of an influence on stress and strain. In summary, SMILE surgery increased the equivalent stress and strain on the human cornea.

Conclusions

The equivalent stress and strain of the anterior and posterior human corneal surfaces increased after SMILE refractive surgery; these increases were particularly noticeable on the posterior surface of the cornea.

Keywords:

Citation: Yinyu Song, Lihua Fang, Qinyue Zhu, Ruirui Du, Binhui Guo, Jiahui Gong, Jixia Huang. Biomechanical responses of the cornea after small incision lenticule extraction (SMILE) refractive surgery based on a finite element model of the human eye[J]. Mathematical Biosciences and Engineering, 2021, 18(4): 4212-4225. doi: 10.3934/mbe.2021211

Related Papers:

[1]	Mustafa Gürbüz, Yakup Taşdan, Erhan Set . Ostrowski type inequalities via the Katugampola fractional integrals. AIMS Mathematics, 2020, 5(1): 42-53. doi: 10.3934/math.2020004
[2]	Shuang-Shuang Zhou, Saima Rashid, Muhammad Aslam Noor, Khalida Inayat Noor, Farhat Safdar, Yu-Ming Chu . New Hermite-Hadamard type inequalities for exponentially convex functions and applications. AIMS Mathematics, 2020, 5(6): 6874-6901. doi: 10.3934/math.2020441
[3]	Hong Yang, Shahid Qaisar, Arslan Munir, Muhammad Naeem . New inequalities via Caputo-Fabrizio integral operator with applications. AIMS Mathematics, 2023, 8(8): 19391-19412. doi: 10.3934/math.2023989
[4]	Dazhao Chen . Weighted boundedness for Toeplitz type operator related to singular integral transform with variable Calderón-Zygmund kernel. AIMS Mathematics, 2021, 6(1): 688-697. doi: 10.3934/math.2021041
[5]	Haoliang Fu, Muhammad Shoaib Saleem, Waqas Nazeer, Mamoona Ghafoor, Peigen Li . On Hermite-Hadamard type inequalities for $n$ -polynomial convex stochastic processes. AIMS Mathematics, 2021, 6(6): 6322-6339. doi: 10.3934/math.2021371
[6]	Tekin Toplu, Mahir Kadakal, İmdat İşcan . On n-Polynomial convexity and some related inequalities. AIMS Mathematics, 2020, 5(2): 1304-1318. doi: 10.3934/math.2020089
[7]	Muhammad Imran Asjad, Waqas Ali Faridi, Mohammed M. Al-Shomrani, Abdullahi Yusuf . The generalization of Hermite-Hadamard type Inequality with exp-convexity involving non-singular fractional operator. AIMS Mathematics, 2022, 7(4): 7040-7055. doi: 10.3934/math.2022392
[8]	Gültekin Tınaztepe, Sevda Sezer, Zeynep Eken, Sinem Sezer Evcan . The Ostrowski inequality for $s$ -convex functions in the third sense. AIMS Mathematics, 2022, 7(4): 5605-5615. doi: 10.3934/math.2022310
[9]	Suriyakamol Thongjob, Kamsing Nonlaopon, Sortiris K. Ntouyas . Some (p, q)-Hardy type inequalities for (p, q)-integrable functions. AIMS Mathematics, 2021, 6(1): 77-89. doi: 10.3934/math.2021006
[10]	Chunyan Luo, Yuping Yu, Tingsong Du . Estimates of bounds on the weighted Simpson type inequality and their applications. AIMS Mathematics, 2020, 5(5): 4644-4661. doi: 10.3934/math.2020298

Abstract

Purpose

To investigate the biomechanical responses of the human cornea after small incision lenticule extraction (SMILE) procedures, especially their effects of SMILE surgery on stress and strain.

Methods

Results

Conclusions

1. Introduction

The term "infectious diseases" refers to illnesses caused by living organisms, usually parasites and microorganisms, that spread from one host to another either directly or indirectly. Infectious diseases have been a significant concern in our society. These illnesses can sometimes pose severe threats and lead to epidemics. Medical research has rapidly advanced in identifying and controlling three primary characteristics: Infectivity, epidemic potential, and uncertainty, which have helped reduce the burden of such diseases. Infectious diseases are modelled mathematically as a dynamic transmission cycle involving interactions between susceptible and diseased hosts. These interactions are usually expressed as a coupled ordinary differential equation. The SIR model takes into account several variables while examining infectious disease systems. Several models that examine significant aspects in the research of infectious diseases include SIS, SIRV, SEIR, SVEIRS, SIQR, SEIRV, SEIQS, and others ^{[1,2,3,4,5,6,7]}. Understanding disease dynamics is complex due to factors like environmental fluctuations, making it challenging to collect precise data. For example, temperature variations can significantly impact the reproduction rate of disease-spreading bacteria, influencing the disease's dynamics. As a result, it's essential to consider how uncertainty in biological parameters affects disease behavior. If incorrect parameter values are used, model outputs may be biased, so sensitivity analysis is necessary for understanding how model outputs change in response to parameter variations. When studying the transmission of infectious diseases, sensitivity analysis is often used to develop strategies for slowing infection spread by targeting key model parameters, such as the primary reproduction number $R_0$ . The literature provides information on dynamic parameter sensitivities in infectious spread for traditional SIR and SEIR models to maximize the usefulness of observed data, with the primary source of uncertainty being the random variability of the included parameters.

Interactions between species can impact the quality of life for all involved. These interactions are crucial for an organism's survival since its primary objective is to stay alive. Without these relationships, life as we know it would be impossible. It's remarkable how fascinating organism interactions can be and can be studied using mathematical models, an exciting field within ecology research ^[5,7]. Different factors, such as environmental conditions, impact the growth and balance of living populations, and modeling can be an invaluable tool for understanding their dynamics. Organisms, typically proto-microorganisms and parasites can cause infectious illnesses through direct or indirect transmission from one host to another. Due to their reason to periodically spread like epidemics and cause serious problems, infectious diseases have recently become a significant worry in our culture. The three major components that lessen the severity of these disorders are infectivity, epidemic potential, and uncertainty. The medical community has made great strides in recognising and managing these components. The dynamic transmission cycle is a mathematical model of an infectious illness that involves interactions between susceptible and infected hosts. Typically, these interactions are expressed as a linked set of ordinary differential equations.

Fractional calculus is a mathematical analysis that extends beyond differentiation and integration conventions. Fractional calculus concerns derivatives and integrals of non-integer orders, such as fractions or decimals, while classical calculus concentrates on the powers of integers. The Caputo operator is a valuable tool for problem-solving and is commonly used to handle real-world challenges. The main difference between the Riemann-Liouville operator and other fractional operators is the singularity quality of the kernels. They are applied in roughly related ways to observe different model dynamics. It follows that using these operations would inevitably result in better outcomes. Researchers attempting to understand better a model's dynamics advise using fractional operators with non-singular kernels. Establishing the C-F derivative has addressed certain standard limits on fractional derivatives ^[8]. The C-F derivative outperforms the most commonly used techniques for assessing human knowledge ^[9]. The operator's performance and its benefits were validated using various approaches. Studies on numerous scientific, engineering, and mathematical models have proved their usefulness. However, experts from real-world cases are a more meaningful example ^[10,11]. Atangana and Baleanu introduced a derivative operator with a non-singular kernel. It utilizes the Mittag-Leffer function and is in a convenient area for modellers of real-world situations because of its non-local and non-singular kernel. Over the past few decades, the Atangana-Baleanu variant has obtained strong research value and applicability in diverse fields, particularly in biological models ^[12]. This fractional derivative introduced non-local and complex dynamic behavior and numerous natural results. In addition, since biological models are expected to have a memory effect and hereditary features, which can be defined more precisely using fractional calculus, the solution of the fractional order system is expected to obtain how to control epidemic diseases ^[13]. The authors of this study investigated the sensitivity of the model solution in delay differential systems using variational and direct methods ^[14]. Corruption and terrorism have become significant issues in many nations worldwide. However, more needs to be written about the relationship between the two. To address this, the author has developed a novel fractional-order mathematical model to explore the coexistence of terrorism and corruption ^[15]. This article examines the authors' mathematical model of COVID-19, which incorporates a fractional-order system and considers the effectiveness of vaccination ^[16].

Fractional-order models have become increasingly popular in control engineering, physics, neural networks, and medicine. They can predict the condition of a system at any future moment by considering its current state and all of its past states. Due to its practicality, researchers are interested in fractional-order calculus. The concept of fractional-order differentiation and integration was first introduced by two prominent figures in mathematics history, Riemann and Liouville. The R-L and Caputo fractional derivatives are examples of fractional differential operators that use a power law kernel. These operators have precise definitions and have been widely researched and applied across various disciplines ^[17]. However, their limited applicability and inability to represent other natural and artificial systems have hindered their development. To overcome these limitations, novel C-F fractional-order derivatives have been developed based on the exponential kernel, successfully describing various real-world events ^[18]. Several mathematical models with C-F fractional order have been discussed in the literature.

Numerous numerical techniques have been developed to solve biological models and fractional-order differential equations ^[19]. The Newton interpolation formula, the Toufik-Atangana approach, the Adam-Bashforth method, the predictor-corrector method, and the Runge-Kutta technique are a few of the numerous numerical algorithms that have been the subject of extensive literature ^{[17,18,19,20,21,22,23]}. To solve non-linear systems numerically, the Toufik-Atangana approach is widely used ^[24,25]. Recently, a novel numerical system called the Toufik-Atangana numerical system has emerged as a potential solution to the problems with the Adams-Bashforth approach ^[26,27]. Combining the fundamental theorem of fractional calculus with the two-step Lagrange polynomial, this new methodology develops a novel and very efficient numerical method ^[26]. Using this approach, problems can be resolved quickly and accurately. A biological model has been subjected to many fractional operators using the Toufik-Atangana numerical framework to get the required findings. Stability analysis is an essential method in numerical analysis that yields effective results. When expressed mathematically, stability is necessary for spectrum analysis of real-world problems ^[28]. Hyers expounded on Ulam's 1940 introduction of the U-H stability hypothesis in 1941. As previously stated, the optimal approximation or exact solution to the issue is the basis for computing the stability. Furthermore, it is simple to implement and evaluate the recommended stability ideas ^[28,29].

The following subjects in this work: The definitions and basic ideas of fractional calculus are covered in Section 2. We examine a non-integer C-F model for infectious diseases in Section 3. The positivity and sensitivity analysis of the model under discussion is covered in detail in Section 4. In Section 5, we focuses on the proposed model solution's existence, uniqueness, and stability in the U-H. In Section 6, we present the numerical method for the C-F operator Adams-Bashforth scheme. Section 7 features numerical findings and a graphic analysis. Finally, in Section 8, we discuss the findings of the investigation.

2. Definitions and basic concepts

The essential definitions and theorems relevant to the C-F have been provided in this section.

Definition 2.1. ^[18,30] Suppose that $\vartheta$ $\in$ $H^1(q_1, q_2),$ $q_1 < q_2,$ and $\varsigma \in (0, 1),$ so the C-F fractional differential operator. Then, we have

$\begin{align} {}^{CF}{}_{q_1}D_{{\mathfrak{k}}}^{\varsigma}[\vartheta({\mathfrak{k}})] = \frac{B(\varsigma)}{1-{\varsigma}}\int_{q_1}^{\mathfrak{k}} \vartheta^{'}(\mathfrak{z}){exp\left[ -\varsigma \frac{{\mathfrak{k}}-\mathfrak{z}}{1-\varsigma}\right] }d{\mathfrak{z}}, \end{align}$

(2.1)

where $B(\varsigma)$ is the normalization function with $B(0) = B(1) = 1.$

However, if $\vartheta$ $\notin$ $H^1(q_1, q_2)$ , then we have

$\begin{align} {}^{CF}{}_{q_1}D_{{\mathfrak{k}}}^{\varsigma}[\vartheta({\mathfrak{k}})] = \frac{B(\varsigma)}{1-{\varsigma}}\int_{q_1}^{\mathfrak{k}} \left( \vartheta({{\mathfrak{k}}}) - \vartheta(\mathfrak{z})\right) {exp\left[ -\varsigma \frac{{\mathfrak{k}}-\mathfrak{z}}{1-\varsigma}\right] }d{\mathfrak{z}}. \end{align}$

(2.2)

Definition 2.2. ^[18,31] Suppose that $\vartheta$ $\in$ $H^1(q_1, q_2),$ $q_1 < q_2,$ and $\varsigma \in (0, 1),$ so the C-F fractional differential operator. Then, we have

$\begin{align} {}^{CF}{}_{q_1}D_{{\mathfrak{k}}}^{\varsigma}[\vartheta({\mathfrak{k}})] = \frac{(2-\varsigma)B(\varsigma)}{2(1-{\varsigma})}\int_{q_1}^{\mathfrak{k}} \vartheta^{'}(\mathfrak{z}){exp\left[ -\varsigma \frac{{\mathfrak{k}}-\mathfrak{z}}{1-\varsigma}\right] }d{\mathfrak{z}}. \end{align}$

(2.3)

Definition 2.3. ^[18,32] Suppose that $\varsigma$ is the order of integral in the C-F integral operator. Then, we have

$\begin{align} {}^{CF}{}_{\; }I^{\varsigma}_{{\mathfrak{k}}}[\vartheta({\mathfrak{k}})] = \frac{2(1-\varsigma)}{(2-\varsigma)B({\varsigma})}\vartheta ({\mathfrak{k}}) + \frac{2\varsigma}{(2-\varsigma)B({\varsigma})}\int_{q_1}^{{\mathfrak{k}}}\vartheta (\mathfrak{z})d\mathfrak{z}. \end{align}$

(2.4)

3. Fractional order of infectious diseases model

This section presents the generalized form of the infectious spread transmission dynamics with specific parameter values for a given population. In this paper, we build upon the fundamental SEIR contagious disease model by expanding it to the SEQIRDV model ^[33], which considers different disease stages and traits. There are seven subpopulations within the overall population $N$ : susceptible $(S)$ , infected $(I)$ , exposed $(E)$ , recovered $(R)$ , quarantined $(Q)$ , dead $(D)$ , and vaccinated $(V)$ . Realistically, the natural death rate caused for each subpopulation is represented by a value $\mu$ . A second parameter, $\Pi$ , indicates the recruitment of susceptible individuals in any infected population at any time $\mathfrak{k}$ . At the disease transmission rate $\beta$ , susceptible people can be exposed and become members of the diseased class. We assume that specific individuals receive vaccinations for a particular infectious disease at a rate of $\nu$ . The rate at which people contract the disease is $\beta$ , and they move on to the affected group after a latent period of $\gamma$ . Based on the vaccination's efficiency measure $\sigma$ , the interaction of those who have had vaccinations falls into this category. With a predetermined period until death $\rho$ , the infected population is confined for a proposed average length of $\delta$ , either entering the dead population with a disease mortality rate of $\tau$ or the recovered class with a recovery rate of $\omega$ . Depending on the availability of vaccines for the particular condition, the class "vaccinated" is added. The flowchart of the infectious dynamical disease system (3.1) is represented by Figure 1. After accounting for these variables, the dynamical system that results is as follows ^[33]:

$\begin{align} \frac{dS}{d{\mathfrak{k}}} & = \Pi - \beta SI - \upsilon S - \mu S, \\ \frac{dE}{d{\mathfrak{k}}} & = \beta SI - \gamma E + \sigma \beta VI - \mu E, \\ \frac{dI}{d{\mathfrak{k}}} & = \gamma E - \delta I - \mu I, \\ \frac{dQ}{d{\mathfrak{k}}} & = \delta I - (1-\tau)\omega Q - \tau \rho Q - \mu Q, \\ \frac{dR}{d{\mathfrak{k}}} & = (1-\tau)\omega Q -\mu R, \\ \frac{dD}{d{\mathfrak{k}}} & = \tau \rho Q, \\ \frac{dV}{d{\mathfrak{k}}} & = \upsilon S - \sigma \beta VI - \mu V, \end{align}$

(3.1)

Figure 1. A diagrammatic graph that illustrates the transmission of the disease.

DownLoad: Full-Size Img PowerPoint

with initial conditions $S(0) \geq 0,$ $E(0) \geq 0,$ $I(0) \geq 0,$ $Q(0) \geq 0,$ $R(0) \geq 0,$ $D(0) \geq 0$ and $V(0) \geq 0.$

Fractional derivatives in mathematical biology offer a versatile framework for modeling various biological phenomena, including vaccination processes and recovery dynamics. By extending beyond traditional epidemic and noninfectious disease models, researchers gain a deeper understanding of the complex behaviors exhibited by biological systems. This enables them to develop innovative solutions to critical vaccination research and healthcare issues by incorporating fractional calculus into their mathematical models. Initially, we transformed the model with integer order into one with fractional order. We apply the C-F derivative in the system (3.1), and then we get

$\begin{align} {}^{CF}{}_{0}D^{\varsigma}_{{\mathfrak{k}}}S & = \Pi - \beta SI - \upsilon S - \mu S, \\ {}^{CF}{}_{0}D^{\varsigma}_{{\mathfrak{k}}}E & = \beta SI - \gamma E + \sigma \beta VI - \mu E, \\ {}^{CF}{}_{0}D^{\varsigma}_{{\mathfrak{k}}}I & = \gamma E - \delta I - \mu I, \\ {}^{CF}{}_{0}D^{\varsigma}_{{\mathfrak{k}}}Q & = \delta I - (1-\tau)\omega Q - \tau \rho Q - \mu Q, \\ {}^{CF}{}_{0}D^{\varsigma}_{{\mathfrak{k}}}R & = (1-\tau)\omega Q -\mu R, \\ {}^{CF}{}_{0}D^{\varsigma}_{{\mathfrak{k}}}D & = \tau \rho Q, \\ {}^{CF}{}_{0}D^{\varsigma}_{{\mathfrak{k}}}V & = \upsilon S - \sigma \beta VI - \mu V. \end{align}$

(3.2)

4. Preliminary analysis of the model

4.1. Positivity and boundedness

In fractional order modeling of biological systems, the population remains positive and bounded over time to represent biological constraints accurately. Positivity and boundedness are essential in fractional order models to ensure stability, realistic behavior, and relevance to real-world situations. These features are necessary due to the complex dynamics introduced by fractional calculus. Positive states play a crucial role in maintaining a system's stability and meaningful behavior. In contrast, positive outputs are vital for applying and understanding the model's results in practical applications. Boundedness ensures that the system does not exhibit unlimited growth, which is essential for stability and control. Furthermore, boundedness guarantees that the behavior of a system remains predictable and controllable, thus facilitating the analysis of system performance and reaction. The proposed model has a built-in feature where its solutions are always positive and bounded. We ensure that all state variables have non-negative values for any time ${\mathfrak{k}} > 0$ . This ^[34] means a trajectory starting with a positive initial condition will stay favorable for ${\mathfrak{k}} > 0$ . Therefore, system (3.2) provides that

$\begin{align} {}^{CF}{}_{0}D^{\varsigma}_{{\mathfrak{k}}}S|_{S = 0} & = \Pi \geq 0, \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; {}^{CF}{}_{0}D^{\varsigma}_{{\mathfrak{k}}}E|_{E = 0} = \sigma \beta VI \geq 0, \\ {}^{CF}{}_{0}D^{\varsigma}_{{\mathfrak{k}}}I|_{I = 0} & = \gamma E \geq 0, \; \; \; \; \; \; \; \; \; \; \; \; \; {}^{CF}{}_{0}D^{\varsigma}_{{\mathfrak{k}}}Q|_{Q = 0} = \delta I \geq 0, \\ {}^{CF}{}_{0}D^{\varsigma}_{{\mathfrak{k}}}R|_{R = 0} & = (1-\tau)\omega Q \geq 0, \; \; {}^{CF}{}_{0}D^{\varsigma}_{{\mathfrak{k}}}D|_{D = 0} = \tau \rho Q \geq 0, \\ {}^{CF}{}_{0}D^{\varsigma}_{{\mathfrak{k}}}V|_{V = 0} & = \upsilon S \geq 0. \end{align}$

(4.1)

Since $\mathbb{N}(t) = S({\mathfrak{k}}) + I({\mathfrak{k}}) + E({\mathfrak{k}}) + Q({\mathfrak{k}}) + D({\mathfrak{k}}) + R({\mathfrak{k}}) + V({\mathfrak{k}})$ is total population, we have

$\begin{align*} {}^{CF}{}_{0}D^{\varsigma}_{{\mathfrak{k}}}\mathbb{N}({\mathfrak{k}}) & = \Pi -\mu S - \mu E - \mu I - \mu Q - \mu R - \mu V \nonumber \\ & \leq \Pi - \mu S, \end{align*}$

then one has

$\begin{align*} \mathbb{N}({\mathfrak{k}}) \leq \left( \mathbb{N}(0) - \frac{\Pi}{\mu}\right)E_{\varsigma}(-\mu {\mathfrak{k}}) + \frac{\Pi}{\mu}. \end{align*}$

Therefore, we have

$\begin{align*} \Omega = \left\lbrace \left( S({\mathfrak{k}}), I({\mathfrak{k}}), E({\mathfrak{k}}), Q({\mathfrak{k}}), D({\mathfrak{k}}), R({\mathfrak{k}}), V({\mathfrak{k}})\right) \in \Re^7_+ : 0 \leq \mathbb{N}({\mathfrak{k}}) \leq \frac{\Pi}{\mu} \right\rbrace, \end{align*}$

which provides the feasible region for the infectious diseases model, and $\Omega$ is positively invariant. Thus, the proposed model (3.2) is well-posed mathematically.

4.2. Sensitivity analysis

This study gives us insights into each parameter's importance in the disease's spread. This information is crucial not only for planning experiments but also for integrating data and simplifying complex models. Sensitivity analysis is often used to assess the resilience of model predictions to changes in parameter values, as there are sometimes errors in data collection and assumed parameter values. This tool helps identify the characteristics that significantly impact the threshold $R_0$ and should be the focus of intervention strategies. Sensitivity indices provide a way to measure the proportional change in a variable resulting from a change in the parameter. We use a variable's normalized forward sensitivity index concerning a specific parameter to achieve this objective. This index is defined as the ratio of the relative change in the variable to the relative change in the parameter. The sensitivity index can be determined using partial derivatives if a variable is differentiable concerning the parameter. We conducted a sensitivity study in this section to determine the impact of each parameter on the $R_0$ . A sensitivity analysis can also determine the critical parameters vital to illness management. This method ascertains the relative contribution of each parameter value to the $R_0$ . Therefore, we have

$\begin{align*} R_0 = \frac{\gamma \beta \Pi (\mu + \sigma \nu )}{\mu (\nu + \mu)(\delta + \mu)(\gamma + \mu)} . \end{align*}$

The system's qualitative characteristics are being considered while organising the relevant parameters based on their impact on the value of $R_0$ , which is quite beneficial. The results of this investigation will help identify the most crucial disease control factors. We can utilize sensitivity indices to determine the effect of a parameter change on the state variable. These indices are calculated using the definition provided in the paper ^[34,35]. To estimate the sensitivity index, we use partial derivatives, as shown below:

$\begin{align} \Lambda^{R_0}_{\amalg} = \frac{\partial R_0}{\partial \amalg}\times \frac{\amalg}{R_0}. \end{align}$

(4.2)

The sensitivity indices for different parameters are present in . These indices have been calculate using the starting values, except for $\sigma$ . shows the sensitivity indices of $R_0$ for the considered parameters of interest. The results indicate that $\Pi$ , $\beta$ , $\delta$ , and $\mu$ are highly significant characteristics, as seen in and . Based on the sensitivity analysis results, it was found that the $R_0$ values increase or decrease in direct proportion to changes in the values of $\Pi,$ $\beta,$ $\sigma,$ $\delta,$ $\mu,$ $\nu$ and $\gamma$ . The reproduction number within the defined limits has been analyzed, considering the following factors. After reviewing the analyses and graphs, it is evident that precautionary measures should be taken to prevent the spread of the disease by improving adverse conditions and reducing factors that contribute to increased reproduction numbers. is the plot of $R_0$ versus the $\beta$ and $\Pi$ . is the plot of $R_0$ versus the $\gamma$ and $\delta$ . is the plot of $R_0$ versus the $\Pi$ and $\mu$ . is the plot of $R_0$ versus the $\sigma$ and $\nu$ . Notably, we have observed that $\mu$ is quite sensitive, and raising this parameter could result in a notable drop in the $R_0$ value. Consequently, limiting these factors can aid in the prevention of the spread of illness. The variables are considered when examining the relevant reproduction number within the specified limits. Based on the evaluations and visual representations, it is concluded that appropriate measures should be implemented to prevent the spread of the disease. This can be achieved by decreasing the factors that contribute to the positive growth rate of reproduction numbers and raising the factors that have a detrimental effect on it. Hence, it can be concluded that increasing awareness about quarantine and vaccination among affected individuals can significantly reduce the spread of infection within the population.

Table 1. Sensitivity indices.

Parameters	Sensitivity index
$\Pi$	$1$
$\beta$	$1$
$\sigma$	$0.1429$
$\delta$	$-0.9444$
$\mu$	$-1.2506$
$\gamma$	$0.0556$
$\nu$	$0.1395$

| Show Table

DownLoad: CSV

Figure 2. Sensitivity plot for the basic reproduction number

$R_0$ .

DownLoad: Full-Size Img PowerPoint

Figure 3. Sensitivity analysis of

$R_0$ according to the model parameters for (a)

$\beta$ and

$\Pi$ , (b)

$\gamma$ and

$\delta$ (c)

$\Pi$ and

$\mu$ , (d)

$\sigma$ and

$\nu$ .

DownLoad: Full-Size Img PowerPoint

5. Essential existence and uniqueness of the solutions

Here, we examine the uniqueness and existence of solutions using the fixed point theorem, which is essential for the proposed model ^[31,36]. The C-F fractional integral operator, when applied to system (3.2), yields

$\begin{align} S({\mathfrak{k}}) - S(0) = &\frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\left[ \Pi - \beta SI - \upsilon S - \mu S \right] \\ & + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\left[ \Pi - \beta S(\mathfrak{Y})I(\mathfrak{Y}) - \upsilon S(\mathfrak{Y}) - \mu S(\mathfrak{Y}) \right] d(\mathfrak{Y}), \\ E({\mathfrak{k}}) - E(0) = & \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\left[ \beta SI - \gamma E + \sigma \beta VI - \mu E \right] \\ & + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\left[ \beta S(\mathfrak{Y})I(\mathfrak{Y}) - \gamma E(\mathfrak{Y}) + \sigma \beta V(\mathfrak{Y})I(\mathfrak{Y}) - \mu E(\mathfrak{Y})\right] d(\mathfrak{Y}), \\ I({\mathfrak{k}}) - I(0) = & \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\left[ \gamma E - \delta I - \mu I \right] + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\left[ \gamma E(\mathfrak{Y}) - \delta I(\mathfrak{Y}) - \mu I(\mathfrak{Y})\right] d(\mathfrak{Y}), \\ Q({\mathfrak{k}}) - Q(0) = & \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\left[ \delta I - (1-\tau)\omega Q - \tau \rho Q - \mu Q \right] \\ & + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\left[ \delta I(\mathfrak{Y}) - (1-\tau)\omega Q(\mathfrak{Y}) - \tau \rho Q(\mathfrak{Y}) - \mu Q(\mathfrak{Y}) \right] d(\mathfrak{Y}), \\ R({\mathfrak{k}}) - R(0) = & \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\left[ \gamma E - \delta I - \mu I \right] + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\left[ \gamma E(\mathfrak{Y}) - \delta I(\mathfrak{Y}) - \mu I(\mathfrak{Y})\right] d(\mathfrak{Y}), \\ D({\mathfrak{k}}) - D(0) = & \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\left[ \tau \rho Q \right] + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\left[ \tau \rho Q(\mathfrak{Y})\right] d(\mathfrak{Y}), \\ V({\mathfrak{k}}) - V(0) = &\frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\left[ \upsilon S - \sigma \beta VI - \mu V \right] \\ & + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\left[ \upsilon S(\mathfrak{Y}) - \sigma \beta V(\mathfrak{Y})I(\mathfrak{Y}) - \mu V(\mathfrak{Y}) \right] d(\mathfrak{Y}). \end{align}$

(5.1)

Define the following kernels:

$\begin{align} M_1({\mathfrak{k}}, S) & = \Pi - \beta SI - \upsilon S - \mu S, \\ M_2({\mathfrak{k}}, E) & = \beta SI - \gamma E + \sigma \beta VI - \mu E, \\ M_3({\mathfrak{k}}, I) & = \gamma E - \delta I - \mu I, \\ M_4({\mathfrak{k}}, Q) & = \delta I - (1-\tau)\omega Q - \tau \rho Q - \mu Q, \\ M_5({\mathfrak{k}}, R) & = (1-\tau)\omega Q -\mu R, \\ M_6({\mathfrak{k}}, D) & = \tau \rho Q, \\ M_7({\mathfrak{k}}, V) & = \upsilon S - \sigma \beta VI - \mu V, \end{align}$

(5.2)

then we have

$\begin{align} S({\mathfrak{k}}) - S(0) & = \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\left[ M_1({\mathfrak{k}}, S) \right] + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\left[ M_1(\mathfrak{Y}, S) \right] d\mathfrak{Y}, \\ E({\mathfrak{k}}) - E(0) & = \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\left[ M_2({\mathfrak{k}}, E) \right] + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\left[ M_2(\mathfrak{Y}, E)\right] d\mathfrak{Y}, \\ I({\mathfrak{k}}) - I(0) & = \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\left[ M_3({\mathfrak{k}}, I) \right] + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\left[ M_3(\mathfrak{Y}, I) \right] d\mathfrak{Y}, \\ Q({\mathfrak{k}}) - Q(0) & = \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\left[ M_4({\mathfrak{k}}, Q) \right] + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\left[ M_4(\mathfrak{Y}, Q) \right] d\mathfrak{Y}, \\ R({\mathfrak{k}}) - R(0) & = \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\left[ M_5({\mathfrak{k}}, R) \right] + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\left[ M_5(\mathfrak{Y}, R)\right] d\mathfrak{Y}, \\ D({\mathfrak{k}}) - D(0) & = \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\left[ M_6({\mathfrak{k}}, D) \right] + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\left[ M_6(\mathfrak{Y}, D)\right] d\mathfrak{Y}, \\ V({\mathfrak{k}}) - V(0) & = \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\left[ M_7({\mathfrak{k}}, V) \right] + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\left[ M_7(\mathfrak{Y}, V) \right] d\mathfrak{Y}. \end{align}$

(5.3)

Theorem 5.1. If the dissimilarity listed below holds:

$\begin{equation*} 0 \leq \beta \varpi_2 + \upsilon + \mu < 1. \end{equation*}$

Afterwards, the contraction mapping and Lipschitz condition are convinced by the kernel $M_1$ .

Proof. Assume that $S$ and $S_1$ are any two functions, then we get

$\begin{align*} \|M_1({\mathfrak{k}}, S) - M_1({\mathfrak{k}}, S_1)\| & = \left\| \left( \Pi - \beta S({\mathfrak{k}})I({\mathfrak{k}}) - \upsilon S({\mathfrak{k}}) - \mu S({\mathfrak{k}}) \right) \right. - \left. \left( \Pi - \beta S_1({\mathfrak{k}})I({\mathfrak{k}}) - \upsilon S_1({\mathfrak{k}}) - \mu S_1({\mathfrak{k}})\right) \right\| \\ & \leq \left( \beta \varpi_2 + \upsilon + \mu \right)\left\|S({\mathfrak{k}}) - S_1({\mathfrak{k}}) \right\| \\ & \leq {\Upsilon}_1\left\|S({\mathfrak{k}}) - S_1({\mathfrak{k}}) \right\| . \end{align*}$

Let ${\Upsilon}_1 = \beta \varpi_2 + \upsilon + \mu$ , we do suppose that $S, E, I, Q, R, D$ and $V$ is bounded functions, i.e., $\|S({\mathfrak{k}})\| \leq \varpi_1,$ $\|E({\mathfrak{k}})\|\leq \varpi_2,$ $\|I({\mathfrak{k}})\|\leq \varpi_3,$ $\|Q({\mathfrak{k}})\|\leq \varpi_4,$ $\|R({\mathfrak{k}})\|\leq \varpi_5,$ $\|D({\mathfrak{k}})\|\leq \varpi_6$ and $\|V({\mathfrak{k}})\|\leq \varpi_7.$

$\begin{align} \|M_1({\mathfrak{k}}, S) - M_1({\mathfrak{k}}, S_1)\| & \leq \Upsilon_1 \left\| S({\mathfrak{k}}) - S_1({\mathfrak{k}})\right\| . \end{align}$

(5.4)

Hence, kernel $M_1$ satisfies the Lipschitz condition and since

$0 \leq \beta \varpi_2 + \upsilon + \mu < 1,$

then it is also a contraction for $M_1$ .

Similarly, kernels $M_2,$ $M_3,$ $M_4,$ $M_5,$ $M_6$ , and $M_7$ satisfy the Lipschitz condition, respectively, as follows:

$\begin{align} \|M_2({\mathfrak{k}}, E) - M_2({\mathfrak{k}}, E_{1})\| & \leq \Upsilon_2 \left\| E({\mathfrak{k}}) - E_{1}({\mathfrak{k}})\right\|, \\ \|M_3({\mathfrak{k}}, I) - M_3({\mathfrak{k}}, I_{1})\| & \leq \Upsilon_3 \left\| I({\mathfrak{k}}) - I_{1}({\mathfrak{k}})\right\|, \\ \|M_4({\mathfrak{k}}, Q) - M_4({\mathfrak{k}}, Q_{1})\| & \leq \Upsilon_4 \left\| Q({\mathfrak{k}}) - Q_{1}({\mathfrak{k}})\right\|, \\ \|M_5({\mathfrak{k}}, R) - M_5({\mathfrak{k}}, R_{1})\| & \leq \Upsilon_5 \left\| R({\mathfrak{k}}) - R_{1}({\mathfrak{k}})\right\|, \\ \|M_6({\mathfrak{k}}, D) - M_6({\mathfrak{k}}, D_{1})\| & \leq \Upsilon_6 \left\| D({\mathfrak{k}}) - D_{1}({\mathfrak{k}})\right\|, \\ \|M_7({\mathfrak{k}}, V) - M_7({\mathfrak{k}}, V_{1})\| & \leq \Upsilon_7 \left\| V({\mathfrak{k}}) - V_{1}({\mathfrak{k}})\right\|. \end{align}$

(5.5)

Equations (5.4) and (5.5) are used to apply the previously described kernels, which transform system (5.1) into

$\begin{align} S({\mathfrak{k}}) & = S(0) + \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}M_1({\mathfrak{k}}, S) + \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}} M_1(\mathfrak{Y}, S)d\mathfrak{Y}, \\ E({\mathfrak{k}}) & = E(0) + \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}M_2({\mathfrak{k}}, E) + \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}} M_2(\mathfrak{Y}, E)d\mathfrak{Y}, \\ I({\mathfrak{k}}) & = I(0) + \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}M_3({\mathfrak{k}}, I) + \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}} M_3(\mathfrak{Y}, I)d\mathfrak{Y}, \\ Q({\mathfrak{k}}) & = Q(0) + \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}M_4({\mathfrak{k}}, S) + \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}} M_4(\mathfrak{Y}, Q)d\mathfrak{Y}, \\ R({\mathfrak{k}}) & = R(0) + \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}M_5({\mathfrak{k}}, R) + \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}} M_5(\mathfrak{Y}, R)d\mathfrak{Y}, \\ D({\mathfrak{k}}) & = D(0) + \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}M_6({\mathfrak{k}}, D) + \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}} M_6(\mathfrak{Y}, D)d\mathfrak{Y}, \\ V({\mathfrak{k}}) & = V(0) + \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}M_7({\mathfrak{k}}, V) + \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}} M_7(\mathfrak{Y}, V)d\mathfrak{Y}. \end{align}$

(5.6)

We now introduce the following recursive formulas:

$\begin{align} S_r({\mathfrak{k}}) & = \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}M_1({\mathfrak{k}}, S_{r-1}) + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}M_1(\mathfrak{Y}, S_{r-1})d\mathfrak{Y}, \\ E_r({\mathfrak{k}}) & = \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}M_2({\mathfrak{k}}, E_{r-1}) + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}M_2(\mathfrak{Y}, E_{r-1})d\mathfrak{Y}, \\ I_r({\mathfrak{k}}) & = \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}M_3({\mathfrak{k}}, I_{r-1}) + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}M_3(\mathfrak{Y}, I_{r-1})d\mathfrak{Y}, \\ Q_r({\mathfrak{k}}) & = \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}M_4({\mathfrak{k}}, Q_{r-1}) + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}M_4(\mathfrak{Y}, Q_{r-1})d\mathfrak{Y}, \\ R_r({\mathfrak{k}}) & = \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}M_5({\mathfrak{k}}, R_{r-1}) + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}M_5(\mathfrak{Y}, R_{r-1})d\mathfrak{Y}, \\ D_r({\mathfrak{k}}) & = \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}M_6({\mathfrak{k}}, D_{r-1}) + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}M_6(\mathfrak{Y}, D_{r-1})d\mathfrak{Y}, \\ V_r({\mathfrak{k}}) & = \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}M_7({\mathfrak{k}}, V_{r-1}) + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}M_7(\mathfrak{Y}, V_{r-1})d\mathfrak{Y}, \end{align}$

(5.7)

where the initial condition are

$\begin{align*} S_0({\mathfrak{k}}) & = S(0), \; E_0({\mathfrak{k}}) = E(0), \; I_0({\mathfrak{k}}) = I(0), \; Q_0({\mathfrak{k}}) = Q(0), \\ R_0({\mathfrak{k}}) & = R(0), \; D_0({\mathfrak{k}}) = D(0), \; V_0({\mathfrak{k}}) = V(0). \nonumber \end{align*}$

With regard to the recursive formulas, the differences between successive terms can be represented as

$\begin{align*} \chi_{1, r}({\mathfrak{k}}) & = S_r({\mathfrak{k}}) - S_{r-1}({\mathfrak{k}}) \nonumber \\ & = \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\left[ M_1({\mathfrak{k}}, S_{r-1}) - M_1({\mathfrak{k}}, S_{r-2})\right] \nonumber \\ &\quad + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\left[ M_1(\mathfrak{Y}, S_{r-1})- M_1(\mathfrak{Y}, S_{r-2})\right] d\mathfrak{Y}, \nonumber \end{align*}$

$\begin{align*} \chi_{2, r}({\mathfrak{k}}) & = E_r({\mathfrak{k}}) - E_{r-1}({\mathfrak{k}}) \nonumber \\ & = \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\left[ M_2({\mathfrak{k}}, E_{r-1}) - M_2({\mathfrak{k}}, E_{r-2})\right] \nonumber \\ & \quad + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\left[ M_2(\mathfrak{Y}, E_{r-1})- M_2(\mathfrak{Y}, E_{r-2})\right] d\mathfrak{Y}, \nonumber \end{align*}$

$\begin{align*} \chi_{3, r}({\mathfrak{k}}) & = I_r({\mathfrak{k}}) - I_{r-1}({\mathfrak{k}}) \nonumber \\ & = \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\left[ M_3({\mathfrak{k}}, I_{r-1}) - M_3({\mathfrak{k}}, I_{r-2})\right] \nonumber \\ &\quad + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\left[ M_3(\mathfrak{Y}, I_{r-1})- M_3(\mathfrak{Y}, I_{r-2})\right] d\mathfrak{Y}, \nonumber \end{align*}$

$\begin{align} \chi_{4, r}({\mathfrak{k}}) & = Q_r({\mathfrak{k}}) - Q_{r-1}({\mathfrak{k}}) \\ & = \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\left[ M_4({\mathfrak{k}}, Q_{r-1}) - M_4({\mathfrak{k}}, Q_{r-2})\right] \\ & \quad+ \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\left[ M_4(\mathfrak{Y}, Q_{r-1})- M_4(\mathfrak{Y}, Q_{r-2})\right] d\mathfrak{Y}, \end{align}$

(5.8)

$\begin{align*} \chi_{5, r}({\mathfrak{k}}) & = R_r({\mathfrak{k}}) - R_{r-1}({\mathfrak{k}}) \nonumber \\ & = \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\left[ M_5({\mathfrak{k}}, R_{r-1}) - M_5({\mathfrak{k}}, R_{r-2})\right] \nonumber \\ &\quad + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\left[ M_5(\mathfrak{Y}, R_{r-1})- M_5(\mathfrak{Y}, R_{r-2})\right] d\mathfrak{Y}, \nonumber \end{align*}$

$\begin{align*} \chi_{6, r}({\mathfrak{k}}) & = D_r({\mathfrak{k}}) - D_{r-1}({\mathfrak{k}}) \nonumber \\ & = \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\left[ M_6({\mathfrak{k}}, D_{r-1}) - M_6({\mathfrak{k}}, D_{r-2})\right] \nonumber \\ &\quad + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\left[ M_6(\mathfrak{Y}, D_{r-1})- M_6(\mathfrak{Y}, D_{r-2})\right] d\mathfrak{Y}, \nonumber \end{align*}$

$\begin{align*} \chi_{7, r}({\mathfrak{k}}) & = V_r({\mathfrak{k}}) - V_{r-1}({\mathfrak{k}}) \nonumber \\ & = \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\left[ M_7({\mathfrak{k}}, V_{r-1}) - M_7({\mathfrak{k}}, V_{r-2})\right] \nonumber \\ & \quad+ \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\left[ M_7(\mathfrak{Y}, V_{r-1})- M_7(\mathfrak{Y}, V_{r-2})\right] d\mathfrak{Y}. \nonumber \end{align*}$

Note that

$\begin{align} S_r({\mathfrak{k}}) & = \sum\limits_{j = 1}^{r} \chi_{1, j}, \; \; E_r({\mathfrak{k}}) = \sum\limits_{j = 1}^{r} \chi_{2, j}, \; \; I_r({\mathfrak{k}}) = \sum\limits_{j = 1}^{r} \chi_{3, j}, \; \; Q_r({\mathfrak{k}}) = \sum\limits_{j = 1}^{r} \chi_{4, j}, \\ R_r({\mathfrak{k}}) & = \sum\limits_{j = 1}^{r} \chi_{5, j}, \; \; D_r({\mathfrak{k}}) = \sum\limits_{j = 1}^{r} \chi_{6, j}, \; \; V_r({\mathfrak{k}}) = \sum\limits_{j = 1}^{r} \chi_{7, j}. \end{align}$

(5.9)

After solving system (5.8), we apply the usual supremum norm to both sides of the first equation of system (5.8), then we obtain

$\begin{align} \left\| \chi_{1, r}({\mathfrak{k}})\right\| & = \left\| S_r({\mathfrak{k}}) - S_{r-1}({\mathfrak{k}})\right\| \\ & = \left\| \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\left[ M_1({\mathfrak{k}}, S_{r-1}) - M_1({\mathfrak{k}}, S_{r-2})\right]\right. \\ &\quad \left. + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\left[ M_1(\mathfrak{Y}, S_{r-1})- M_1(\mathfrak{Y}, S_{r-2})\right] d\mathfrak{Y}\right\| . \end{align}$

(5.10)

Using Eq (5.10) and the triangle inequality, we get

$\begin{align} \left\| S_r({\mathfrak{k}}) - S_{r-1}({\mathfrak{k}})\right\| & \leq \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\left\| \left[ M_1({\mathfrak{k}}, S_{r-1}) - M_1({\mathfrak{k}}, S_{r-2})\right]\right\| \\ & \quad + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\left\| \int_{0}^{{\mathfrak{k}}}\left[ M_1(\mathfrak{Y}, S_{r-1})- M_1(\mathfrak{Y}, S_{r-2})\right] d\omega\right\| . \end{align}$

(5.11)

Thus, using the Lipschitz constant $\Upsilon_1$ to propitiate the Lipschitz condition, the kernel $M_1$ enables us to determine

$\begin{align} \|S_r({\mathfrak{k}}) - S_{r-1}({\mathfrak{k}})\| \leq \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\Upsilon_1 \left\| S_{r-1} - S_{r-2}\right\| + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\Upsilon_1 \int_{0}^{{\mathfrak{k}}} \left\|S_{r-1} - S_{r-2}\right\| d\mathfrak{Y} . \end{align}$

(5.12)

Thus, we obtain

$\begin{align} \|\chi_{1, r}({\mathfrak{k}})\| & \leq \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\Upsilon_1\|\chi_{1, r-1}({\mathfrak{k}})\| + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\Upsilon_1 \int_{0}^{{\mathfrak{k}}}\|\chi_{1, r-1}(\mathfrak{Y})\|d\mathfrak{Y}. \end{align}$

(5.13)

In a similar manner, we get

$\begin{align} \|\chi_{2, r}({\mathfrak{k}})\| & \leq \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\Upsilon_2\|\chi_{2, r-1}({\mathfrak{k}})\| + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\Upsilon_2 \int_{0}^{{\mathfrak{k}}}\|\chi_{2, r-1}(\mathfrak{Y})\|d\mathfrak{Y}, \\ \|\chi_{3, r}({\mathfrak{k}})\| & \leq \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\Upsilon_3\|\chi_{3, r-1}({\mathfrak{k}})\| + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\Upsilon_3 \int_{0}^{{\mathfrak{k}}}\|\chi_{3, r-1}(\mathfrak{Y})\|d\mathfrak{Y}, \\ \|\chi_{4, r}({\mathfrak{k}})\| & \leq \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\Upsilon_4\|\chi_{4, r-1}({\mathfrak{k}})\| + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\Upsilon_4 \int_{0}^{{\mathfrak{k}}}\|\chi_{4, r-1}(\mathfrak{Y})\|d\mathfrak{Y}, \\ \|\chi_{5, r}({\mathfrak{k}})\| & \leq \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\Upsilon_5\|\chi_{5, r-1}({\mathfrak{k}})\| + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\Upsilon_5 \int_{0}^{{\mathfrak{k}}}\|\chi_{5, r-1}(\mathfrak{Y})\|d\mathfrak{Y}, \\ \|\chi_{6, r}({\mathfrak{k}})\| & \leq \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\Upsilon_6\|\chi_{6, r-1}({\mathfrak{k}})\| + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\Upsilon_6 \int_{0}^{{\mathfrak{k}}}\|\chi_{6, r-1}(\mathfrak{Y})\|d\mathfrak{Y}, \\ \|\chi_{7, r}({\mathfrak{k}})\| & \leq \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\Upsilon_7\|\chi_{7, r-1}({\mathfrak{k}})\| + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\Upsilon_7 \int_{0}^{{\mathfrak{k}}}\|\chi_{7, r-1}(\mathfrak{Y})\|d\mathfrak{Y}. \end{align}$

(5.14)

Theorem 5.2. If there $\exists$ a time ${\mathfrak{k}}_0 > 0$ , then the ensuing disparities are valid:

$\begin{equation} \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\Upsilon_1 + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\Upsilon_1 {\mathfrak{k}}_0 < 1, \end{equation}$

(5.15)

solutions exist for the infectious disease system.

Proof. Let $S({\mathfrak{k}}),$ $E({\mathfrak{k}}),$ $I({\mathfrak{k}}),$ $Q({\mathfrak{k}}),$ $R(t),$ $D({\mathfrak{k}})$ , and $V({\mathfrak{k}})$ be bounded functions and use the Lipschitz condition. Now, in Eqs (5.13) and (5.14), using the recursive method, we obtain

$\begin{align} \left\| \chi_{1, r}({\mathfrak{k}})\right\| & \leq \left\| S_r(0)\right\| \left[ \frac{2(1-{\varsigma})}{(2-\varsigma)B(\varsigma)}\Upsilon_1 + \frac{2{\varsigma}}{(2-\varsigma)B(\varsigma)}\Upsilon_1{\mathfrak{k}}\right]^r , \\ \left\| \chi_{2, r}({\mathfrak{k}})\right\| & \leq \left\| E_r(0)\right\| \left[ \frac{2(1-{\varsigma})}{(2-\varsigma)B(\varsigma)}\Upsilon_2 + \frac{2{\varsigma}}{(2-\varsigma)B(\varsigma)}\Upsilon_2 {\mathfrak{k}}\right]^r , \\ \left\| \chi_{3, r}({\mathfrak{k}})\right\| & \leq \left\| I_r(0)\right\| \left[ \frac{2(1-{\varsigma})}{(2-\varsigma)B(\varsigma)}\Upsilon_3 + \frac{2{\varsigma}}{(2-\varsigma)B(\varsigma)}\Upsilon_3 {\mathfrak{k}}\right]^r , \\ \left\| \chi_{4, r}({\mathfrak{k}})\right\| & \leq \left\| Q_r(0)\right\| \left[ \frac{2(1-{\varsigma})}{(2-\varsigma)B(\varsigma)}\Upsilon_4 + \frac{2{\varsigma}}{(2-\varsigma)B(\varsigma)}\Upsilon_4 {\mathfrak{k}}\right]^r , \\ \left\| \chi_{5, r}({\mathfrak{k}})\right\| & \leq \left\| R_r(0)\right\| \left[ \frac{2(1-{\varsigma})}{(2-\varsigma)B(\varsigma)}\Upsilon_5 + \frac{2{\varsigma}}{(2-\varsigma)B(\varsigma)}\Upsilon_5 {\mathfrak{k}}\right]^r , \\ \left\| \chi_{6, r}({\mathfrak{k}})\right\| & \leq \left\| D_r(0)\right\| \left[ \frac{2(1-{\varsigma})}{(2-\varsigma)B(\varsigma)}\Upsilon_6 + \frac{2{\varsigma}}{(2-\varsigma)B(\varsigma)}\Upsilon_6 {\mathfrak{k}}\right]^r , \\ \left\| \chi_{7, r}({\mathfrak{k}})\right\| & \leq \left\| V_r(0)\right\| \left[ \frac{2(1-{\varsigma})}{(2-\varsigma)B(\varsigma)}\Upsilon_7 + \frac{2{\varsigma}}{(2-\varsigma)B(\varsigma)}\Upsilon_7 {\mathfrak{k}}\right]^r . \end{align}$

(5.16)

Thus, it is proven that the aforementioned solutions exist and continue to exist. To demonstrate the function is a solution of system (3.2), we make the following assumptions:

$\begin{align} S({\mathfrak{k}}) - S(0) & = S_r({\mathfrak{k}}) - \mathfrak{A}_{1, r}({\mathfrak{k}}), \\ E({\mathfrak{k}}) - E(0) & = E_r({\mathfrak{k}}) - \mathfrak{A}_{2, r}({\mathfrak{k}}), \\ I({\mathfrak{k}}) - I(0) & = I_r({\mathfrak{k}}) - \mathfrak{A}_{3, r}({\mathfrak{k}}), \\ Q({\mathfrak{k}}) - Q(0) & = Q_r({\mathfrak{k}}) - \mathfrak{A}_{4, r}({\mathfrak{k}}), \\ R({\mathfrak{k}}) - R(0) & = R_r({\mathfrak{k}}) - \mathfrak{A}_{5, r}({\mathfrak{k}}), \\ D({\mathfrak{k}}) - D(0) & = D_r({\mathfrak{k}}) - \mathfrak{A}_{6, r}({\mathfrak{k}}), \\ V({\mathfrak{k}}) - V(0) & = V_r({\mathfrak{k}}) - \mathfrak{A}_{7, r}({\mathfrak{k}}). \end{align}$

(5.17)

Then, we have

$\begin{align} \left\| \mathfrak{A}_{1, r}({\mathfrak{k}})\right\| & = \left\| \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\left[ M_1({\mathfrak{k}}, S) - M_1({\mathfrak{k}}, S_{r-1})\right]\right. \\ &\quad \left. + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\left[ M_1(\mathfrak{Y}, S)- M_1(\mathfrak{Y}, S_{r-1})\right] d\mathfrak{Y}\right\| \\ & \leq \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\left\| M_1({\mathfrak{k}}, S) - M_1({\mathfrak{k}}, S_{r-1})\right\| \\ & \quad + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\left\| M_1(\mathfrak{Y}, S)- M_1(\mathfrak{Y}, S_{r-1})\right\| d\mathfrak{Y} \\ & \leq \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\Upsilon_1 \left\| S - S_{r-1}\right\| + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\Upsilon_1 \left\| S - S_{r-1}\right\| {\mathfrak{k}}. \end{align}$

(5.18)

Repeating this process recursively, it becomes

$\begin{align} \left\| \mathfrak{A}_{1, r}({\mathfrak{k}})\right\| \leq \left[ \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)} + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)} {\mathfrak{k}} \right]^{r+1} {\Upsilon}^{r+1}_1 a . \end{align}$

(5.19)

At the point ${\mathfrak{k}}_0$ , we get

(5.20)

As $r$ $\rightarrow \infty$ in Eq (5.20), then we get

$\begin{equation*} \|\mathfrak{A}_{1, r}({\mathfrak{k}})\| \rightarrow 0. \end{equation*}$

Similarly, we get

$\begin{equation*} \|\mathfrak{A}_{2, r}({\mathfrak{k}})\| \rightarrow 0, \; \|\mathfrak{A}_{3, r}({\mathfrak{k}})\| \rightarrow 0, \; \|\mathfrak{A}_{4, r}({\mathfrak{k}})\| \rightarrow 0, \; \|\mathfrak{A}_{5, r}({\mathfrak{k}})\| \rightarrow 0, \; \|\mathfrak{A}_{6, r}({\mathfrak{k}})\| \rightarrow 0, \; \|\mathfrak{A}_{7, r}({\mathfrak{k}})\| \rightarrow 0 . \end{equation*}$

We establish that a system of system (3.2) solutions is unique. Consider the possibility that there is another set of model (3.2) to determine the uniqueness of the solution. We have $S_1({\mathfrak{k}}),$ $E_1({\mathfrak{k}}),$ $I_1({\mathfrak{k}}),$ $Q_1({\mathfrak{k}}),$ $R_1({\mathfrak{k}}),$ $D_1({\mathfrak{k}}),$ and $V_1({\mathfrak{k}});$ then

$\begin{align} S({\mathfrak{k}}) - S_1({\mathfrak{k}}) = &\frac{2(1-{\varsigma})}{(2-{\varsigma})B(\varsigma)}\left[ M_1({\mathfrak{k}}, S) - M_1({\mathfrak{k}}, S_{1})\right] \\ & + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\left[ M_1(\mathfrak{Y}, S) - M_1(\mathfrak{Y}, S_{1})\right] d\mathfrak{Y}. \end{align}$

(5.21)

When we apply the norm to Eq (5.21), we obtain

$\begin{align} \|S({\mathfrak{k}}) - S_1({\mathfrak{k}})\| & \leq \frac{2(1-{\varsigma})}{(2-\varsigma) B(\varsigma)}\left\| \left[ M_1({\mathfrak{k}}, S) - M_1({\mathfrak{k}}, S_{1})\right] \right\| \\ & \quad+ \frac{2\varsigma}{(2- \varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\left\| M_1(\mathfrak{Y}, S) - M_1(\mathfrak{Y}, S_{1})\right\| d\mathfrak{Y}. \end{align}$

(5.22)

Using the kernel's Lipschitz condition, one can achieve

$\begin{align} \|S({\mathfrak{k}}) - S_1({\mathfrak{k}})\| & \leq \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\Upsilon_1 \left\| S({\mathfrak{k}}) - S_{1}({\mathfrak{k}}) \right\| + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\Upsilon_1 {\mathfrak{k}} \left\| S({\mathfrak{k}}) - S_{1}({\mathfrak{k}}) \right\|. \end{align}$

(5.23)

By reducing the complexity of Eq (5.23), we get

$\begin{align} \|S({\mathfrak{k}}) - S_1({\mathfrak{k}})\|\left[ 1 - \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\Upsilon_1 - \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\Upsilon_1 {\mathfrak{k}} \right] \leq 0. \end{align}$

(5.24)

Theorem 5.3. The proposed model has a unique solution if the following conditions are hold:

$\begin{equation} \|S({\mathfrak{k}}) - S_1({\mathfrak{k}})\|\left[ 1 - \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\Upsilon_1 - \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\Upsilon_1 {\mathfrak{k}} \right] > 0. \end{equation}$

(5.25)

Proof. Using Eq (5.24), we have

(5.26)

Implying that

$\begin{align} \left\| S({\mathfrak{k}}) - S_1({\mathfrak{k}}) \right\| = 0, \end{align}$

(5.27)

we obtain

$\begin{align} S({\mathfrak{k}}) = S_1({\mathfrak{k}}). \end{align}$

(5.28)

Continuing in the same manner, we have

$\begin{equation} \begin{array}{ll} E({\mathfrak{k}}) & = E_{1}({\mathfrak{k}}), \; \; I({\mathfrak{k}}) = I_{1}({\mathfrak{k}}), \; \; Q({\mathfrak{k}}) = Q_{1}({\mathfrak{k}}), \\ R({\mathfrak{k}}) & = R_{1}({\mathfrak{k}}), \; \; D({\mathfrak{k}}) = D_{1}({\mathfrak{k}}), \; \; V({\mathfrak{k}}) = V_{1}({\mathfrak{k}}). \end{array} \end{equation}$

(5.29)

As a result, we proved that system (3.2)'s system of solutions is unique.

Stability analysis

We employ nonlinear functional analysis to investigate the Ulam-Hyers (U-H) stability ^[28,37] of the proposed fractional model (3.2). Eighty-four years have passed since Professor Ulam presented the stability problem to the University of Wisconsin Mathematics Club. Ulam asked whether a proposition remains valid or approximately valid when the hypothesis is slightly modified ^[28,38]. This subject is intriguing and significant in multiple scientific disciplines, motivating numerous individuals to pursue its investigation. The current name for this topic is the U-H stability problem. Initially, the focus was on the stability of group homomorphisms. One year later, Hyers resolved the inquiry of the additive mappings across Banach spaces using the "contraction mapping theorem". After this significant advancement, numerous investigations were conducted on Ulam's challenge using diverse approaches and variations. Rassias extended the findings of Hyers ^[39]. Instead of using a positive constant, he employed a dominating function to regulate the estimate. This problem is commonly referred to as the Ulam-Hyers-Rassias stability problem or the generalized U-H stability problem.

In recent decades, numerous research articles have been published on U-H stability, focusing on ordinary differential equations (ODEs) ^[40,41]. ODEs result in productive outcomes, encompassing both linear and nonlinear equations ^[42,43]. Many recent studies have focused on the U-H stability concerning fractional systems ^[44]. The concepts of existence, uniqueness, and U-H stability have been established in various functional systems. U-H stability measures the difference in solutions of differential equations with fractional order. This stability ensures that controlled systems remain stable and perform predictably, even with slight variations in the system dynamics. When modeling physical processes with fractional-order dynamics, U-H stability guarantees that the model stays accurate even when minor changes or approximations are made to the governing equations. U-H stability is used to analyze the sensitivity of solutions to small changes in boundary data for fractional differential equations with boundary conditions. This analysis ensures the reliability of the solutions in practical applications. The concept of U-H stability is beneficial for analyzing and implementing fractional-order systems. It provides a framework to ensure that solutions remain stable even when subjected to minor disturbances, which is essential in numerous practical situations. For simplicity, we consider the proposed model (3.2) to be as follows:

$\begin{align} {}^{CF}{}_{0}D^{\varsigma}_{\mathfrak{k}} \mathfrak{U}({{\mathfrak{k}}}) & = \mathfrak{V}({{{\mathfrak{k}}}}, \mathfrak{U}({{{\mathfrak{k}}}})), \\ \mathfrak{U}({0}) & = \mathfrak{U}_0 \geq 0, \end{align}$

(5.30)

where

$\begin{align*} &\mathfrak{U}({{\mathfrak{k}}}) = \left( S({\mathfrak{k}}), E({\mathfrak{k}}), I({\mathfrak{k}}), Q({\mathfrak{k}}), R({\mathfrak{k}}), D({\mathfrak{k}}), V({\mathfrak{k}})\right)^T, \\ &\mathfrak{V}({{{\mathfrak{k}}}}, \mathfrak{U}({{{\mathfrak{k}}}})) = \left( M_1, M_2, M_3, M_4, M_5, M_6, M_7\right)^T. \end{align*}$

Applying the fractional integral on (5.30), we obtain

$\begin{align} \mathfrak{U}({\mathfrak{k}}) - \mathfrak{U}(0) = \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\mathfrak{V}({{{\mathfrak{k}}}}, \mathfrak{U}({{{\mathfrak{k}}}})) + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\mathfrak{V}(\mathfrak{Y}, \mathfrak{U}(\mathfrak{Y}))d\mathfrak{Y} . \end{align}$

(5.31)

Definition 5.1. The proposed model (3.2) is UH stable if there exists $\varsigma > 0$ and along with the condition for any $\epsilon > 0$ and $\bar{\mathfrak{U}} \in \mathbb{B},$ if

$\begin{align} \left| {}^{CF}{}_{0}D^{\varsigma}_{\mathfrak{k}} \mathfrak{U}({{\mathfrak{k}}}) - \mathfrak{V}({{{\mathfrak{k}}}}, \mathfrak{U}({{{\mathfrak{k}}}}))\right| \leq \epsilon, \end{align}$

(5.32)

then $\exists$ $\mathfrak{U} \in \mathbb{B}$ , and $\mathbb{B}$ in the model (3.2),

$\begin{align*} \mathfrak{U}(0) = \bar{\mathfrak{U}}(0) = \bar{\mathfrak{U}_0}, \end{align*}$

such that

$\begin{align*} \left\| \bar{\mathfrak{U}} - \mathfrak{U} \right\| \leq \varsigma \epsilon, \end{align*}$

where

$\begin{align*} \begin{cases} & \bar{\mathfrak{U}}({{\mathfrak{k}}}) = \left( \bar{S}({\mathfrak{k}}), \bar{E}({\mathfrak{k}}), \bar{I}({\mathfrak{k}}), \bar{Q}({\mathfrak{k}}), \bar{R}({\mathfrak{k}}), \bar{D}({\mathfrak{k}}), \bar{V}({\mathfrak{k}})\right)^T, \\ & \mathfrak{V}({{{\mathfrak{k}}}}, \bar{\mathfrak{U}}({{\mathfrak{k}}})) = \left( \bar{M_1}, \bar{M_2}, \bar{M_3}, \bar{M_4}, \bar{M_5}, \bar{M_6}, \bar{M_7}\right)^T, \\ & \epsilon = \max(\epsilon_i)^T; \; i = 1, 2, ..., 7, \\ & \varsigma = \max(\varsigma_i)^T; \; i = 1, 2, ..., 7. \end{cases} \end{align*}$

Remark 5.1. Let a small perturbation $\kappa \in \mathbb{C}[0, \mathbb{T}],$ such that $\kappa(0) = 0$ with the following condition: $\| \kappa({\mathfrak{k}}) \| \leq \bar{\epsilon},$ for ${\mathfrak{k}} \in [0, \mathbb{T}]$ and $\bar{\epsilon} > 0.$

Lemma 5.1. Let the solution $\bar{\mathfrak{U}}_{\kappa}({{\mathfrak{k}}})$ of the perturbed system

$\begin{align} {}^{CF}{}_{0}D^{\varsigma}_{\mathfrak{k}} \bar{\mathfrak{U}}({{\mathfrak{k}}}) & = \mathfrak{V}({{{\mathfrak{k}}}}, \bar{\mathfrak{U}}({{{\mathfrak{k}}}})) + \kappa({\mathfrak{k}}), \\ \bar{\mathfrak{U}}(0) & = \bar{\mathfrak{U}}_0 \end{align}$

(5.33)

hold the condition

$\begin{align*} \left\| \bar{\mathfrak{U}}_\kappa({{\mathfrak{k}}}) - \bar{\mathfrak{U}}({{\mathfrak{k}}})\right\| \leq \Phi \bar{\epsilon}, \end{align*}$

where

$\begin{align*} \Phi & = \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)} + \frac{2 \varsigma}{(2-\varsigma)B(\varsigma)} \mathbb{T}, \\ \kappa({\mathfrak{k}}) & = \left( \kappa_1({\mathfrak{k}}), \kappa_2({\mathfrak{k}}), \kappa_3({\mathfrak{k}}), \kappa_4({\mathfrak{k}}), \kappa_5({\mathfrak{k}}), \kappa_6({\mathfrak{k}}), \kappa_7({\mathfrak{k}})\right)^T. \end{align*}$

Proof. Applying the fractional integral in Eq (5.33), we obtain

$\begin{align} \bar{\mathfrak{U}}_{\kappa}({\mathfrak{k}}) - \bar{\mathfrak{U}}(0) & = \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\mathfrak{V}({{{\mathfrak{k}}}}, \bar{\mathfrak{U}}({{\mathfrak{k}}})) + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\mathfrak{V}(\mathfrak{Y}, \bar{\mathfrak{U}}(\mathfrak{Y}))d\mathfrak{Y} \\ & \quad+ \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\kappa({\mathfrak{k}}) + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\kappa(\mathfrak{Y})d\mathfrak{Y}. \end{align}$

(5.34)

Also,

$\begin{align} \bar{\mathfrak{U}}({\mathfrak{k}}) & = \bar{\mathfrak{U}}(0) + \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\mathfrak{V}({{{\mathfrak{k}}}}, \bar{\mathfrak{U}}({{\mathfrak{k}}})) + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\mathfrak{V}(\mathfrak{Y}, \bar{\mathfrak{U}}(\mathfrak{Y}))d\mathfrak{Y}. \end{align}$

(5.35)

Using Remark 5.1, we obtain

$\begin{align*} \left\| \bar{\mathfrak{U}}_{\kappa}({\mathfrak{k}}) - \bar{\mathfrak{U}}({\mathfrak{k}})\right\| & \leq \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\left| \kappa({\mathfrak{k}})\right| + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\left| \kappa(\mathfrak{Y})\right| d\mathfrak{Y} \nonumber \\ & \leq \left( \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)} + \frac{2 \varsigma}{(2-\varsigma)B(\varsigma)} \mathbb{T}\right) \bar{\epsilon} \\ & \leq \Phi \bar{\epsilon}. \end{align*}$

This completes the proof.

Theorem 5.4. The proposed model (3.2) is UM stable if

$\begin{align*} \left\| \bar{\mathfrak{U}}({\mathfrak{k}}) - \mathfrak{U}({\mathfrak{k}})\right\| \leq \bar{\varsigma \epsilon}. \end{align*}$

Proof. Let $\bar{\mathfrak{U}}$ be the solution of Eq (5.32) and with help of uniqueness, $\mathfrak{U}$ be a unique solution of the system (5.30), then we have

$\begin{align*} \left\| \bar{\mathfrak{U}}({\mathfrak{k}}) - \mathfrak{U}({\mathfrak{k}})\right\| & \leq \Phi \bar{\epsilon} + \frac{2(1-\varsigma)}{(2-\varsigma)B(\varsigma)}\left\| \mathfrak{V}({{{\mathfrak{k}}}}, \bar{\mathfrak{U}}({{\mathfrak{k}}})) - \mathfrak{V}({{{\mathfrak{k}}}}, {\mathfrak{U}}({{\mathfrak{k}}})) \right\| \\ & \; \; \; + \frac{2\varsigma}{(2-\varsigma)B(\varsigma)}\int_{0}^{{\mathfrak{k}}}\left\| \mathfrak{V}({\mathfrak{Y}}, \bar{\mathfrak{U}}({\mathfrak{Y}})) - \mathfrak{V}({\mathfrak{Y}}, {\mathfrak{U}}({\mathfrak{Y}})) \right\|d\mathfrak{Y} + \Phi \bar{\epsilon} \\ & \leq 2 \Phi \bar{\epsilon} + \Phi \bar{\delta} \left\| \bar{\mathfrak{U}}({\mathfrak{k}}) - \mathfrak{U}({\mathfrak{k}})\right\|. \end{align*}$

Simplifying the above equation, we obtain

$\begin{align*} \left\| \bar{\mathfrak{U}}({\mathfrak{k}}) - \mathfrak{U}({\mathfrak{k}})\right\| & \frac{2 \Phi \bar{\epsilon}}{1 - \Phi \bar{\delta}} = \bar{\varsigma \epsilon}, \end{align*}$

where

$\begin{align*} \varsigma = \frac{2 \Phi}{1 - \Phi \bar{\delta}}. \end{align*}$

Hence, the proposed model (3.2) is UM stable.

6. Numerical scheme

We present a numerical solution for the model (3.2) using the Adams-Bashforth (A-B) technique. Owolabi and Atangana et al. ^[45] introduced a three-step A-B technique with the C-F fractional derivative, which we used to determine the numerical scheme for the fractional-order system (3.2). Take into consideration the C-F derivative of the fractional differential equation

$\begin{equation} \begin{array}{ll} {}^{CF}{}_{0}D^{\varsigma}_{\mathfrak{k}} \mathfrak{E}({{\mathfrak{k}}}) = \mathfrak{F}({{{\mathfrak{k}}}}, \mathfrak{E}({{{\mathfrak{k}}}})). \end{array} \end{equation}$

(6.1)

By utilising the C-F fractional integral on both sides of Eq (6.1), we obtain

$\begin{align} \mathfrak{E}({\mathfrak{k}}) - \mathfrak{E}(0) = \frac{(1-\varsigma)}{B(\varsigma)}{\mathfrak{F}}({{{\mathfrak{k}}}}, \mathfrak{E}({{{\mathfrak{k}}}})) + \frac{\varsigma}{B(\varsigma)}\int_{0}^{{\mathfrak{k}}}{\mathfrak{F}}(\mathfrak{Y}, \mathfrak{E}(\mathfrak{Y}))d\mathfrak{Y}. \end{align}$

(6.2)

We have divided the time interval into smaller intervals with steps of h, then we get ${\mathfrak{k}}_{e+1} = {\mathfrak{k}}_e + h$ , ${\mathfrak{k}}_0 = 0$ , $e = 0, 1, 2, ..., e-1.$ We put ${\mathfrak{k}} = {\mathfrak{k}}_{e+1}$ and ${\mathfrak{k}} = {\mathfrak{k}}_{e}$ in Eq (6.2). We have found the difference in the resulting equations by performing a computation, then we get

$\begin{align} \mathfrak{E}({\mathfrak{k}}_{e+1}) - \mathfrak{E}({\mathfrak{k}}_{e}) & = \frac{(1-\varsigma)}{B(\varsigma)}\left[ {\mathfrak{F}}({{{\mathfrak{k}}}_{e+1}}, \mathfrak{E}({{{\mathfrak{k}}}_{e+1}})) - {\mathfrak{F}}({{{\mathfrak{k}}}_{e}}, \mathfrak{E}({{{\mathfrak{k}}}_{e}}))\right] + \frac{\varsigma}{B(\varsigma)} \int_{{\mathfrak{k}}_e}^{{\mathfrak{k}}_{e+1}}{\mathfrak{F}}(\mathfrak{Y}, \mathfrak{E}(\mathfrak{Y}))d\mathfrak{Y} . \end{align}$

(6.3)

We are estimated the integral $\int_{{\mathfrak{k}}_e}^{{\mathfrak{k}}_{e+1}}{\mathfrak{F}}({\mathfrak{Y}}, \mathfrak{E}({\mathfrak{Y}}))d\mathfrak{Y}$ with the approximation of $\int_{{\mathfrak{k}}_e}^{{\mathfrak{k}}_{e+1}}\mathfrak{Q}_2({\mathfrak{Y}})d\mathfrak{Y},$ where $\mathfrak{Q}_2({\mathfrak{Y}})$ is the Lagrange polynomial of points $({\mathfrak{k}}_{e-2}, {\mathfrak{F}}({{{\mathfrak{k}}}_{e-2}}, \mathfrak{E}({{{\mathfrak{k}}}_{e-2}}))),$ $({\mathfrak{k}}_{e-1}, {\mathfrak{F}}({{{\mathfrak{k}}}_{e-1}}, \mathfrak{E}({{{\mathfrak{k}}}_{e-1}})))$ and $({\mathfrak{k}}_{e}, {\mathfrak{F}}({{{\mathfrak{k}}}_{e}}, \mathfrak{E}({{{\mathfrak{k}}}_{e}})))$ . Thus,

$\begin{align} \mathfrak{Q}_2(\mathfrak{Y}) = \sum\limits_{j = 0}^{j = 2}{\mathfrak{F}}({{{\mathfrak{k}}}_{e-j}}, \mathfrak{E}({{{\mathfrak{k}}}_{e-j}}))\mathfrak{L}_{j}(\mathfrak{Y}), \end{align}$

(6.4)

where $\mathfrak{L}_{j}(\mathfrak{Y})$ is the Lagrange basis polynomials on the point ${\mathfrak{k}}_e, {\mathfrak{k}}_{e-1}, {\mathfrak{k}}_{e-2}.$

Let $\mathfrak{E}_{e} = \mathfrak{E}({{{\mathfrak{k}}}_{e}}),$ and $v = \frac{{\mathfrak{k}}_{e+1}-\mathfrak{Y}}{h},$ after replacing the Lagrange basis polynomials and performing integration, we obtain the following result:

$\begin{align} \int_{{\mathfrak{k}}_e}^{{\mathfrak{k}}_{e+1}}{\mathfrak{F}}(\mathfrak{Y}, \mathfrak{E}(\mathfrak{Y}))d\mathfrak{Y} & = h \int_{0}^{1}\left( {\mathfrak{F}}({\mathfrak{k}}_e, \mathfrak{E}_{e})\frac{(v-2)(v-3)}{(1-2)(1-3)}\right. \\ & \quad\left. + {\mathfrak{F}}({\mathfrak{k}}_{e-1}, \mathfrak{E}_{e-1})\frac{(v-1)(v-3)}{(2-1)(2-3)} + {\mathfrak{F}}({\mathfrak{k}}_{e-2}, \mathfrak{E}_{e-2})\frac{(v-2)(v-1)}{(3-2)(3-1)}\right) dv \\ & = \frac{23h}{12}{\mathfrak{F}}({\mathfrak{k}}_{e}, \mathfrak{E}_{e}) - \frac{16h}{12}{\mathfrak{F}}({\mathfrak{k}}_{e-1}, \mathfrak{E}_{e-1}) + \frac{5h}{12}{\mathfrak{F}}({\mathfrak{k}}_{e-2}, \mathfrak{E}_{e-2}). \end{align}$

(6.5)

Using Eq (6.5) in Eq (6.3), we have

$\begin{align} \mathfrak{E}({\mathfrak{k}}_{e+1}) - \mathfrak{E}({\mathfrak{k}}_{e}) = & \left( \frac{(1-\varsigma)}{B(\varsigma)} + \frac{23h\varsigma}{12 B(\varsigma)}\right) {\mathfrak{F}}({{{\mathfrak{k}}}_{e}}, \mathfrak{E}({{{\mathfrak{k}}}_{e}})) \\ & - \left( \frac{(1-\varsigma)}{B(\varsigma)} + \frac{16h\varsigma}{12 B(\varsigma)}\right) {\mathfrak{F}}({{{\mathfrak{k}}}_{e-1}}, \mathfrak{E}({{{\mathfrak{k}}}_{e-1}})) + \left( \frac{5h\varsigma}{12 B(\varsigma)}\right) {\mathfrak{F}}({{{\mathfrak{k}}}_{e-2}}, \mathfrak{E}({{{\mathfrak{k}}}_{e-2}})). \end{align}$

(6.6)

The error with this technique is

$\begin{align} \mathfrak{R}^{\varsigma}_{e}({\mathfrak{k}}) = \frac{\varsigma}{B(\varsigma)}\int_{{\mathfrak{k}}_e}^{{\mathfrak{k}}_{e+1}}\frac{3h^3}{8}\mathfrak{F}^4(s)ds = \frac{3\varsigma h^3}{8B(\varsigma)}{\mathfrak{F}^3}(x_e, {\mathfrak{E}(x_e)}), x_e \in ({\mathfrak{k}}_e, {\mathfrak{k}}_{e+1}). \end{align}$

(6.7)

The proposed model's numerical simulations are generated using the three-step A-B approach for the C-F fractional derivative in Eq (6.6).

Let the vectors

$\mathfrak{E}({\mathfrak{k}}) = \left[ S({\mathfrak{k}}), E({\mathfrak{k}}), I({\mathfrak{k}}) Q({\mathfrak{k}}), R({\mathfrak{k}}), D({\mathfrak{k}}), V({\mathfrak{k}})\right],$

and

$\mathfrak{F}({\mathfrak{k}}, \mathfrak{E}({\mathfrak{k}})) = \left[ \mathfrak{F}_1({\mathfrak{k}}, \mathfrak{E}({\mathfrak{k}})), \mathfrak{F}_2({\mathfrak{k}}, \mathfrak{E}({\mathfrak{k}})), \mathfrak{F}_3({\mathfrak{k}}, \mathfrak{E}({\mathfrak{k}})), \mathfrak{F}_4({\mathfrak{k}}, \mathfrak{E}({\mathfrak{k}})), \mathfrak{F}_5({\mathfrak{k}}, \mathfrak{E}({\mathfrak{k}})), \mathfrak{F}_6({\mathfrak{k}}, \mathfrak{E}({\mathfrak{k}})), \mathfrak{F}_7({\mathfrak{k}}, \mathfrak{E}({\mathfrak{k}}))\right],$

where the following are the functions specified by system (3.2), we have

$\begin{align} \begin{cases} \mathfrak{F}_1({\mathfrak{k}}, \mathfrak{E}({\mathfrak{k}})) & = \; \Pi - \beta SI - \upsilon S - \mu S, \\ \mathfrak{F}_2({\mathfrak{k}}, \mathfrak{E}({\mathfrak{k}})) & = \; \beta SI - \gamma E + \sigma \beta VI - \mu E, \\ \mathfrak{F}_3({\mathfrak{k}}, \mathfrak{E}({\mathfrak{k}})) & = \; \gamma E - \delta I - \mu I, \\ \mathfrak{F}_4({\mathfrak{k}}, \mathfrak{E}({\mathfrak{k}})) & = \; \delta I - (1-\tau)\omega Q - \tau \rho Q - \mu Q, \\ \mathfrak{F}_5({\mathfrak{k}}, \mathfrak{E}({\mathfrak{k}})) & = \; (1-\tau)\omega Q -\mu R, \\ \mathfrak{F}_6({\mathfrak{k}}, \mathfrak{E}({\mathfrak{k}})) & = \; \tau \rho Q, \\ \mathfrak{F}_7({\mathfrak{k}}, \mathfrak{E}({\mathfrak{k}})) & = \; \upsilon S - \sigma \beta VI - \mu V. \end{cases} \end{align}$

(6.8)

We can express model (3.2) in vector form as shown below:

(6.9)

Equation (6.6) is utilized to obtain the solution of model (3.2), which is expressed through the iterative formula given below:

$\begin{align} \mathfrak{E}({\mathfrak{k}}_{e+1}) = & \mathfrak{E}({\mathfrak{k}}_{e}) + \left( \frac{(1-\varsigma)}{B(\varsigma)} + \frac{23h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}({{{\mathfrak{k}}}_{e}}, \mathfrak{E}_{e}) \\ & - \left( \frac{(1-\varsigma)}{B(\varsigma)} + \frac{16h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}({{{\mathfrak{k}}}_{e-1}}, \mathfrak{E}_{e-1}) + \left( \frac{5h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}({{{\mathfrak{k}}}_{e-2}}, \mathfrak{E}_{e-2}). \end{align}$

(6.10)

Let $\mathfrak{E}_0 = \mathfrak{E}({\mathfrak{k}}_0) = \left[S({\mathfrak{k}}_0), E({\mathfrak{k}}_0), I({\mathfrak{k}}_0), Q({\mathfrak{k}}_0), R({\mathfrak{k}}_0), D({\mathfrak{k}}_0), V({\mathfrak{k}}_0) \right]^{T},$ $\mathfrak{E}_{e-2} = \mathfrak{E}({\mathfrak{k}}_{e-2}),$ $\mathfrak{E}_{e-1} = \mathfrak{E}({\mathfrak{k}}_{e-1}),$ $\mathfrak{E}_{e} = \mathfrak{E}({\mathfrak{k}}_{e})$ and $\mathfrak{E}_{e+1} = \mathfrak{E}({\mathfrak{k}}_{e+1}),$ then

$\begin{align} \mathfrak{E}_{e+1} = & \mathfrak{E}_{e} + \left( \frac{(1-\varsigma)}{B(\varsigma)} + \frac{23h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}({{{\mathfrak{k}}}_{e}}, \mathfrak{E}_{e}) \\ & - \left( \frac{(1-\varsigma)}{B(\varsigma)} + \frac{16h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}({{{\mathfrak{k}}}_{e-1}}, \mathfrak{E}_{e-1}) + \left( \frac{5h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}({{{\mathfrak{k}}}_{e-2}}, \mathfrak{E}_{e-2}). \end{align}$

(6.11)

We have obtained the iterative formula:

$\begin{align*} S_{e+1} = & S_{e} + \left( \frac{(1-\varsigma)}{B(\varsigma)} + \frac{23h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}_1(S_{e}, E_{e}, I_{e}, Q_{e}, R_{e}, D_{e}, V_{e}) \nonumber \\ & - \left( \frac{(1-\varsigma)}{B(\varsigma)} + \frac{16h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}_1(S_{e-1}, E_{e-1}, I_{e-1}, Q_{e-1}, R_{e-1}, D_{e-1}, V_{e-1}) \nonumber \\ & + \left( \frac{5h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}_1(S_{e-2}, E_{e-2}, I_{e-2}, Q_{e-2}, R_{e-2}, D_{e-2}, V_{e-2}), \end{align*}$

$\begin{align*} E_{e+1} = & E_{e} + \left( \frac{(1-\varsigma)}{B(\varsigma)} + \frac{23h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}_2(S_{e}, E_{e}, I_{e}, Q_{e}, R_{e}, D_{e}, V_{e}) \nonumber \\ & - \left( \frac{(1-\varsigma)}{B(\varsigma)} + \frac{16h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}_2(S_{e-1}, E_{e-1}, I_{e-1}, Q_{e-1}, R_{e-1}, D_{e-1}, V_{e-1}) \nonumber \\ & + \left( \frac{5h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}_2(S_{e-2}, E_{e-2}, I_{e-2}, Q_{e-2}, R_{e-2}, D_{e-2}, V_{e-2}), \end{align*}$

$\begin{align*} I_{e+1} = & I_{e} + \left( \frac{(1-\varsigma)}{B(\varsigma)} + \frac{23h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}_3(S_{e}, E_{e}, I_{e}, Q_{e}, R_{e}, D_{e}, V_{e}) \nonumber \\ & - \left( \frac{(1-\varsigma)}{B(\varsigma)} + \frac{16h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}_3(S_{e-1}, E_{e-1}, I_{e-1}, Q_{e-1}, R_{e-1}, D_{e-1}, V_{e-1}) \nonumber \\ & + \left( \frac{5h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}_3(S_{e-2}, E_{e-2}, I_{e-2}, Q_{e-2}, R_{e-2}, D_{e-2}, V_{e-2}), \end{align*}$

$\begin{align} Q_{e+1} = & Q_{e} + \left( \frac{(1-\varsigma)}{B(\varsigma)} + \frac{23h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}_4(S_{e}, E_{e}, I_{e}, Q_{e}, R_{e}, D_{e}, V_{e}) \\ & - \left( \frac{(1-\varsigma)}{B(\varsigma)} + \frac{16h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}_4(S_{e-1}, E_{e-1}, I_{e-1}, Q_{e-1}, R_{e-1}, D_{e-1}, V_{e-1}) \\ & + \left( \frac{5h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}_4(S_{e-2}, E_{e-2}, I_{e-2}, Q_{e-2}, R_{e-2}, D_{e-2}, V_{e-2}), \end{align}$

(6.12)

$\begin{align*} R_{e+1} = & R_{e} + \left( \frac{(1-\varsigma)}{B(\varsigma)} + \frac{23h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}_5(S_{e}, E_{e}, I_{e}, Q_{e}, R_{e}, D_{e}, V_{e}) \nonumber \\ & - \left( \frac{(1-\varsigma)}{B(\varsigma)} + \frac{16h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}_5(S_{e-1}, E_{e-1}, I_{e-1}, Q_{e-1}, R_{e-1}, D_{e-1}, V_{e-1}) \nonumber \\ & + \left( \frac{5h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}_5(S_{e-2}, E_{e-2}, I_{e-2}, Q_{e-2}, R_{e-2}, D_{e-2}, V_{e-2}), \end{align*}$

$\begin{align*} D_{e+1} = & D_{e} + \left( \frac{(1-\varsigma)}{B(\varsigma)} + \frac{23h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}_6(S_{e}, E_{e}, I_{e}, Q_{e}, R_{e}, D_{e}, V_{e}) \nonumber \\ & - \left( \frac{(1-\varsigma)}{B(\varsigma)} + \frac{16h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}_6(S_{e-1}, E_{e-1}, I_{e-1}, Q_{e-1}, R_{e-1}, D_{e-1}, V_{e-1}) \nonumber \\ & + \left( \frac{5h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}_6(S_{e-2}, E_{e-2}, I_{e-2}, Q_{e-2}, R_{e-2}, D_{e-2}, V_{e-2}), \end{align*}$

$\begin{align*} V_{e+1} = & V_{e} + \left( \frac{(1-\varsigma)}{B(\varsigma)} + \frac{23h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}_7(S_{e}, E_{e}, I_{e}, Q_{e}, R_{e}, D_{e}, V_{e}) \nonumber \\ & - \left( \frac{(1-\varsigma)}{B(\varsigma)} + \frac{16h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}_7(S_{e-1}, E_{e-1}, I_{e-1}, Q_{e-1}, R_{e-1}, D_{e-1}, V_{e-1}) \nonumber \\ & + \left( \frac{5h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}_7(S_{e-2}, E_{e-2}, I_{e-2}, Q_{e-2}, R_{e-2}, D_{e-2}, V_{e-2}). \end{align*}$

Stability analysis of numerical scheme

Here, we are analysing the stability of the fractional A-B scheme

(6.13)

Now, rewrite Eq (6.11), then we have

$\begin{align} \mathfrak{E}_{e+1} & = \mathfrak{E}_{e} + \left( \frac{(1-\varsigma)}{B(\varsigma)} + \frac{23h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}({{{\mathfrak{k}}}_{e}}, \mathfrak{E}_{e}) \\ &\quad - \left( \frac{(1-\varsigma)}{B(\varsigma)} + \frac{16h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}({{{\mathfrak{k}}}_{e-1}}, \mathfrak{E}_{e-1}) + \left( \frac{5h\varsigma}{12 B(\varsigma)}\right) \mathfrak{F}({{{\mathfrak{k}}}_{e-2}}, \mathfrak{E}_{e-2}) \\ & = \mathfrak{E}_{e} + A \mathfrak{F}({{{\mathfrak{k}}}_{e}}, \mathfrak{E}_{e}) - B \mathfrak{F}({{{\mathfrak{k}}}_{e-1}}, \mathfrak{E}_{e-1}) + C \mathfrak{F}({{{\mathfrak{k}}}_{e-2}}, \mathfrak{E}_{e-2}). \end{align}$

(6.14)

Using Eq (5.30), we have

$\begin{align} \mathfrak{E}_{e+1} = (1+A) \mathfrak{E}_{e} - B \mathfrak{E}_{e-1} + C \mathfrak{E}_{e-2}. \end{align}$

(6.15)

Afterward, we apply the Von-Neumann stability analysis for the terms in the aforementioned equation.

$\begin{align} \mathfrak{E}_{e+1} & = \bar{\mathfrak{E}_{e+1}}\mathfrak{e}^{(e+1)i \Delta \mathfrak{k}}, \; \; \mathfrak{E}_{e} = \bar{\mathfrak{E}_{e}}\mathfrak{e}^{(e)i \Delta \mathfrak{k}}, \\ \mathfrak{E}_{e-1} & = \bar{\mathfrak{E}_{e-1}}\mathfrak{e}^{(e-1)i \Delta \mathfrak{k}}, \; \; \mathfrak{E}_{e-2} = \bar{\mathfrak{E}_{e-2}}\mathfrak{e}^{(e-2)i \Delta \mathfrak{k}}. \end{align}$

(6.16)

So that

$\begin{align} \bar{\mathfrak{E}_{e+1}}\mathfrak{e}^{(e+1)i \Delta \mathfrak{k}} = (1 + A)\bar{\mathfrak{E}_{e}}\mathfrak{e}^{e i \Delta \mathfrak{k}} - B \bar{\mathfrak{E}_{e-1}}\mathfrak{e}^{(e-1)i \Delta \mathfrak{k}} + C \bar{\mathfrak{E}_{e-2}}\mathfrak{e}^{(e-2)i \Delta \mathfrak{k}}, \end{align}$

(6.17)

which reduces to

$\begin{align} \bar{\mathfrak{E}_{e+1}}\mathfrak{e}^{i \Delta \mathfrak{k}} = (1 + A)\bar{\mathfrak{E}_{e}} - B \bar{\mathfrak{E}_{e-1}}\mathfrak{e}^{- i \Delta \mathfrak{k}} + C \bar{\mathfrak{E}_{e-2}}\mathfrak{e}^{-2i \Delta \mathfrak{k}}. \end{align}$

(6.18)

A recursive formula is applied, when $e = 0,$ we have

$\begin{align} \bar{\mathfrak{E}_{1}}\mathfrak{e}^{i \Delta \mathfrak{k}} = (1 + A)\bar{\mathfrak{E}_{0}}. \end{align}$

(6.19)

Simplifying the previous equation (6.15), we have

$\begin{align*} \left\| \bar{\mathfrak{E}_{1}}\mathfrak{e}^{i \Delta \mathfrak{k}}\right\| & = \left\| (1 + A)\bar{\mathfrak{E}_{0}} \right\|, \\ \implies \left\| \bar{\mathfrak{E}_{1}}\right\| & = (1 + A) \left\| \bar{\mathfrak{E}_{0}} \right\|. \end{align*}$

$\begin{align*} h \leq 0.55 \; \; \text{and} \; \; \varsigma > \frac{1}{1 - 1.9167h}, \end{align*}$

then, we suppose that

$\begin{align*} \forall e \geq \bar{\mathfrak{E}_{e}} < \bar{\mathfrak{E}_{0}}, \end{align*}$

we have

$\begin{align} \left| \bar{\mathfrak{E}_{e+1}}\right| & \leq (1 + A)\left| \bar{\mathfrak{E}_{e}}\right|\left| \mathfrak{e}^{- i \Delta \mathfrak{k}}\right| - B \left| \bar{\mathfrak{E}_{e-1}}\right| \left| \mathfrak{e}^{- 2i \Delta \mathfrak{k}}\right| + C \left| \bar{\mathfrak{E}_{e-2}}\right| \left| \mathfrak{e}^{-3i \Delta \mathfrak{k}}\right| \\ & < (1 + A)\left| \bar{\mathfrak{E}_{e}}\right| - B \left| \bar{\mathfrak{E}_{e-1}}\right| + C \left| \bar{\mathfrak{E}_{e-2}}\right|. \end{align}$

(6.20)

Simplifying the previous equation (6.15), we have

$\begin{align} \left| \bar{\mathfrak{E}_{e+1}}\right| & < (1 + A)\left| \bar{\mathfrak{E}_{0}}\right| - B \left| \bar{\mathfrak{E}_{0}}\right| + C \left| \bar{\mathfrak{E}_{0}}\right| \\ & < (1 + A - B + C)\left| \bar{\mathfrak{E}_{0}}\right| \\ & < (1 + A - B + C) < 1. \end{align}$

(6.21)

Hence, when applied to Eq (5.30), the three-step A-B method with the C-F derivative is conditionally stable.

7. Numerical results and graphical analysis

In this section, we simulate the non-linear fractional proposed model using a novel numerical technique that uses the C-F non-integer operator. We displayed graphical findings for various fractional orders and infection rates during the simulation. To comprehend the system's behavior and operation, tracking and controlling the infection rate is imperative. Utilizing the numerical approach outlined in Section 6, we followed the instructions and finished the simulation results. We examined numerical findings for various fractional orders and time intervals using MATLAB R2016a. As such, we can see how adjusting the parameters and starting conditions affects the model's predictions through the simulations and improves our grasp of the dynamics of the model. The fractional order analysis's graphical results have been more illuminating and broadly applicable than other related efforts. The following parameter values have been applied to the simulation: $S(0) = 100,$ $E(0) = 50,$ $I(0) = 30,$ $Q(0) = 20,$ $R(0) = 18,$ $D(0) = 2,$ $V(0) = 60,$ $\Pi = 20,$ $\beta = 2 \times 10^{-4},$ $\nu = 1 \times 10^{-3},$ $\mu = 0.3,$ $\gamma = 5.1,$ $\sigma = 0.05,$ $\delta = 5.1,$ $\tau = 0.04,$ $\omega = 0.03$ and $\rho = 15.$ For this study, we have chosen specific orders to represent the complex nature of the system under examination, as arbitrary order derivatives offer greater degrees of freedom and provide a wide range of geometries. depicts the state variable plots of time graphically in the suggested model using the ABM method for different orders $(\varsigma)$ : $0.99, 0.96, 0.93$ and $0.90$ . The fractional suggested model (3.2) for various values of $\varsigma$ has been numerically simulated and shown in Figures 5 and 6.

Figure 4. Solution trajectories for the proposed model (3.2) when the derivative order varies.

DownLoad: Full-Size Img PowerPoint

Figure 5. The plots for the proposed model (3.2) when varies values of

$\varsigma$ .

DownLoad: Full-Size Img PowerPoint

Figure 6. The plots for the proposed model (3.2) when varies values of

$\varsigma$ .

DownLoad: Full-Size Img PowerPoint

When we place a classical order, the infected population is higher, but when we place an order with $\varsigma = 0.85,$ the infected population is smaller, as illustrated in . When placing a classical order, the recovered and quarantined population is smaller. However, when placing an order with $\varsigma = 0.85,$ the recovered and quarantined population is more significant.

We consider a small order, leading to a more significant impact on the dead population than a more substantial higher order. depicts the 3D phase trajectory of the system (3.2) as the order of derivative = 1 and $\mu$ varies. depicts the time series graph of the system (3.2) when the order of derivative $(\varsigma) = 0.99$ and $\mu$ varies. When we place the value of $\mu = 0.45$ , the susceptible, quarantined, and vaccinated population is smaller. However, when we set $\mu = 0.30$ , the susceptible, quarantined, and vaccinated population is more significant, as shown in . If we choose the smaller value for $\mu$ , we have large populations of these populations. The dead population becomes more extensive if we choose the smaller value for $\mu$ . depicts the time series graph of the system (3.2) when the order of derivative $(\varsigma) = 0.99$ , $\delta$ varies and $\nu$ varies. The infected and quarantined populations are sensitive to parameter $\delta$ increasing order. The exposed and infected populations are sensitive to parameter $\nu$ . depicts the 3D phase trajectory of the system (3.2) as the order of derivative = 1 and $\beta$ varies. depicts the time series graph of the system (3.2) when the order of derivative $(\varsigma) = 0.99$ and $\beta$ varies. When we place the value of $\beta$ is small, we get the susceptible, infected, and recovered population is more significant compared to the more considerable value of $\beta$ . The dead population becomes more extensive if we choose the smaller value for $\beta$ . depicts the 3D phase trajectory of the system (3.2) as the order of derivative = 1 and $\tau$ varies. depicts the time series graph of the system (3.2) when the order of derivative $(\varsigma) = 0.99$ and $\tau$ varies. The exposed and recovered populations are not very sensitive to parameter $\tau$ . We consider a small value for $\tau$ , leading to a more significant impact on the quarantined population than a more substantial value of $\tau$ .

Figure 7. Three-dimensional numerical results for proposed model (3.2) when

$\mu$ varies and

$\varsigma = 1$ .

DownLoad: Full-Size Img PowerPoint

Figure 8. The plots for the proposed model (3.2) when

$\mu$ varies and

$\varsigma = 0.99$ .

DownLoad: Full-Size Img PowerPoint

Figure 9. The plots for the proposed model (3.2) when

$\delta$ and

$\nu$ are varies and

$\varsigma = 0.99$ .

DownLoad: Full-Size Img PowerPoint

Figure 10. Three-dimensional numerical results for proposed model (3.2) when

$\beta$ varies and

$\varsigma = 1$ .

DownLoad: Full-Size Img PowerPoint

Figure 11. The plots for the proposed model (3.2) when

$\beta$ varies and

$\varsigma = 0.99$ .

DownLoad: Full-Size Img PowerPoint

Figure 12. Three-dimensional numerical results for proposed model (3.2) when

$\tau$ varies and

$\varsigma = 1$ .

DownLoad: Full-Size Img PowerPoint

Figure 13. The plots for the proposed model (3.2) when

$\tau$ varies and

$\varsigma = 0.99$ .

DownLoad: Full-Size Img PowerPoint

The numerical simulation suggests we gain more reliable insights into behavior within specific ranges with more precise parameter information. We can explore the model's behavior under different input parameter values. Therefore, this study highlights the advantages of using numerical techniques in real-world disease models.

8. Conclusions

In this paper, we discuss the SEIQRDV model, which focuses on the spread of infectious diseases and considers positivity, boundedness, and sensitivity. We have created a fractional order model based on a described integer order contagious disease model. Using fixed-point theory, we have proven the uniqueness and existence of solutions to this model under specific circumstances. We have also used a novel numerical technique supported by the A-B scheme to find the numerical solution. By plotting the model's output at different fractional orders and parameter values, we have observed how the population curves respond to these variations. However, it is important to note that the uncertainty resulting from random fluctuations in the model parameters cannot be eliminated. Additionally, determining which parameters are essential in the model can be challenging, causing sensitivity analysis for such a model to be difficult. The sensitivity analysis of the basic reproduction number has shown that the model's parameters have a significant impact. The results indicate that the recruitment and transmission rates are the most influential factors contributing to a substantial increase in the $R_0$ . These findings provide valuable insights into the mechanisms that affect the dynamics of infectious disease transmission and offer strategies for reducing disease transmission.

We observe that even slight changes in parameters such as $\beta$ , $\delta$ , $\mu$ , $\nu$ , and $\tau$ result in different population dynamics. These efforts will enhance our understanding of this paradigm and enable us to obtain more information successfully. Moreover, the results indicate that as the number of infectious individuals in the community decreases, the number of recovered individuals in the system increases. This suggests a correlation between the decrease in the number of sick individuals and the increase in the number of those who have effectively recuperated from the illness. The proposed model provides valuable insights into disease transmission dynamics using fractional-order derivatives. The results of our investigation are crucial for improving the accuracy of the infectious disease model and formulating successful strategies. Furthermore, our investigation indicated that meaningful parameters can be found and examined even in uncertainty. These research results suggest expanding this work to encompass control theory and other fractional operators.

Author contributions

Parveen Kumar: Formal analysis, Writing–original draft, Methodology, Software, Validation; Sunil Kumar: Supervision, Conceptualization, Investigation, Writing–review and editing; Badr Saad T Alkahtani: Validation, Writing–review and editing, Conceptualization, Visualization; Sara S Alzaid: Investigation, Formal analysis, Visualization, Writing–original draft. All authors have read and agreed to the published version of the manuscript.

Acknowledgement

This research was funded by Researchers Supporting Project number (RSPD2024R526), King Saud University, Riyadh, Saudi Arabia.

Funding

This research was funded by Researchers Supporting Project number (RSPD2024R526), King Saud University, Riyadh, Saudi Arabia.

Conflict of interest

All authors declare no conflicts of interest.

References

[1]	J. Wang, X. Xun, D. O. Ophthalmology, S. G. Hospital, Environmental risk factors for myopia and advances in its control and prevention, Shanghai Med. Pharm. J., 38 (2017), 3-7.
[2]	M. Kohlhaas, Corneal sensation after cataract and refractive surgery, J. Cataract Refract. Surg., 24 (1998), 1399-1409. doi: 10.1016/S0886-3350(98)80237-X
[3]	D. Z. Reinstein, T. J. Archer, M. Gobbe, The Key characteristics of corneal refractive surgery: biomechanics, spherical aberration, and corneal sensitivity after SMILE, in Small Incision Lenticule Extraction (SMILE), Springer, Cham, (2015), 123-142.
[4]	A. Gyldenkerne, A. Ivarsen, J. Hjortdal, Comparison of corneal shape changes and aberrations induced By FS-LASIK and SMILE for myopia, J. Refractive Surg., 31 (2015), 223-229. doi: 10.3928/1081597X-20150303-01
[5]	P. Wei, G. P. Cheng, J. Zhang, A.L. Ng, Y. Wang, Changes in corneal volume at different areas and its correlation with corneal biomechanics after SMILE and FS-LASIK surgery, J. Ophthalmol., 2020 (2020), 1-7.
[6]	G. Scarcelli, S. H. Yun, Brillouin scanning microscopy in keratoconus, in Keratoconus: Recent Advances in Diagnosis and Treatment, Springer International Publishing, (2017), 167-173.
[7]	N. Cartwright, J. R. Tyrer, P. D. Jaycock, J. Marshall, Effects of variation in depth and side cut angulations in LASIK and thin-flap LASIK using a femtosecond laser: a biomechanical study, J. Refractive Surg., 28 (2012), 419-425. doi: 10.3928/1081597X-20120518-07
[8]	D. V. Franus, Change in the stress-strain state of the cornea after refractive surgery, in 2015 International Conference on Mechanics-Seventh Polyakhov's Reading, IEEE, (2015), 1-4.
[9]	J. C. He, J. Gwiazda, F. Thorn, R. Held, Wave-front aberrations in the anterior corneal surface and the whole eye, J. Opt. Soc. Am. A, 20 (2003), 1155-1163. doi: 10.1364/JOSAA.20.001155
[10]	A. S. Roy, W. J. Dupps, Effects of altered corneal stiffness on native and postoperative LASIK corneal biomechanical behavior: a whole-eye finite element analysis, J. Refractive Surge., 25 (2009), 875-887. doi: 10.3928/1081597X-20090917-09
[11]	K. J. Bathe, D. Chapelle, Computational Fluid And Solid Mechanics: the Finite Element Analysis of Shells- Fundamentals, Springer, 2003.
[12]	M. Á. Ariza-Gracia, J. F. Zurita, D. P. Pinero, J. F. Rodriguez, Coupled biomechanical response of the cornea assessed by non-contact tonometry, a simulation study, PLOS ONE, 10 (2015), e0121486.
[13]	M. Kaliske, A formulation of elasticity and viscoelasticity for fibre reinforced material at small and finite strains, Comput. Methods Appl. Mech. Eng., 185 (2000), 225-243. doi: 10.1016/S0045-7825(99)00261-3
[14]	K. Anderson, A. El-Sheikh, T. Newson, Application of structural analysis to the mechanical behaviour of the cornea, J. R. Soc. Interface, 1 (2004), 1742-5662.
[15]	P. M. Pinsky, D. V. Datye, A microstructurally-based finite element model of the incised human cornea, J. Biomech., 24 (1991), 907-922. doi: 10.1016/0021-9290(91)90169-N
[16]	R. Shah, S. Shah, S. Sengupta, Results of small incision lenticule extraction: all-in-one femtosecond laser refractive surgery, J. Cataract Refractive Surg., 37 (2011), 127-137. doi: 10.1016/j.jcrs.2010.07.033
[17]	J. B. Randleman, B. Russell, M. A. Ward, K. P. Thompson, R. D. Stulting, Risk factors and prognosis for corneal ectasia after LASIK, Ophthalmology, 110 (2003), 267-275. doi: 10.1016/S0161-6420(02)01727-X
[18]	K. A. John, Comparison of corneal biomechanics after myopic small-incision lenticule extraction compared to LASIK: an ex vivo study, Clin. Ophthalmol., 12 (2018), 237-245. doi: 10.2147/OPTH.S153509
[19]	I. M. Osman, H. A. Helaly, M. Abdalla, M. A. Shousha, Corneal biomechanical changes in eyes with small incision lenticule extraction and laser assisted in situ keratomileusis, BMC Ophthalmol., 16 (2016), 123. doi: 10.1186/s12886-016-0304-3
[20]	T. F. Peinado, D. P. Piñ ero, I. A. López, J. L. Alio, Correlation of both corneal surfaces in corneal ectasia after myopic LASIK, Optom. Vision Sci., 88 (2011), E539-E542.
[21]	R. Grytz, K. Krishnan, R. Whitley, V. Libertiaux, J. C. Downs, A mesh-free approach to incorporate complex anisotropic and heterogeneous material properties into eye-specific finite element models, Comput. Methods Applied Mech. Eng., 358 (2020), 112654. doi: 10.1016/j.cma.2019.112654
[22]	J. G. Yu, F. J. Bao, Y. F. Feng, C. Whitford, A. Elsheikh, Assessment of corneal biomechanical behavior under posterior and anterior pressure, J. Refractive Surg., 29 (2013), 64-70. doi: 10.3928/1081597X-20121228-05
[23]	L. Fang, The influence of the aspheric profiles for transition zone on optical performance of human eye after conventional ablation, J. Eur. Opt. Soc. Rapid Publ., 9 (1990).
[24]	C. J. Roberts, W. J. Dupps, Biomechanics of corneal ectasia and biomechanical treatments, J. Cataract Refractive Surg., 40 (2014), 991-998. doi: 10.1016/j.jcrs.2014.04.013
[25]	M. Girard, W. J. Dupps, M. Baskaran, G. Scarcelli, S. H. Yun, H. A. Quigley, et al., Translating ocular biomechanics into clinical practice: current state and future prospects, Curr. Eye Res., 40 (2015), 1-18. doi: 10.3109/02713683.2014.914543
[26]	S. Ganesh, S. Brar, R. R. Arra, Refractive lenticule extraction small incision lenticule extraction: a new refractive surgery paradigm, Indian J. Ophthalmol., 66 (2018), 10-19. doi: 10.4103/ijo.IJO_761_17
[27]	W. Di, Y. Wang, L. Zhang, S. Wei, T. Xin, Corneal biomechanical effects: Small-incision lenticule extraction versus femtosecond laser-assisted laser in situ keratomileusis, J. Cataract Refractive Surg., 40 (2014), 954-962. doi: 10.1016/j.jcrs.2013.07.056
[28]	M. Balidis, Biomechanical profile of refractive surgery procedures, Acta Ophthalmol., 97 (2019).
[29]	P. J. Shih, I. J. Wang, W. F. Cai, J. Y. Yen, Biomechanical simulation of stress concentration and intraocular pressure in corneas subjected to myopic refractive surgical procedures, Sci. Rep., 7 (2017), 13906. doi: 10.1038/s41598-017-14293-0
[30]	M. I. Cordero-Mendieta, E. Pinos-Vélez, R. Coronel-Berrezueta, Study of corneal biomechanics and modeling of Young's module in healthy and pathological corneas, in International Conference on Applied Human Factors and Ergonomics, Springer, Cham, (2020), 99-105.
[31]	N. T. Mohammad, F. Craig, G. Dipika, Finite element modelling of cornea mechanics: a review, Arq. Bras. Oftalmol., 77 (2014), 60-65.

Reader Comments

Your name:*

Email:*
© 2021 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)