On Theoretical and Numerical Results of Serum Hepatitis Disease Using Piecewise Fractal–Fractional Perspectives

Abstract

In the present research, we consider a biological model of serum hepatitis disease. We carry out a detailed analysis of the mentioned model along with a class with asymptomatic carriers to explore its theoretical and numerical aspects. We initiate the study by using the piecewise fractal–fractional derivative (FFD) by which the crossover effects within the model are examined. We split the time interval into subintervals. In one subinterval, FFD with a power law kernel is taken, while in the second one, FFD with an exponential decay kernel of the proposed model is considered. This model is then studied for its disease-free equilibrium point, existence, and Hyers–Ulam (H-U) stability results. For numerical results of the model and a visual presentation, we apply the Lagrange interpolation method and an extended Adams–Bashforth–Moulton (ABM) method, respectively.

1. Introduction

Serum hepatitis, also known as Hepatitis B, is a serious disease that causes liver infection. This disease is caused by the Hepatitis B virus (HBV) and poses a significant health challenge worldwide. It can lead to both acute (short-term) and chronic (long-term) illnesses. Initially, it may be transmitted to a child from an infected mother before his/her birth (pregnancy period) and during delivery or childbirth. Its spread also takes place from an infected individual via sexual contact, through contact with the blood of an infected individual, or through unsafe injections. Its transmission may take place via contact with poisonous or polluted medical equipment or other objects. Similarly, injection practices can be a reason for its transmission. As claimed by the World Health Organization (WHO), about 296 million individuals globally are living with chronic infection. It is reported that in 2019, around 820,000 people died from this infection, with the majority of deaths attributed to cirrhosis and hepatocellular carcinoma (primary liver cancer). Moreover, 22% of the diagnosed population, which is almost 6.6 million individuals, were medically treated. This number makes up almost 10% of infected people.

The WHO reported a significant decline in the spread of chronic infection among children under five years of age. During the pre-vaccination period from the 1980s to the 2000s, the estimated rate was about 5%, whereas in 2019, it was less than 1%. However, we should note that the WHO predicts almost 1.5 million new cases of this serious disease every year. Preventive measures include antiviral prophylaxis during pregnancy, as well as safe and effective vaccinations. These interventions are of great importance in the prevention and control of Hepatitis B. The HBV models studied in the sense of classical derivatives are unable to capture memory and genetic characteristics, which can be observed in non-integer-order models.

Researchers and mathematicians have focused on fractional calculus due to its effective use in the modeling of various mathematical problems in many real-world situations. More specifically, researchers working in applied mathematics have taken an interest in non-integer-order derivatives to achieve higher accuracy in mathematical modeling (see [1,2,3,4,5,6,7]). Fractional calculus [8,9] is continuously being developed to more deeply and accurately understand the properties of real-world problems. It allows researchers in this field to analyze and explore phenomena that cannot be described accurately by integer-order calculus. Complex systems and phenomena can be better understood by using fractional derivatives and integrals. The advancements in fractional order calculus offer opportunities to enhance our understanding of real-world problems. In recent literature, various types of non-integer-order differential operators, along with their associated integral operators, have been explored. Caputo and Fabrizio [10] investigated non-integer-order derivatives with the exponential kernel. In [11], authors investigated properties of the new fractional derivative without singular kernel. Atangana and Baleanu [12] formulated fractional derivatives and integrals with Mittag-Leffler kernels and extended them to higher arbitrary orders.

In some real-world phenomena, transitions are observed when they shift from power law to exponential decay, or from deterministic to stochastic randomness. Atangana and Araz [13] proposed a new concept of piecewise differential and integral operators, which are effective for modeling and solving such phenomena. This approach is not like conventional methods, which do not exhibit crossover behavior due to their lack of abrupt changes. Additionally, the time interval is split into two subintervals at the point of discontinuity. In contrast to the usual fractional derivatives, the piecewise concept is useful in describing these crossover effects among various forms (see [14]).

On the other hand, researchers have focused on fractal–fractional calculus due to its significant use in real-world problems. The idea of fractal–fractional was formally started a few years ago when the author of reference [15] published the first work about this. The concept of fractals is as old as fractional calculus. Researchers have used the applications of fractal–fractional to investigate various problems in physical sciences. For instance, see [16,17]. For more applications, see [18,19].

Gul et al. [20] considered the HBV model. They studied the dynamics of the Caputo fractional order HBV model with asymptomatic carriers. They studied the existence, stability results, and numerical solution of the considered model and simulated the results. In [21], the authors investigated this model under the piecewise Atangana Baleanu derivative and derived the aforementioned results. In both research studies, the authors made significant efforts toward accurate modeling of the disease; however, these models lack a memory effect. As with many real-world problems, it is necessary to know how much information a system carries. In such a situation, the concept of a fractal–fractional derivative is beneficial. In our research paper, we reformulate the HBV model in the sense of the piecewise fractal–fractional order derivative, which incorporates the additional property of memory and manages the crossover effect. In this sense, our proposed model is more informative. Moreover, typically in a piecewise derivative, a classical derivative is taken in the first subinterval of the time domain with the concerned fractional derivative in the second subinterval. In our paper, we consider an FFD with a power law kernel in the first subinterval and an FFD with an exponential decay kernel in the second subinterval. This makes the results more interesting and different.

The innovations and objectives of this research are as follows:

• To define a piecewise Caputo FFD by combining the FFD with the power law kernel and FFD with the exponential decay-type kernel.
• To reformulate serum hepatitis disease in the sense of piecewise Caputo FFD.
• To study the existence and H-U type stability results for the proposed model under piecewise Caputo FFD.
• To find the numerical solution of the proposed model under piecewise Caputo FFD by applying the Lagrange interpolation method and the extended ABM method.
• To visually present our results.

The rest of the manuscript is structured as follows: In Section 2, we present the basic definitions and preliminary results. In Section 3, the proposed mathematical model and its formulation are presented. Section 4 is allocated to the derivation of the basic reproduction number of the concerned model. In Section 5, we analyze the fractional order system, derive the main results of existence, and present the uniqueness of the solution along with their stability results based on H-U’s concept of stability. In Section 6, the computational results along with graphical representations of the proposed model are presented. In Section 7, we present the conclusion of the main work with some future work directions.

2. Basic Results

Definitions 1 and 2 are adopted from [13]. In the upcoming definitions, $\varsigma$ denotes the fractional order derivative while $\eta$ denotes the fractal dimension.

Definition 1.

For $y\in C[0,T]$ and $0<\varsigma ,\eta \le 1,$ the piecewise Caputo-fractal–fractional derivative is defined by

