*Article* **Numerical Simulation for a Hybrid Variable-Order Multi-Vaccination COVID-19 Mathematical Model**

**Nasser Sweilam 1,\*,†, Seham M. Al-Mekhlafi 2,3,†, Reem G. Salama <sup>4</sup> and Tagreed A. Assiri <sup>5</sup>**


**Abstract:** In this paper, a hybrid variable-order mathematical model for multi-vaccination COVID-19 is analyzed. The hybrid variable-order derivative is defined as a linear combination of the variableorder integral of Riemann–Liouville and the variable-order Caputo derivative. A symmetry parameter *σ* is presented in order to be consistent with the physical model problem. The existence, uniqueness, boundedness and positivity of the proposed model are given. Moreover, the stability of the proposed model is discussed. The theta finite difference method with the discretization of the hybrid variableorder operator is developed for solving numerically the model problem. This method can be explicit or fully implicit with a large stability region depending on values of the factor Θ. The convergence and stability analysis of the proposed method are proved. Moreover, the fourth order generalized Runge–Kutta method is also used to study the proposed model. Comparative studies and numerical examples are presented. We found that the proposed model is also more general than the model in the previous study; the results obtained by the proposed method are more stable than previous research in this area.

**Keywords:** variable-order hybrid operator; Pfizer vaccine; Moderna vaccine; Janssen vaccine; theta finite difference method; generalized fourth order Runge–Kutta method

**MSC:** 65L05; 37N30; 65M06

#### **1. Introduction**

Coronaviruses are a large family of viruses known to cause illnesses ranging from the common cold to more serious illnesses such as severe acute respiratory syndrome [1]. The World Health Organization has designated this variant as a variant of serious concern. The United States Centers for Disease Control and Prevention has granted Emergency Use Authorization to the following vaccines: Pfizer-BioNTech with 95% efficacy against symptomatic COVID-19, Moderna vaccine with 94.5% efficacy and Janssen vaccine manufactured by Johnson & Johnson, which has an efficacy rating of 67%, as well as many others [1,2]. SARS-CoV-2 vaccinations have been shown to be effective against infections, including both silent and symptomatic cases, of severe COVID-19 illness and deaths [2]. Mathematical modeling is a valuable tool to study disease spread and control very effectively. Several mathematical models have been proposed in the literature to study and understand the novel complex transmission pattern of the COVID-19 pandemic; see, for example, [3–8].

In the meantime, there are now extensive articles explaining the advantage of fractional order models for studying real mathematical models in various fields [9]. The variableorder fractional derivatives (VOFDs) can describe the effects of the long variable memory of a time-dependent system. In [10], Samko et al. proposed this interesting extension of the

**Citation:** Sweilam, N.; Al-Mekhlafi, S.M.; Salama, R.G.; Assiri, T.A. Numerical Simulation for a Hybrid Variable-Order Multi-Vaccination COVID-19 Mathematical Model. *Symmetry* **2023**, *15*, 869. https:// doi.org/10.3390/sym15040869

Academic Editors: Francisco Martínez González, Mohammed K. A. Kaabar and Juan Luis García Guirao

Received: 4 February 2023 Revised: 26 March 2023 Accepted: 30 March 2023 Published: 6 April 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

classical calculation of fractions. In the concept of fractional derivative with variable order, the order may vary either as a function of the independent differentiation variable (*t*) or as a function of another (possibly spatial) variable (*x*), or both. Therefore, the derivative models described using variable-order fractional derivatives are useful and appropriate for the epidemic models. We can obtain the results of fractional order and integer order as a special case from variable-order mathematical models [11–18].

In this article, we will present the theta finite different method with the discretization of new hybrid fractional operator. This operator is called the constant proportional Caputo variable-order fractional derivative (CPC-Θ FDM) and is used to study the proposed model numerically. In the literature, the theta finite differences method (ΘFDM) method, also called the weighted average finite differences method (WAFDM), is one of the finite difference methods [19,20]. This method could be an explicit method or an implicit method (more stable and efficient), depending on the weight factor Θ ∈ [0, 1]. Using Caputo and Riesz–Feller derivatives, this method was developed for a nonstandard finite difference method [21,22].

The goal of this work is to present and analyze a hybrid variable-order fractional model of multi-vaccination for COVID-19. The new variable-order hybrid derivatives are defined as the linear combination of the variable-order Riemann–Liouville integral and the variable-order derivative of Caputo. This is one of the most effective and reliable of these operators; it is more general than the Caputo fractional operator. Positivity, boundedness and stability will be proved in the current model.

Moreover, one of the aims of this article is developing CPC-Θ FDM for solving the variable-order fractional differential equations numerically and we will compare the obtained results with the results obtained with the fourth order generalized Runge–Kutta method (GRK4M) [23] and the method in [24]. Moreover, we extended the method in [24] to variable order. The analysis of stability and convergence of the proposed method will be studied. Numerical simulations will be given to confirm the efficiency and wide applicability of the proposed method.

To our knowledge, no numerical investigations of a hybrid variable-order fractional for multi-vaccination for a COVID-19 mathematical model utilizing CPC-Θ FDM have been conducted.

This paper is organized as follows: Some notations and definitions of variable-order fractional derivatives are introduced in Section 2. In Section 3, the model with a hybrid variable order is presented; moreover, the positivity, boundedness, existence and uniqueness of the solutions and the stability of the present model are discussed. In Section 4, the numerical methods GRK4M and CPC-Θ SFDM are studied; moreover, stability analyses for these methods are proved. In Section 5, numerical simulations are presented. The conclusions are ultimately outlined in Section 6.

#### **2. Notations and Preliminaries**

In this section, we review several key definitions of variable-order calculus that will be utilized throughout the remainder of this article.

**Definition 1.** *Caputo's derivatives (right–left side variable-order fractional α*(*t*)*) are defined, respectively, as follows [25]:*

$$(^{\mathbb{C}}D\_{b-}^{a(t)}f)(\mathbf{x}) = (^{\mathbb{C}}\_{l}D\_{b}^{a(t)}f)(t) = \frac{(-1)^{n}}{\Gamma(n-a(t))} \int\_{t}^{b} \frac{f^{(n)}(\mathbf{s})}{(s-t)^{-n+a(t)+1}}d\mathbf{s}, \quad b > t,\tag{1}$$

$$(\,^\mathbb{C}D\_{a+}^{a(t)}f\rangle(t) = (\,^\mathbb{C}D\_t^{a(t)}f\rangle(t) = \frac{1}{\Gamma(n-a(t))} \int\_a^t \frac{f^{(n)}(s)}{(t-\mathfrak{f})^{-n+a(t)+1}}ds, \; \; t > a,\tag{2}$$

$$f(t) \in AC^n[a, b], \ n = 1 + [\mathfrak{R}(\alpha(t))], \ \mathfrak{R}(\alpha(t)) \notin \mathbb{N}\_0.$$

**Definition 2.** *Let* 1 > *α*(*t*) > 0, −∞ < *a* < *b* < +∞*; the right–left side variable-order fractional Riemann–Liouville's integral and f*(*t*) <sup>∈</sup> *ACn*[*a*, *<sup>b</sup>*] *are given as follows [25]:*

$$\,\_{t}I\_{b}^{a(t)}f(t) = \left[\int\_{t}^{b} f(s)(t-s)^{a(t)-1}ds\right] \frac{1}{\Gamma(a(t))}, \; t < b,\tag{3}$$

$$\,\_aI\_t^{a(t)}f(t) = \left[\int\_a^t f(s)(t-s)^{a(t)-1}ds\right] \frac{1}{\Gamma(a(t))'}, \quad t > a. \tag{4}$$

*<sup>α</sup>*(*t*) <sup>∈</sup> <sup>C</sup>.

**Definition 3** ([26])**.** *The variable-order fractional Caputo proportional operator (CP) is given as follows:*

$$\begin{split} \, \_0^\mathrm{CP} D\_t^{a(t)} y(t) &= \int\_0^t (\Gamma(1-a(t)))^{-1} (t-s)^{-a(t)} (y'(s)K\_0(s,a(t)) + y(s)K\_1(s,a(t))) ds \\ &= \left( \frac{\Gamma(1-a(t))^{-1}}{t^{a(t)}} \right) (y'(t)K\_0(t,a(t)) + y(t)K\_1(t,a(t))). \end{split} \tag{5}$$

*K*1(*α*(*t*), *t*)=(−*α*(*t*) + 1)*t <sup>α</sup>*(*t*), *K*0(*α*(*t*), *t*) = *t* (1−*α*(*t*))*α*(*t*), 1 > *α*(*t*) > 0.

Alternatively, the constant proportional Caputo (CPC) variable-order fractional hybrid operator can be formulated as follows [26]:

$$\begin{aligned} \, \_0^\text{CPC} D\_t^{a(t)} y(t) &= \left( \int\_0^t (t-s)^{-a(t)} \frac{1}{\Gamma(1-a(t))} (K\_1(a(t))y(s) + y'(s)K\_0(a(t))) ds \right) \\ &= K\_1(a(t))\_0^{RL} I\_t^{1-a(t)} y(t) + K\_0(a(t))\_0^\mathbb{C} D\_t^{a(t)} y(t), \end{aligned} \tag{6}$$

*<sup>K</sup>*0(*α*(*t*)) = *<sup>Q</sup>*(−*α*(*t*)+1)*α*(*t*), *<sup>K</sup>*1(*α*(*t*)) = *<sup>Q</sup>α*(*t*)(−*α*(*t*) + <sup>1</sup>), where *<sup>Q</sup>* is a constant.

**Definition 4.** *Moreover, its inverse operator is [26]:*

$$\,\_{0}^{\text{CPC}}I\_{t}^{a(t)}y(t) = \left(\int\_{0}^{t} \exp\left[\frac{\mathcal{K}\_{1}(a(t))}{\mathcal{K}\_{0}(a(t))}(t-s)\right]^{RL} \boldsymbol{D}\_{t}^{1-a(t)}y(s)ds\right)\frac{1}{\mathcal{K}\_{0}(a(t))}.\tag{7}$$

#### **3. A Hybrid Variable-Order Mathematical Model**

A variable-order multiple vaccination model for COVID-19 is presented below; it is an extension of the model given in [24]. To satisfy the dimensional fit between the two sides of the resulting variable-order fraction equations, the variable-order operator is modified by an auxiliary parameter *σ*. As a result, the dimension of the left side is (day)−<sup>1</sup> [27]. The following is the updated variable-order nonlinear fractional mathematical model:

$$\begin{split} \frac{1}{\sigma^{1-a(t)}t} \mathbb{C}^{\text{PC}}D\_{t}^{a(t)}S &= \Lambda - \nu\_{1}S - \nu\_{2}S - \nu\_{3}S - \lambda S - \mu S\_{\prime} \\ \frac{1}{\sigma^{1-a(t)}t} \mathbb{C}^{\text{PC}}D\_{t}^{a(t)}V\_{1} &= \nu\_{1}S - (1 - \xi\_{1})\lambda V\_{1} - \mu V\_{1}, \\ \frac{1}{\sigma^{1-a(t)}t} \mathbb{C}^{\text{PC}}D\_{t}^{a(t)}V\_{2} &= \nu\_{2}S - (1 - \xi\_{2})\lambda V\_{2} - \mu V\_{2}, \\ \frac{1}{\sigma^{1-a(t)}t} \mathbb{C}^{\text{PC}}D\_{t}^{a(t)}V\_{3} &= \nu\_{3}S - (1 - \xi\_{3})\lambda V\_{3} - \mu V\_{3}, \\ \frac{1}{\sigma^{1-a(t)}t} \mathbb{C}^{\text{PC}}D\_{t}^{a(t)}A &= f\_{3}(1 - \xi\_{3})\lambda V\_{3} + f\_{2}(1 - \xi\_{2})\lambda V\_{2} + f\_{1}(1 - \xi\_{1})\lambda V\_{1} - (\gamma\_{A} + \mu)A + p\lambda S\_{A}, \end{split}$$

1 *σ*1−*α*(*t*) *CPC* <sup>0</sup> *<sup>D</sup>α*(*t*) *<sup>t</sup> IU* =(1 − *p*)*λS* − (*γIU* + *dIU* + *α*1*μ*)*IU*, 1 *σ*1−*α*(*t*) *CPC* <sup>0</sup> *<sup>D</sup>α*(*t*) *<sup>t</sup> IV* =(1 − *f*2)(1 − *ξ*2)*λV*<sup>2</sup> + (1 − *f*3)(1 − *ξ*3)+(1 − *f*1)(1 − *ξ*1)*λV*1*λV*<sup>3</sup> − (*γIV* + (1 − *φ*)*αμ* + *dIV*)*IV*, 1 *σ*1−*α*(*t*) *CPC* <sup>0</sup> *<sup>D</sup>α*(*t*) *<sup>t</sup> IS* =*α*1(1 − *φ*)*IV* − (*dIS* + *μ* + *γIS*)*IS* + *α*<sup>1</sup> *IU*, 1 *σ*1−*α*(*t*) *CPC* <sup>0</sup> *<sup>D</sup>α*(*t*) *<sup>t</sup> R* =*γ<sup>A</sup> A* + *γIU IU* + *γIV IV* + *γIS IS* − *μR*. (8) *λ* = *βN*−<sup>1</sup> *H IU* + *θA* + *η<sup>v</sup> Iv* , *S* + *V*<sup>1</sup> + *V*<sup>2</sup> + *V*<sup>3</sup> + *A* + *IU* + *IV* + *IS* + *R* = *NH*(*t*),

with the initial conditions

$$S(0) = s\_0 \ge 0, \ V\_1(0) = v\_{10} \ge 0, \ V\_2 = v\_{20} \ge 0, \ V\_3 = v\_{30} \ge 0, \ A = a\_0 \ge 0, \ I\_{\mathcal{U}} = i\_{\mathcal{U}0} \ge 0,\tag{9}$$

$$I\_V = i\_{v0} \ge 0, \ I\_S = i\_{s0}, \ R(0) = r\_0 \ge 0.\tag{9}$$

Figure 1 shows the flowchart of the model (8). Table 1 shows the definitions of variables for system (8). The hypotheses of the model for the rate of each type of vaccination are the same as in [24], as follows:


We can verify the boundedness of the solution for the suggested model (8) as follows:

$$\frac{1}{\sigma^{1-a(t)}}\text{C}\_{0}^{\text{CPC}}D\_{t}^{a(t)}S + \stackrel{\text{CPC}}{\rightarrow}^{\text{CPC}}D\_{t}^{a(t)}R + \stackrel{\text{CPC}}{\rightarrow}^{\text{CPC}}D\_{t}^{a(t)}V\_{3} + \stackrel{\text{CPC}}{\rightarrow}^{\text{CPC}}D\_{t}^{a(t)}V\_{2} + \stackrel{\text{CPC}}{\rightarrow}^{\text{LP}} +$$

$$\frac{1}{\sigma^{1-a(t)}}\stackrel{\text{CPC}}{D}\_{t}^{a(t)}A + \stackrel{\text{CPC}}{\rightarrow}^{\text{PCC}}D\_{t}^{a(t)}S + \stackrel{\text{CPC}}{\rightarrow}^{\text{CPC}}D\_{t}^{a(t)}I\_{5} + \stackrel{\text{CPC}}{\rightarrow}^{\text{CPC}}D\_{t}^{a(t)}I\_{II} = \sigma^{-1+a(t)}\stackrel{\text{CPC}}{\rightarrow}^{\text{PCC}}D\_{t}^{a(t)}N\_{H}(t)$$

$$\frac{1}{\sigma^{1-a(t)}}\stackrel{\text{CPC}}{D}\_{t}^{a(t)}N\_{H}(t) = \Lambda - \mu N\_{H}(t) - [d\_{IV}I\_{V} + d\_{III}I\_{U} + d\_{IS}I\_{S}], \quad N\_{H}(0) = A \ge 0,\tag{10}$$

$$
\Lambda - (\mu + 3\delta)N\_H \le \frac{dN\_H}{dt} < \Lambda - \mu N\_{H\prime} \quad \delta = \min\{d\_{IV\prime}, d\_{III\prime}d\_{IS}\}.
$$

Therefore, we have *NH*(*t*) <sup>≤</sup> <sup>Λ</sup>*μ*−1, at *<sup>t</sup>* −→ <sup>∞</sup>. The feasible region <sup>Ω</sup> <sup>=</sup> {*S*, *<sup>A</sup>*, *IU*, *IV*, *IS*, *<sup>R</sup>*, *<sup>V</sup>*3, *<sup>V</sup>*1, *<sup>V</sup>*<sup>2</sup> <sup>∈</sup> <sup>R</sup>9, *NH*(*t*) <sup>≤</sup> <sup>Λ</sup>*μ*−1}.

System (8) has a solution in Ω. This verifies the boundedness of the solution.

**Table 1.** Variables of system (8).


**Figure 1.** Flowchart for system (8).

**Theorem 1.** *Using (9), for t* ≥ 0 *solutions of (8) are still nonnegative.*

**Proof.** Using (9), we obtain [28]:

1 *σ*1−*α*(*t*) *CPC* <sup>0</sup> *<sup>D</sup>α*(*t*) *<sup>t</sup> S* |*S*=<sup>0</sup> =Λ ≥ 0, 1 *σ*1−*α*(*t*) *CPC* <sup>0</sup> *<sup>D</sup>α*(*t*) *<sup>t</sup> V*<sup>1</sup> |*V*1=<sup>0</sup> =*V*1*S* ≥ 0, 1 *σ*1−*α*(*t*) *CPC* <sup>0</sup> *<sup>D</sup>α*(*t*) *<sup>t</sup> V*<sup>2</sup> |*V*2=<sup>0</sup> =*V*2*S* ≥ 0, 1 *σ*1−*α*(*t*) *CPC* <sup>0</sup> *<sup>D</sup>α*(*t*) *<sup>t</sup> V*<sup>3</sup> |*V*3=<sup>0</sup> =*V*3*S* ≥ 0, 1 *σ*1−*α*(*t*) *CPC* <sup>0</sup> *<sup>D</sup>α*(*t*) *<sup>t</sup> A* |*A*=<sup>0</sup> =(1 − *ξ*2)*f*2*λV*<sup>2</sup> + (1 − *ξ*3)*f*3*λV*<sup>3</sup> + *pλS* + (1 − *ξ*1)*f*1*λV*<sup>1</sup> ≥ 0, 1 *σ*1−*α*(*t*) *CPC* <sup>0</sup> *<sup>D</sup>α*(*t*) *<sup>t</sup> IU* |*IU*=<sup>0</sup> =(1 − *p*)*λS* ≥ 0, 1 *σ*1−*α*(*t*) *CPC* <sup>0</sup> *<sup>D</sup>α*(*t*) *<sup>t</sup> IV* |*IV*=<sup>0</sup> =(1 − *ξ*2)*λ*(1 − *f*2)*V*<sup>2</sup> + (1 − *ξ*3)*λV*3(1 − *f*3)+(1 − *ξ*1)*λV*1(1 − *f*1) ≥ 0, 1 *σ*1−*α*(*t*) *CPC* <sup>0</sup> *<sup>D</sup>α*(*t*) *<sup>t</sup> IS* |*IS*=<sup>0</sup> =*α*<sup>1</sup> *IU* + (1 − *φ*)*α*<sup>1</sup> *IV* ≥ 0, 1 *σ*−*α*(*t*)+<sup>1</sup> *CPC* <sup>0</sup> *<sup>D</sup>α*(*t*) *<sup>t</sup> R* |*R*=<sup>0</sup> =*γ<sup>A</sup> A* + *γIU IU* + *γIV IV* + *γIS IS* ≥ 0. (11)

#### *3.1. Uniqueness and Existence*

The existence and uniqueness of the solutions of the proposed model will be established using Banach fixed point theorem. Let system (8) be written as follows [4]:

$$\,\_0^{\text{CPC}}D\_t^{a(t)}\varepsilon(t) = \alpha(\varepsilon(t), t), \; \varepsilon(0) = \varepsilon\_0 \ge 0,\tag{12}$$

*<sup>ε</sup>*(*t*) = *S*, *A*, *IU*, *IV*, *IS*, *R*, *V*3, *V*1, *V*<sup>2</sup> *<sup>T</sup>* represents the variables of the proposed system (8)

and is a vector that represents the equations in the right of the system (8).

⎛ ⎜⎜⎜⎜⎜⎜⎜⎜⎝ 1 <sup>2</sup> <sup>3</sup> 4 <sup>5</sup> <sup>6</sup> <sup>7</sup> <sup>8</sup> <sup>9</sup> ⎞ ⎟⎟⎟⎟⎟⎟⎟⎟⎠ = ⎛ ⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝ *<sup>σ</sup>*1−*α*(*t*)(<sup>Λ</sup> <sup>−</sup> *<sup>ν</sup>*1*<sup>S</sup>* <sup>−</sup> *<sup>ν</sup>*2*<sup>S</sup>* <sup>−</sup> *<sup>ν</sup>*3*<sup>S</sup>* <sup>−</sup> *<sup>λ</sup><sup>S</sup>* <sup>−</sup> *<sup>μ</sup>S*) *<sup>σ</sup>*1−*α*(*t*)(*ν*1*<sup>S</sup>* <sup>−</sup> (<sup>1</sup> <sup>−</sup> *<sup>ξ</sup>*1)*λV*<sup>1</sup> <sup>−</sup> *<sup>μ</sup>V*1) *<sup>σ</sup>*1−*α*(*t*)(*ν*1*<sup>S</sup>* <sup>−</sup> (<sup>1</sup> <sup>−</sup> *<sup>ξ</sup>*2)*λV*<sup>2</sup> <sup>−</sup> *<sup>μ</sup>V*2) *<sup>σ</sup>*1−*α*(*t*)(*ν*1*<sup>S</sup>* <sup>−</sup> (<sup>1</sup> <sup>−</sup> *<sup>ξ</sup>*3)*λV*<sup>3</sup> <sup>−</sup> *<sup>μ</sup>V*3) *<sup>σ</sup>*1−*α*(*t*)((<sup>1</sup> <sup>−</sup> *<sup>ξ</sup>*2)*<sup>λ</sup> <sup>f</sup>*2*V*<sup>2</sup> + (<sup>1</sup> <sup>−</sup> *<sup>ξ</sup>*1)*<sup>λ</sup> <sup>f</sup>*1*V*<sup>1</sup> + (<sup>1</sup> <sup>−</sup> *<sup>ξ</sup>*3)*<sup>λ</sup> <sup>f</sup>*3*V*<sup>3</sup> <sup>−</sup> (*γ<sup>A</sup>* <sup>+</sup> *<sup>μ</sup>*)*A*) + *<sup>p</sup>λ<sup>S</sup> <sup>σ</sup>*1−*α*(*t*)((<sup>1</sup> <sup>−</sup> *<sup>p</sup>*)*λ<sup>S</sup>* <sup>−</sup> (*γIU* <sup>+</sup> *dIU* <sup>+</sup> *<sup>α</sup>*1*μ*)*IU* ) *<sup>σ</sup>*1−*α*(*t*)((<sup>1</sup> <sup>−</sup> *<sup>ξ</sup>*2)*λ*(<sup>1</sup> <sup>−</sup> *<sup>f</sup>*2)*V*<sup>2</sup> + (<sup>1</sup> <sup>−</sup> *<sup>ξ</sup>*3)*λ*(<sup>1</sup> <sup>−</sup> *<sup>f</sup>*3)*V*<sup>3</sup> + (<sup>1</sup> <sup>−</sup> *<sup>ξ</sup>*1)*λ*(<sup>1</sup> <sup>−</sup> *<sup>f</sup>*1)*V*<sup>1</sup> <sup>−</sup> (*dIV* <sup>+</sup> *<sup>α</sup>*(<sup>1</sup> <sup>−</sup> *<sup>φ</sup>*)*μ*)*IV* <sup>+</sup> *<sup>γ</sup>IV* ) *<sup>σ</sup>*1−*α*(*t*)(*α*<sup>1</sup> *IU* <sup>−</sup> (*dIS* <sup>+</sup> *<sup>μ</sup>* <sup>+</sup> *<sup>γ</sup>IS*)*IS*) + *<sup>α</sup>*1(<sup>1</sup> <sup>−</sup> *<sup>φ</sup>*)*IV <sup>σ</sup>*1−*α*(*t*)(*γIU IU* <sup>+</sup> *<sup>γ</sup>IV IV* <sup>+</sup> *<sup>γ</sup>IS IS* <sup>−</sup> *<sup>μ</sup><sup>R</sup>* <sup>+</sup> *<sup>γ</sup><sup>A</sup> <sup>A</sup>*) ⎞ ⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠ ,

with an initial condition *ε*0. Furthermore, Lipschitz requirements as in [4] are satisfied:

$$\left\|\left|\mathcal{O}\left(\varepsilon\_{1}(t),t\right)-\mathcal{O}\left(\varepsilon\_{2}(t),t\right)\right|\right\| \leq \mathcal{W}^{0} \left\|\varepsilon\_{1}(t)-\varepsilon\_{2}(t)\right\|, \ \mathcal{W}^{0} \in \mathbb{R}.\tag{13}$$

**Theorem 2.** *If the following conditions are met:*

$$\frac{W^0 F\_{\text{max}}^{a(t)} X\_{\text{max}}^{a(t)}}{\Gamma(a(t) - 1) K\_0(a(t))} < 1,\tag{14}$$

*the hybrid variable-order fractional model (8) has a unique solution.*

**Proof.** Applying (6) in (12), we have:

$$\varepsilon(t) = \varepsilon(t\_0) + \frac{1}{K\_0(a(t))} \int\_0^t \exp(-\frac{K\_1(a(t))}{K\_0(a(t))} (t-s))\_0^{RL} D\_t^{1-a(t)} \mathcal{O}(\varepsilon(s), s) ds. \tag{15}$$

$$\begin{aligned} \text{Let } B: \mathbb{C}(\mathbb{K}, \mathbb{R}^9) &\longrightarrow \mathbb{C}(\mathbb{K}, \mathbb{R}^9) \text{ and } \newline K = (0, T); \text{ then:}\\ B[\varepsilon(t)] &= \varepsilon(t\_0) + \frac{1}{K\_0(\boldsymbol{a}(t))} \int\_0^t \exp(-\frac{K\_1(\boldsymbol{a}(t))}{K\_0(\boldsymbol{a}(t))} (t - s))\_0^{\text{RL}} D\_t^{1 - \boldsymbol{a}(t)} \boldsymbol{\mathcal{o}}(\varepsilon(s), s) ds. \end{aligned} \tag{16}$$

We have:

$$B[\varepsilon(t)] = \varepsilon(t).$$

The supremum norm on *K* is represented by .*K*. Thus *ε*(*t*)*<sup>K</sup>* = sup *t*∈*K ε*(*t*), *<sup>ε</sup>*(*t*) <sup>∈</sup> *<sup>C</sup>*(*K*, <sup>R</sup>9).

So, .*<sup>K</sup>* with *<sup>C</sup>*(*K*, <sup>R</sup>9) is a Banach space. Then, the following relation holds:

$$\|\Lambda\| \|\varrho(s,t)\|\_{K} \|\varepsilon(s)\|\_{K} \geq \|\int\_{0}^{t} \varrho(s,t)\varepsilon(s)ds\|, \quad 0 < t < \Lambda < \infty$$

with *<sup>ϕ</sup>*(*s*, *<sup>t</sup>*) <sup>∈</sup> *<sup>C</sup>*(*K*2, <sup>R</sup>9) *<sup>ε</sup>*(*t*) <sup>∈</sup> *<sup>C</sup>*(*K*, <sup>R</sup>9), then sup*t*,*s*∈*<sup>K</sup>* <sup>|</sup>*ϕ*(*s*, *<sup>t</sup>*)<sup>|</sup> <sup>=</sup> *ϕ*(*s*, *<sup>t</sup>*)*K*. Relation (16) can be written as:

*B*[*ε*1(*t*)] <sup>−</sup> *<sup>B</sup>*[*ε*2(*t*)]*<sup>K</sup>* ≤ <sup>1</sup> *<sup>K</sup>*0(*α*(*t*)) *<sup>t</sup>* 0 *exp*(−*K*1(*α*(*t*)) *<sup>K</sup>*0(*α*(*t*))(*<sup>t</sup>* <sup>−</sup> *<sup>s</sup>*))(*RL* <sup>0</sup> *<sup>D</sup>*1−*α*(*t*) *<sup>t</sup>* (*ε*1(*s*),*s*) <sup>−</sup>*RL* <sup>0</sup> *<sup>D</sup>*1−*α*(*t*) *<sup>t</sup>* (*ε*2(*s*),*s*))*dsK*. <sup>≤</sup> *α*(*t*) *max K*0(*α*(*t*))Γ(*α*(*t*) − 1) *t* 0 (*<sup>t</sup>* <sup>−</sup> *<sup>s</sup>*)*α*(*t*)−2((*ε*1(*s*),*s*) <sup>−</sup> (*ε*2(*s*),*s*))*dsK*, <sup>≤</sup> *<sup>α</sup>*(*t*) *maxX<sup>α</sup>*(*t*) *max K*0(*α*(*t*))Γ(*α*(*t*) − 1) (*ε*1(*t*), *t*) − (*ε*2(*t*), *t*)*K*, ≤ *W*0*<sup>α</sup>*(*t*) *maxX<sup>α</sup>*(*t*) *max K*0(*α*(*t*))Γ(*α*(*t*) − 1) *ε*1(*t*) − *ε*2(*t*)*K*. (17)

Then

$$\|\|B[\varepsilon\_1(t)] - B[\varepsilon\_2(t)]\|\_{K} \le L \|\|\varepsilon\_1(t) - \varepsilon\_2(t)\|\_{K} \tag{18}$$

.

where

$$L = \frac{W^0 F\_{\max}^{\alpha(t)} X\_{\max}^{\alpha(t)}}{K\_0(\alpha(t)) \Gamma(\alpha(t) - 1)}$$

*B* is a contraction operator if 1 > *L*. So (8) has a unique solution.

#### **Table 2.** The definition of all parameters of system (8).


#### *3.2. Local Stability*

The basic reproduction number is calculated in this section. The next generation operator method is used to investigate the local stability of the disease-free equilibrium (DFE), which is given by solving <sup>1</sup> *σ*1−*α*(*t*) *CPC* <sup>0</sup> *<sup>D</sup>α*(*t*) *<sup>t</sup>* (.) = 0 of model (8) and considering *IU* = *IV* = *IS* = 0. Then, we obtained *D*0, where *D*<sup>0</sup> is the DFE and is given by [33]:

$$D\_0 = (\tilde{S}\_\nu \tilde{V}\_3 \tilde{V}\_2 \tilde{V}\_1 \tilde{A}\_\nu \tilde{I}\_V \tilde{I}\_V \tilde{I}\_\nu \tilde{I}\_{\mathrm{LI}} \tilde{R}) = $$

$$(\frac{\Lambda}{(\nu\_3 + \nu\_2 + \nu\_1 + \mu)}, \frac{\nu\_3 \Lambda}{(\nu\_3 + \nu\_2 + \nu\_1 + \mu)}, \frac{\nu\_2 \Lambda}{(\nu\_3 + \nu\_2 + \nu\_1 + \mu)}, \frac{\nu\_1 \Lambda}{(\nu\_1 + \nu\_2 + \nu\_3 + \mu)})$$

$$\frac{\Lambda}{(\nu\_1 + \nu\_2 \Lambda + \nu\_3 + \mu)}, 0, 0, 0, 0).$$

,

As a result, the matrix *V* of the transfer of individuals between compartments and the matrix *F* of new infection terms are provided by

$$F = \sigma^{1-\alpha(t)} \begin{pmatrix} \frac{\beta\theta\bar{Q}}{\bar{N}\_H} & \frac{\beta\bar{Q}}{\bar{N}\_H} & \frac{\beta\eta\_V\bar{Q}}{\bar{N}\_H} & 0\\ (1-p)\frac{\beta\theta\bar{S}}{\bar{N}\_H} & (1-p)\frac{\beta\bar{S}}{\bar{N}\_H} & (1-p)\frac{\beta\eta\_V\bar{S}}{\bar{N}\_H} & 0\\ \frac{\beta\theta\bar{\theta}}{\bar{N}\_H} & \frac{\beta\bar{\theta}}{\bar{N}\_H} & \frac{\beta\eta\_V\bar{\theta}}{\bar{N}\_H} & 0\\ 0 & 0 & 0 & 0 \end{pmatrix},$$
 
$$\text{with } \tilde{v} = (1-\tilde{\xi}\_3)(1-f\_3)\tilde{V}\_3 + (1-\tilde{\xi}\_2)(1-f\_2)\tilde{V}\_2 + (1-\tilde{\xi}\_1)(1-f\_1)\tilde{V}\_1.$$

with *<sup>ν</sup>*˜ = (<sup>1</sup> <sup>−</sup> *<sup>ξ</sup>*3)(<sup>1</sup> <sup>−</sup> *<sup>f</sup>*3)*V*˜ <sup>3</sup> + (<sup>1</sup> <sup>−</sup> *<sup>ξ</sup>*2)(<sup>1</sup> <sup>−</sup> *<sup>f</sup>*2)*V*˜ <sup>2</sup> + (<sup>1</sup> <sup>−</sup> *<sup>ξ</sup>*1)(<sup>1</sup> <sup>−</sup> *<sup>f</sup>*1)*V*˜ 1, *<sup>Q</sup>*˜ = (<sup>1</sup> <sup>−</sup> *<sup>ξ</sup>*3)*f*3*V*˜ <sup>3</sup> + (<sup>1</sup> <sup>−</sup> *<sup>ξ</sup>*1)*f*1*V*˜ <sup>1</sup> + (<sup>1</sup> <sup>−</sup> *<sup>ξ</sup>*2)*f*2*V*˜ <sup>2</sup> + *pS*˜.

$$V = \sigma^{1-a(t)} \begin{pmatrix} \mu + \gamma\_A & 0 & 0 & 0 \\ 0 & \gamma\_{IU} + d\_{IU} + a\_1 + \mu & 0 & 0 \\ 0 & 0 & \gamma\_{IV} + d\_{IV} + (1 - \phi)a\_1 + \mu & 0 \\ 0 & -a\_1 & -(1 - \phi)a\_1 & \gamma\_{IU} + d\_{IU} + \mu \\ \tau u & \ddots \ddots \ddots & \ddots & \ddots & 1 & 1 & 1 \\ \tau u & \ddots & \ddots & \ddots & \ddots & 1 & 1 & 1 \end{pmatrix} \cdot \begin{pmatrix} 1 \\ 0 \\ 0 \\ \tau \end{pmatrix}$$

The model's basic reproduction number, denoted by *R*0, is given by [34,35]:

$$\rho(FV^{-1}) = R\_0 = \sigma^{1-a(t)} \beta \left( \frac{(1-p)E\_1 E\_3 \mu + E\_1 E\_2 \eta\_V Y\_1 + E\_2 E\_3 \eta\_A \theta Y\_2}{\mu(\nu\_1 + \nu\_2 + \nu\_3)E\_1 E\_2 E\_3} \right). \tag{19}$$

with *E*<sup>1</sup> = (*γ<sup>A</sup>* + *μ*), *E*<sup>2</sup> = (*γIU* + *dIU* + *α*<sup>1</sup> + *μ*), *E*<sup>3</sup> = (*α*1(1 − *φ*) + *dIV* + *μ* + *γIV*), *Y*<sup>1</sup> = (1 − *ξ*3)(1 − *f*3)*ν*<sup>3</sup> + (1 − *ξ*1)(1 − *f*1)*ν*<sup>1</sup> + (1 − *ξ*2)(1 − *f*2)*ν*2, *Y*<sup>2</sup> = (1 − *ξ*3)*f*3*ν*<sup>3</sup> + (1 − *ξ*2)*f*2*ν*<sup>2</sup> + *μp*(1 − *ξ*1)*f*1*ν*1.

**Theorem 3.** *The disease-free equilibrium point D*<sup>0</sup> *of model* (8) *is locally asymptotically stable (LAS) if R*<sup>0</sup> < 1 *and unstable if R*<sup>0</sup> > 1*.*

**Proof.** The Jacobian matrix of the system (8) at the DFE is used to investigate the local stability of model (8) [33,36].

$$f(D\_0) = \sigma^{1-a(t)} \begin{pmatrix} X & 0 & 0 & 0 & A\_1 & A\_2 & A\_3 & 0 & 0 \\ \nu\_1 & -\mu & 0 & 0 & B\_1 & B\_2 & B\_3 & 0 & 0 \\ \nu\_2 & 0 & -\mu & 0 & F\_1 & F\_2 & F\_3 & 0 & 0 \\ \nu\_3 & 0 & 0 & -\mu & G\_1 & G\_2 & G\_3 & 0 & 0 \\ 0 & 0 & 0 & 0 & M\_1 & M\_2 & M\_3 & 0 & 0 \\ 0 & 0 & 0 & 0 & N\_1 & N\_2 & N\_3 & 0 & 0 \\ 0 & 0 & 0 & 0 & Z\_1 & Z\_2 & Z\_3 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & \alpha\_1 & (1-\phi)\alpha\_1 & -E\_4 & 0 \\ 0 & 0 & 0 & 0 & \gamma\_A & \gamma\_{III} & \gamma\_{IV} & \gamma\_{IS} & -\mu \end{pmatrix}.$$

where *<sup>X</sup>* <sup>=</sup> <sup>−</sup>(*ν*<sup>1</sup> <sup>+</sup> *<sup>ν</sup>*<sup>2</sup> <sup>+</sup> *<sup>ν</sup>*<sup>3</sup> <sup>+</sup> *<sup>μ</sup>*), *<sup>A</sup>*<sup>1</sup> <sup>=</sup> <sup>−</sup>*βθS*˜ *H N*˜ *<sup>H</sup>* , *<sup>A</sup>*<sup>2</sup> <sup>=</sup> <sup>−</sup>*βS*˜ *H N*˜ *<sup>H</sup>* , *<sup>A</sup>*<sup>3</sup> <sup>=</sup> <sup>−</sup>*βη<sup>V</sup> <sup>S</sup>*˜ *H N*˜ *<sup>H</sup>* , *<sup>B</sup>*<sup>1</sup> <sup>=</sup> <sup>−</sup>(<sup>1</sup> <sup>−</sup> *<sup>ξ</sup>*1) *βθV*˜ 1 *N*˜ *<sup>H</sup>* , *<sup>B</sup>*<sup>2</sup> <sup>=</sup> <sup>−</sup>(<sup>1</sup> <sup>−</sup> *<sup>ξ</sup>*1) *<sup>β</sup>V*˜ 1 *N*˜ *<sup>H</sup>* , *<sup>B</sup>*<sup>3</sup> <sup>=</sup> <sup>−</sup>(<sup>1</sup> <sup>−</sup> *<sup>ξ</sup>*1) *βηVV*˜ 1 *N*˜ *<sup>H</sup>* , *<sup>F</sup>*<sup>1</sup> <sup>=</sup> <sup>−</sup>(<sup>1</sup> <sup>−</sup> *<sup>ξ</sup>*2) *βθV*˜ 2 *N*˜ *<sup>H</sup>* , *<sup>F</sup>*<sup>2</sup> <sup>=</sup> <sup>−</sup>(<sup>1</sup> <sup>−</sup> *<sup>ξ</sup>*2) *<sup>β</sup>V*˜ 2 *N*˜ *<sup>H</sup>* , *<sup>F</sup>*<sup>3</sup> <sup>=</sup> <sup>−</sup>(<sup>1</sup> <sup>−</sup> *<sup>ξ</sup>*2) *βηVV*˜ 2 *N*˜ *<sup>H</sup>* , *<sup>G</sup>*<sup>1</sup> <sup>=</sup> <sup>−</sup>(<sup>1</sup> <sup>−</sup> *<sup>ξ</sup>*3) *βθV*˜ 3 *N*˜ *<sup>H</sup>* , *<sup>G</sup>*<sup>2</sup> <sup>=</sup> <sup>−</sup>(<sup>1</sup> <sup>−</sup> *<sup>ξ</sup>*3) *<sup>β</sup>V*˜ 3 *N*˜ *<sup>H</sup>* , *<sup>G</sup>*<sup>3</sup> <sup>=</sup> <sup>−</sup>(<sup>1</sup> <sup>−</sup> *<sup>ξ</sup>*3) *βηVV*˜ 3 *N*˜ *<sup>H</sup>* , *<sup>M</sup>*<sup>1</sup> <sup>=</sup> *βθQ*˜ *N*˜ *<sup>H</sup>* <sup>−</sup> *<sup>E</sup>*1, *<sup>M</sup>*<sup>2</sup> <sup>=</sup> *<sup>β</sup>Q*˜ *N*˜ *<sup>H</sup>* , *<sup>M</sup>*<sup>3</sup> <sup>=</sup> *βηVQ*˜ *N*˜ *<sup>H</sup>* , *<sup>N</sup>*<sup>1</sup> <sup>=</sup> *βθ*(1−*p*)*S*˜ *H N*˜ *<sup>H</sup>* , *<sup>N</sup>*<sup>2</sup> <sup>=</sup> *<sup>β</sup>*(1−*p*)*S*˜ *H N*˜ *<sup>H</sup>* <sup>−</sup> *<sup>E</sup>*2, *<sup>N</sup>*<sup>3</sup> <sup>=</sup> *<sup>β</sup>*(1−*p*)*η<sup>V</sup> <sup>S</sup>*˜ *H N*˜ *<sup>H</sup>* , *<sup>Z</sup>*<sup>1</sup> = *βθν*˜ *N*˜ *<sup>H</sup>* , *<sup>Z</sup>*<sup>2</sup> = *βν*˜ *N*˜ *<sup>H</sup>* , *<sup>Z</sup>*<sup>3</sup> = *βη<sup>V</sup> <sup>ν</sup>*˜ *N*˜ *<sup>H</sup>* − *E*3, *E*<sup>4</sup> = (*γIS* + *dIS* + *μ*). The characteristic equation: (*ν*<sup>3</sup> <sup>+</sup> *<sup>ν</sup>*<sup>1</sup> <sup>+</sup> *<sup>ν</sup>*<sup>2</sup> <sup>+</sup> *<sup>μ</sup>* <sup>+</sup> *<sup>λ</sup>*)(*λ*<sup>3</sup> + (*E*<sup>1</sup> <sup>+</sup> *<sup>E</sup>*<sup>2</sup> <sup>+</sup> *<sup>E</sup>*<sup>3</sup> <sup>−</sup> (1−*p*)*S*˜+*η<sup>V</sup> <sup>ν</sup>*˜+*η<sup>V</sup> <sup>θ</sup>Q*˜ *N*˜ *H β*)*λ*<sup>2</sup> + (*E*1*E*<sup>2</sup> <sup>+</sup> *<sup>E</sup>*1*E*<sup>3</sup> <sup>+</sup> *<sup>E</sup>*2*E*<sup>3</sup> <sup>−</sup> *<sup>β</sup>*[(*E*<sup>1</sup> <sup>+</sup> *<sup>E</sup>*3)(<sup>1</sup> <sup>−</sup> *<sup>p</sup>*)*S*˜ + (*E*<sup>1</sup> <sup>+</sup> *<sup>E</sup>*2)*ηVν*˜ + (*E*<sup>2</sup> <sup>+</sup> *<sup>E</sup>*3)*ηAθQ*˜])*<sup>λ</sup>* <sup>+</sup> *<sup>E</sup>*1*E*2*E*3(<sup>1</sup> <sup>−</sup> *<sup>R</sup>*0))(*<sup>μ</sup>* <sup>+</sup> *<sup>λ</sup>*)4(*<sup>λ</sup>* <sup>+</sup> *<sup>E</sup>*4) = 0. Then, we have (*λ* + *μ*) = 0,(*λ* + *E*4) = 0,(*λ* + *ν*<sup>1</sup> + *ν*<sup>2</sup> + *ν*<sup>3</sup> + *μ*) = 0; the arguments are *arg*(*λk*) > *<sup>π</sup> <sup>a</sup>* <sup>&</sup>gt; *<sup>k</sup>* <sup>2</sup>*<sup>π</sup> <sup>a</sup>* <sup>&</sup>gt; *<sup>π</sup> <sup>M</sup>* <sup>&</sup>gt; *<sup>π</sup>* <sup>2</sup>*<sup>M</sup>* , where *k* = 0, 1, 2, 3, ..., *a* − 1.

$$\begin{aligned} (\lambda^3 + (E\_1 + E\_2 + E\_3 - \beta \frac{(1 - p)\mathcal{S} + \eta\_V \vartheta + \eta\_V \vartheta \mathcal{Q}}{N})\lambda^2 + (E\_1 E\_2 + E\_1 E\_3 + E\_2 E\_3 - \\ \beta[(E\_1 + E\_3)(1 - p)\mathcal{S} + (E\_1 + E\_2)\eta\_V \mathcal{v} + (E\_2 + E\_3)\eta\_A \theta \mathcal{Q}])\lambda + E\_1 E\_2 E\_3 (1 - R\_0)) = 0. \end{aligned}$$

We can rewrite the above equation as:

$$
\lambda^3 + a\lambda^2 + b\lambda + c = 0,\tag{20}
$$

where

$$a = (E\_1 + E\_2 + E\_3 - \beta \frac{(1 - p)\tilde{S} + \eta\_V \tilde{v} + \eta\_V \theta \tilde{Q}}{\tilde{N}})\_\prime$$

$$b = (E\_1 E\_2 + E\_1 E\_3 + E\_2 E\_3 - \beta [(E\_1 + E\_3)(1 - p)\tilde{S} + (E\_1 + E\_2)\eta\_V \tilde{v} + (E\_2 + E\_3)\eta\_A \theta \tilde{Q}])\_\prime$$

$$c = E\_1 E\_2 E\_3 (1 - R\_0).$$

$$
\lambda^3 + a\lambda^2 + b\lambda + c = 0,\tag{21}
$$

We obtain

$$
\lambda^3 + a\lambda^2 + b\lambda + c = (\lambda - \zeta\_{11})(\lambda^2 - \tau\lambda + \zeta\_{11})\,\tag{22}
$$

$$
\pi = -(a + \zeta\_{11})\_\prime \tag{23}
$$

$$\mathcal{J}\_{11} = b + \mathcal{J}\_{11}(a + \mathcal{J}\_{11})\_\prime \tag{24}$$

$$
\mathcal{L} = -\mathcal{J}\_{11}\delta\_{11\prime} \tag{25}
$$

Hence, the other two roots are given by

$$
\zeta\_{11,2,3} = \frac{1}{2} (\tau \pm \sqrt{\triangle}) , \tag{26}
$$

$$
\triangle = \pi^2 - 4\delta\_{11} = a^2 - 2a\zeta\_{11} - (3\zeta\_{11}^2 + 4b). \tag{27}
$$

These two roots are complex conjugate when < 0, real and distinct when > 0, and real and conincident when = 0.

Considering that = 0 occurs *a* = *ζ*<sup>11</sup> ± 2 *ζ*2 <sup>11</sup> + *<sup>b</sup>*, we have that if *<sup>ζ</sup>*<sup>2</sup> <sup>11</sup> + *b* < 0, then > 0 and two distinct real roots given by

$$
\zeta\_{11,2,3} = \frac{1}{2} (\pi \pm \sqrt{\triangle}).
$$

If *ζ*<sup>2</sup> <sup>11</sup> <sup>+</sup> *<sup>b</sup>* <sup>=</sup> 0 then = (*<sup>a</sup>* <sup>−</sup> *<sup>ζ</sup>*11)<sup>2</sup> and two distinct real roots exist given by *<sup>ζ</sup>*11,2,3 <sup>=</sup> <sup>1</sup> 2 (*τ* ± |*a* − *ζ*11|).

So that *<sup>ζ</sup>*11,2 <sup>=</sup> <sup>−</sup>*ζ*11,1 and *<sup>ζ</sup>*11,3 <sup>=</sup> <sup>−</sup>*a*, if *<sup>ζ</sup>*<sup>2</sup> 11,1 + *b* > 0 and (*ζ*<sup>11</sup> − 2 *ζ*2 <sup>11</sup> + *b*) < *a* < (*ζ*<sup>11</sup> + 2 *ζ*2 <sup>11</sup> + *b*), then < 0 and two complex conjugate roots exist, given by *ζ*11,2,3 = *<sup>α</sup>*<sup>11</sup> <sup>±</sup> *iB*<sup>11</sup> where *<sup>α</sup>*<sup>22</sup> <sup>=</sup> *<sup>τ</sup>* <sup>2</sup> , *B*<sup>11</sup> = <sup>√</sup>4*δ*11−*τ*<sup>2</sup> <sup>2</sup> = *<sup>δ</sup>*<sup>11</sup> − *<sup>α</sup>*<sup>2</sup> <sup>11</sup>. *a* = (*ζ*<sup>11</sup> − 2 (*ζ*<sup>2</sup> <sup>11</sup>) + *b*) or *a* = (*ζ*<sup>11</sup> − 2 (*ζ*<sup>2</sup> <sup>11</sup>) + *a*2), then ( = 0) and two concident real roots exist given by *ζ*11,2 = *ζ*11,3 = *<sup>τ</sup>* <sup>2</sup> <sup>=</sup> *<sup>a</sup>*+*ζ*<sup>11</sup> <sup>2</sup> *a* < (*ζ*<sup>11</sup> − 2 (*ζ*<sup>2</sup> <sup>11</sup>) + *a*2) or *a*<sup>1</sup> > (*ζ*<sup>11</sup> − 2 (*ζ*<sup>2</sup> <sup>11</sup>) + *b*). Then, = 0 and two distinct real roots exists given by

$$
\zeta\_{11,2,3} = \frac{1}{2} (\pi \pm \sqrt{\triangle}).
$$

Applying the Routh–Hurwitz criterion [37], Equation (27) has roots with negative real parts if and only if *R*<sup>0</sup> < 1. Thus, the DFE is locally asymptotically stable.

#### **4. Numerical Methods for Solving the Proposed Model**

#### *4.1. GRK4M*

Consider the fractional derivatives with variable order given by the following equation: *C* <sup>0</sup> *<sup>D</sup>α*(*t*) *<sup>t</sup> ε*(*t*) = *f*(*t*,*ε*(*t*)), *Tf* ≥ *t* > 0, 1 ≥ *α*(*t*) > 0, (28)

$$
\varepsilon(0) = \varepsilon\_o.
$$

Using GRK4M [23], the approximate solution of (28) is:

$$\varepsilon\_{n+1} = \varepsilon\_n + \frac{1}{6}(K\_1 + 2K\_2 + 2K\_3 + K\_4),\tag{29}$$

$$K\_1 = \Upsilon f(t\_n, \varepsilon\_n),$$

$$K\_2 = \Upsilon f(t\_n + \frac{1}{2}\Upsilon, \varepsilon\_n + \frac{1}{2}K\_1),$$

$$K\_3 = \Upsilon f(t\_n + \frac{1}{2}\Upsilon, \varepsilon\_n + \frac{1}{2}K\_2),$$

$$K\_4 = \Upsilon f(t\_n + \Upsilon, \varepsilon\_n + K\_3),$$

where <sup>Υ</sup> <sup>=</sup> *<sup>τ</sup>α*(*tn*) Γ(*α*(*tn*) + 1) .

#### *4.2. Stability of GRK4M*

To investigate the stability of GRK4M, we shall utilize the following test problem of variable-order linear differential equation for simplicity:

$$\,\_{0}^{\mathbb{C}}D\_{t}^{a(t)}\varepsilon(t) = \varepsilon(t)\upsilon,\;\,T\_{f} \ge t > 0,\;\,\upsilon < 0,\;\,1 \ge a(t) > 0,\tag{30}$$

$$\varepsilon(0) = \varepsilon\_{\upsilon}.$$

As in [23], Equation (30) is written as follows:

$$
\varepsilon(t\_{i+1}) = \varepsilon(t\_i) + \frac{1}{6} \frac{\upsilon \tau^{a(t\_i)}}{\Gamma(1 + a(t\_i))} \varepsilon(t\_i), \ i = 0, 1, \ldots, n - 1. \tag{31}
$$

Then, we have the following equation [38]:

$$
\varepsilon(t\_{i+1}) = (1 + \frac{1}{6} \frac{\tau^{a(t\_i)} \upsilon}{\Gamma(1 + \alpha(t\_i))})^i \varepsilon\_0. \tag{32}
$$

The condition of stability [38]:

$$-1 < \left(\frac{1}{6} \frac{\tau^{\alpha(t\_i)} \upsilon}{\Gamma(1 + \alpha(t\_i))} + 1\right) < 1.$$

#### *4.3. CPC-*Θ*FDM*

Consider:

$$\prescript{\text{CPC}}{0}{D}\_{t}^{\mathfrak{a}(t)}\varepsilon(t) = \mathfrak{f}(t,\varepsilon(t)), \quad \varepsilon(0) = \varepsilon\_{0}, \quad 1 \ge \mathfrak{a}(t) > 0. \tag{33}$$

Relationship (6) can be expressed as follows:

$$\begin{split}{}^{\mathbb{C}\mathbb{CP}}D\_{t}^{\boldsymbol{a}(t)}\boldsymbol{\varepsilon}(t) &= \frac{1}{\Gamma(1-\boldsymbol{a}(t))} \int\_{0}^{t} (t-s)^{-\boldsymbol{a}(t)} (K\_{1}(\boldsymbol{a}(t))\boldsymbol{\varepsilon}(s) + K\_{0}(\boldsymbol{a}(t))\boldsymbol{\varepsilon}'(s)) ds, \\ &= K\_{1}(\boldsymbol{a})\_{0}^{\text{RL}}I\_{t}^{1-\boldsymbol{a}(t)}\boldsymbol{\varepsilon}(t) + K\_{0}(\boldsymbol{a}(t))\_{0}^{\mathbb{C}}D\_{t}^{\boldsymbol{a}(t)}\boldsymbol{\varepsilon}(t), \\ &= K\_{1}(\boldsymbol{a})\_{0}^{\text{RL}}D\_{t}^{\boldsymbol{a}(t)-1}\boldsymbol{\varepsilon}(t) + K\_{0}(\boldsymbol{a}(t))\_{0}^{\mathbb{C}}D\_{t}^{\boldsymbol{a}(t)}\boldsymbol{\varepsilon}(t), \end{split} \tag{34}$$

Using ΘFDM and GL-approximation, we can discretize (34) as shown below:

$$\begin{split} \left. \mathcal{C}\_{0}^{\text{CP}} D\_{t}^{\alpha(t)} \varepsilon(t) \right|\_{t=l^{n}} &= \frac{K\_{1}(\boldsymbol{a}(t\_{n}))}{\mathcal{T}^{\boldsymbol{a}(t\_{n})-1}} \left( \varepsilon\_{n+1} + \sum\_{i=1}^{n+1} \omega\_{i} \varepsilon\_{n+1-i} \right) \\ &+ \frac{K\_{0}(\boldsymbol{a}(t\_{n}))}{\mathcal{T}^{\boldsymbol{a}\_{n}}} \left( \varepsilon\_{n+1} - \sum\_{i=1}^{n+1} \varrho\_{i} \varepsilon\_{n+1-i} - \zeta\_{n+1} \varepsilon\_{0} \right), \end{split} \tag{35}$$

$$\frac{K\_1(\mathfrak{a}(t\_n))}{\tau^{\mathfrak{a}(t\_n)-1}} \left(\varepsilon\_{n+1} + \sum\_{i=1}^{n+1} \omega\_i \varepsilon\_{n+1-i}\right) + \frac{K\_0(\mathfrak{a}(t\_n))}{\tau^{\mathfrak{a}(t\_n)}} \left(\varepsilon\_{n+1} - \sum\_{i=1}^{n+1} \varrho\_i \varepsilon\_{n+1-i} - \mathfrak{g}\_{n+1} \varepsilon\_0\right)$$

$$\widetilde{\varepsilon} = (\Theta)\widetilde{\varepsilon}(\varepsilon(t\_n), t\_n) + (1 - \Theta)\widetilde{\varepsilon}(\varepsilon(t\_{n+1}), t\_{n+1}), \tag{36}$$

where, *<sup>ω</sup>*<sup>0</sup> <sup>=</sup> 1, *<sup>ω</sup><sup>i</sup>* = (<sup>1</sup> <sup>−</sup> *<sup>α</sup>*(*tn*) *<sup>i</sup>* )*ωi*−1, *<sup>t</sup> <sup>n</sup>* = *<sup>n</sup>τ*, *<sup>τ</sup>* = *Tf <sup>N</sup>* , *N* is a natural number, *<sup>i</sup>* = (−1)*i*−<sup>1</sup> *α*(*tn*) *i* , <sup>1</sup> = *α*(*tn*), *ς<sup>i</sup>* = *<sup>i</sup> α*(*tn*) <sup>Γ</sup>(1−*α*(*tn*)). Moreover, consider that [39]: 0 < *i*+<sup>1</sup> < *<sup>i</sup>* < ... < <sup>1</sup> = *α*(*tn*) < 1,

$$0 < \varsigma\_{i+1} < \varsigma\_i < \dots < \varsigma\_1 = \frac{1}{\Gamma(-\alpha(t\_n) + 1)}, \quad i = 1, 2, \dots, n+1.$$

**Remark 1.** *If K*1(*α*(*t*)) = 0 *and K*0(*α*(*t*)) = 1 *in (36), we can obtain the discretization of Caputo operator with theta finite difference technique (C-*Θ *FDM).*

#### *4.4. CPC-*Θ*FDM Stability Analysis*

The stability of method (36) will be considered here. We shall utilize the test problem of variable-order linear differential equation, for simplicity:

$$(\,^{\mathbb{C}P\mathbb{C}}D\_t^{\mathfrak{a}(t)})\varepsilon(t) = A\varepsilon(t), \quad t>0, \quad A<0, \quad 0<\mathfrak{a}(t)\le 1. \tag{37}$$

By (34) and GL-approximation, we can discretize (37) as shown below:

$$\begin{split} \frac{K\_{1}(a(t\_{n}))}{\tau^{a(t\_{n})-1}} \left( \varepsilon\_{n+1} + \sum\_{i=1}^{n+1} \omega\_{i} \varepsilon\_{n+1-i} \right) &+ \frac{K\_{0}(a(t\_{n}))}{\tau^{a(t\_{n})}} \left( \varepsilon\_{n+1} - \sum\_{i=1}^{n+1} \varrho\_{i} \varepsilon\_{n+1-i} - \varsigma\_{n+1} \varepsilon\_{0} \right) \\ = \Theta A \varepsilon\_{n} + (1 - \Theta) A \varepsilon\_{n+1}; \end{split} \tag{38}$$

$$\text{output } \mathbb{C} = \frac{K\_1(a(t\_n))}{\tau^{a(t\_n) - 1}}, \quad B = \frac{K\_0(a(t\_n))}{\tau^{a(t\_n)}}. \text{ Then, from boundedness theorem [40], we have:}$$

$$\varepsilon\_{n+1} = \frac{1}{\mathbb{C} + B} \left( A \varepsilon\_n - \mathbb{C} \sum\_{i=1}^{n+1} \omega\_i \varepsilon\_{n+1-i} + B \left( \sum\_{i=1}^{n+1} \varrho\_i \varepsilon\_{n+1-i} + \varrho\_{n+1} \varepsilon\_0 \right) \right) \le \varepsilon\_{n\prime} \tag{39}$$

This means *ε*<sup>0</sup> ≥ *ε*<sup>1</sup> ≥ ... ≥ *εn*−<sup>1</sup> ≥ *ε<sup>n</sup>* ≥ *εn*+1. Then, method (36) is stable.

#### *4.5. Convergence of the Method*

Equation (34) can be discretized as shown below:

$$\begin{split} \, \_0^{\mathbb{C}P\mathbb{C}} D\_t^{\alpha(t)} \varepsilon(t)|\_{t=t^n} &= \frac{K\_1(\mathfrak{a}(t\_n))}{\tau^{\mathfrak{a}(t\_n)-1}} \left( \varepsilon\_{n+1} + \sum\_{i=1}^{n+1} \omega\_i \varepsilon\_{n+1-i} \right) \\ &+ \frac{K\_0(\mathfrak{a}(t\_n))}{\tau^{\mathfrak{a}\_n}} \left( \varepsilon\_{n+1} - \sum\_{i=1}^{n+1} \varrho\_i \varepsilon\_{n+1-i} - \varsigma\_{n+1} \varepsilon\_0 \right), \end{split} \tag{40}$$

$$\frac{\mathcal{K}\_{1}(\boldsymbol{a}(t\_{n}))}{\tau^{\mathbf{a}(t\_{n})-1}} \left(\boldsymbol{\varepsilon}\_{n+1} + \sum\_{i=1}^{n+1} \omega\_{i} \boldsymbol{\varepsilon}\_{n+1-i}\right) + \frac{\mathcal{K}\_{0}(\boldsymbol{a}(t\_{n}))}{\tau^{\mathbf{a}(t\_{n})}} \left(\boldsymbol{\varepsilon}\_{n+1} - \sum\_{i=1}^{n+1} \varrho\_{i} \boldsymbol{\varepsilon}\_{n+1-i} - \boldsymbol{\zeta}\_{n+1} \boldsymbol{\varepsilon}\_{0}\right)$$

$$-\Theta \xi(\boldsymbol{\varepsilon}(t\_{n}), t\_{n}) - (1 - \Theta)\xi(\boldsymbol{\varepsilon}(t\_{n+1}), t\_{n+1}) = T\_{\mathrm{Rtu}}\tag{41}$$

where

$$||T\_{Rn}||\_{\infty} < W, \quad W = \mathbb{C} \max\_{0 \le i \le n+1} |\varepsilon\_{i+1}|\_{\prime}.$$

$$\mathcal{C} = (\tau^{a(t\_i)-1} + \tau^{a(t\_i)}).$$

The proposed method is convergent because it is stable and consistent [41], then (41) is convergent.

#### **5. Numerical Results**

In the following, we solved (8) numerically using GRK4M (29) and CPC-ΘFDM (36). Using CPC-Θ FDM for solving (8), we obtained (9*N* + 9) of the nonlinear algebraic system with (9*N* + 9) unknown.

$$\left(S\_{\prime} \, V\_{1\prime} \, V\_{2\prime} \, V\_{3\prime} \, A \, I\_{\prime\prime} \, I\_{V\prime} \, I\_{S\prime} \, R\right)$$

can be solved using an appropriate iterative method based on the assumed beginning conditions. For the real data, we use [24]; the authors in this reference used the literature to obtain some parameter values and the remaining values were fitted to the data for the state of Texas, USA. They fitted the data of (8) solutions with the data for the state of Texas from 13 March to 29 June 2021 [29,42]. The model was fitted with three datasets, Moderna, Janssen, and Pfizer, with immunization data for Texas state. The three vaccination rates *υ*1, *υ*<sup>2</sup> and *υ*<sup>3</sup> corresponding to each vaccine as well as the effective contact rate for COVID-19 transmission, *β*, are estimated. According to publicly available data, the total population of the state of Texas, USA, for the year 2021 was 29,200,000 [1]. Let *R*(0) = 5000, *V*2(0) = 4,016,005, *A*(0) = 50,000, *V*3(0) = 129,859, *S*(0) = 24,000,000, *IU*(0) = 17,000, *IV*(0) = 15,000, *V*1(0) = 4,115,127 and *IS*(0) = 10,000. The parameter values are given in Table 2. To show that the proposed scheme is efficient, we compare the results that we obtained in this paper with the results that were found in reference [24], which are given in Figure 2 in constant fractional order. Figure 3 shows the behavior of the approximate

solution of (8) (using the method in [24]) with different values of *α*(*t*). As can be seen from this figure, when the value of the fractional derivative changes over time, the results are different and this can dramatically affect the behavior of the model. This confirms the generality of the variable-order derivatives. Unfortunately, this method gives us unstable solutions, as in Figure 4a, when the value of the step size equals one. Moreover, we obtained the stable solutions using the proposed method CPC-ΘFDM and Θ = 0, in the fully implicit case given in this paper. This confirms that the method in [24] is stable only when the step size is very small, while our used method is stable regardless of the value of the step size. Figure 5 shows the behavior of the approximate solution of (8) (using CPC-ΘFDM and Θ = 0.5, *Q* = 0.00025) with different values of *α*(*t*). The approximate solution behavior of (8) is shown in Figure 6 (Θ = 1 and using CPC-ΘFDM) with different values of *α*(*t*), *Q* = 0.00025. The approximate solution behavior of (8) (using GRK4M with different values of *α*(*t*)) is shown in Figure 7. Figure 8 shows the behavior of the approximate solution of (8) (using CPC-Θ FDM when *K*0(*α*(*t*)) = 1, *K*1(*α*(*t*)) = 0 and Θ = 0) with different values of *α*(*t*). We noted that by comparing our results with different variable orders and constant orders as given in [24] and Figure 5, the result in the case of constant order is agreement. Moreover, by compering the results given in Figures 7 and 8, the result given using CPC-ΘFDM (fully implicit case) is convergent, better than the results given using GRK4 when we use nonlinear *α*(*t*). Figure 9 shows the relation between *R* and *Iv*, *Iu*, *Is*, *A* using CPC-ΘFDM (fully implicit case) and nonlinear *α*(*t*). Furthermore, we found that the variable-order derivative order model is a more general model than the fractional order model given in [24] and integer order; a new behavior of the solution appears by using different values of *α*(*t*). Moreover, we can obtain the fractional Caputo operator as a special case from the CPC operator when *K*0(*α*(*t*)) = 1, *K*1(*α*(*t*)) = 0. Moreover, we can obtain the fractional Caputo operator as a special case from the CPC operator if *K*0(*α*(*t*)) = 1, *K*1(*α*(*t*)) = 0. The solutions obtained using the new method CPC-ΘFDM can be explicit (Θ = 1) or implicit (0 ≤ Θ ≤ 1,) and fully implicit with accurate solution when (Θ = 0).

**Figure 2.** Real data [24] versus fitting model (8).

**Figure 3.** The solution behavior using the method in [24] with different values of *α*(*t*).

**Figure 4.** The solution behavior using the method [24] in (**a**) and using CPC-ΘFDM and Θ = 0, in (**b**).

**Figure 5.** The solution behavior acquired via CPC-ΘFDM and Θ = 0.5, of (8).

**Figure 6.** The solution behavior acquired via CPC-ΘFDM and Θ = 1 of (8).

**Figure 7.** The solution behavior acquired via GRK4M of (8).

**Figure 8.** The solution behavior acquired via CPC-ΘFDM and Θ = 0 of (8).

**Figure 9.** The relation between the variables concerning nonlinear *α*(*t*) using CPC-ΘFDM and Θ = 0.

#### **6. Conclusions**

A novel hybrid variable-order fractional multi-vaccination model for COVID-19 is presented in this paper in order to further explore the spread of COVID-19. The main advantage of the hybrid variable-order fractional operator is that it can be defined as a linear combination of the variable-order integral of Riemann–Liouville and the variable-order Caputo derivative; it is one of the most effective and reliable operators and it is more general than the Caputo fractional operator. The proposed model's dynamics are improved and its complexity is increased by employing variable-order fractional derivatives. Furthermore, the variable-order fractional Caputo operator can be derived as a special case from the CPC operator. Existence, boundedness, uniqueness, positivity and stability of the proposed model are established for the model. To be compatible with the physical model, a new parameter *σ* is added. The proposed model is numerically studied using CPC-ΘFDM and GRK4M. CPC-ΘFDM depends on the values of the factor Θ. It can be explicit (Θ = 1) or fully implicit (Θ = 0) with a large stability region. We compared our results with the real data from the state of Texas in the United States. Moreover, the results obtained from the CPC-ΘFDM are more stable than the results obtained from the proposed method in [24]. As a result, some graphs are provided for various linear and non-linear variable-order derivatives. In the future, the presented study can be extended to optimal control and to examine the impact of multiple vaccination strategies on the dynamics of COVID-19 in a population.

**Author Contributions:** Conceptualization, T.A.A.; Methodology, N.S. and S.M.A.-M.; Software, S.M.A.-M. and R.G.S.; Formal analysis, N.S.; Investigation, R.G.S. and T.A.A.; Resources, S.M.A.-M.; Writing—original draft, T.A.A.; Writing —review and editing, N.S. and R.G.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