$${}^{PCFF}{D}^{\varsigma ,\eta }y\left(\mathrm{t}\right)=\left\{\begin{array}{c}{}^{CFF}{D}_p^{\varsigma ,\eta }y\left(\mathrm{t}\right),\mathrm{t}\in [0,{\mathrm{t}}_1],\\ {}^{CFF}{D}_e^{\varsigma ,\eta }y\left(\mathrm{t}\right),\mathrm{t}\in ({\mathrm{t}}_1,T],\end{array}\right.$$

(1)

where the notion ${}^{CFF}{D}_p^{\varsigma ,\eta }y\left(\mathrm{t}\right)$ represents FFD with the power law kernel while the notion ${}^{CFF}{D}_e^{\varsigma ,\eta }y\left(\mathrm{t}\right)$ represents FFD with the exponential decay-type kernel.

Definition 2.

Let $y\in C(0,T)$. Then, we define the ςth order piecewise fractal–fractional integral (FFI) of $y\left(\mathrm{t}\right)$ as follows:

$${}^{PFF}{I}^{\varsigma ,\eta }y\left(\mathrm{t}\right)=\left\{\begin{array}{c}{}^{FF}{I}_p^{\varsigma ,\eta }y\left(\mathrm{t}\right),\mathrm{t}\in [0,{\mathrm{t}}_1],\\ {}^{FF}{I}_e^{\varsigma ,\eta }y\left(\mathrm{t}\right),\mathrm{t}\in ({\mathrm{t}}_1,T];\end{array}\right.$$

(2)

${}^{FF}{I}_p^{\varsigma ,\eta }y\left(\mathrm{t}\right)$ and ${}^{FF}{I}_e^{\varsigma ,\eta }y\left(\mathrm{t}\right)$ are defined bellow.

Definition 3

([15]). Let $y\in C(0,T)$. Then, the FFI of y associated with the power law kernel is defined as follows:

$${}^{FF}{I}_p^{\varsigma ,\eta }y\left(\mathrm{t}\right)=\frac{\eta }{\mathrm{\Gamma }\left(\varsigma \right)}{\int }_0^{\mathrm{t}}{v}^{\eta -1}{(\mathrm{t}-v)}^{\varsigma -1}y\left(v\right)dv.$$

(3)

Definition 4

([15]). Let $y\in C(0,T)$. Then, the FFI of $y\left(\mathrm{t}\right)$ associated with the exponential decay kernel is defined as follows:

$${}_{0 }{}^{FF}{I}_e^{\varsigma ,\eta }y\left(\mathrm{t}\right)=\frac{\eta (1-\varsigma ){\mathrm{t}}^{\eta -1}y\left(\mathrm{t}\right)}{\mathbf{B}\left(\varsigma \right)}+\frac{\varsigma \eta }{\mathbf{B}\left(\varsigma \right)}{\int }_a^{\mathrm{t}}{v}^{\varsigma -1}y\left(v\right)dv.$$

(4)

Definition 5

([15]). For $y\in C[0,T]$ and $0<\varsigma ,\eta \le 1$, if y is fractal-differentiable on $(a,b)$ with the order of $\eta ,$ then FFD with the power law kernel is defined as follows:

$${}^{CFF}{D}_p^{\varsigma ,\eta }y\left(\mathrm{t}\right)=\frac{1}{\mathrm{\Gamma }(1-\varsigma )}\frac{d}{d{t}^{\eta }}{\int }_0^{\mathrm{t}}{(\mathrm{t}-v)}^{-\varsigma }dv.$$

(5)

Definition 6

([15]). For $0\le \varsigma ,\eta \le 1,$ the FFD of y with the exponential decay kernel is given by the following:

$${}_{0 }{}^{CFF}{D}_e^{\varsigma ,\eta }y\left(\mathrm{t}\right)=\frac{\mathbf{B}\left(\varsigma \right)}{1-\varsigma }\frac{d}{d{t}^{\eta }}{\int }_a^{\mathrm{t}}exp\left(\frac{-\varsigma {(\mathrm{t}-v)}^{\alpha -\varsigma -1}}{1-\varsigma }\right)y\left(z\right)dv,$$

where $M\left(0\right)=M\left(1\right)=1.$

Theorem 1.

If ${\mathbb{Q}}_1,{\mathbb{Q}}_2$ are two operators such that the first is a contraction and the second is completely continuous over a closed bounded subset $\mathbf{H}$ of a Banach space J, then the operator equation $\mathbb{Q}\mathcal{F}+{\mathbb{Q}}_2\mathcal{F}=\mathcal{F}$ has at least one solution.

Definition 7.

The formula used for the Adams–Bashforth method of ordinary problems is given by the following:

$$\begin{array}{cc}\hfill y_{m+1}& =y_m+h\sum _{i=1}^r{b}_{i}f(x_m+ih,y_m+ih),\end{array}$$

where

• $y_m$ is the approximate solution at time $x_m$.
• h is the step size.
• $f(y,x)$ is the ODE.
• $b_i$ are coefficients that depend on the order of the method.

3. Mathematical Model and Its Formulation

This section is subdivided into two subsections, where the mathematical model and its formulation are presented, respectively.

3.1. Mathematical Model

We generalize the HBV model [20] in the sense of the piecewise FFD with the power law kernel and the exponential law kernel as follows:

$$\left\{\begin{array}{c}{}^{PCFF}{D}_{\mathrm{t}}^{\varsigma ,\eta }\mathbf{X}\left(\mathrm{t}\right)=\left\{\begin{array}{c}\varrho -\vartheta \left(\mathbf{A}+{\omega }_1{\mathbf{A}}_c+{\wp }_1\mathbf{C}\right)\mathbf{X}-\xi \mathbf{X},\\ \mathbf{X}\left(0\right)>0,\end{array}\right.\\ {}^{PCFF}{D}_{\mathrm{t}}^{\varsigma ,\eta }\mathbf{E}\left(\mathrm{t}\right)=\left\{\begin{array}{c}\vartheta \left(\mathbf{A}+{\omega }_1{\mathbf{A}}_c+{\wp }_1\mathbf{C}\right)\mathbf{X}-\left(\xi +{\psi }_1\right)\mathbf{E},\\ \mathbf{E}\left(0\right)>0,\end{array}\right.\\ {}^{PCFF}{D}_{\mathrm{t}}^{\varsigma ,\eta }\mathbf{A}\left(\mathrm{t}\right)=\left\{\begin{array}{c}{\psi }_{1}n\mathbf{E}-\left(\xi +\mu +{\rho }_1+{\eta }_1\right)\mathbf{A},\\ \mathbf{A}\left(0\right)>0,\end{array}\right.\\ {}^{PCFF}{D}_{\mathrm{t}}^{\varsigma ,\eta }{\mathbf{A}}_c\left(\mathrm{t}\right)=\left\{\begin{array}{c}{\psi }_1\left(1-n\right)\mathbf{E}-\left(\xi +x_1+\lambda \right){\mathbf{A}}_c,\\ {\mathbf{A}}_c\left(0\right)>0,\end{array}\right.\\ {}^{PCFF}{D}_{\mathrm{t}}^{\varsigma ,\eta }\mathbf{C}\left(\mathrm{t}\right)=\left\{\begin{array}{c}{\rho }_1\mathbf{A}+x_1{\mathbf{A}}_c-\left(\xi +\delta +v_1\right)\mathbf{C},\\ \mathbf{C}\left(0\right)>0,\end{array}\right.\\ {}^{PCFF}{D}_{\mathrm{t}}^{\varsigma ,\eta }{\mathbf{R}}_p\left(\mathrm{t}\right)=\left\{\begin{array}{c}{\eta }_1\mathbf{A}+v_1\mathbf{C}+\lambda {\mathbf{A}}_c-\xi {\mathbf{R}}_p,\\ {\mathbf{R}}_p\left(0\right)>0;\end{array}\right.\end{array}\right.$$

(6)

$\varrho$ is the susceptible birth rate, and $\vartheta$ and $\xi$ denote the effective contact rate and natural fatality rate, respectively. ${\psi }_1\left(1-n\right)$ is the infection rate of the exposed population, with a portion of ${\psi }_1\left(1-n\right)$, shifting to class $\mathbf{A}$ at a rate of ${\psi }_{1}n$. Class ${\mathbf{A}}_c$ is used for symptomatically infected individuals. ${\rho }_1$ and $x_1$ denote the rates at which individuals of acute and asymptomatic classes become carriers, respectively. ${\eta }_1$, $\lambda$, and $v_1$ denote the recovery rates for acute, asymptomatic, and carrier individuals, respectively. $\mu$ and $\delta$ denote the death rates in the acute and chronic classes due to disease, respectively. ${\omega }_1$ and ${\wp }_1$ denote the coefficients of asymptomatic and carrier individuals, respectively.

3.2. Model Formulation

In deriving the integer-order model of serum hepatitis disease, we take the total human population, $\mathbf{N}\left(\mathrm{t}\right)$ as a sum of six classes; that is, $\mathbf{N}\left(\mathrm{t}\right)=\mathbf{X}\left(\mathrm{t}\right)+\mathbf{E}\left(\mathrm{t}\right)+\mathbf{A}\left(\mathrm{t}\right)+{\mathbf{A}}_c\left(\mathrm{t}\right)+\mathbf{C}\left(\mathrm{t}\right)+{\mathbf{R}}_p\left(\mathrm{t}\right),$ where the parameters are defined in

Table 1. Parameters for system (6).

Parameters	Parameter Definition
$\mathbf{X}\left(\mathrm{t}\right)$	Susceptible class of individuals.
$\mathbf{E}\left(\mathrm{t}\right)$	Exposed class of population.
$\mathbf{A}\left(\mathrm{t}\right)$	Acute class of infected individuals.
${\mathbf{A}}_c\left(\mathrm{t}\right)$	Asymptomatic carrier.
$\mathbf{C}\left(\mathrm{t}\right)$	Chronic class of infected individuals.
${\mathbf{R}}_p\left(\mathrm{t}\right)$	Recovered class of individuals.

.

Susceptible class: The class $\mathbf{X}\left(\mathrm{t}\right)$ is recruited at the rate $\varrho$ and decreased due to the natural fatality (death) rate $\xi .$ Also, this class diminishes as a proportion of the population becomes infected after contact with either Class $\mathbf{A}$ of acutely infected individuals or the asymptomatic carrier and chronic class $\mathbf{C}\left(\mathrm{t}\right)$ of infected individuals, at rates of ${\omega }_1$ and ${\wp }_1,$ respectively. Therefore, we formulated the dynamics of the susceptible class of individuals, as follows:

$$\frac{d\mathbf{X}\left(\mathrm{t}\right)}{dt}=\varrho -\vartheta \left(\mathbf{A}+{\omega }_1{\mathbf{A}}_c+{\wp }_1\mathbf{C}\right)\mathbf{X}-\xi \mathbf{X}.$$

Exposed class: Class $\mathbf{E}\left(\mathrm{t}\right)$ increases at the rate $\vartheta$, which represents the proportion of the population becoming newly infected through previously described contacts. This class decreases due to the natural fatality rate, $\xi$, and by the rate at which the exposed population becomes acutely infected. Therefore, we formulate the dynamics of the exposed class of the population as follows:

$$\frac{d\mathbf{E}\left(\mathrm{t}\right)}{dt}=\vartheta \left(\mathbf{A}+{\omega }_1{\mathbf{A}}_c+{\wp }_1\mathbf{C}\right)\mathbf{X}-\left(\xi +{\psi }_1\right)\mathbf{E}.$$

Acute class: Class $\mathbf{A}\left(\mathrm{t}\right)$ is recruited at the rate ${\psi }_{1}n$ of the exposed population. This class decreases at the natural fatality rate, $\xi$, by the death rate, $\mu$, due to disease in the acute class. This class also decreases by the rate ${\rho }_1$, at which the individuals of the class go to the asymptomatic carrier, and by the rate ${\eta }_1$, at which the individuals of the class recover. Therefore, we formulate the dynamics of the acute class of the population as follows:

$$\frac{d\mathbf{A}\left(\mathrm{t}\right)}{dt}={\psi }_{1}n\mathbf{E}-\left(\xi +\mu +{\rho }_1+{\eta }_1\right)\mathbf{A}.$$

Asymptomatic carrier: Class ${\mathbf{A}}_c\left(\mathrm{t}\right)$ is recruited from the exposed class at the rate ${\psi }_1$ and decreases at the rate ${\psi }_{1}n$, which represents the transition of individuals from the exposed to the acute class. This class decreases due to the natural fatality rate $\xi$, the rate, $x_1$, at which members transition to the chronic class of infected individuals. Also, it decreases at the rate $\lambda$, at which members of the class recover. Therefore, we formulate the dynamics of the asymptomatic carrier as follows:

$$\frac{d{\mathbf{A}}_c\left(\mathrm{t}\right)}{dt}={\psi }_1\left(1-n\right)\mathbf{E}-\left(\xi +x_1+\lambda \right){\mathbf{A}}_c.$$

Chronic class of infected individuals: Class $\mathbf{C}\left(\mathrm{t}\right)$ is recruited at the rates of ${\rho }_1$ and $x_1$ of the acute class of infected individuals and asymptomatic carriers, respectively. The class is reduced by the natural fatality rate, $\xi$, as well as the $\delta$ fatality rate, due to disease. The class also reduces by the rate, $v_1$, at which the individuals in the class recover. Therefore, we formulate the dynamics of the chronic class of infected individuals as follows:

$$\frac{d\mathbf{C}\left(\mathrm{t}\right)}{dt}={\rho }_1\mathbf{A}+x_1{\mathbf{A}}_c-\left(\xi +\delta +v_1\right)\mathbf{C}.$$

Recovered class of individuals: Class ${\mathbf{R}}_p\left(\mathrm{t}\right)$ grows at rates ${\eta }_1$, $v_1$, and $\lambda$, where the individuals of the aforementioned classes recover and decrease at a natural fatality rate, $\xi$. Therefore, we formulate the dynamics of a recovered class of individuals as follows:

$$\frac{d{\mathbf{R}}_p\left(\mathrm{t}\right)}{dt}={\eta }_1\mathbf{A}+v_1\mathbf{C}+\lambda {\mathbf{A}}_c-\xi {\mathbf{R}}_p.$$

As a result of the above derivations of the aforementioned classes, the integer-order model of serum hepatitis disease with a class of asymptomatic carriers is presented as follows:

$$\left\{\begin{array}{c}\frac{d\mathbf{X}\left(\mathrm{t}\right)}{dt}=\left\{\begin{array}{c}\varrho -\vartheta \left(\mathbf{A}+{\omega }_1{\mathbf{A}}_c+{\wp }_1\mathbf{C}\right)\mathbf{X}-\xi \mathbf{X},\\ \mathbf{X}\left(0\right)>0,\end{array}\right.\\ \frac{d\mathbf{E}\left(\mathrm{t}\right)}{dt}=\left\{\begin{array}{c}\vartheta \left(\mathbf{A}+{\omega }_1{\mathbf{A}}_c+{\wp }_1\mathbf{C}\right)\mathbf{X}-\left(\xi +{\psi }_1\right)\mathbf{E},\\ \mathbf{E}\left(0\right)>0,\end{array}\right.\\ \frac{d\mathbf{A}\left(\mathrm{t}\right)}{dt}=\left\{\begin{array}{c}{\psi }_{1}n\mathbf{E}-\left(\xi +\mu +{\rho }_1+{\eta }_1\right)\mathbf{A},\\ \mathbf{A}\left(0\right)>0,\end{array}\right.\\ \frac{d{\mathbf{A}}_c\left(\mathrm{t}\right)}{dt}=\left\{\begin{array}{c}{\psi }_1\left(1-n\right)\mathbf{E}-\left(\xi +x_1+\lambda \right){\mathbf{A}}_c,\\ {\mathbf{A}}_c\left(0\right)>0,\end{array}\right.\\ \frac{d\mathbf{C}\left(\mathrm{t}\right)}{dt}=\left\{\begin{array}{c}{\rho }_1\mathbf{A}+x_1{\mathbf{A}}_c-\left(\xi +\delta +v_1\right)\mathbf{C},\\ \mathbf{C}\left(0\right)>0,\end{array}\right.\\ \frac{d{\mathbf{R}}_p\left(\mathrm{t}\right)}{dt}=\left\{\begin{array}{c}{\eta }_1\mathbf{A}+v_1\mathbf{C}+\lambda {\mathbf{A}}_c-\xi {\mathbf{R}}_p,\\ {\mathbf{R}}_p\left(0\right)>0.\end{array}\right.\end{array}\right.$$

(7)

We extend the integer-order model presented in (7) to a piecewise fractal–fractional order model with power law and exponential law kernels as given in (6).

4. Equilibrium Point and Basic Reproduction Number

For disease-free equilibrium point ${\mathrm{\Theta }}_0$, we have the following:

$$\left\{\begin{array}{c}\varrho -\vartheta \left(\mathbf{A}+{\omega }_1{\mathbf{A}}_c+{\wp }_1\mathbf{C}\right)\mathbf{X}-\xi \mathbf{X}=0\\ \vartheta \left(\mathbf{A}+{\omega }_1{\mathbf{A}}_c+{\wp }_1\mathbf{C}\right)\mathbf{X}-\left(\xi +{\psi }_1\right)\mathbf{E}=0\\ {\psi }_{1}n\mathbf{E}-\left(\xi +\mu +{\rho }_1+{\eta }_1\right)\mathbf{A}=0\\ {\psi }_1\left(1-n\right)\mathbf{E}-\left(\xi +x_1+\lambda \right){\mathbf{A}}_c=0\\ {\rho }_1\mathbf{A}+x_1{\mathbf{A}}_c-\left(\xi +\delta +v_1\right)\mathbf{C}=0\\ {\eta }_1\mathbf{A}+v_1\mathbf{C}+\lambda {\mathbf{A}}_c-\xi {\mathbf{R}}_p=0.\end{array}\right.$$

And is provided as follows:

$${\mathrm{\Theta }}_0=\left(\frac{\varrho }{\xi },0,0,0,0,0\right);$$

$\varrho$ is the birth rate of the affected population and $\xi$ denotes the natural dying rate, respectively. From [22], the basic reproduction number, $R_0$, for model (6) is given by the following:

$$R_0=\frac{\vartheta \varrho n{\psi }_1}{\xi {\theta }_2{\theta }_1}+\frac{\vartheta {\omega }_1\varrho {\psi }_1\left(1-n\right)}{\xi {\theta }_3{\theta }_1}+\frac{\vartheta {\wp }_1\varrho {\rho }_{1}n{\psi }_1}{\xi {\theta }_4{\theta }_1{\theta }_2}+\frac{\vartheta x_1\varrho {\psi }_1{\wp }_1\left(1-n\right)}{\xi {\theta }_4{\theta }_3{\theta }_1},$$

where

$$\begin{array}{ccc}\hfill {\theta }_1& =& \left(\xi +{\psi }_1\right),{\theta }_2=\left(\xi +\mu +{\rho }_1+{\eta }_1\right),\\ \hfill {\theta }_3& =& \left(\xi +x_1+\lambda \right),{\theta }_4=\left(\xi +\delta +v_1\right).\end{array}$$

Significance of the Basic Reproduction Number

The basic reproduction number is the average number of infected contacts per sick individual. It is a dimensionless number rather than a rate. Therefore, $R_0$ is the number of secondary infections produced by a single typical infection in a rarefied population. It plays an important role in predicting the infection in society. If there are no environmental changes or outside factors that intervene, a virus will continue to spread among susceptible hosts at the population level if $R_0$ is greater than 1.

5. Existence and Stability Analysis of the Piecewise Fractal–Fractional Model 6

In this section, the main results, such as the existence, uniqueness, and stability of the solution of the proposed model (6), are established. Let $\mathbf{I}=\left[0,T\right].$ We define a Banach space as follows: $\mathbf{J}=C\left(\mathbf{I},{\mathbf{R}}^{+}\right)\times C\left(\mathbf{I},{\mathbb{R}}^{+}\right)\times C\left(\mathbf{I},{\mathbb{R}}^{+}\right)\times C\left(\mathbf{I},{\mathbb{R}}^{+}\right)\times C\left(\mathbf{I},{\mathbb{R}}^{+}\right)\times C\left(\mathbf{I},{\mathbb{R}}^{+}\right)$ under the given norm, as follows:

$$∥s∥=sup\left\{\left|\mathbf{X}\left(\mathrm{t}\right)\right|+\left|\mathbf{E}\left(\mathrm{t}\right)\right|+\left|\mathbf{A}\left(\mathrm{t}\right)\right|+\left|{\mathbf{A}}_c\left(\mathrm{t}\right)\right|+\left|\mathbf{C}\left(\mathrm{t}\right)\right|+\left|{\mathbf{R}}_P\left(\mathrm{t}\right)\right|\right\};$$

$\mathbf{X},\mathbf{E},\mathbf{A},{\mathbf{A}}_c,\mathbf{C},{\mathbf{R}}_P\in \mathbf{J}.$

Lemma 1.

The piecewise-fractal–fractional problem is as follows:

$$\left\{\begin{array}{c}{}^{CFF}{D}_p^{\varsigma ,\eta }s\left(\mathrm{t}\right)=\left\{\begin{array}{c}\phi \left(\mathrm{t}\right),0<\varsigma ,\eta \le 1,if \mathrm{t}\in \left[0,{\mathrm{t}}_1\right],\\ s\left(0\right)=s_0,\end{array}\right.\\ {}^{CFF}{D}_e^{\varsigma ,\eta }s\left(\mathrm{t}\right)=\left\{\begin{array}{c}\phi \left(\mathrm{t}\right),0<\varsigma ,\eta \le 1,if \mathrm{t}\in \left({\mathrm{t}}_1,T\right],\\ s\left({\mathrm{t}}_1\right)=s_1,\end{array}\right.\end{array}\right.$$

(8)

has the following solution:

$$s\left(\mathrm{t}\right)=\left\{\begin{array}{c}{s}_0+\frac{\eta }{\mathrm{\Gamma }\left(\varsigma \right)}{\int }_0^{{\mathrm{t}}_1}{v}^{\eta -1}{(x-v)}^{\varsigma -1}\phi \left(v\right)dv,if \mathrm{t}\in \left[0,{\mathrm{t}}_1\right],\\ s\left({\mathrm{t}}_1\right)+\frac{(1-\varsigma )}{\mathbf{B}\left(\varsigma \right)}\eta {\mathrm{t}}^{\eta -1}\phi \left(\mathrm{t}\right)+\frac{\varsigma \eta }{\mathbf{B}\left(\varsigma \right)}{\int }_{{\mathrm{t}}_1}^{\mathrm{t}}{v}^{\eta -1}\phi \left(v\right)dv,if \mathrm{t}\in \left({\mathrm{t}}_1,T\right].\end{array}\right.$$

We reformulate model (6) as follows:

$$\left\{\begin{array}{c}{}^{PCFF}{D}^{\varsigma ,\eta }\mathbf{X}\left(\mathrm{t}\right)={\mathrm{\Phi }}_1(\mathrm{t},s\left(\mathrm{t}\right)),\\ {}^{PCFF}{D}^{\varsigma ,\eta }\mathbf{E}\left(\mathrm{t}\right)={\mathrm{\Phi }}_2(\mathrm{t},s\left(\mathrm{t}\right)),\\ {}^{PCFF}{D}^{\varsigma ,\eta }\mathbf{A}\left(\mathrm{t}\right)={\mathrm{\Phi }}_3(\mathrm{t},s\left(\mathrm{t}\right)),\\ {}^{PCFF}{D}^{\varsigma ,\eta }{\mathbf{A}}_{\mathbf{C}}\left(\mathrm{t}\right)={\mathrm{\Phi }}_4(\mathrm{t},s\left(\mathrm{t}\right)),\\ {}^{PCFF}{D}^{\varsigma ,\eta }\mathbf{C}\left(\mathrm{t}\right)={\mathrm{\Phi }}_5(\mathrm{t},s\left(\mathrm{t}\right)),\\ {}^{PCFF}{D}^{\varsigma ,\eta }{\mathbf{R}}_p\left(\mathrm{t}\right)={\mathrm{\Phi }}_6(\mathrm{t},s\left(\mathrm{t}\right)),\end{array}\right.$$

(9)

where

$$\left\{\begin{array}{c}{\mathrm{\Phi }}_1(\mathrm{t},s\left(\mathrm{t}\right))=\varrho -\vartheta \left(\mathbf{A}+{\omega }_1{\mathbf{A}}_c+{\wp }_1\mathbf{C}\right)\mathbf{X}-\xi \mathbf{X},\\ {\mathrm{\Phi }}_2(\mathrm{t},s\left(\mathrm{t}\right))=\vartheta \left(\mathbf{A}+{\omega }_1{\mathbf{A}}_c+{\wp }_1\mathbf{C}\right)\mathbf{X}-\left(\xi +{\psi }_1\right)\mathbf{E},\\ {\mathrm{\Phi }}_3(\mathrm{t},s\left(\mathrm{t}\right))={\psi }_{1}n\mathbf{E}-\left(\xi +\mu +{\rho }_1+{\eta }_1\right)\mathbf{A},\\ {\mathrm{\Phi }}_4(\mathrm{t},s\left(\mathrm{t}\right))={\psi }_1\left(1-n\right)\mathbf{E}-\left(\xi +x_1+\lambda \right){\mathbf{A}}_c,\\ {\mathrm{\Phi }}_5(\mathrm{t},s\left(\mathrm{t}\right))={\rho }_1\mathbf{A}+x_1{\mathbf{A}}_c-\left(\xi +\delta +v_1\right)\mathbf{C},\\ {\mathrm{\Phi }}_6(\mathrm{t},s\left(\mathrm{t}\right))={\eta }_1\mathbf{A}+v_1\mathbf{C}+\lambda {\mathbf{A}}_c-\xi {\mathbf{R}}_p.\end{array}\right.$$

We take our model as follows:

$$\left\{\begin{array}{c}{}^{PCFF}{D}^{\varsigma ,\eta }s\left(\mathrm{t}\right)=\mathcal{F}(\mathrm{t},s\left(\mathrm{t}\right)),\\ s\left(0\right)=s_0>0.\end{array}\right.$$

(10)

In view of Definitions 2 and 3 and Lemma 1, the equivalent form of model (10) is given by the following:

$$s\left(\mathrm{t}\right)=\left\{\begin{array}{c}{s}_0+\frac{\eta }{\mathrm{\Gamma }\left(\varsigma \right)}{\int }_0^{\mathrm{t}}{v}^{\eta -1}{(\mathrm{t}-v)}^{\varsigma -1}\mathcal{F}(v,s\left(v\right))dv,\mathrm{i}\mathrm{f} \mathrm{t}\in \left[0,{\mathrm{t}}_1\right],\\ s\left({\mathrm{t}}_1\right)+\frac{(1-\varsigma )}{\mathbf{B}\left(\varsigma \right)}\eta {\mathrm{t}}^{\eta -1}\mathcal{F}(\mathrm{t},s\left(\mathrm{t}\right))+\frac{\varsigma \eta }{\mathbf{B}\left(\varsigma \right)}{\int }_{{\mathrm{t}}_1}^{\mathrm{t}}{v}^{\eta -1}\mathcal{F}(v,s\left(v\right))dv,\mathrm{i}\mathrm{f} \mathrm{t}\in \left({\mathrm{t}}_1,T\right].\end{array}\right.$$

(11)

where

$$s\left(\mathrm{t}\right)=\left(\begin{array}{c}\mathbf{X}\left(\mathrm{t}\right)\\ \mathbf{E}\left(\mathrm{t}\right)\\ \mathbf{A}\left(\mathrm{t}\right)\\ {\mathbf{A}}_c\left(\mathrm{t}\right)\\ \mathbf{C}\left(\mathrm{t}\right)\\ {\mathbf{R}}_p\left(\mathrm{t}\right)\end{array}\right),s\left(0\right)=\left(\begin{array}{c}{\mathbf{X}}_0\\ {\mathbf{E}}_0\\ {\mathbf{A}}_0\\ {\mathbf{A}}_c\\ {\mathbf{C}}_0\\ {\mathbf{R}}_{p0}\end{array}\right),s\left({\mathrm{t}}_1\right)=\left(\begin{array}{c}{\mathbf{X}}_{{\mathrm{t}}_1}\\ {\mathbf{E}}_{{\mathrm{t}}_1}\\ {\mathbf{A}}_{{\mathrm{t}}_1}\\ {\mathbf{A}}_c\\ {\mathbf{C}}_{{\mathrm{t}}_1}\\ {\mathbf{R}}_{p{\mathrm{t}}_1}\end{array}\right).$$

and

$$\mathcal{F}(\mathrm{t},s\left(\mathrm{t}\right))=\left(\begin{array}{c}{\mathrm{\Phi }}_1(\mathrm{t},s\left(\mathrm{t}\right))\\ {\mathrm{\Phi }}_2(\mathrm{t},s\left(\mathrm{t}\right))\\ {\mathrm{\Phi }}_3(\mathrm{t},s\left(\mathrm{t}\right))\\ {\mathrm{\Phi }}_4(\mathrm{t},s\left(\mathrm{t}\right))\\ {\mathrm{\Phi }}_5(\mathrm{t},s\left(\mathrm{t}\right))\\ {\mathrm{\Phi }}_6(\mathrm{t},s\left(\mathrm{t}\right))\end{array}\right).$$

Now, we define an operator $W:\mathbf{J}\to \mathbf{J}$ by the following:

$$W\left(s\left(\mathrm{t}\right)\right)=\left\{\begin{array}{c}{s}_0+\frac{\eta }{\mathrm{\Gamma }\left(\varsigma \right)}{\int }_0^{\mathrm{t}}{v}^{\eta -1}{(\mathrm{t}-v)}^{\varsigma -1}\mathcal{F}(v,s\left(v\right))dv,\mathrm{i}\mathrm{f} \mathrm{t}\in \left[0,{\mathrm{t}}_1\right],\\ s\left({\mathrm{t}}_1\right)+\frac{(1-\varsigma )}{\mathbf{B}\left(\varsigma \right)}\eta {\mathrm{t}}^{\eta -1}\mathcal{F}(\mathrm{t},s\left(\mathrm{t}\right))+\frac{\varsigma \eta }{\mathbf{B}\left(\varsigma \right)}{\int }_{{\mathrm{t}}_1}^{\mathrm{t}}{v}^{\eta -1}\mathcal{F}(v,s\left(v\right))dv,\mathrm{i}\mathrm{f} \mathrm{t}\in \left({\mathrm{t}}_1,T\right].\end{array}\right.$$

The following assumptions are necessary for the analysis of existence and uniqueness.

Hypothesis 1

(H1). $\mathcal{F}:J\times \mathbf{J}\to \mathbf{R}$ is continuous and there exist two constants $\mathbb{k},\mathbf{q}>0$, such that we have the following:

$$\left|\mathcal{F}(\mathrm{t},s(\mathrm{t}\left)\right)\right|\le \mathbb{k}+\left|s\left(\mathrm{t}\right)\right|\mathbf{q},for v\in Jand \mathcal{Y}\in \mathbf{J}.$$

Hypothesis 2

(H2). Assume that the real number $\mathcal{L}>0$ satisfies the following:

$$\left|\mathcal{F}(\mathrm{t},s_1\left(\mathrm{t}\right))-\mathcal{F}(\mathrm{t},s_2\left(\mathrm{t}\right))\right|\le \mathcal{L}\left|s_1\left(\mathrm{t}\right)-s_2\left(\mathrm{t}\right)\right|,for \mathrm{t}\in Jand s_1,s_2\in \mathbf{J}.$$

Theorem 2.

Under assumptions (H1)–(H2), system (9) has at least one solution with the following condition:

$$0<max\left\{\frac{\eta \mathcal{L}\beta (\eta ,\varsigma ){{\mathrm{t}}_1}^{\eta +\varsigma -1}}{\mathrm{\Gamma }\left(\varsigma \right)},\left(1-\frac{(1-\varsigma )}{\mathbf{B}\left(\varsigma \right)}\eta {\mathrm{t}}^{\eta -1}\mathcal{L}\right)\right\}<1.$$

(12)

Proof. 

We are transforming the fractional system (9) into a fixed point problem as the following equation:

$$s=W\left(s\left(\mathrm{t}\right)\right),s\in \mathbf{J}.$$

where the operator $W:\mathbf{J}\to \mathbf{J}$ is defined by the following:

$$W\left(s\left(\mathrm{t}\right)\right)=\left\{\begin{array}{c}{s}_0+\frac{\eta }{\mathrm{\Gamma }\left(\varsigma \right)}{\int }_0^{\mathrm{t}}{v}^{\eta -1}{(\mathrm{t}-v)}^{\varsigma -1}\mathcal{F}(v,s\left(v\right))dv,\mathrm{i}\mathrm{f} \mathrm{t}\in \left[0,{\mathrm{t}}_1\right],\\ s\left({\mathrm{t}}_1\right)+\frac{(1-\varsigma )}{\mathbf{B}\left(\varsigma \right)}\eta {\mathrm{t}}^{\eta -1}\mathcal{F}(\mathrm{t},s\left(\mathrm{t}\right))+\frac{\varsigma \eta }{\mathbf{B}\left(\varsigma \right)}{\int }_{{\mathrm{t}}_1}^{\mathrm{t}}{v}^{\eta -1}\mathcal{F}(v,s\left(v\right))dv,\mathrm{i}\mathrm{f} \mathrm{t}\in \left({\mathrm{t}}_1,T\right].\end{array}\right.$$

(13)

Let ${\mathrm{\Lambda }}_{\zeta }=\left\{s\in \mathbf{J}:∥s∥\le \zeta \right\}$ be a closed ball with the following:

$$\zeta \ge \left\{\begin{array}{c}\frac{∥s\left(0\right)∥+\frac{\eta \mathbb{k}{{\mathrm{t}}_1}^{\eta +\varsigma -1}}{\mathrm{\Gamma }\left(\varsigma \right)}\beta (\eta ,\varsigma )}{1-\frac{\eta \mathbf{q}{{\mathrm{t}}_1}^{\eta +\varsigma -1}}{\mathrm{\Gamma }\left(\varsigma \right)}\beta (\eta ,\varsigma )},\mathrm{i}\mathrm{f} \mathrm{t}\in \left[0,{\mathrm{t}}_1\right],\\ \frac{\left|s\left({\mathrm{t}}_1\right)\right|+\left(\frac{(1-\varsigma )}{\mathbf{B}\left(\varsigma \right)}\eta {\mathrm{t}}^{\eta -1}+\frac{\varsigma T^{\varsigma }}{\mathbf{B}\left(\varsigma \right)}\right)\mathbb{k}}{1-\left(\frac{(1-\varsigma )}{\mathbf{B}\left(\varsigma \right)}\eta {\mathrm{t}}^{\eta -1}+\frac{\varsigma T^{\varsigma }}{\mathbf{B}\left(\varsigma \right)}\right)\mathbf{q}},\mathrm{i}\mathrm{f} \mathrm{t}\in \left({\mathrm{t}}_1,T\right].\end{array}\right.$$

Define the operators $W_1$ and $W_2$, such that $W=W_1+W_2$, as follows:

$$W_{1}s\left(\mathrm{t}\right)=\left\{\begin{array}{c}{s}_0+\frac{\eta }{\mathrm{\Gamma }\left(\varsigma \right)}{\int }_0^{\mathrm{t}}{v}^{\eta -1}{(\mathrm{t}-v)}^{\varsigma -1}\mathcal{F}(v,s\left(v\right))dv,\mathrm{i}\mathrm{f}\mathrm{t}\in \left[0,{\mathrm{t}}_1\right],\\ s\left({\mathrm{t}}_1\right)+\frac{(1-\varsigma )}{\mathbf{B}\left(\varsigma \right)}\eta {\mathrm{t}}^{\eta -1}\mathcal{F}(\mathrm{t},s\left(\mathrm{t}\right)),\mathrm{i}\mathrm{f} \mathrm{t}\in \left({\mathrm{t}}_1,T\right].\end{array}\right.$$

and

$$W_{2}s\left(\mathrm{t}\right)=\left\{\begin{array}{l}0,\mathrm{i}\mathrm{f}\mathrm{t}\in \left[0,{\mathrm{t}}_1\right],\\ \frac{\varsigma \eta }{\mathbf{B}\left(\varsigma \right)}{\int }_{{\mathrm{t}}_1}^{\mathrm{t}}{v}^{\eta -1}\mathcal{F}(v,s\left(v\right))dv,\mathrm{i}\mathrm{f} \mathrm{t}\in \left({\mathrm{t}}_1,T\right].\end{array}\right.$$

Now, we will divide the proof into several steps, as follows:

Step 1: $W_{1}s\left(\mathrm{t}\right)+W_{2}s\left(\mathrm{t}\right)\in {\mathrm{\Lambda }}_{\zeta }.$ If $\mathrm{t}\in \left[0,{\mathrm{t}}_1\right],$ $s\in {\mathrm{\Lambda }}_{\zeta },$ with (H1), we have the following:

$$\begin{array}{l}\left|W_{1}s\left(\mathrm{t}\right)+W_{2}s\left(\mathrm{t}\right)\right|=\left|s_0+\frac{\eta }{\mathrm{\Gamma }\left(\varsigma \right)}{\int }_0^{\mathrm{t}}{v}^{\eta -1}{(\mathrm{t}-v)}^{\varsigma -1}\mathcal{F}(v,s\left(v\right))dv\right|\\ \le \left|s_0\right|+\frac{\eta }{\mathrm{\Gamma }\left(\varsigma \right)}{\int }_0^{\mathrm{t}}{v}^{\eta -1}{(\mathrm{t}-v)}^{\varsigma -1}\left|\mathcal{F}(v,s(v\left)\right)\right|dv.\end{array}$$

(14)

Consider the integral ${\int }_0^{\mathrm{t}}{v}^{\eta -1}{(\mathrm{t}-v)}^{\varsigma -1}dv$. Let $v={\mathrm{t}}_{1}u.$ This implies that $dv={\mathrm{t}}_{1}du.$ If $v=0$ then $u=0$ and if $v={\mathrm{t}}_1$ then $u=1$. Thus,

$${\int }_0^{\mathrm{t}}{v}^{\eta -1}{(\mathrm{t}-v)}^{\varsigma -1}dv={{\mathrm{t}}_1}^{\eta +\varsigma -1}{\int }_0^1{\rho }^{\eta -1}{(1-\rho )}^{\varsigma -1}dv.$$

(15)

Hence, from (14), we obtain the following:

$$\begin{array}{l}∥W_{1}s+W_{2}s∥\le ∥s\left(0\right)∥+\frac{\eta \left(\mathbb{k}+\zeta \mathbf{q}\right){{\mathrm{t}}_1}^{\eta +\varsigma -1}}{\mathrm{\Gamma }\left(\varsigma \right)}{\int }_0^1{\rho }^{\eta -1}{(1-\rho )}^{\varsigma -1}dv\\ \le ∥s\left(0\right)∥+\frac{\eta \left(\mathbb{k}+\zeta \mathbf{q}\right){{\mathrm{t}}_1}^{\eta +\varsigma -1}}{\mathrm{\Gamma }\left(\varsigma \right)}\beta (\eta ,\varsigma )\le \zeta ,\end{array}$$

(16)

where $\beta (\eta ,\varsigma )$ is the well-known beta function, which is defined by $\beta (\eta ,\varsigma )={\int }_0^1{\rho }^{\eta -1}{(1-\rho )}^{\varsigma -1}dv.$

For $\mathrm{t}\in \left({\mathrm{t}}_1,T\right],$ $s\in {\mathrm{\Lambda }}_{\zeta },$ with (H1), we have the following:

$$\begin{array}{ccc}\hfill \left|W_{1}s\left(\mathrm{t}\right)+W_{2}s\left(\mathrm{t}\right)\right|& =& \left|s\left({\mathrm{t}}_1\right)+\frac{(1-\varsigma )}{\mathbf{B}\left(\varsigma \right)}\eta {\mathrm{t}}^{\eta -1}\mathcal{F}(\mathrm{t},s\left(\mathrm{t}\right))+\frac{\varsigma \eta }{\mathbf{B}\left(\varsigma \right)}{\int }_{{\mathrm{t}}_1}^{\mathrm{t}}{v}^{\eta -1}\mathcal{F}(v,s\left(v\right))dv\right|\hfill \\ & \le & \left|s\left({\mathrm{t}}_1\right)\right|+\frac{(1-\varsigma )}{\mathbf{B}\left(\varsigma \right)}\eta {\mathrm{t}}^{\eta -1}\left|\mathcal{F}(\mathrm{t},s(\mathrm{t}\left)\right)\right|+\frac{\varsigma \eta }{\mathbf{B}\left(\varsigma \right)}{\int }_{{\mathrm{t}}_1}^{\mathrm{t}}{v}^{\eta -1}\left|\mathcal{F}(v,s(v\left)\right)\right|dv\hfill \\ & \le & \left|s\left({\mathrm{t}}_1\right)\right|+\frac{(1-\varsigma )}{\mathbf{B}\left(\varsigma \right)}\eta {\mathrm{t}}^{\eta -1}\left(\mathbb{k}+\left|s\left(\mathrm{t}\right)\right|\mathbf{q}\right)+\frac{\varsigma \eta }{\mathbf{B}\left(\varsigma \right)}{\int }_{{\mathrm{t}}_1}^{\mathrm{t}}{v}^{\eta -1}\left(\mathbb{k}+\left|s\left(v\right)\right|\mathbf{q}\right)dv.\hfill \end{array}$$

Hence,

$$\begin{array}{lll}\hfill ∥W_{1}s+W_{2}s∥& \le & \left|s\left({\mathrm{t}}_1\right)\right|+\left(\frac{(1-\varsigma )}{\mathbf{B}\left(\varsigma \right)}\eta {\mathrm{t}}^{\eta -1}+\frac{\varsigma T^{\varsigma }}{\mathbf{B}\left(\varsigma \right)}\right)\left(\mathbb{k}+∥s∥\mathbf{q}\right)\hfill \\ & \le & \left|s\left({\mathrm{t}}_1\right)\right|+\left(\frac{(1-\varsigma )}{\mathbf{B}\left(\varsigma \right)}\eta {\mathrm{t}}^{\eta -1}+\frac{\varsigma T^{\varsigma }}{\mathbf{B}\left(\varsigma \right)}\right)\mathbb{k}\hfill \\ & & +\left(\frac{(1-\varsigma )}{\mathbf{B}\left(\varsigma \right)}\eta {\mathrm{t}}^{\eta -1}+\frac{\varsigma T^{\varsigma }}{\mathbf{B}\left(\varsigma \right)}\right)\mathbf{q}\zeta \hfill \\ & \le & \zeta .\hfill \end{array}$$

This demonstrates that $W_{1}s\left(\mathrm{t}\right)+W_{2}s\left(\mathrm{t}\right)\in {\mathrm{\Lambda }}_{\zeta }.$

Step 2: $W_1$ is the contraction.

For $\mathrm{t}\in \left[0,{\mathrm{t}}_1\right],$ $s_1,s_2\in {\mathrm{\Lambda }}_{\zeta }$. Then,

$$\begin{array}{ccc}\left|W_1{s}_1\left(\mathrm{t}\right)-W_1{s}_2\left(\mathrm{t}\right)\right|& \le & {\displaystyle \underset{\mathrm{t}\in \left[0,{\mathrm{t}}_1\right]}{sup}}\frac{\eta }{\mathrm{\Gamma }\left(\varsigma \right)}{\int }_0^{\mathrm{t}}{v}^{\eta -1}{(\mathrm{t}-v)}^{\varsigma -1}\left|\mathcal{F}(v,s_1\left(v\right))-\mathcal{F}(v,s_2\left(v\right))\right|dv.\end{array}$$

Using the transformation given in (15) and assumption (H2), we have the following:

$$\begin{array}{ccc}\hfill \left|W_1{s}_1\left(\mathrm{t}\right)-W_1{s}_2\left(\mathrm{t}\right)\right|& \le & {\displaystyle \underset{\mathrm{t}\in \left[0,{\mathrm{t}}_1\right]}{sup}}\frac{\eta }{\mathrm{\Gamma }\left(\varsigma \right)}{{\mathrm{t}}_1}^{\eta +\varsigma -1}{\int }_0^1{\rho }^{\eta -1}{(1-\rho )}^{\varsigma -1}\left|\mathcal{F}(v,s_1\left(v\right))-\mathcal{F}(v,s_2\left(v\right))\right|dv\hfill \\ & \le & \frac{\eta \mathcal{L}\beta (\eta ,\varsigma ){{\mathrm{t}}_1}^{\eta +\varsigma -1}}{\mathrm{\Gamma }\left(\varsigma \right)}∥s_1-s_2∥.\hfill \end{array}$$

For $\mathrm{t}\in \left({\mathrm{t}}_1,T\right],$ $s_1,s_2\in {\mathrm{\Lambda }}_{\zeta }$. Then, via (H2), we obtain the following:

$$\begin{array}{ccc}\hfill \left|W_1{s}_1\left(\mathrm{t}\right)-W_1{s}_2\left(\mathrm{t}\right)\right|& \le & {\displaystyle \underset{\mathrm{t}\in \left({\mathrm{t}}_1,T\right]}{sup}}\left(\left|s_1\left({\mathrm{t}}_1\right)-s_2\left({\mathrm{t}}_1\right)\right|+\left|\frac{(1-\varsigma )}{\mathbf{B}\left(\varsigma \right)}\eta {\mathrm{t}}^{\eta -1}\left|\mathcal{F}(\mathrm{t},s_1\left(\mathrm{t}\right))-\mathcal{F}(\mathrm{t},s_2\left(\mathrm{t}\right))\right|\right|\right)\hfill \\ & \le & ∥s_1-s_2∥+\frac{(1-\varsigma )}{\mathbf{B}\left(\varsigma \right)}\eta {\mathrm{t}}^{\eta -1}\mathcal{L}∥s_1-s_2∥\hfill \\ & =& \left(1-\frac{(1-\varsigma )}{\mathbf{B}\left(\varsigma \right)}\eta {\mathrm{t}}^{\eta -1}\mathcal{L}\right)∥s_1-s_2∥.\hfill \end{array}$$

Hence,

$$∥W_1{s}_1-W_1{s}_2∥\le \left(1-\frac{(1-\varsigma )}{\mathbf{B}\left(\varsigma \right)}\eta {\mathrm{t}}^{\eta -1}\mathcal{L}\right)∥s_1-s_2∥.$$

From (12), we see that $W_1$ is the contraction.

Step 3: Relative compactness of $W_2$.

Part (1): $W_2$ is continuous.

Since $\mathcal{F}(\mathrm{t},s(\mathrm{t}\left)\right)$ is continuous, then $W_2$ is continuous.

Part (2): $W_2$ is uniformly bounded on ${\mathrm{\Lambda }}_{\zeta }.$

For $\mathrm{t}\in \left[0,{\mathrm{t}}_1\right],$$s\in {\mathrm{\Lambda }}_{\zeta }$, the result can be obtained immediately.

For $\mathrm{t}\in \left({\mathrm{t}}_1,T\right],$$s\in {\mathrm{\Lambda }}_{\zeta }$, we have the following:

$$\begin{array}{ccc}\hfill \left|W_{2}s\left(\mathrm{t}\right)\right|& \le & {\displaystyle \underset{\mathrm{t}\in \left({\mathrm{t}}_1,T\right]}{sup}}\frac{\varsigma \eta }{\mathbf{B}\left(\varsigma \right)}{\int }_{{\mathrm{t}}_1}^{\mathrm{t}}{v}^{\eta -1}\left|\mathcal{F}(v,s(v\left)\right)\right|dv\hfill \\ & \le & {\displaystyle \underset{\mathrm{t}\in \left({\mathrm{t}}_1,T\right]}{sup}}\frac{\varsigma \eta }{\mathbf{B}\left(\varsigma \right)}{\int }_{{\mathrm{t}}_1}^{\mathrm{t}}{v}^{\eta -1}\left(\mathbb{k}+\left|s\left(v\right)\right|\mathbf{q}\right)dv.\hfill \end{array}$$

Hence,

$$∥W_{2}s∥\le \frac{\varsigma T^{\eta }}{\mathbf{B}\left(\varsigma \right)}\left(\mathbb{k}+\zeta \mathbf{q}\right).$$

Hence, $W_2$ is uniformly bounded on ${\mathrm{\Lambda }}_{\zeta }$.

Part (3): $W_2$ is equicontinuous. We discuss two cases as follows:

Case (1). For $s\in {\mathrm{\Lambda }}_{\zeta }$ and ${\mathrm{t}}_a,{\mathrm{t}}_b\in \left(0,{\mathrm{t}}_1\right]$ with the condition ${\mathrm{t}}_a<{\mathrm{t}}_b,$ we have the following:

$$∥W_{2}s\left({\mathrm{t}}_b\right)-W_{2}s\left({\mathrm{t}}_a\right)∥=0.$$

Case (2) For any ${\mathrm{t}}_a,{\mathrm{t}}_b\in \left({\mathrm{t}}_1,T\right],{\mathrm{t}}_a<{\mathrm{t}}_b$ and $s\in {\mathrm{\Lambda }}_{\zeta }$, we have $\left(\mathbb{k}+\left|s\left(\mathrm{t}\right)\right|\mathbf{q})\right)$

$$\begin{array}{ccc}\hfill ∥W_{2}s\left({\mathrm{t}}_b\right)-W_{2}s\left({\mathrm{t}}_a\right)∥& \le & \frac{\varsigma \eta }{\mathbf{B}\left(\varsigma \right)}{\int }_{{\mathrm{t}}_1}^{{\mathrm{t}}_b}{v}^{\eta -1}\left|\mathcal{F}(v,s(v\left)\right)\right|dv\hfill \\ & & -\frac{\varsigma \eta }{\mathbf{B}\left(\varsigma \right)}{\int }_{{\mathrm{t}}_1}^{{\mathrm{t}}_a}{v}^{\eta -1}\left|\mathcal{F}(v,s(v\left)\right)\right|dv\hfill \\ & \le & \frac{\varsigma \eta }{\mathbf{B}\left(\varsigma \right)}\left({\int }_{{\mathrm{t}}_1}^{{\mathrm{t}}_b}{v}^{\eta -1}\left|\mathcal{F}(v,s(v\left)\right)\right|dv-{\int }_{{\mathrm{t}}_1}^{{\mathrm{t}}_a}{v}^{\eta -1}\left|\mathcal{F}(v,s(v\left)\right)\right|dv\right)\hfill \\ & =& \frac{\varsigma }{\mathbf{B}\left(\varsigma \right)}({({\mathrm{t}}_b-{\mathrm{t}}_1)}^{\eta }-{({\mathrm{t}}_a-{\mathrm{t}}_1)}^{\eta })\left(\mathbb{k}+\zeta \mathbf{q}\right)\hfill \\ & \to & 0 as {\mathrm{t}}_b\to {\mathrm{t}}_a.\hfill \end{array}$$

Thus, $W_2$ is equicontinuous. In view of the Arzelá–Ascoli theorem, together with the above steps, W turns relatively compact and, hence, it is completely continuous. Therefore, Theorem 1 guarantees at least one solution for problem (9). □

Theorem 3.

Assuming that $\left(H_2\right)$, together with the condition

$$0<max\left[\frac{\eta \mathcal{L}\beta (\eta ,\varsigma ){{\mathrm{t}}_1}^{\eta +\varsigma -1}}{\mathrm{\Gamma }\left(\varsigma \right)},\left(1+\frac{\mathcal{L}}{\mathbf{B}\left(\varsigma \right)}\left(\eta \left(1-\varsigma \right){{\mathrm{t}}_1}^{\eta -1}+\varsigma (T^{\eta }-{\mathrm{t}}_1^{\eta })\right)\right)\right]<1$$

hold, then problem (9) has a unique result.

Proof. 

For $\mathrm{t}\in \left[0,{\mathrm{t}}_1\right],$ $s_1,s_2\in {\mathrm{\Lambda }}_{\zeta }$ with (H2), we have the following:

$$\begin{array}{ccc}\hfill \left|W{s}_1\left(\mathrm{t}\right)-W{s}_2\left(\mathrm{t}\right)\right|& \le & \underset{\mathrm{t}\in \left[0,{\mathrm{t}}_1\right]}{sup}\frac{\eta }{\mathrm{\Gamma }\left(\varsigma \right)}{\int }_0^{\mathrm{t}}{v}^{\eta -1}{(\mathrm{t}-v)}^{\varsigma -1}\left|\mathcal{F}(v,s_1\left(v\right))-\mathcal{F}(v,s_2\left(v\right))\right|dv.\end{array}$$

Using the transformation given in (15) and assumption (H2), we have the following:

$$\begin{array}{ccc}\hfill \left|W{s}_1\left(\mathrm{t}\right)-W{s}_2\left(\mathrm{t}\right)\right|& \le & {\displaystyle \underset{\mathrm{t}\in \left[0,{\mathrm{t}}_1\right]}{sup}}\frac{\eta }{\mathrm{\Gamma }\left(\varsigma \right)}{{\mathrm{t}}_1}^{\eta +\varsigma -1}{\int }_0^1{\rho }^{\eta -1}{(1-\rho )}^{\varsigma -1}\left|\mathcal{F}(v,s_1\left(v\right))-\mathcal{F}(v,s_2\left(v\right))\right|dv\hfill \\ & \le & \frac{\eta \mathcal{L}\beta (\eta ,\varsigma ){{\mathrm{t}}_1}^{\eta +\varsigma -1}}{\mathrm{\Gamma }\left(\varsigma \right)}∥s_1-s_2∥.\hfill \end{array}$$

Thus,

$$∥W{s}_1-W{s}_2∥\le \frac{\eta \mathcal{L}\beta (\eta ,\varsigma ){{\mathrm{t}}_1}^{\eta +\varsigma -1}}{\mathrm{\Gamma }\left(\varsigma \right)}∥s_1-s_2∥.$$

For $\mathrm{t}\in \left({\mathrm{t}}_1,T\right],$ $s_1,s_2\in {\mathrm{\Lambda }}_{\zeta }$ with (H2), we have the following:

$$\begin{array}{ccc}\hfill \left|W{s}_1\left(\mathrm{t}\right)-W{s}_2\left(\mathrm{t}\right)\right|& \le & \left|s_1\left({\mathrm{t}}_1\right)-s_2\left({\mathrm{t}}_1\right)\right|+{\displaystyle \underset{\mathrm{t}\in \left({\mathrm{t}}_1,T\right]}{sup}}[\frac{(1-\varsigma )}{\mathbf{B}\left(\varsigma \right)}\eta {\mathrm{t}}^{\eta -1}\left|\mathcal{F}(\mathrm{t},s_1\left(\mathrm{t}\right))-\mathcal{F}(\mathrm{t},s_2\left(\mathrm{t}\right))\right|\hfill \\ & +& \frac{\varsigma \eta }{\mathbf{B}\left(\varsigma \right)}{\int }_{{\mathrm{t}}_1}^{\mathrm{t}}{v}^{\eta -1}\left|\mathcal{F}(v,s_1\left(v\right))-\mathcal{F}(v,s_2\left(v\right))\right|dv]\hfill \\ & \le & \left|s_1\left(\mathrm{t}\right)-s_2\left(\mathrm{t}\right)\right|+{\displaystyle \underset{\mathrm{t}\in \left({\mathrm{t}}_1,T\right]}{sup}}[\frac{(1-\varsigma )\mathcal{L}}{\mathbf{B}\left(\varsigma \right)}\eta {\mathrm{t}}^{\eta -1}\left|s_1\left(\mathrm{t}\right)-s_2\left(\mathrm{t}\right)\right|\hfill \\ & +& \frac{\varsigma \mathcal{L}(T^{\eta }-{\mathrm{t}}_1^{\eta })}{\mathbf{B}\left(\varsigma \right)}\left|s_1\left(\mathrm{t}\right)-s_2\left(\mathrm{t}\right)\right|].\hfill \end{array}$$

Hence,

$$∥W{s}_1-W{s}_2∥\le \left(1+\frac{\mathcal{L}}{\mathbf{B}\left(\varsigma \right)})\eta \left(1-\varsigma \right){{\mathrm{t}}_1}^{\eta -1}+\varsigma (T^{\eta }-{\mathrm{t}}_1^{\eta }))\right)∥s_1-s_2∥.$$

This shows that W is a contraction. Which guarantees the unique solution of model (9). □

Hyers–Ulam (H-U) Stability

Definition 8.

Model (9) is H-U stable if there exists a real number:

$$\mathfrak{M}=max\left\{{\mathfrak{M}}_1,{\mathfrak{M}}_2,{\mathfrak{M}}_3,{\mathfrak{M}}_4,{\mathfrak{M}}_5,{\mathfrak{M}}_6\right\}>0,$$

such that for each $ϵ=max\left\{{ϵ}_1,{ϵ}_2,{ϵ}_3,{ϵ}_4,{ϵ}_5,{ϵ}_6\right\}>0,$ and for any solution $\tilde{s}\in \mathbf{J}$ of inequality, we have the following:

$$\left|{}^{PCFF}{D}^{\varsigma ,\eta }\tilde{s}\left(\mathrm{t}\right)-\mathcal{F}(\mathrm{t},\tilde{s}\left(\mathrm{t}\right))\right|\le ϵ,\mathrm{t}\in \mathbf{I},$$

there exists a unique solution $s\in \mathbf{J}$ of model (9), which satisfies the inequality given as follows:

$$∥\tilde{s}-s∥\le \mathfrak{M}ϵ,\mathrm{t}\in \mathbf{I},$$

where

$$\widehat{s}\left(\mathrm{t}\right)=\left(\begin{array}{c}\widehat{\mathbf{X}}\left(\mathrm{t}\right)\\ \widehat{\mathbf{E}}\left(\mathrm{t}\right)\\ \widehat{\mathbf{A}}\left(\mathrm{t}\right)\\ \widehat{{\mathbf{A}}_c}\left(\mathrm{t}\right)\\ \widehat{\mathbf{C}}\left(\mathrm{t}\right)\\ \widehat{{\mathbf{R}}_p}\left(\mathrm{t}\right)\end{array}\right),\widehat{s}\left(0\right)=\left(\begin{array}{c}\widehat{\mathbf{X}}\left(0\right)\\ \widehat{\mathbf{E}}\left(0\right)\\ \widehat{\mathbf{A}}\left(0\right)\\ \widehat{{\mathbf{A}}_c}\left(0\right)\\ \widehat{\mathbf{C}}\left(0\right)\\ \widehat{{\mathbf{R}}_p}\left(0\right)\end{array}\right),\mathcal{F}(\mathrm{t},\widehat{s}\left(\mathrm{t}\right))=\left(\begin{array}{c}{\mathrm{\Phi }}_1(\mathrm{t},\widehat{s}\left(\mathrm{t}\right))\\ {\mathrm{\Phi }}_2(\mathrm{t},\widehat{s}\left(\mathrm{t}\right))\\ {\mathrm{\Phi }}_3(\mathrm{t},\widehat{s}\left(\mathrm{t}\right))\\ {\mathrm{\Phi }}_4(\mathrm{t},\widehat{s}\left(\mathrm{t}\right))\\ {\mathrm{\Phi }}_5(\mathrm{t},\widehat{s}\left(\mathrm{t}\right))\\ {\mathrm{\Phi }}_6(\mathrm{t},\widehat{s}\left(\mathrm{t}\right))\end{array}\right).$$

Remark 1.

Let there exist a small perturbation $\mathrm{\Psi }\in \mathbf{J},$ such that

(i) $\left|\mathrm{\Psi }\left(\mathrm{t}\right)\right|\le ϵ,\mathrm{t}\in J;$

(ii) ${}^{PCFF}{D}^{\varsigma ,\eta }\widehat{s}\left(\mathrm{t}\right)=\mathcal{F}(\mathrm{t},\widehat{s}\left(\mathrm{t}\right))+\mathrm{\Psi }\left(\mathrm{t}\right),\mathrm{t}\in J,$ where

$$\mathrm{\Psi }\left(\mathrm{t}\right)=\left(\mathrm{t},{\mathrm{\Psi }}_1\left(\mathrm{t}\right),{\mathrm{\Psi }}_2\left(\mathrm{t}\right),{\mathrm{\Psi }}_3\left(\mathrm{t}\right),{\mathrm{\Psi }}_4\left(\mathrm{t}\right),{\mathrm{\Psi }}_5\left(\mathrm{t}\right),{\mathrm{\Psi }}_6\left(\mathrm{t}\right)\right).$$

A problem with a small perturbation function is obtained by Remark 1:

$$\left\{\begin{array}{c}{}^{PCFF}{D}^{\varsigma ,\eta }\widehat{s}\left(\mathrm{t}\right)=\mathcal{F}(\mathrm{t},\widehat{s}\left(\mathrm{t}\right))+\mathrm{\Psi }\left(\mathrm{t}\right),\\ \widehat{s}\left(0\right)={\widehat{s}}_0>0.\end{array}\right.$$

(17)

Lemma 2.

The solution of the above problem with a given small perturbation function is given by the following:

$$\widehat{s}\left(\mathrm{t}\right)=\left\{\begin{array}{c}{\widehat{s}}_0+\frac{\eta }{\mathrm{\Gamma }\left(\varsigma \right)}{\int }_0^{\mathrm{t}}{v}^{\eta -1}{(\mathrm{t}-v)}^{\varsigma -1}(\mathcal{F}(v,\widehat{s}\left(v\right))+\mathrm{\Psi }\left(v\right))dv,if \mathrm{t}\in \left[0,{\mathrm{t}}_1\right],\\ \widehat{s}\left({\mathrm{t}}_1\right)+\frac{(1-\varsigma )}{\mathbf{B}\left(\varsigma \right)}\eta {\mathrm{t}}^{\eta -1}(\mathcal{F}(v,\widehat{s}\left(v\right))+\mathrm{\Psi }\left(v\right))+\frac{\varsigma \eta }{\mathbf{B}\left(\varsigma \right)}{\int }_{{\mathrm{t}}_1}^{\mathrm{t}}{v}^{\eta -1}(\mathcal{F}(v,\widehat{s}\left(v\right))+\mathrm{\Psi }\left(v\right))dv,\\ if \mathrm{t}\in \left({\mathrm{t}}_1,T\right].\end{array}\right.$$

(18)

Proof. 

The proof is adopted from Equation (11). □

Theorem 4.

Let conditions of Theorem 3 hold. Model (9) is H-U stable with the following condition: $0<max\left[\frac{\eta \mathcal{L}\beta (\eta ,\varsigma ){{\mathrm{t}}_1}^{\eta +\varsigma -1}}{\mathrm{\Gamma }\left(\varsigma \right)},\left(1+\frac{\mathcal{L}}{\mathbf{B}\left(\varsigma \right)}\left(\eta \left(1-\varsigma \right){{\mathrm{t}}_1}^{\eta -1}+\varsigma (T^{\eta }-{\mathrm{t}}_1^{\eta })\right)\right)\right]<1.$

Proof. 

Let $ϵ=max\{{ϵ}_1,{ϵ}_2,{ϵ}_3,{ϵ}_4,{ϵ}_5,{ϵ}_6\}>0$ and $\widehat{s}\in \mathbf{J}$ be a function satisfying the following inequality:

$$\left|{}^{PCFF}{D}^{\varsigma ,\eta }\widehat{s}\left(\mathrm{t}\right)-\mathcal{F}(\mathrm{t},\widehat{s}\left(\mathrm{t}\right))\right|\le ϵ,$$

and let $s\in \mathbf{J}$ be the unique solution of model (9). For $\mathrm{t}\in \left[0,{\mathrm{t}}_1\right],$ we have the following:

$$\begin{array}{ccc}\hfill \left|\widehat{s}\left(\mathrm{t}\right)-s\left(\mathrm{t}\right)\right|& \le & {\displaystyle \underset{\mathrm{t}\in \left({\mathrm{t}}_1,T\right]}{sup}}[\frac{\eta }{\mathrm{\Gamma }\left(\varsigma \right)}{\int }_0^{\mathrm{t}}{v}^{\eta -1}{(\mathrm{t}-v)}^{\varsigma -1}\left|\mathcal{F}(v,\widehat{s}\left(v\right))-\mathcal{F}(v,s\left(v\right))\right|dv\hfill \\ & +& \frac{\eta }{\mathrm{\Gamma }\left(\varsigma \right)}{\int }_0^{\mathrm{t}}{v}^{\eta -1}{(\mathrm{t}-v)}^{\varsigma -1}\left|\mathrm{\Psi }\left(v\right))\right|dv]\hfill \\ & \le & {\displaystyle \underset{\mathrm{t}\in \left[0,{\mathrm{t}}_1\right]}{sup}}\left[\frac{\eta \mathcal{L}\beta (\eta ,\varsigma ){\mathrm{t}}^{\eta +\varsigma -1}}{\mathrm{\Gamma }\left(\varsigma \right)}\left|\widehat{s}\left(\mathrm{t}\right)-s\left(\mathrm{t}\right)\right|+\frac{\eta \beta (\eta ,\varsigma ){\mathrm{t}}^{\eta +\varsigma -1}}{\mathrm{\Gamma }\left(\varsigma \right)}ϵ\right].\hfill \end{array}$$

Hence,

$$∥\widehat{s}-s∥\le \frac{\frac{\eta \beta (\eta ,\varsigma ){{\mathrm{t}}_1}^{\eta +\varsigma -1}}{\mathrm{\Gamma }\left(\varsigma \right)}}{1-\frac{\eta \mathcal{L}\beta (\eta ,\varsigma ){{\mathrm{t}}_1}^{\eta +\varsigma -1}}{\mathrm{\Gamma }\left(\varsigma \right)}}ϵ.$$

For $\mathrm{t}\in \left({\mathrm{t}}_1,T\right],$ consider

$$\begin{array}{ccc}\hfill \left|\widehat{s}\left(\mathrm{t}\right)-s\left(\mathrm{t}\right)\right|& \le & \left|\widehat{s}\left({\mathrm{t}}_1\right)-s\left({\mathrm{t}}_1\right)\right|+\underset{\mathrm{t}\in \left({\mathrm{t}}_1,T\right]}{sup}[\frac{(1-\varsigma )}{\mathbf{B}\left(\varsigma \right)}\eta {\mathrm{t}}^{\eta -1}\left|\mathcal{F}(v,\widehat{s}\left(v\right))-\mathcal{F}(v,s\left(v\right))\right|\hfill \\ & +& \frac{(1-\varsigma )}{\mathbf{B}\left(\varsigma \right)}\eta {\mathrm{t}}^{\eta -1}\left|\mathrm{\Psi }\left(v\right)\right|\hfill \\ & +& \frac{\varsigma \eta }{\mathbf{B}\left(\varsigma \right)}{\int }_{{\mathrm{t}}_1}^{\mathrm{t}}{v}^{\eta -1}|\mathcal{F}(v,\widehat{s}\left(v\right))-\mathcal{F}(v,s\left(v\right))|+\frac{\varsigma \eta }{\mathbf{B}\left(\varsigma \right)}{\int }_{{\mathrm{t}}_1}^{\mathrm{t}}{v}^{\eta -1}\left|\mathrm{\Psi }\left(v\right)\right|dv]\hfill \\ & \le & \left|s_1\left(\mathrm{t}\right)-s_2\left(\mathrm{t}\right)\right|+\underset{\mathrm{t}\in \left({\mathrm{t}}_1,T\right]}{sup}[\frac{(1-\varsigma )\mathcal{L}}{\mathbf{B}\left(\varsigma \right)}\eta {\mathrm{t}}^{\eta -1}\left|s_1\left(\mathrm{t}\right)-s_2\left(\mathrm{t}\right)\right|\hfill \\ & +& \frac{\varsigma \mathcal{L}(T^{\eta }-{\mathrm{t}}_1^{\eta })}{\mathbf{B}\left(\varsigma \right)}\left|s_1\left(\mathrm{t}\right)-s_2\left(\mathrm{t}\right)\right|+\frac{(1-\varsigma )}{\mathbf{B}\left(\varsigma \right)}\eta {\mathrm{t}}^{\eta -1}\left|\mathrm{\Psi }\left(v\right)\right|+\frac{\varsigma \eta }{\mathbf{B}\left(\varsigma \right)}{\int }_{{\mathrm{t}}_1}^{\mathrm{t}}{v}^{\eta -1}\left|\mathrm{\Psi }\left(v\right)\right|dv]\hfill \\ & \le & \left(1+\frac{\mathcal{L}}{\mathbf{B}\left(\varsigma \right)}\left(\eta \left(1-\varsigma \right){{\mathrm{t}}_1}^{\eta -1}+\varsigma (T^{\eta }-{\mathrm{t}}_1^{\eta })\right)\right)∥\widehat{s}-s∥\hfill \\ & +& \frac{1}{\mathbf{B}\left(\varsigma \right)}\left(\eta \left(1-\varsigma \right){{\mathrm{t}}_1}^{\eta -1}+\varsigma (T^{\eta }-{\mathrm{t}}_1^{\eta })\right)ϵ.\hfill \end{array}$$

Thus,

$$∥\widehat{s}-s∥\le \frac{\frac{1}{\mathbf{B}\left(\varsigma \right)}\left(\eta \left(1-\varsigma \right){{\mathrm{t}}_1}^{\eta -1}+\varsigma (T^{\eta }-{\mathrm{t}}_1^{\eta })\right)}{1-\left(1+\frac{\mathcal{L}}{\mathbf{B}\left(\varsigma \right)}\left(\eta \left(1-\varsigma \right){{\mathrm{t}}_1}^{\eta -1}+\varsigma (T^{\eta }-{\mathrm{t}}_1^{\eta })\right)\right)}ϵ.$$

(19)

Since $0<\mathcal{L}{\mathrm{t}}_1<1,$ and $0<\frac{\mathcal{L}}{\mathbf{B}\left(\varsigma \right)}\left(\left(1-\varsigma \right)+\frac{{(T-{\mathrm{t}}_1)}^{\varsigma }}{\mathrm{\Gamma }\left(\varsigma \right)}\right)<1.$ Then, by choosing $\mathfrak{M}>0,$ such that

$$\mathfrak{M}=\left\{\begin{array}{c}\frac{\frac{\eta \beta (\eta ,\varsigma ){{\mathrm{t}}_1}^{\eta +\varsigma -1}}{\mathrm{\Gamma }\left(\varsigma \right)}}{1-\frac{\eta \mathcal{L}\beta (\eta ,\varsigma ){{\mathrm{t}}_1}^{\eta +\varsigma -1}}{\mathrm{\Gamma }\left(\varsigma \right)}},\mathrm{i}\mathrm{f} \mathrm{t}\in \left[0,{\mathrm{t}}_1\right],\\ \frac{\frac{1}{\mathbf{B}\left(\varsigma \right)}\left(\eta \left(1-\varsigma \right){{\mathrm{t}}_1}^{\eta -1}+\varsigma (T^{\eta }-{\mathrm{t}}_1^{\eta })\right)}{1-\left(1+\frac{\mathcal{L}}{\mathbf{B}\left(\varsigma \right)}\left(\eta \left(1-\varsigma \right){{\mathrm{t}}_1}^{\eta -1}+\varsigma (T^{\eta }-{\mathrm{t}}_1^{\eta })\right)\right)},\mathrm{i}\mathrm{f} \mathrm{t}\in \left({\mathrm{t}}_1,T\right],\end{array}\right.$$

from (19), we have the following:

$$∥\widehat{s}-s∥\le \mathfrak{M}ϵ.$$

Therefore, the solution of problem (9) is U-H stable. Consequently, the solution of the proposed model is U-H stable. □

6. Computational Results

From the solution of (10), we have the following integral equation at $t=t_{m+1}$ by using (11) with $\mathrm{h}={\mathrm{t}}_i-{\mathrm{t}}_{i-1}$

$$\begin{array}{l}s\left({\mathrm{t}}_{m+1}\right)=\left\{\begin{array}{l}{s}_0+\frac{\eta }{\mathrm{\Gamma }\left(\varsigma \right)}{\int }_0^{{\mathrm{t}}_{m+1}}{v}^{\eta -1}{({\mathrm{t}}_{m+1}-v)}^{\varsigma -1}\mathcal{F}(v,s\left(v\right))dv,\mathrm{i}\mathrm{f}{\mathrm{t}}_{m+1}\in \left[0,{\mathrm{t}}_1\right],\\ s\left({\mathrm{t}}_1\right)+\frac{(1-\varsigma )}{\mathbf{B}\left(\varsigma \right)}\eta {\mathrm{t}}_{m+1}^{\eta -1}\mathcal{F}({\mathrm{t}}_{m+1},s\left({\mathrm{t}}_{m+1}\right))-\frac{(1-\varsigma )}{\mathbf{B}\left(\varsigma \right)}\eta {\mathrm{t}}_{m-1}^{\eta -1}\mathcal{F}({\mathrm{t}}_{m-1},s\left({\mathrm{t}}_{m-1}\right))\\ +\frac{\varsigma \eta }{\mathbf{B}\left(\varsigma \right)}{\int }_{{\mathrm{t}}_1}^{{\mathrm{t}}_{m+1}}{v}^{\eta -1}\mathcal{F}(v,s\left(v\right))dv,\mathrm{i}\mathrm{f} \mathrm{t}\in \left({\mathrm{t}}_1,T\right].\end{array}\right.\end{array}$$

(20)

We can further write (20) as follows:

$$\begin{array}{c}\hfill s\left({\mathrm{t}}_{m+1}\right)=\left\{\begin{array}{l}{s}_0+\frac{\eta }{\mathrm{\Gamma }\left(\varsigma \right)}{\displaystyle \sum _{\ell =0}^m}{\int }_{{\mathrm{t}}_{\ell }}^{{\mathrm{t}}_{\ell +1}}{v}^{\eta -1}{({\mathrm{t}}_{\ell +1}-v)}^{\varsigma -1}\mathcal{F}(v,s\left(v\right))dv,\mathrm{i}\mathbf{f}{\mathrm{t}}_{\ell +1}\in \left[0,{\mathrm{t}}_1\right],\\ s\left({\mathrm{t}}_1\right)+\frac{(1-\varsigma )}{\mathbf{B}\left(\varsigma \right)}\eta {\mathrm{t}}_{m+1}^{\eta -1}\mathcal{F}({\mathrm{t}}_{m+1},s\left({\mathrm{t}}_{m+1}\right))-\frac{(1-\varsigma )}{\mathbf{B}\left(\varsigma \right)}\eta {\mathrm{t}}_{m-1}^{\eta -1}\mathcal{F}({\mathrm{t}}_{m-1},s\left({\mathrm{t}}_{m-1}\right))\\ +\frac{\varsigma \eta }{\mathbf{B}\left(\varsigma \right)}{\int }_{{\mathrm{t}}_1}^{{\mathrm{t}}_{m+1}}{v}^{\eta -1}\mathcal{F}(v,s\left(v\right))dv,\mathrm{i}\mathrm{f} \mathrm{t}\in \left({\mathrm{t}}_1,T\right].\end{array}\right.\end{array}$$

(21)

We write (21) on using the Lagrange interpolation as follows:

$$\begin{array}{c}s({\mathrm{t}}_{m+1})=\left\{\begin{array}{c}{s}_0+\frac{\eta h^{\varsigma }}{\mathrm{\Gamma }(\varsigma +2)}\sum _{\ell =1}^m[{\mathrm{t}}_{\ell }^{\eta -1}\mathcal{F}({\mathrm{t}}_{\ell },s({\mathrm{t}}_{\ell }))({(m+1-\ell )}^{\varsigma }(m-\ell +2+\varsigma )\hfill \\ -{(m-\ell )}^{\varsigma }(m-\ell +2+2\varsigma ))-{\mathrm{t}}_{\ell -1}^{\eta -1}\mathcal{F}(t_{\ell -1},s({\mathrm{t}}_{{\ell }_1}))({(m-\ell +1)}^{\varsigma -1}\hfill \\ -{(m-\ell )}^{\varsigma }(m-\ell +1+\varsigma ))],\mathrm{i}\mathrm{f}{\mathrm{t}}_{\ell +1}\in \left[0,{\mathrm{t}}_1\right],\hfill \\ s({\mathrm{t}}_1)+\frac{(1-\varsigma )}{\mathbf{B}(\varsigma )}\eta {\mathrm{t}}_{m+1}^{\eta -1}\mathcal{F}({\mathrm{t}}_{m+1},s({\mathrm{t}}_{m+1}))-\frac{(1-\varsigma )}{\mathbf{B}(\varsigma )}\eta {\mathrm{t}}_{m-1}^{\eta -1}\mathcal{F}({\mathrm{t}}_{m-1},s({\mathrm{t}}_{m-1}))\hfill \\ +\frac{\varsigma \eta }{\mathbf{B}(\varsigma )}[\frac{3h{\mathrm{t}}_m^{\eta -1}}{2}\mathcal{F}({\mathrm{t}}_m,s({\mathrm{t}}_m))-\frac{h{\mathrm{t}}_{m-1}^{\eta -1}}{2}\mathcal{F}({\mathrm{t}}_{m-1},s({\mathrm{t}}_{m-1})),\mathrm{i}\mathrm{f}\mathrm{t}\in \left({\mathrm{t}}_1,T\right],\hfill \end{array}\right.\end{array}$$

which, in the most compact form, can be written as follows:

$$\begin{array}{c}\hfill s({\mathrm{t}}_{m+1})=\left\{\begin{array}{c}{s}_0+\frac{\eta h^{\varsigma }}{\mathrm{\Gamma }(\varsigma +2)}\sum _{\ell =1}^m[{\mathrm{t}}_{\ell }^{\eta -1}\mathcal{F}({\mathrm{t}}_{\ell },s({\mathrm{t}}_{\ell }))({(m+1-\ell )}^{\varsigma }(m-\ell +2+\varsigma )\hfill \\ -{(m-\ell )}^{\varsigma }(m-\ell +2+2\varsigma ))-{\mathrm{t}}_{\ell -1}^{\eta -1}\mathcal{F}(t_{\ell -1},s({\mathrm{t}}_{{\ell }_1}))({(m-\ell +1)}^{\varsigma -1}\hfill \\ -{(m-\ell )}^{\varsigma }(m-\ell +1+\varsigma ))],\mathrm{i}\mathrm{f}{\mathrm{t}}_{\ell +1}\in \left[0,{\mathrm{t}}_1\right],\hfill \\ s({\mathrm{t}}_1)+\eta {\mathrm{t}}_{m+1}^{\eta -1}[\frac{(1-\varsigma )}{U(\varsigma )}+\frac{3\varsigma h}{U(\varsigma )}]\mathcal{F}({\mathrm{t}}_m,s({\mathrm{t}}_m))\hfill \\ -\eta {\mathrm{t}}_{m-1}^{\eta -1}[\frac{(1-\varsigma )}{U(\varsigma )}+\frac{\varsigma h}{U(\varsigma )}]\mathcal{F}({\mathrm{t}}_{m-1},s({\mathrm{t}}_{m-1})),\mathrm{i}\mathrm{f}\mathrm{t}\in \left({\mathrm{t}}_1,T\right].\hfill \end{array}\right.\end{array}$$

(22)

To deduce the numerical scheme for the proposed model (6), consider $s=\left(\mathbf{X},\mathbf{E},\mathbf{A},{\mathbf{A}}_c,\mathbf{C},{\mathbf{R}}_p\right)$, and we write the right side of the proposed model (6) as follows:

$$\left\{\begin{array}{c}{}^{PCFF}{D}_{\mathrm{t}}^{\varsigma ,\eta }\mathbf{X}\left(\mathrm{t}\right)={\mathrm{\Phi }}_1(\mathrm{t},s\left(\mathrm{t}\right)),\\ {}^{PCFF}{D}_{\mathrm{t}}^{\varsigma ,\eta }\mathbf{E}\left(\mathrm{t}\right)={\mathrm{\Phi }}_2(\mathrm{t},s\left(\mathrm{t}\right)),\\ {}^{PCFF}{D}_{\mathrm{t}}^{\varsigma ,\eta }\mathbf{A}\left(\mathrm{t}\right)={\mathrm{\Phi }}_3(\mathrm{t},s\left(\mathrm{t}\right)),\\ {}^{PCFF}{D}_{\mathrm{t}}^{\varsigma ,\eta }{\mathbf{A}}_c\left(\mathrm{t}\right)={\mathrm{\Phi }}_4(\mathrm{t},s\left(\mathrm{t}\right)),\\ {}^{PCFF}{D}_{\mathrm{t}}^{\varsigma ,\eta }\mathbf{C}\left(\mathrm{t}\right)={\mathrm{\Phi }}_5(\mathrm{t},s\left(\mathrm{t}\right)),\\ {}^{PCFF}{D}_{\mathrm{t}}^{\varsigma ,\eta }{\mathbf{R}}_p\left(\mathrm{t}\right)={\mathrm{\Phi }}_6(\mathrm{t},s\left(\mathrm{t}\right)).\end{array}\right.$$

(23)

Now, in view of (22), the numerical scheme for system (23) can be written as follows:

$$\mathbf{X}({\mathrm{t}}_{m+1})=\left\{\begin{array}{c}{\mathbf{X}}_0+\frac{\eta h^{\varsigma }}{\mathrm{\Gamma }(\varsigma +2)}\sum _{\ell =1}^m[{\mathrm{t}}_{\ell }^{\eta -1}{\mathrm{\Phi }}_1({\mathrm{t}}_{\ell },s({\mathrm{t}}_{\ell }))({(m+1-\ell )}^{\varsigma }(m-\ell +2+\varsigma )\hfill \\ -{(m-\ell )}^{\varsigma }(m-\ell +2+2\varsigma ))-{\mathrm{t}}_{\ell -1}^{\eta -1}{\mathrm{\Phi }}_1(t_{\ell -1},s({\mathrm{t}}_{{\ell }_1}))({(m-\ell +1)}^{\varsigma -1}\hfill \\ -{(m-\ell )}^{\varsigma }(m-\ell +1+\varsigma ))],\mathrm{i}\mathrm{f}{\mathrm{t}}_{\ell +1}\in \left[0,{\mathrm{t}}_1\right],\hfill \\ \mathbf{X}({\mathrm{t}}_1)+\eta {\mathrm{t}}_{m+1}^{\eta -1}[\frac{(1-\varsigma )}{U(\varsigma )}+\frac{3\varsigma h}{U(\varsigma )}]{\mathrm{\Phi }}_1({\mathrm{t}}_m,s({\mathrm{t}}_m))\hfill \\ -\eta {\mathrm{t}}_{m-1}^{\eta -1}[\frac{(1-\varsigma )}{U(\varsigma )}+\frac{\varsigma h}{U(\varsigma )}]{\mathrm{\Phi }}_1({\mathrm{t}}_{m-1},s({\mathrm{t}}_{m-1})),\mathrm{i}\mathrm{f}\mathrm{t}\in \left({\mathrm{t}}_1,T\right],\hfill \end{array}\right.$$

$$\mathbf{E}({\mathrm{t}}_{m+1})=\left\{\begin{array}{c}{\mathbf{E}}_0+\frac{\eta h^{\varsigma }}{\mathrm{\Gamma }(\varsigma +2)}\sum _{\ell =1}^m[{\mathrm{t}}_{\ell }^{\eta -1}{\mathrm{\Phi }}_2({\mathrm{t}}_{\ell },s({\mathrm{t}}_{\ell }))({(m+1-\ell )}^{\varsigma }(m-\ell +2+\varsigma )\hfill \\ -{(m-\ell )}^{\varsigma }(m-\ell +2+2\varsigma ))-{\mathrm{t}}_{\ell -1}^{\eta -1}{\mathrm{\Phi }}_2(t_{\ell -1},s({\mathrm{t}}_{{\ell }_1}))({(m-\ell +1)}^{\varsigma -1}\hfill \\ -{(m-\ell )}^{\varsigma }(m-\ell +1+\varsigma ))],\mathrm{i}\mathrm{f}{\mathrm{t}}_{\ell +1}\in \left[0,{\mathrm{t}}_1\right],\hfill \\ \mathbf{E}({\mathrm{t}}_1)+\eta {\mathrm{t}}_{m+1}^{\eta -1}[\frac{(1-\varsigma )}{U(\varsigma )}+\frac{3\varsigma h}{U(\varsigma )}]{\mathrm{\Phi }}_2({\mathrm{t}}_m,s({\mathrm{t}}_m))\hfill \\ -\eta {\mathrm{t}}_{m-1}^{\eta -1}[\frac{(1-\varsigma )}{U(\varsigma )}+\frac{\varsigma h}{U(\varsigma )}]{\mathrm{\Phi }}_2({\mathrm{t}}_{m-1},s({\mathrm{t}}_{m-1})),\mathrm{i}\mathrm{f}\mathrm{t}\in \left({\mathrm{t}}_1,T\right],\hfill \end{array}\right.$$

$$\mathbf{A}({\mathrm{t}}_{m+1})=\left\{\begin{array}{c}{\mathbf{A}}_0+\frac{\eta h^{\varsigma }}{\mathrm{\Gamma }(\varsigma +2)}\sum _{\ell =1}^m[{\mathrm{t}}_{\ell }^{\eta -1}{\mathrm{\Phi }}_3({\mathrm{t}}_{\ell },s({\mathrm{t}}_{\ell }))({(m+1-\ell )}^{\varsigma }(m-\ell +2+\varsigma )\hfill \\ -{(m-\ell )}^{\varsigma }(m-\ell +2+2\varsigma ))-{\mathrm{t}}_{\ell -1}^{\eta -1}{\mathrm{\Phi }}_3(t_{\ell -1},s({\mathrm{t}}_{{\ell }_1}))({(m-\ell +1)}^{\varsigma -1}\hfill \\ -{(m-\ell )}^{\varsigma }(m-\ell +1+\varsigma ))],\mathrm{i}\mathrm{f}{\mathrm{t}}_{\ell +1}\in \left[0,{\mathrm{t}}_1\right],\hfill \\ \mathbf{A}({\mathrm{t}}_1)+\eta {\mathrm{t}}_{m+1}^{\eta -1}[\frac{(1-\varsigma )}{U(\varsigma )}+\frac{3\varsigma h}{U(\varsigma )}]{\mathrm{\Phi }}_3({\mathrm{t}}_m,s({\mathrm{t}}_m))\hfill \\ -\eta {\mathrm{t}}_{m-1}^{\eta -1}[\frac{(1-\varsigma )}{U(\varsigma )}+\frac{\varsigma h}{U(\varsigma )}]{\mathrm{\Phi }}_3({\mathrm{t}}_{m-1},s({\mathrm{t}}_{m-1})),\mathrm{i}\mathrm{f}\mathrm{t}\in \left({\mathrm{t}}_1,T\right],\hfill \end{array}\right.$$

$${\mathbf{A}}_c({\mathrm{t}}_{m+1})=\left\{\begin{array}{c}{\mathbf{A}}_{c,0}+\frac{\eta h^{\varsigma }}{\mathrm{\Gamma }(\varsigma +2)}\sum _{\ell =1}^m[{\mathrm{t}}_{\ell }^{\eta -1}{\mathrm{\Phi }}_4({\mathrm{t}}_{\ell },s({\mathrm{t}}_{\ell }))({(m+1-\ell )}^{\varsigma }(m-\ell +2+\varsigma )\hfill \\ -{(m-\ell )}^{\varsigma }(m-\ell +2+2\varsigma ))-{\mathrm{t}}_{\ell -1}^{\eta -1}{\mathrm{\Phi }}_4(t_{\ell -1},s({\mathrm{t}}_{{\ell }_1}))({(m-\ell +1)}^{\varsigma -1}\hfill \\ -{(m-\ell )}^{\varsigma }(m-\ell +1+\varsigma ))],\mathrm{i}\mathrm{f}{\mathrm{t}}_{\ell +1}\in \left[0,{\mathrm{t}}_1\right],\hfill \\ {\mathbf{A}}_c({\mathrm{t}}_1)+\eta {\mathrm{t}}_{m+1}^{\eta -1}[\frac{(1-\varsigma )}{U(\varsigma )}+\frac{3\varsigma h}{U(\varsigma )}]{\mathrm{\Phi }}_4({\mathrm{t}}_m,s({\mathrm{t}}_m))\hfill \\ -\eta {\mathrm{t}}_{m-1}^{\eta -1}[\frac{(1-\varsigma )}{U(\varsigma )}+\frac{\varsigma h}{U(\varsigma )}]{\mathrm{\Phi }}_4({\mathrm{t}}_{m-1},s({\mathrm{t}}_{m-1})),\mathrm{i}\mathrm{f}\mathrm{t}\in \left({\mathrm{t}}_1,T\right],\hfill \end{array}\right.$$

$$\mathbf{C}({\mathrm{t}}_{m+1})=\left\{\begin{array}{c}{\mathbf{C}}_0+\frac{\eta h^{\varsigma }}{\mathrm{\Gamma }(\varsigma +2)}\sum _{\ell =1}^m[{\mathrm{t}}_{\ell }^{\eta -1}{\mathrm{\Phi }}_5({\mathrm{t}}_{\ell },s({\mathrm{t}}_{\ell }))({(m+1-\ell )}^{\varsigma }(m-\ell +2+\varsigma )\hfill \\ -{(m-\ell )}^{\varsigma }(m-\ell +2+2\varsigma ))-{\mathrm{t}}_{\ell -1}^{\eta -1}{\mathrm{\Phi }}_5(t_{\ell -1},s({\mathrm{t}}_{{\ell }_1}))({(m-\ell +1)}^{\varsigma -1}\hfill \\ -{(m-\ell )}^{\varsigma }(m-\ell +1+\varsigma ))],\mathrm{i}\mathrm{f}{\mathrm{t}}_{\ell +1}\in \left[0,{\mathrm{t}}_1\right],\hfill \\ \mathbf{C}({\mathrm{t}}_1)+\eta {\mathrm{t}}_{m+1}^{\eta -1}[\frac{(1-\varsigma )}{U(\varsigma )}+\frac{3\varsigma h}{U(\varsigma )}]{\mathrm{\Phi }}_5({\mathrm{t}}_m,s({\mathrm{t}}_m))\hfill \\ -\eta {\mathrm{t}}_{m-1}^{\eta -1}[\frac{(1-\varsigma )}{U(\varsigma )}+\frac{\varsigma h}{U(\varsigma )}]{\mathrm{\Phi }}_5({\mathrm{t}}_{m-1},s({\mathrm{t}}_{m-1})),\mathrm{i}\mathrm{f}\mathrm{t}\in \left({\mathrm{t}}_1,T\right],\hfill \end{array}\right.$$

and

$${\mathbf{R}}_p({\mathrm{t}}_{m+1})=\left\{\begin{array}{c}{\mathbf{R}}_{p,0}+\frac{\eta h^{\varsigma }}{\mathrm{\Gamma }(\varsigma +2)}\sum _{\ell =1}^m[{\mathrm{t}}_{\ell }^{\eta -1}{\mathrm{\Phi }}_6({\mathrm{t}}_{\ell },s({\mathrm{t}}_{\ell }))({(m+1-\ell )}^{\varsigma }(m-\ell +2+\varsigma )\hfill \\ -{(m-\ell )}^{\varsigma }(m-\ell +2+2\varsigma ))-{\mathrm{t}}_{\ell -1}^{\eta -1}{\mathrm{\Phi }}_6(t_{\ell -1},s({\mathrm{t}}_{{\ell }_1}))({(m-\ell +1)}^{\varsigma -1}\hfill \\ -{(m-\ell )}^{\varsigma }(m-\ell +1+\varsigma ))],\mathrm{i}\mathrm{f}{\mathrm{t}}_{\ell +1}\in \left[0,{\mathrm{t}}_1\right],\hfill \\ {\mathbf{R}}_p({\mathrm{t}}_1)+\eta {\mathrm{t}}_{m+1}^{\eta -1}[\frac{(1-\varsigma )}{U(\varsigma )}+\frac{3\varsigma h}{U(\varsigma )}]{\mathrm{\Phi }}_6({\mathrm{t}}_m,s({\mathrm{t}}_m))\hfill \\ -\eta {\mathrm{t}}_{m-1}^{\eta -1}[\frac{(1-\varsigma )}{U(\varsigma )}+\frac{\varsigma h}{U(\varsigma )}]{\mathrm{\Phi }}_6({\mathrm{t}}_{m-1},s({\mathrm{t}}_{m-1})),\mathrm{i}\mathrm{f}\mathrm{t}\in \left({\mathrm{t}}_1,T\right].\hfill \end{array}\right.$$

Simulations and Explanation

In view of the numerical scheme given above, we simulate our results by taking the parameter values as used by [20] in the numerical analysis: $\varrho =2,\xi =\frac{1}{67.7},\vartheta =0.042,{\omega }_1={\wp }_1=0.002,{\psi }_1=0.004,n=0.6,\mu =0.001,{\rho }_1={\eta }_1=x_1=0.02,\lambda =0.1,\delta =0.003,z1=0.01,$ and $v_1=0.2.$ Also, the initial data are taken as follows:

$$\left({\mathbf{X}}_0,{\mathbf{E}}_0,{\mathbf{A}}_0,{\mathbf{A}}_c,{\mathbf{C}}_0,{\mathbf{R}}_{p0}\right)=(60,40,3,0.25,0.1,0).$$

Consider three sets of fractal–fractional orders as ${\mathbf{X}}_1=[0,0.60],{\mathbf{X}}_2=[0.60,0.95],{\mathbf{X}}_3=[0.90,1.0]$. Different graphical presentations are given here to simulate our results. First, we consider the fractal–fractional order values in the first set ${\mathbf{X}}_1$, as presented in Figure 1, Figure 2, Figure 3, Figure 4, Figure 5 and Figure 6.

In addition, we present the results for another set of fractal–fractional orders in Figure 7, Figure 8, Figure 9, Figure 10 and Figure 11, respectively.

In addition, we graphically present the results for another set of fractal–fractional orders in Figure 12, Figure 13, Figure 14, Figure 15, Figure 16 and Figure 17, respectively.

We used various fractal–fractional order values to interpret the results graphically in Figure 1, Figure 2, Figure 3, Figure 4, Figure 5, Figure 6, Figure 7, Figure 8, Figure 9, Figure 10, Figure 11, Figure 12, Figure 13, Figure 14, Figure 15, Figure 16 and Figure 17, respectively. From the numerical presentation, we observe the crossover behaviors of each compartment near the point ${\mathrm{t}}_1=50$. The corresponding decline in the affected class is shown in Figure 1 and Figure 12 for different fractal–fractional orders. The exposed class initially grows but then declines, which is different at different values of fractal–fractional orders. The aforementioned dynamics for the exposed class are shown in Figure 2, Figure 7, and Figure 13, respectively. From Figure 3, Figure 8, and Figure 14, we observe a decline in the population with a crossover effect near the point ${\mathrm{t}}_1=50$ using different fractal–fractional order values. Also, the class of asymptomatic carriers is presented graphically using different fractal–fractional order values in Figure 4, Figure 9, and Figure 15, respectively, for numerous values of fractal–fractional orders. This class exhibits growth for the first 50 days and then shows a decline. The population of the model classes, representing exposed, asymptomatic carriers, and chronically infected individuals, grows and peaks at ${\mathrm{t}}_1=25;$ however, it begins to decline in the second sub-interval. These graphical presentations clearly demonstrate the crossover behavior in each class near the point ${\mathrm{t}}_1=25.$ The class of chronically infected individuals shows growth for the first 100 days and then starts to decline, as depicted in Figure 5, Figure 10, and Figure 16. The recovered class shows growth for the initial 100 days, as presented in Figure 6, Figure 11, and Figure 17, respectively, for different values of fractal–fractional orders.

7. Conclusions

In this research work, we conducted a detailed and comprehensive analysis of the Hepatitis B mathematical model using piecewise fractal–fractional analysis. We introduced piecewise fractal–fractional differential operators, which have not previously been considered in the literature. While the available literature contains many results on various piecewise fractional derivatives, models with piecewise fractal-fractional derivatives are rarely studied. We considered our piecewise derivative in a different pattern than the previous concept of piecewise derivatives and integrals. We divided the time interval into two subintervals. In the first one, we considered an FFD with a power law kernel, and in the second one, we considered an FFD with an exponential decay kernel. The disease-free equilibrium point of the proposed model was presented as well. In the main results, we examined the existence and stability of the proposed model. For the numerical results of the model and visual presentation, we used the Lagrange interpolation method and extended the ABM method, respectively. Upon using the said numerical scheme, we presented our results for different values of fractal–fractional orders. Also, we presented results for three different sets of fractal–fractional order values. The concerned dynamics with the crossover effect were presented graphically for each component of the proposed model. In the future, the piecewise fractal–fractional concept could be extended to other dynamical problems.