# Mathematical Modeling and Hyers-Ulam Stability for a Nonlinear
Epidemiological Model with $\Phi_{p}$ Operator and Mittag-Leffler Kernel

Achraf Zinihi<EMAIL_ADDRESS>Moulay Rchid Sidi Ammi
<EMAIL_ADDRESS>Matthias Ehrhardt<EMAIL_ADDRESS>Department of Mathematics, MAIS Laboratory, AMNEA Group, Faculty of Sciences
and Technics,
Moulay Ismail University of Meknes, Errachidia 52000, Morocco  University of
Wuppertal, Chair of Applied and Computational Mathematics,
Gaußstrasse 20, 42119 Wuppertal, Germany

(February 19, 2024)

###### Abstract

This paper investigates a novel nonlinear singular fractional SI model with
the $\Phi_{p}$ operator and the Mittag-Leffler kernel. The initial
investigation includes the existence, uniqueness, boundedness, and non-
negativity of the solution. We then establish Hyers-Ulam stability for the
proposed model in Banach space. Optimal control analysis is performed to
minimize the spread of infection and maximize the population of susceptible
individuals. Finally, the theoretical results are supported by numerical
simulations.

###### keywords:

epidemiological model , fractional differential equations , $p$-Laplacian
operator , numerical approximations , optimal control.
2020 Mathematics Subject Classification. 92C60, 34A08, 47H20, 33F05, 49J20.

††journal:

##  1 Introduction

Fractional-order models have attracted considerable interest from researchers
in a wide variety of disciplines. Over the past two decades, these models have
found applications in a wide variety of scientific and engineering fields,
including modern physics, signal theory, control theory, hydrodynamics,
viscoelastic theory, fluid dynamics, set theory, computer networks, biology,
etc. Relevant literature on these topics can be found in the works [15, 17,
20, 27, 32, 34, 38].

Recently, several researchers have studied fractional differential equations
(FDEs) with singularities using various mathematical methods. For example, Bai
and Qiu [7] established the existence and uniqueness (EU) of the solution to a
nonlinear singular boundary value problem (BVP) of FDEs using the
Krasnoselskii and Leray-Schauder fixed point theorems. They also demonstrated
applications to underscore their results. Agarwal, O’Regan and Staněk [2]
studied the EU for a singular fractional BVP using the Riemann-Liouville
fractional derivative. Bai and Fang [6] studied a singular nonlinear coupled
system of FDEs, using Leray-Schauder and Krasnoselskii fixed point techniques
for the $\mathrm{EU}$ of the solution. Vong [37] studied FDEs with singularity
and non-local boundary conditions using the Schauder fixed-point approach and
upper-lower solution techniques. Pu et al. [33] studied positive solutions of
a multipoint BVP with singularity and applied their results to a specific
example. Khan, Chen and Sun [25] studied nonlinear FDEs with singularity and
$p$-Laplacian to establish the EU of the solution and performed stability
analysis.

Mathematical models have long been indispensable tools for understanding and
predicting the dynamics of infectious diseases. By quantifying the complex
interactions between pathogens and populations, these models enable
researchers and policymakers to gain insight into the spread of disease and
evaluate potential control strategies. Some of the pioneering work in epidemic
modeling can be attributed to Kermack and McKendrick [23], who introduced the
SIR model in 1927. This model divided the population into three compartments:
susceptible ($S$), infectious ($I$), and recovered ($R$). Using differential
equations, the model captured the transitions between these compartments and
laid the foundation for subsequent advances in epidemic modeling.

Over time, mathematical models of epidemics have evolved and expanded to
incorporate additional complexities. Researchers recognized the importance of
accounting for factors such as age structure, spatial heterogeneity, and
varying transmission rates. This led to the development of more sophisticated
compartmental models, such as the SEIR model, which introduced an exposed
compartment [4]. For a recent review on epidemiological models we refer the
interested reader to [18] and the references therein.

In addition, spatial epidemic models and network-based models emerged to
capture the influence of geographic location and social connectedness on
disease spread [22]. In recent years, advanced mathematical techniques have
further enhanced the capabilities of epidemic models. Network theory has
provided insights into the role of social connections in disease transmission,
allowing the exploration of targeted intervention strategies. Nonlinear
dynamics and chaos theory have shed light on complex epidemic behavior,
including the emergence of periodic outbreaks and bifurcations.

The main goal of this paper is to minimize the number of infected individuals
for the fractional SI model, which describes the evolution of two
compartments: susceptible individuals ($S$) and infected individuals ($I$),
considering the Atangana-Baleanu fractional derivative in the Caputo sense
($\mathcal{ABC}$ fractional derivative, for short) and the $\Phi_{p}$
operator. The inclusion of these operators in our epidemic model offers
several compelling motivations. First, the use of the $\mathcal{ABC}$
fractional derivative allows for the inclusion of memory effects in the model
and long-range interactions in disease transmission. Traditional derivative
operators assume instantaneous changes, which may not accurately capture the
dynamics of infectious diseases. By introducing fractional calculus, we can
account for the persistence of past infection rates, allowing for a more
realistic representation of disease transmission.

Second, the $\Phi_{p}$ operator introduces nonlinearity and non-local
interactions, reflecting the impact of infection on susceptible and infected
populations within the SI model. This extension allows us to analyze the
impact of localized outbreaks and potential hotspots within the epidemic
dynamics.

This SI model has the potential to provide insight into the spread of
infectious diseases under fractional order dynamics. The $\mathcal{ABC}$
fractional derivatives and the $\Phi_{p}$ operator contribute to a more
nuanced understanding of the interplay between susceptible and infected
individuals, providing a valuable tool for designing effective disease control
and mitigation strategies. The fractional nature of the system allows for the
incorporation of memory effects, enhancing the model’s realism and
applicability to real-world scenarios. Further exploration of this coupled
operator in this context provides an opportunity to delve into the intricate
fractional calculus aspects of the model, contributing to a broader
understanding of fractional-order epidemiological systems.

Nevertheless, we encountered challenges in establishing the positivity and
boundedness of the solutions associated with the proposed model. These
challenges arise from the nonlinearity introduced by the operator
${}^{\mathcal{ABC}}\mathcal{D}^{\alpha}_{t}[\Phi_{p}(\cdot)]$. In addition,
difficulties were encountered in determining the adjoint system, which is
crucial for establishing the necessary optimality conditions. At the same
time, difficulties were encountered in numerical approximations aimed at
determining the forward-backward $\mathcal{ABC}$ derivative.

Due to the distinct characteristics of these two operators, numerous experts
and scholars have made significant contributions to the scientific literature
in various fields. For example, but not exhaustively, the work [35] delves
into the analysis of a reaction-diffusion SIR biological model formulated as a
parabolic system of PDEs incorporating the $p$-Laplacian operator. This study
emphasizes the critical role of vaccine distribution in the induction of
immunity. The primary goal of this work is to develop an optimal control
strategy that is carefully designed to limit both the spread of infection and
the associated vaccination costs.

Jena, Chakraverty, and Baleanu [21] used the homotopy perturbation Elzaki
transformation method to derive solutions to an epidemic model of childhood
diseases involving the $\mathcal{ABC}$ fractional derivatives. Their primary
goal was to protect children from diseases preventable by vaccination. They
proposed the homotopy perturbation method (HPM) transform to address the
problem at hand because of the occasional challenges the Elzaki transform
faces in handling nonlinear terms within the FDEs. The study by [39] focused
primarily on nonlinear coupling $(p_{1},p_{2})$ Laplacian systems with the
nonsingular $\mathcal{ABC}$ fractional derivative. He established sufficient
criteria for the existence and uniqueness of solutions based on the parameter
values $p_{1}$ and $p_{2}$. This interest and wide application of these two
operators in various scientific disciplines underscores their broad and
versatile influence on the research landscape. Further insights and relevant
literature on these topics can be explored in works such as [8, 9, 16, 19,
29].

This paper is organized as follows. In Section 2, several key and crucial
definitions are given. In Section 3 we present the fractional optimal control
SI model with $\Phi_{p}$ operator and Mittag-Leffler kernel. In Section 4 we
prove the existence, uniqueness, non-negativity, and boundedness of the
solution. Section 5 deals with the Hyers-Ulam stability of the proposed
problem. Furthermore, Section 6 is devoted to the determination of the
necessary optimality conditions. Before concluding the present study,
interesting numerical approximations are explained in Section 7. Finally, we
conclude our study in Section 8.

##  2 Notations and preliminaries

In this section, we will recall some definitions and properties that we will
need in the next sections on the Mittag-Leffler function, which plays an
important role in the solution of fractional order differential and integral
equations. To do this, we will introduce the $\mathcal{ABC}$ fractional
derivative, the $\mathcal{AB}$ fractional integral, and the $\Phi_{p}$
operator. For the rest of this paper, we will choose $\alpha\in(0,1)$ and
$2\leq p<\infty$.

First, the two-parameter and one-parameter Mittag-Leffler functions are
defined as follows.

######  Definition 1 (Mittag-Leffler function, [5]).

Let $\omega\in\mathbb{C}$ such that $\operatorname{Re}(\omega)>0$, the Mittag-
Leffler function $E_{\alpha,\beta}$ is defined by

$E_{\alpha,\beta}(\omega)=\sum_{k=0}^{\infty}\frac{\omega^{k}}{\Gamma(\alpha
k+\beta)},\quad\beta>0,$ (1)

where

$\Gamma(\omega)=\int_{0}^{+\infty}e^{-t}t^{\omega-1}\,dt.$

If $\beta=1$, then

$E_{\alpha}(\omega)=E_{\alpha,1}(\omega)=\sum_{k=0}^{\infty}\frac{\omega^{k}}{\Gamma(\alpha
k+1)}.$

We note that (1) was introduced by Magnus Gösta Mittag-Leffler [28] in 1903
with $\beta=1$. Next, we provide basic facts about the Gamma function and a
link to the exponential function.

######  Remark 1.

  * i.

$\Gamma(\omega+1)=\omega\,\Gamma(\omega)$.

  * ii.

By definition we have $\Gamma(1)=1$, and for any $n\in\mathbb{N}$ we find
$\Gamma(n+1)=n!$.

  * iii.

$E_{\alpha,\beta}$ is a generalization of the exponential function, and we
have $E_{1}(\omega)=E_{1,1}(\omega)=e^{\omega}$.

We now define the $\mathcal{ABC}$ fractional derivative and the $\mathcal{AB}$
fractional integral. Let $T>0$, $f\in H^{1}(0,T)$ and $t\in(0,T)$. We set
$B(\alpha)=(1-\alpha)+\frac{\alpha}{\Gamma(\alpha)}$ and
$\gamma=\frac{\alpha}{1-\alpha}$.

######  Definition 2 ($\mathcal{ABC}$ fractional derivative, [5]).

  * a.

The $\mathcal{ABC}$ derivative of order $\alpha$ of $f$ with base point $0$ is
defined at point $t$ by

${}^{\mathcal{ABC}}\mathcal{D}^{\alpha}_{t}f(t)=\frac{B(\alpha)}{1-\alpha}\int_{0}^{t}f^{\prime}(y)E_{\alpha}[-\gamma(t-y)^{\alpha}]\,dy.$
(2)

  * b.

The backward $\mathcal{ABC}$ derivative with base point $T$, is given by

${}^{\mathcal{ABC}}_{T}\mathcal{D}^{\alpha}_{t}f(t)=-\frac{B(\alpha)}{1-\alpha}\int_{t}^{T}f^{\prime}(y)E_{\alpha}[-\gamma(y-t)^{\alpha}]\,dy.$
(3)

Note that if we let $\alpha\to 1$ in (2), then we get the usual derivative
$\partial_{t}$.

######  Definition 3 ($\mathcal{AB}$ fractional integral, [5]).

The $\mathcal{AB}$ fractional integral operator with base point $0$, is
written as

${}^{\mathcal{AB}}I^{\alpha}f(t)=\frac{1-\alpha}{B(\alpha)}f(t)+\frac{\alpha}{B(\alpha)\Gamma(\alpha)}\int_{0}^{t}f(s)(t-s)^{\alpha-1}\,ds.$
(4)

The following Lemma 1 describes the basic relation of the $\mathcal{ABC}$
fractional derivative and the $\mathcal{AB}$ fractional integral. It will help
us in the later Sections 4 and 6 to prove the positivity and boundedness of
the solution and to derive the adjoint system associated with the state model.

######  Lemma 1 ([1, Proposition 3.4.]).

Under the previous assumptions, we have

${}^{\mathcal{AB}}I^{\alpha}\bigl{(}{}^{\mathcal{ABC}}\mathcal{D}^{\alpha}_{t}f(t)\bigr{)}=f(t)-f(0).$
(5)

The last part of this section deals with the $\Phi_{p}$ operator. Let
$\Phi_{p}$ be the function defined on $\mathbb{R}\to\mathbb{R}$ by

$\Phi_{p}(w)=|w|^{p-2}w,\quad\forall w\in\mathbb{R}.$

In addition, we have,

$\Phi_{p}^{-1}(\cdot)=\Phi_{p^{*}}(\cdot),\quad\text{where}\quad\frac{1}{p}+\frac{1}{p^{*}}=1.$

The two lemmas below will help us to prove that the proposed problem has a
unique solution and to prove the associated Hyers-Ulam stability.

######  Lemma 2 ([25, Lemma 1.3]).

  * i.

Let $x,y\in\mathbb{R}$ such that $|x|,|y|\leq k$, then

$|\Phi_{p}(x)-\Phi_{p}(y)|\leq(p-1)k^{p-2}|x-y|.$ (6)

  * ii.

If $1<p\leq 2$ and $x,y>0$ such that $|x|,|y|\geq k>0$, then

$|\Phi_{p}(x)-\Phi_{p}(y)|\leq(p-1)k^{p-2}|x-y|.$ (7)

######  Lemma 3.

Let $\varphi$ be a continuous bounded function on $[0,T]$ into $\mathbb{R}$.
Then

${}^{\mathcal{ABC}}\mathcal{D}^{\alpha}_{t}\bigl{[}\Phi_{p}(\varphi(t))\bigr{]}\cdot\varphi(t)\geq\frac{p-1}{p}\,{}^{\mathcal{ABC}}\mathcal{D}^{\alpha}_{t}\bigl{(}\Phi_{p}(\varphi(t))\varphi(t)\bigr{)}.$
(8)

###### Proof.

The inequality (8) can be rewritten as

$\mathcal{X}:={}^{\mathcal{ABC}}\mathcal{D}^{\alpha}_{t}\bigl{[}\Phi_{p}(\varphi(t))\bigr{]}\cdot\varphi(t)-\frac{p-1}{p}{}^{\mathcal{ABC}}\mathcal{D}^{\alpha}_{t}\bigl{(}\Phi_{p}(\varphi(t))\varphi(t)\bigr{)}\geq
0.$

Since

$\bigl{[}\Phi_{p}(\varphi(t))\bigr{]}^{\prime}=\bigl{[}|\varphi(t)|^{p-2}\varphi(t)\bigr{]}^{\prime}=\Bigl{[}e^{(p-2)\ln|\varphi(t)|}\varphi(t)\Bigr{]}^{\prime}=(p-1)|\varphi(t)|^{p-2}\varphi^{\prime}(t),$
(9)

we have

$\bigl{[}\Phi_{p}(\varphi(t))\varphi(t)\bigr{]}^{\prime}=\bigl{[}\Phi_{p}(\varphi(t))\bigr{]}^{\prime}\varphi(t)+\Phi_{p}(\varphi(t))\varphi^{\prime}(t)=\frac{p}{p-1}\bigl{[}\Phi_{p}(\varphi(t))\bigr{]}^{\prime}\varphi(t).$

Consequently,

$\begin{split}\mathcal{X}&={}^{\mathcal{ABC}}\mathcal{D}^{\alpha}_{t}\bigl{[}\Phi_{p}(\varphi(t))\bigr{]}\cdot\varphi(t)-\frac{p-1}{p}{}^{\mathcal{ABC}}\mathcal{D}^{\alpha}_{t}\bigl{(}\Phi_{p}(\varphi(t))\varphi(t)\bigr{)}\\\
&=\frac{B(\alpha)}{1-\alpha}\biggl{(}\varphi(t)\int_{0}^{t}\bigl{[}\Phi_{p}(\varphi(z))\bigr{]}^{\prime}E_{\alpha}\bigl{[}-\gamma(t-z)^{\alpha}\bigr{]}\,dz-\frac{p-1}{p}\int_{0}^{t}\bigl{[}\Phi_{p}(\varphi(z))\varphi(z)\bigr{]}^{\prime}E_{\alpha}\bigl{[}-\gamma(t-z)^{\alpha}\bigr{]}\,dz\biggr{)}\\\
&=\frac{B(\alpha)}{1-\alpha}\biggl{(}\varphi(t)\int_{0}^{t}\bigl{[}\Phi_{p}(\varphi(z))\bigr{]}^{\prime}E_{\alpha}\bigl{[}-\gamma(t-z)^{\alpha}\bigr{]}\,dz-\int_{0}^{t}\bigl{[}\Phi_{p}(\varphi(z))\bigr{]}^{\prime}\varphi(z)E_{\alpha}\bigl{[}-\gamma(t-z)^{\alpha}\bigr{]}\,dz\biggr{)}.\end{split}$

Thus,

$\begin{split}\mathcal{X}&=\frac{B(\alpha)}{1-\alpha}\int_{0}^{t}\bigl{(}\varphi(t)-\varphi(z)\bigr{)}\bigl{[}\Phi_{p}(\varphi(z))\bigr{]}^{\prime}E_{\alpha}\bigl{[}-\gamma(t-z)^{\alpha}\bigr{]}\,dz\\\
&=\frac{B(\alpha)}{1-\alpha}\int_{0}^{t}\Bigl{(}\int_{z}^{t}\varphi^{\prime}(s)\,ds\Bigr{)}\bigl{[}\Phi_{p}(\varphi(z))\bigr{]}^{\prime}E_{\alpha}\bigl{[}-\gamma(t-z)^{\alpha}\bigr{]}\,dz\\\
&=\frac{B(\alpha)}{1-\alpha}\int_{0}^{t}\varphi^{\prime}(s)\Bigl{(}\int_{0}^{s}\bigl{[}\Phi_{p}(\varphi(z))\bigr{]}^{\prime}E_{\alpha}\bigl{[}-\gamma(t-z)^{\alpha}\bigr{]}\,dz\Bigr{)}\,ds\\\
&\stackrel{{\scriptstyle\eqref{E2.8}}}{{=}}\frac{B(\alpha)}{2(1-\alpha)(p-1)}\int_{0}^{t}|\varphi(s)|^{2-p}\bigl{(}E_{\alpha}[-\gamma(t-s)^{\alpha}]\bigr{)}^{-1}\frac{\partial}{\partial
s}\biggl{(}\Bigl{(}\int_{0}^{s}\bigl{[}\Phi_{p}(\varphi(z))\bigr{]}^{\prime}E_{\alpha}\bigl{[}-\gamma(t-z)^{\alpha}\bigr{]}\,dz\Bigr{)}^{2}\biggr{)}\,ds.\end{split}$

Since $p\geq 2$ and $\varphi$ is bounded, then $|\varphi(s)|^{p-2}\leq C$ and
that means $|\varphi(s)|^{2-p}\geq\sigma:=\frac{1}{C}$. Therefore,

$\mathcal{X}\geq\frac{\sigma
B(\alpha)}{2(1-\alpha)(p-1)}\int_{0}^{t}\bigl{(}E_{\alpha}[-\gamma(t-s)^{\alpha}]\bigr{)}^{-1}\frac{\partial}{\partial
s}\biggl{(}\Bigl{(}\int_{0}^{s}\bigl{[}\Phi_{p}(\varphi(z))\bigr{]}^{\prime}E_{\alpha}\bigl{[}-\gamma(t-z)^{\alpha}\bigr{]}\,dz\Bigr{)}^{2}\biggr{)}\,ds.$

Knowing that

$\frac{\partial}{\partial
s}\Bigl{(}\frac{1}{E_{\alpha}[-\gamma(t-s)^{\alpha}]}\Bigr{)}=-\frac{\gamma(t-s)^{\alpha-1}E_{\alpha,\alpha}[-\gamma(t-s)^{\alpha}]}{E_{\alpha}[-\gamma(t-s)^{\alpha}]^{2}}\leq
0,\ \text{ and }\ E_{\alpha}(0)=1.$

Then, an integration by part gives

$\begin{split}\mathcal{X}\geq&\frac{\sigma
B(\alpha)}{2(1-\alpha)(p-1)}\biggl{[}\bigl{(}E_{\alpha}[-\gamma(t-s)^{\alpha}]\bigr{)}^{-1}\Bigl{(}\int_{0}^{s}\bigl{[}\Phi_{p}(\varphi(z))\bigr{]}^{\prime}E_{\alpha}\bigl{[}-\gamma(t-z)^{\alpha}\bigr{]}\,dz\Bigr{)}^{2}\biggr{]}_{s=0}^{s=t}\\\
&+\frac{\sigma\gamma
B(\alpha)}{2(1-\alpha)(p-1)}\int_{0}^{t}\frac{(t-s)^{\alpha-1}E_{\alpha,\alpha}[-\gamma(t-s)^{\alpha}]}{E_{\alpha}[-\gamma(t-s)^{\alpha}]^{2}}\Bigl{(}\int_{0}^{s}\bigl{[}\Phi_{p}(\varphi(z))\bigr{]}^{\prime}E_{\alpha}\bigl{[}-\gamma(t-z)^{\alpha}\bigr{]}\,dz\Bigr{)}^{2}\,ds.\end{split}$

Thus,

$\begin{split}\mathcal{X}&\geq\frac{\sigma
B(\alpha)}{2(1-\alpha)(p-1)}\bigg{[}\Bigl{(}\int_{0}^{s}\bigl{[}\Phi_{p}(\varphi(z))\bigr{]}^{\prime}E_{\alpha}\bigl{[}-\gamma(t-z)^{\alpha}\bigr{]}\,dz\Bigr{)}^{2}\\\
&\quad+\int_{0}^{t}\frac{(t-s)^{\alpha-1}E_{\alpha,\alpha}[-\gamma(t-s)^{\alpha}]}{E_{\alpha}[-\gamma(t-s)^{\alpha}]^{2}}\Bigl{(}\int_{0}^{s}\bigl{[}\Phi_{p}(\varphi(z))\bigr{]}^{\prime}E_{\alpha}\bigl{[}-\gamma(t-z)^{\alpha}\bigr{]}dz\Bigr{)}^{2}\,ds\bigg{]}.\end{split}$

This concludes the proof. ∎

##  3 Mathematical model with vital dynamics

The SI model is a mathematical framework used to study the dynamics of an
epidemic within a population. It divides the population ($N$) into different
compartments based on their disease status. Let’s look at the compartments and
how the epidemic spreads from one to the other:

Susceptible $S$: This compartment represents individuals who are susceptible
to the disease and can become infected if they come in contact with an
infectious individual.

Infected $I$: This compartment includes individuals who are currently infected
with the disease and can transmit it to susceptible individuals through
various modes of transmission, such as respiratory droplets, physical contact,
contaminated surfaces, etc.

The following Table 1 contains the transmission coefficients for the proposed
SI model, and Figure 1 summarizes how the epidemic spreads from one
compartment to another.

Table 1: Vital dynamics parameters of the SI model. Symbol | Description
---|---
$\mu$ ($\geq 0$) | Birth rate
$\beta$ ($\geq 0$) | Effective contact rate
$\eta$ ($\geq 0$) | Natural mortality rate

SI$\mu N$$\eta S$$\beta SI$$\eta I$$uS$ Figure 1: Conversion scheme for the
proposed SI model.

Let $\mathcal{I}_{T}:=[0,T]$. Given the assumptions explained above, our SI
model is formulated as follows

$\begin{cases}&{}^{\mathcal{ABC}}\mathcal{D}^{\alpha}_{t}\bigl{[}\Phi_{p}(S(t))\bigr{]}=\mu
N-\beta SI-uS-\eta S,\\\
&{}^{\mathcal{ABC}}\mathcal{D}^{\alpha}_{t}\bigl{[}\Phi_{p}(I(t))\bigr{]}=\beta
SI+uS-\eta I,\end{cases}\quad t\in\mathcal{I}_{T},$ (10)

with the total population $N(t)=S(t)+I(t)$ and supplied with the positive
initial data

$S(0)=S_{0}\geq 0,\text{ and }\ I(0)=I_{0}\geq 0.$ (11)

The term $\mu N-\beta SI-uS-\eta S$ in the first equation in (10) represents
the change in the susceptible population over time, influenced by infection
(quadratic contact term) and external factors, where

$u\in\mathcal{U}_{ad}=\bigl{\\{}u\text{ is Lebesgue integrable in
}\mathcal{I}_{T}\ |\ 0\leq u(t)<1\bigr{\\}}.$ (12)

The variable $u$ in (10) represents the optimal control (vaccination), while
the term $uS$ denotes the vaccinated susceptible individuals. By this variable
$u$, we want to minimize the spread of infection and maximize the population
of susceptible individuals, cf. Section 6.

Let $\vartheta=(x,y)=(S,I)$ such that
$\vartheta_{0}=(x_{0},y_{0})=(S_{0},I_{0})$. We consider the vector function
$\Psi$ defined by

$\Psi(\vartheta(t))=\bigl{(}\Psi_{1}(\vartheta(t)),\Psi_{2}(\vartheta(t))\bigr{)},$

where

$\begin{cases}\Psi_{1}(\vartheta(t))=\mu(x(t)+y(t))-\beta
x(t)y(t)-u(t)x(t)-\eta x(t),\\\ \Psi_{2}(\vartheta(t))=\beta
x(t)y(t)+u(t)x(t)-\eta y(t).\end{cases}\ t\in\mathcal{I}_{T},$ (13)

with the initial condition $\vartheta_{0}=(x_{0},y_{0})$. Furthermore, if the
solutions of (10) are bounded, then $\Psi$ satisfies a Lipschitz condition

$\bigl{\|}\Psi(\tilde{\vartheta}(t))-\Psi(\bar{\vartheta}(t))\bigr{\|}\leq\xi\bigl{\|}\tilde{\vartheta}(t)-\bar{\vartheta}(t)\bigr{\|}.$
(14)

The supremum norm of $\vartheta$, denoted by $\|\cdot\|_{\mathcal{I}_{T}}$, is
defined as
$\|\vartheta\|_{\mathcal{I}_{T}}:=\sup_{t\in\mathcal{I}_{T}}\|\vartheta(t)\|$.
It is obvious that $\mathcal{C}(\mathcal{I}_{T},\mathbb{R}^{2})$, equipped
with the norm $\|\vartheta\|_{\mathcal{I}_{T}}$, forms a Banach space.
Consequently, the problem (10)–(11) can be reformulated in
$\mathcal{C}(\mathcal{I}_{T},\mathbb{R}^{2})$ as

$\begin{cases}{}^{\mathcal{ABC}}\mathcal{D}^{\alpha}_{t}\bigl{[}\Phi_{p}(\vartheta(t))\bigr{]}=\Psi(\vartheta(t)),\\\
\vartheta(0)=\vartheta_{0},\end{cases}\quad t\in\mathcal{I}_{T},$ (15)

where

$\Phi_{p}(\vartheta(t))=\bigl{(}\Phi_{p}(x(t)),\Phi_{p}(y(t))\bigr{)}.$

##  4 Existence and uniqueness of the solution

In this section we will prove the existence and uniqueness of solutions of the
nonlinear fractional system (10)–(11) using a fixed point theorem.

######  Theorem 1 (Boundedness).

All solutions of the SI model (10)–(11) are bounded if

$|\mu-\eta|<\frac{2^{2-p}B(\alpha)}{1-\alpha}=\frac{\alpha+(1-\alpha)\Gamma(\alpha)}{(1-\alpha)\Gamma(\alpha)2^{p-2}},$
(16)

i.e. if $\Lambda:=2^{2-p}-\frac{(1-\alpha)|\mu-\eta|}{B(\alpha)}>0$ holds.

###### Proof.

Let $(x,y)$ be a solution of (10)–(11). We have

${}^{\mathcal{ABC}}\mathcal{D}^{\alpha}_{t}\bigl{[}\Phi_{p}(x)+\Phi_{p}(y)\bigr{]}={}^{\mathcal{ABC}}\mathcal{D}^{\alpha}_{t}\bigl{[}\Phi_{p}(x)\bigr{]}+{}^{\mathcal{ABC}}\mathcal{D}^{\alpha}_{t}\bigl{[}\Phi_{p}(y)\bigr{]}=(\mu-\eta)N,$

with $N(t)=x(t)+y(t)$. Applying the $\mathcal{AB}$ fractional integral (4), we
get

$\bigl{(}\Phi_{p}(x)+\Phi_{p}(y)\bigr{)}-\bigl{(}\Phi_{p}(x_{0})+\Phi_{p}(y_{0})\bigr{)}=(\mu-\eta){}^{\mathcal{AB}}I^{\alpha}N.$

Using Tartar’s inequality [30, Lemma 2.2.], we obtain

$\begin{split}\bigl{|}\Phi_{p}(x)+\Phi_{p}(y)\bigr{|}&=\bigl{|}|x|^{p-2}x+|y|^{p-2}y\bigr{|}\\\
&=\bigl{|}|x|^{p-2}x-|y|^{p-2}(-y)\bigr{|}=\bigl{|}|x|^{p-2}x-|(-y)|^{p-2}(-y)\bigr{|}\\\
&\geq 2^{2-p}|x-(-y)|^{p-1}=2^{2-p}|x+y|^{p-1}=2^{2-p}|N|^{p-1}.\end{split}$

Since we assume positive initial data $x_{0}\geq 0$ and $y_{0}\geq 0$, it
follows

$\displaystyle\bigl{|}\Phi_{p}x_{0})+\Phi_{p}(y_{0})\bigr{|}$
$\displaystyle=\bigl{|}|x_{0}|^{p-2}x_{0}+|y_{0}|^{p-2}y_{0}\bigr{|}=|x_{0}|^{p-2}x_{0}+|y_{0}|^{p-2}y_{0}$
$\displaystyle\leq|x_{0}+y_{0}|^{p-2}x_{0}+|x_{0}+y_{0}|^{p-2}y_{0}\leq|N_{0}|^{p-2}N_{0}.$

By the inverse triangular inequality, we have

$2^{2-p}|N|^{p-1}-|N_{0}|^{p-2}N_{0}\leq|\mu-\eta|{}^{\mathcal{AB}}I^{\alpha}|N|.$

If $N\in[-1,1]$, then $N$ is bounded. Else, we get $|N|\leq|N|^{p-1}$ (since
$p\geq 2$). After that,

$\begin{split}2^{2-p}|N|-|N_{0}|^{p-2}N_{0}&\leq|\mu-\eta|{}^{\mathcal{AB}}I^{\alpha}|N|\\\
&\leq\frac{(1-\alpha)|\mu-\eta|}{B(\alpha)}|N|+\frac{\alpha|\mu-\eta|}{B(\alpha)\Gamma(\alpha)}\int_{0}^{t}|N(s)|(t-s)^{\alpha-1}\,ds,\end{split}$

which gives

$\Lambda|N|\leq|N_{0}|^{p-2}N_{0}+\frac{\alpha(\mu+\eta)}{B(\alpha)\Gamma(\alpha)}\int_{0}^{t}|N(s)|(t-s)^{\alpha-1}\,ds.$

Finally, Gronwall’s inequality yields

$\begin{split}|N(t)|&\leq\frac{|N_{0}|^{p-2}N_{0}}{\Lambda}\exp\Bigl{(}\frac{\alpha(\mu+\eta)}{\Lambda
B(\alpha)\Gamma(\alpha)}\int_{0}^{t}(t-s)^{\alpha-1}\,ds\Bigr{)}\\\
&\leq\frac{|N_{0}|^{p-2}N_{0}}{\Lambda}\exp\Bigl{(}\frac{\alpha(\mu+\eta)t^{\alpha}}{\Lambda
B(\alpha)\Gamma(\alpha)\alpha}\Bigr{)}\\\
&=\frac{|N_{0}|^{p-2}N_{0}}{\Lambda}\exp\Bigl{(}\frac{(\mu+\eta)t^{\alpha}}{\Lambda
B(\alpha)\Gamma(\alpha)}\Bigr{)}.\end{split}$

This establishes the boundedness of solutions. ∎

In the rest of this work, we assume that the boundedness condition (16) is
satisfied.

######  Theorem 2 (Positivity).

All solutions of the SI model (10)–(11) are positive for positive data.

###### Proof.

Let $(x,y)$ be a solution of the SI model (10)–(11).

First, suppose $x<0$. Applying (5) to the first equation of (10), we have

$\Phi_{p}(x)-\Phi_{p}(x_{0})={}^{\mathcal{AB}}I^{\alpha}(\mu N-\beta xy-
ux-\eta x).$

Since $x_{0}\geq 0$ and by the hypothesis, we get

${}^{\mathcal{AB}}I^{\alpha}(\mu N-\beta xy-ux-\eta x)<0.$

Given that ${}^{\mathcal{AB}}I^{\alpha}$ is increasing, we get

$\mu N-\beta xy<ux+\eta x<0,\quad\text{i.e.}\quad$ (17)

From (10) and (17), we get

${}^{\mathcal{ABC}}\mathcal{D}^{\alpha}_{t}\bigl{[}\Phi_{p}(x)\bigr{]}=\mu
N-\beta xy-ux-\eta x<-ux-\eta x.$ (18)

Knowing that $x<0$, so its negative part $x^{-}=\max(0,-x)>0$. Knowing that
$x$ is bounded, so by multiplying (18) by $x^{-}$, using (8) and (17), we get

$\frac{p-1}{p}{}^{\mathcal{ABC}}\mathcal{D}^{\alpha}_{t}\bigl{[}|x|^{p-2}(x^{-})^{2}\bigr{]}\leq{}^{\mathcal{ABC}}\mathcal{D}^{\alpha}_{t}\bigl{[}\Phi_{p}(x)\bigr{]}x^{-}<-(u+\eta)(x^{-})^{2}<0,$

i.e.

${}^{\mathcal{ABC}}\mathcal{D}^{\alpha}_{t}\bigl{[}|x|^{p-2}(x^{-})^{2}\bigr{]}<0.$

Applying the Laplace transform to this equation (cf. [10]), we have

$\frac{B(\alpha)s^{\alpha}}{\alpha+(1-\alpha)s^{\alpha}}\mathscr{L}\bigl{(}|x|^{p-2}(x^{-})^{2}\bigr{)}-\frac{B(\alpha)s^{\alpha-1}}{\alpha+(1-\alpha)s^{\alpha}}|x_{0}|^{p-2}(x_{0}^{-})^{2}<0.$

Since $x_{0}\geq 0$, then $x_{0}^{-}=0$. Accordingly, we obtain

$\mathscr{L}\bigl{(}|x|^{p-2}(x^{-})^{2}\bigr{)}<0.$

Applying the inverse Laplace transform to this inequality, we have that
$|x|^{p-2}(x^{-})^{2}<0$, which contradicts the hypothesis. So we get $x\geq
0$.

Next, to establish the positivity of $y$, we use the same methodology as in
[3, 31]. We assume that there exists a $0<\tau\in\mathcal{I}_{T}$ such that

$\Phi_{p}(y(\tau))=0\ \text{ and }\ \Phi_{p}(y(t))<0,\ \ \forall
t\in(\tau,\ell],$

where $\ell\in\mathcal{I}_{T}$ is sufficiently close to $\tau$. Using (10), we
get

${}^{\mathcal{ABC}}\mathcal{D}^{\alpha}_{t}\bigl{[}\Phi_{p}(y(\tau))\bigr{]}=u(\tau)x(\tau)\geq
0.$

Applying the mean value theorem for the $\mathcal{ABC}$ derivative [12,
Theorem 2.2.], we find

$\Phi_{p}(y(t))-\Phi_{p}(y(\tau))=\Phi_{p}(y(t))\geq 0,$

which contradicts the assumption $\Phi_{p}(y(t))<0$. Consequently,
$\Phi_{p}(y(t))\geq 0$, and thus $y\geq 0$. ∎

We focus our attention on establishing that the SI model (10)–(11) has a
unique solution.

######  Theorem 3 (Characterization of the solution).

$\vartheta$ is a solution of (10)–(11) if and only if

$\begin{split}\vartheta(t)&=\Phi_{p^{*}}\Bigl{[}\Phi_{p}(\vartheta_{0})+{}^{\mathcal{AB}}I^{\alpha}\bigl{(}\Psi(\vartheta(t))\bigr{)}\Bigr{]}\\\
&=\Phi_{p^{*}}\Bigl{[}\Phi_{p}(\vartheta_{0})+\frac{1-\alpha}{B(\alpha)}\Psi(\vartheta(t))+\frac{\alpha}{B(\alpha)\Gamma(\alpha)}\int_{0}^{t}\Psi(\vartheta(s))(t-s)^{\alpha-1}\,ds\Bigr{]}=:\mathcal{T}(\vartheta(t)).\end{split}$
(19)

###### Proof.

By applying the $\mathcal{AB}$ fractional integral to the reformulated problem
(15), we obtain

$\Phi_{p}(\vartheta(t))=\Phi_{p}(\vartheta_{0})+\frac{1-\alpha}{B(\alpha)}\Psi(\vartheta(t))+\frac{\alpha}{B(\alpha)\Gamma(\alpha)}\int_{0}^{t}\Psi(\vartheta(s))(t-s)^{\alpha-1}\,ds.$

Consequently,

$\vartheta(t)=\mathcal{T}(\vartheta(t)):=\Phi_{p^{*}}\biggl{[}\Phi_{p}(\vartheta_{0})+\frac{1-\alpha}{B(\alpha)}\Psi(\vartheta(t))+\frac{\alpha}{B(\alpha)\Gamma(\alpha)}\int_{0}^{t}\Psi(\vartheta(s))(t-s)^{\alpha-1}\,ds\biggr{]}.$

∎

To study the existence and uniqueness of the solution of the proposed problem
with the fractional derivative $\mathcal{ABC}$, on the basis of the fixed
point theorem and the equation (19), the iterative fixed point formula is
given by

$\begin{cases}\vartheta^{n+1}(t)&=\Phi_{p^{*}}\Bigl{[}\Phi_{p}(\vartheta_{0})+\frac{1-\alpha}{B(\alpha)}\Psi(\vartheta^{n}(t))+\frac{\alpha}{B(\alpha)\Gamma(\alpha)}\int_{0}^{t}\Psi(\vartheta^{n}(s))(t-s)^{\alpha-1}\,ds\Bigr{]},\\\
\vartheta^{0}(t)&=\vartheta_{0}.\end{cases}\quad t\in[0,T],$ (20)

The difference between successive iterations is considered as follows

$e^{n+1}(t)=\vartheta^{n+1}(t)-\vartheta^{n}(t)=\mathcal{T}(\vartheta^{n+1}(t))-\mathcal{T}(\vartheta^{n}(t)),\quad
n\geq 0.$ (21)

This leads us to note that

$\vartheta^{n}(t)=\sum_{k=0}^{n}e^{k}(t),\ \text{ with }\
e^{0}(t)=\vartheta^{0}(t).$

######  Lemma 4.

Assuming the Lipschitz condition (14), we have

$\|e^{n}\|_{\mathcal{I}_{T}}\leq\biggl{[}\frac{(p^{*}-1)k^{p^{*}-2}\xi}{B(\alpha)}\Bigl{(}1-\alpha+\frac{t^{\alpha}}{\Gamma(\alpha)}\Bigr{)}\biggr{]}^{n}\|\vartheta^{0}\|_{\mathcal{I}_{T}}.$

###### Proof by induction.

Let $n\in\mathbb{N}$, for $n=0$ we get

$\|e^{0}(t)\|_{\mathcal{I}_{T}}=\|\vartheta^{0}(t)\|_{\mathcal{I}_{T}}\leq\biggl{[}\frac{(p^{*}-1)k^{p^{*}-2}\xi}{B(\alpha)}\Bigl{(}1-\alpha+\frac{t^{\alpha}}{\Gamma(\alpha)}\Bigr{)}\biggr{]}^{0}\|\vartheta^{0}(t)\|_{\mathcal{I}_{T}}.$

For $n=1$, by considering the equations (20), (21), Lemma 2, and the
triangular inequality of norms, we find

$\begin{split}\|e^{1}(t)\|_{\mathcal{I}_{T}}&=\|\vartheta^{1}(t)-\vartheta^{0}(t)\|_{\mathcal{I}_{T}}\\\
&=\Bigl{\|}\Phi_{p^{*}}\Bigl{[}\Phi_{p}(\vartheta_{0})+{}^{\mathcal{AB}}I^{\alpha}\bigl{(}\Psi(\vartheta^{0}(t))\bigr{)}\Bigr{]}-\Phi_{p^{*}}\Bigl{[}\Phi_{p}(\vartheta_{0})+{}^{\mathcal{AB}}I^{\alpha}\bigl{(}\Psi(\vartheta^{-1}(t))\bigr{)}\Bigr{]}\Bigr{\|}_{\mathcal{I}_{T}}\\\
&\leq(p^{*}-1)k^{p^{*}-2}\Bigl{\|}{}^{\mathcal{AB}}I^{\alpha}\bigl{(}\Psi(\vartheta^{0}(t))\bigr{)}-{}^{\mathcal{AB}}I^{\alpha}\bigl{(}\Psi(\vartheta^{-1}(t))\bigr{)}\Bigr{\|}_{\mathcal{I}_{T}}\\\
&\leq(p^{*}-1)k^{p^{*}-2}\Big{\|}\frac{1-\alpha}{B(\alpha)}\bigl{(}\Psi(\vartheta^{0}(t))-\Psi(\vartheta^{-1}(t))\bigr{)}\\\
&\qquad+\frac{\alpha}{B(\alpha)\Gamma(\alpha)}\int_{0}^{t}\bigl{(}\Psi(\vartheta^{0}(s))-\Psi(\vartheta^{-1}(s))\bigr{)}(t-s)^{\alpha-1}\,ds\Big{\|}_{\mathcal{I}_{T}}\\\
&\leq(p^{*}-1)k^{p^{*}-2}\frac{1-\alpha}{B(\alpha)}\bigl{\|}\Psi(\vartheta^{0}(t))-\Psi(\vartheta^{-1}(t))\bigr{\|}_{\mathcal{I}_{T}}\\\
&\qquad+(p^{*}-1)k^{p^{*}-2}\frac{\alpha}{B(\alpha)\Gamma(\alpha)}\int_{0}^{t}\bigl{\|}\Psi(\vartheta^{0}(s))-\Psi(\vartheta^{-1}(s))\bigr{\|}_{\mathcal{I}_{T}}(t-s)^{\alpha-1}\,ds.\end{split}$

Afterwards,

$\begin{split}\|e^{1}(t)\|_{\mathcal{I}_{T}}&\leq(p^{*}-1)k^{p^{*}-2}\frac{1-\alpha}{B(\alpha)}\xi\|\vartheta^{0}(t)-\vartheta^{-1}(t)\|_{\mathcal{I}_{T}}\\\
&\qquad+(p^{*}-1)k^{p^{*}-2}\frac{\alpha}{B(\alpha)\Gamma(\alpha)}\xi\int_{0}^{t}\|\vartheta^{0}(s)-\vartheta^{-1}(s)\|_{\mathcal{I}_{T}}(t-s)^{\alpha-1}\,ds\\\
&\leq(p^{*}-1)k^{p^{*}-2}\frac{1-\alpha}{B(\alpha)}\xi\|e^{0}(t)\|_{\mathcal{I}_{T}}+(p^{*}-1)k^{p^{*}-2}\frac{\alpha}{B(\alpha)\Gamma(\alpha)}\xi\int_{0}^{t}\|e^{0}(s)\|_{\mathcal{I}_{T}}(t-s)^{\alpha-1}\,ds\\\
&\leq(p^{*}-1)k^{p^{*}-2}\frac{1-\alpha}{B(\alpha)}\xi\|\vartheta^{0}\|_{\mathcal{I}_{T}}+(p^{*}-1)k^{p^{*}-2}\frac{\alpha}{B(\alpha)\Gamma(\alpha)}\xi\|\vartheta^{0}\|_{\mathcal{I}_{T}}\int_{0}^{t}(t-s)^{\alpha-1}\,ds\\\
&\leq(p^{*}-1)k^{p^{*}-2}\frac{1-\alpha}{B(\alpha)}\xi\|\vartheta^{0}\|_{\mathcal{I}_{T}}+(p^{*}-1)k^{p^{*}-2}\frac{\alpha}{B(\alpha)\Gamma(\alpha)\alpha}\xi\|\vartheta^{0}\|_{\mathcal{I}_{T}}t^{\alpha}\\\
&\leq\biggl{[}\frac{(p^{*}-1)k^{p^{*}-2}\xi}{B(\alpha)}\Bigl{(}1-\alpha+\frac{t^{\alpha}}{\Gamma(\alpha)}\Bigr{)}\biggr{]}^{1}\|\vartheta^{0}\|_{\mathcal{I}_{T}}.\end{split}$

We then assume that the property is true for order $n$, and show that it is
true for $n+1$. Using the same methodology as before, we obtain

$\begin{split}\|e^{n+1}(t)\|_{\mathcal{I}_{T}}&=\|\vartheta^{n+1}(t)-\vartheta^{n}(t)\|_{\mathcal{I}_{T}}\\\
&\leq(p^{*}-1)k^{p^{*}-2}\frac{1-\alpha}{B(\alpha)}\xi\|e^{n}(t)\|_{\mathcal{I}_{T}}+(p^{*}-1)k^{p^{*}-2}\frac{\alpha}{B(\alpha)\Gamma(\alpha)}\xi\int_{0}^{t}\|e^{n}(s)\|_{\mathcal{I}_{T}}(t-s)^{\alpha-1}\,ds.\end{split}$

Thus, we have

$\begin{split}\|e^{n+1}(t)\|_{\mathcal{I}_{T}}&\leq(p^{*}-1)k^{p^{*}-2}\frac{1-\alpha}{B(\alpha)}\xi\Biggl{(}\biggl{[}\frac{(p^{*}-1)k^{p^{*}-2}\xi}{B(\alpha)}\Bigl{(}1-\alpha+\frac{t^{\alpha}}{\Gamma(\alpha)}\Bigr{)}\biggr{]}^{n}\|\vartheta^{0}\|_{\mathcal{I}_{T}}\Biggr{)}\\\
&\qquad+\frac{(p^{*}-1)k^{p^{*}-2}\alpha}{B(\alpha)\Gamma(\alpha)}\xi\int_{0}^{t}\Biggl{(}\biggl{[}\frac{(p^{*}-1)k^{p^{*}-2}\xi}{B(\alpha)}\Bigl{(}1-\alpha+\frac{t^{\alpha}}{\Gamma(\alpha)}\Bigr{)}\biggr{]}^{n}\|\vartheta^{0}\|_{\mathcal{I}_{T}}\Biggr{)}(t-s)^{\alpha-1}\,ds\\\
&\leq\biggl{[}\frac{(p^{*}-1)k^{p^{*}-2}\xi}{B(\alpha)}\Bigl{(}1-\alpha+\frac{t^{\alpha}}{\Gamma(\alpha)}\Bigr{)}\biggr{]}^{n}\|\vartheta^{0}\|_{\mathcal{I}_{T}}\bigg{(}(p^{*}-1)k^{p^{*}-2}\frac{1-\alpha}{B(\alpha)}\xi\\\
&\qquad+(p^{*}-1)k^{p^{*}-2}\frac{\alpha}{B(\alpha)\Gamma(\alpha)}\xi\int_{0}^{t}(t-s)^{\alpha-1}\,ds\bigg{)}\\\
&\leq\biggl{[}\frac{(p^{*}-1)k^{p^{*}-2}\xi}{B(\alpha)}\Bigl{(}1-\alpha+\frac{t^{\alpha}}{\Gamma(\alpha)}\Bigr{)}\biggr{]}^{n+1}\|\vartheta^{0}\|_{\mathcal{I}_{T}}.\end{split}$

∎

######  Theorem 4 (Uniqueness).

If there is a $t_{0}$ which satisfies

$\frac{(p^{*}-1)k^{p^{*}-2}\xi}{B(\alpha)}\Bigl{(}1-\alpha+\frac{t_{0}^{\alpha}}{\Gamma(\alpha)}\Bigr{)}<1,$
(22)

then the problem (15) admits one and only one solution.

###### Proof.

Recall from (21) that $\vartheta^{n}=\sum_{j=0}^{n}e^{j}$. Then we have,

$\|\vartheta^{n}\|_{\mathcal{I}_{T}}\leq\|\vartheta^{0}\|_{\mathcal{I}_{T}}\sum_{j=0}^{n}\biggl{[}\frac{(p^{*}-1)k^{p^{*}-2}\xi}{B(\alpha)}\Bigl{(}1-\alpha+\frac{t^{\alpha}}{\Gamma(\alpha)}\Bigr{)}\biggr{]}^{j}.$

For $t=t_{0}$, the above relation reads

$\|\vartheta^{n}\|_{\mathcal{I}_{T}}\leq\|\vartheta^{0}\|_{\mathcal{I}_{T}}\sum_{j=0}^{n}\biggl{[}\frac{(p^{*}-1)k^{p^{*}-2}\xi}{B(\alpha)}\Bigl{(}1-\alpha+\frac{t_{0}^{\alpha}}{\Gamma(\alpha)}\Bigr{)}\biggr{]}^{j}.$
(23)

Since we have

$\frac{(p^{*}-1)k^{p^{*}-2}\xi}{B(\alpha)}\Bigl{(}1-\alpha+\frac{t_{0}^{\alpha}}{\Gamma(\alpha)}\Bigr{)}<1,$

the geometric series

$\sum_{j=0}^{n}\biggl{[}\frac{(p^{*}-1)k^{p^{*}-2}\xi}{B(\alpha)}\Bigl{(}1-\alpha+\frac{t_{0}^{\alpha}}{\Gamma(\alpha)}\Bigr{)}\biggr{]}^{j}$

is convergent. Consequently, the series $(\vartheta^{n})$ exists and is
bounded for any $n\in\mathbb{N}$, and we have

$\lim_{n\to\infty}\|\vartheta^{n}\|_{\mathcal{I}_{T}}<\infty.$

We also consider the relation

$\mathcal{R}^{n}=\vartheta-\vartheta^{n}.$

Using the same technique as in the proof of Lemma 3, we can show the following
inequality

$\begin{split}\|\mathcal{R}^{n}\|_{\mathcal{I}_{T}}&\leq\frac{(p^{*}-1)k^{p^{*}-2}\xi}{B(\alpha)}\Bigl{(}1-\alpha+\frac{t^{\alpha}}{\Gamma(\alpha)}\Bigr{)}\|\mathcal{R}^{n-1}\|_{\mathcal{I}_{T}}\\\
&\leq\biggl{[}\frac{(p^{*}-1)k^{p^{*}-2}\xi}{B(\alpha)}\Bigl{(}1-\alpha+\frac{t^{\alpha}}{\Gamma(\alpha)}\Bigr{)}\biggr{]}^{n}\|\mathcal{R}^{0}\|_{\mathcal{I}_{T}}.\end{split}$

Now, taking the limit $n\to\infty$ in the above relation we find

$\lim_{n\to\infty}\|\vartheta-\vartheta^{n}\|_{\mathcal{I}_{T}}=0,\ \text{
i.e. }\ \lim_{n\to\infty}\vartheta^{n}=\vartheta.$

Therefore, the existence of a solution is proven.

Now we need to prove uniqueness. Let us assume that there are two solutions to
the reformulated system (15), namely by $\tilde{\vartheta}$ and
$\bar{\vartheta}$, then

$\begin{split}\|\tilde{\vartheta}(t)-\bar{\vartheta}(t)\|_{\mathcal{I}_{T}}&\leq(p^{*}-1)k^{p^{*}-2}\Big{\|}\frac{1-\alpha}{B(\alpha)}\bigl{(}\Psi(\tilde{\vartheta}(t))-\Psi(\bar{\vartheta}(t))\bigr{)}\\\
&\qquad+\frac{\alpha}{B(\alpha)\Gamma(\alpha)}\int_{0}^{t}\bigl{(}\Psi(\tilde{\vartheta}(s))-\Psi(\bar{\vartheta}(s))\bigr{)}(t-s)^{\alpha-1}\,ds\Big{\|}_{\mathcal{I}_{T}}\\\
&\leq(p^{*}-1)k^{p^{*}-2}\frac{1-\alpha}{B(\alpha)}\bigl{\|}\Psi(\tilde{\vartheta}(t))-\Psi(\bar{\vartheta}(t))\bigr{\|}_{\mathcal{I}_{T}}\\\
&\qquad+(p^{*}-1)k^{p^{*}-2}\frac{\alpha}{B(\alpha)\Gamma(\alpha)}\int_{0}^{t}\bigl{\|}\Psi(\tilde{\vartheta}(s))-\Psi(\bar{\vartheta}(s))\bigr{\|}_{\mathcal{I}_{T}}(t-s)^{\alpha-1}\,ds\\\
&\leq(p^{*}-1)k^{p^{*}-2}\frac{1-\alpha}{B(\alpha)}\xi\bigl{\|}\tilde{\vartheta}(t)-\bar{\vartheta}(t)\bigr{\|}_{\mathcal{I}_{T}}\\\
&\qquad+(p^{*}-1)k^{p^{*}-2}\frac{\alpha}{B(\alpha)\Gamma(\alpha)}\xi\int_{0}^{t}\bigl{\|}\tilde{\vartheta}(s)-\bar{\vartheta}(s)\bigr{\|}_{\mathcal{I}_{T}}(t-s)^{\alpha-1}\,ds\\\
&\leq(p^{*}-1)k^{p^{*}-2}\frac{1-\alpha}{B(\alpha)}\xi\|\tilde{\vartheta}-\bar{\vartheta}\|_{\mathcal{I}_{T}}+(p^{*}-1)k^{p^{*}-2}\frac{t^{\alpha}}{B(\alpha)\Gamma(\alpha)}\xi\|\tilde{\vartheta}-\bar{\vartheta}\|_{\mathcal{I}_{T}}.\end{split}$

This estimate shows that

$\|\tilde{\vartheta}-\bar{\vartheta}\|_{\mathcal{I}_{T}}\leq
0,\quad\text{i.e.}\quad\tilde{\vartheta}=\bar{\vartheta}.$

∎

##  5 Stability analysis

In this section we give an analysis of the Hyers-Ulam stability for the
reformulated problem (15). The study of stability of functional equations
originated from an open question in 1964 (cf. [13, 36]), which focused on the
stability of a group homomorphism. Given a group $\mathcal{G}$ and a metric
group $(\mathcal{G}^{\prime},d)$. Let $\varepsilon>0$, is there a $a>0$ such
that, if $f\colon\mathcal{G}\to\mathcal{G}^{\prime}$ satisfies the condition

$d\bigl{(}f(xy),f(x)f(y)\bigr{)}<a,\quad\text{for all }x,y\in\mathcal{G},$

then there exists a homomorphism $g\colon\mathcal{G}\to\mathcal{G}^{\prime}$
such that

$d\bigl{(}f(x),g(x)\bigr{)}<\varepsilon,\quad\text{for all }x\in\mathcal{G}.$

Essentially, if we have an almost homomorphism, then there is a very near
homomorphism with a small error. Next, we define the stability notion
mentioned above and prove the stability of the reformulated system (15).

######  Definition 4 (Hyers-Ulam stability).

The $\mathcal{ABC}$-FDE problem (15) is Hyers-Ulam stable if, for any given
$\mathscr{S}>0$, there exists a constant $\delta>0$ such that if

$\biggl{|}\vartheta_{i}(t)-\Phi_{p^{*}}\Bigl{[}\Phi_{p}(\vartheta_{0}^{i})+\frac{1-\alpha}{B(\alpha)}\Psi(\vartheta_{i}(t))+\frac{\alpha}{B(\alpha)\Gamma(\alpha)}\int_{0}^{t}\Psi(\vartheta_{i}(s))(t-s)^{\alpha-1}\,ds\Bigr{]}\biggr{|}\leq\mathscr{S},$
(24)

there exists $\varphi(t)=(\varphi_{1}(t),\varphi_{2}(t))$ satisfying

$\varphi_{i}(t)=\Phi_{p^{*}}\Bigl{[}\Phi_{p}(\vartheta_{0}^{i})+\frac{1-\alpha}{B(\alpha)}\Psi(\varphi_{i}(t))+\frac{\alpha}{B(\alpha)\Gamma(\alpha)}\int_{0}^{t}\Psi(\varphi_{i}(s))(t-s)^{\alpha-1}\,ds\Bigr{]},$
(25)

such that

$|\vartheta_{i}(t)-\varphi_{i}(t)|\leq\delta\mathscr{S}.$

The following theorem shows the stability of the problem in the sense defined
above.

######  Theorem 5.

The reformulated problem (15) is Hyers-Ulam stable.

###### Proof.

Let $\vartheta(t)=(\vartheta_{1}(t),\vartheta_{2}(t))=(x(t),y(t))$ be a
solution of (15) and $\varphi(t)=(\varphi_{1}(t),\varphi_{2}(t))$ be an
approximate solution and satisfying (19). Then, for $i=1,2,$ we have

$\begin{split}\bigl{|}\vartheta_{i}(t)-\varphi_{i}(t)\bigr{|}&=\bigg{|}\Phi_{p^{*}}\Bigl{[}\Phi_{p}(\vartheta_{0}^{i})+\frac{1-\alpha}{B(\alpha)}\Psi(\vartheta_{i}(t))+\frac{\alpha}{B(\alpha)\Gamma(\alpha)}\int_{0}^{t}\Psi(\vartheta_{i}(s))(t-s)^{\alpha-1}\,ds\Bigr{]}\\\
&\qquad-\Phi_{p^{*}}\Bigl{[}\Phi_{p}(\vartheta_{0}^{i})+\frac{1-\alpha}{B(\alpha)}\Psi(\varphi_{i}(t))+\frac{\alpha}{B(\alpha)\Gamma(\alpha)}\int_{0}^{t}\Psi(\varphi_{i}(s))(t-s)^{\alpha-1}\,ds\Bigr{]}\bigg{|}.\end{split}$

Using Lemma 2 leads to

$\begin{split}\bigl{|}\vartheta_{i}(t)-\varphi_{i}(t)\bigr{|}&\leq(p-1)k^{p^{*}-2}\Big{|}\frac{1-\alpha}{B(\alpha)}\Psi(\vartheta_{i}(t))+\frac{\alpha}{B(\alpha)\Gamma(\alpha)}\int_{0}^{t}\Psi(\vartheta_{i}(s))(t-s)^{\alpha-1}\,ds\\\
&\qquad-\frac{1-\alpha}{B(\alpha)}\Psi(\varphi_{i}(t))-\frac{\alpha}{B(\alpha)\Gamma(\alpha)}\int_{0}^{t}\Psi(\varphi_{i}(s))(t-s)^{\alpha-1}\,ds\Big{|}\\\
&\leq(p-1)k^{p^{*}-2}\frac{1-\alpha}{B(\alpha)}\bigl{|}\Psi(\vartheta_{i}(t))-\Psi(\varphi_{i}(t))\bigr{|}\\\
&\qquad+(p-1)k^{p^{*}-2}\frac{\alpha}{B(\alpha)\Gamma(\alpha)}\int_{0}^{t}\bigl{|}\Psi(\vartheta_{i}(s))-\Psi(\varphi_{i}(s))\bigr{|}(t-s)^{\alpha-1}\,ds.\end{split}$

Afterwards,

$\begin{split}\bigl{|}\vartheta_{i}(t)-\varphi_{i}(t)\bigr{|}&\leq(p-1)k^{p^{*}-2}\xi\frac{1-\alpha}{B(\alpha)}\bigl{\|}\vartheta_{i}-\varphi_{i}\bigr{\|}_{\mathcal{I}_{T}}\\\
&\qquad+(p-1)k^{p^{*}-2}\xi\frac{\alpha}{B(\alpha)\Gamma(\alpha)}\int_{0}^{t}\bigl{\|}\vartheta_{i}-\varphi_{i}\bigr{\|}_{\mathcal{I}_{T}}(t-s)^{\alpha-1}\,ds.\end{split}$

Then,

$\bigl{|}\vartheta_{i}(t)-\varphi_{i}(t)\bigr{|}\leq\Bigl{[}(p-1)k^{p^{*}-2}\xi\frac{1-\alpha}{B(\alpha)}+(p-1)k^{p^{*}-2}\xi\frac{t^{\alpha}}{B(\alpha)\Gamma(\alpha)}\Bigr{]}\varepsilon\leq\delta\mathscr{S},$

where
$\delta=\frac{(p-1)k^{p^{*}-2}}{B(\alpha)}\Bigl{[}1-\alpha+\frac{T^{\alpha}}{\Gamma(\alpha)}\Bigr{]}\xi$.
∎

After having shown the stability of the problem, we now want to analyze the
optimal control.

##  6 Optimal control analysis

In this section, we will discuss the optimality analysis for the SI model
(10)–(11). We define our objective functional as follows

$\mathscr{T}(u)=\int_{0}^{T}\Bigl{(}\mathscr{G}I^{2}(t)+\frac{\varrho}{2}u^{2}(t)\Bigr{)}\,dt,$
(26)

where $\varrho$ is a weight constant for the vaccination rate
$u\in\mathcal{U}_{ad}$, while $\mathscr{G}$ are the proportional weights
assigned to $S$ and $I$, respectively. In the following, we will first show
the existence of an optimal solution and then investigate the optimality
conditions.

###  6.1 Existence of an optimal solution

Using minimizing sequences, we prove the existence of optimal control. The
proof is based on the following lemma, which we introduce below.

######  Lemma 5 ([11, Proposition 2.1/Corollary 2.1]).

Let $\nu,\varphi\in C^{\infty}(\mathcal{I}_{T})$. Then

$\begin{split}\int_{0}^{T}({}^{\mathcal{ABC}}\mathcal{D}^{\alpha}_{t}\nu)\
\varphi\,dt&=-\int_{0}^{T}\nu\
({}^{\mathcal{ABC}}_{T}\mathcal{D}^{\alpha}_{t}\varphi)\,dt-\frac{B(\alpha)}{1-\alpha}\nu(0)\int_{0}^{T}\varphi(t)E_{\alpha,\alpha}[-\gamma
t^{\alpha}]\,dt\\\
&\quad+\frac{B(\alpha)}{1-\alpha}\varphi(T)\int_{0}^{T}\nu(t)E_{\alpha,\alpha}[-\gamma(T-t)^{\alpha}]\,dt.\end{split}$
(27)

######  Theorem 6.

The reformulated system (15) admits an optimal solution
$\vartheta^{*}=\vartheta(u^{*})\in L^{\infty}(\mathcal{I}_{T})$ that minimizes
the objective functional (26).

###### Proof.

Let $((\vartheta^{n},u^{n}))_{n}$ such as

$\mathscr{T}(u^{*})\leq\mathscr{T}(u^{n})\leq\mathscr{T}(u^{*})+\frac{1}{n},\quad\forall
n\in\mathbb{N}^{*},$

where $u^{n}\in U_{ad}$ and $\vartheta^{n}=\vartheta(u^{n})$ satisfying the
system

$\begin{cases}{}^{\mathcal{ABC}}\mathcal{D}^{\alpha}_{t}\bigl{[}\Phi_{p}(\vartheta^{n})\bigr{]}=\Psi(\vartheta^{n}),&\text{
in }\mathcal{I}_{T},\\\ \vartheta^{n}(0)=\vartheta_{0}.&\end{cases}$

Based on Theorem 1, we can say that $\vartheta^{n}$ is bounded. Consequently,
the second term $\Psi(\vartheta^{n})$ is also bounded. Thus, there exists a
positive constant $c$ such that

$\Bigl{\|}{}^{\mathcal{ABC}}\mathcal{D}^{\alpha}_{t}\bigl{[}\Phi_{p}(\vartheta^{n})\bigr{]}\Bigr{\|}_{2}\leq
c.$

Then there exists a subsequence of ($\vartheta^{n}$), again denoted by
($\vartheta^{n}$), such that

$\displaystyle{}^{\mathcal{ABC}}\mathcal{D}^{\alpha}_{t}\bigl{[}\Phi_{p}(\vartheta^{n})\bigr{]}\rightharpoonup\phi\text{
weakly in }\bigl{(}L^{2}(\mathcal{I}_{T})\bigr{)}^{2},$
$\displaystyle\vartheta^{n}\rightharpoonup\vartheta^{*}\text{ weakly-star in
}\bigl{(}L^{\infty}(\mathcal{I}_{T})\bigr{)}^{2},$ (28)
$\displaystyle\vartheta^{n}\rightharpoonup\vartheta^{*}\text{ weakly in
}\bigl{(}L^{2}(\mathcal{I}_{T})\bigr{)}^{2}.$

It’s worth noting that $D^{\prime}(\mathcal{I}_{T})$ is the dual of
$C_{0}^{\infty}(\mathcal{I}_{T})$. Then, for all $\varphi\in
C_{0}^{\infty}(\mathcal{I}_{T})$, we get

$\displaystyle\int_{0}^{T}\Phi_{p}(x^{n})\
\bigl{(}{}^{\mathcal{ABC}}_{T}\mathcal{D}^{\alpha}_{t}\varphi\bigr{)}\,dt$
$\displaystyle\longrightarrow\int_{0}^{T}\Phi_{p}(x^{*})\
\bigl{(}{}^{\mathcal{ABC}}_{T}\mathcal{D}^{\alpha}_{t}\varphi\bigr{)}\,dt,$
$\displaystyle\varphi(T)\int_{0}^{T}\Phi_{p}(x^{n})E_{\alpha,\alpha}[-\gamma(T-t)^{\alpha}]\,dt$
$\displaystyle\longrightarrow\varphi(T)\int_{0}^{T}\Phi_{p}(x^{*})E_{\alpha,\alpha}[-\gamma(T-t)^{\alpha}]\,dt,$

and

$\displaystyle\int_{0}^{T}\Phi_{p}(y^{n})\
\bigl{(}{}^{\mathcal{ABC}}_{T}\mathcal{D}^{\alpha}_{t}\varphi\bigr{)}\,dt$
$\displaystyle\longrightarrow\int_{0}^{T}\Phi_{p}(y^{*})\
\bigl{(}{}^{\mathcal{ABC}}_{T}\mathcal{D}^{\alpha}_{t}\varphi\bigr{)}\,dt,$
$\displaystyle\varphi(T)\int_{0}^{T}\Phi_{p}(y^{n})E_{\alpha,\alpha}[-\gamma(T-t)^{\alpha}]\,dt$
$\displaystyle\longrightarrow\varphi(T)\int_{0}^{T}\Phi_{p}(y^{*})E_{\alpha,\alpha}[-\gamma(T-t)^{\alpha}]\,dt.$

By Lemma 5, we obtain

$\displaystyle{}^{\mathcal{ABC}}\mathcal{D}^{\alpha}_{t}\bigl{[}\Phi_{p}(x^{n})\bigr{]}\rightharpoonup{}^{\mathcal{ABC}}\mathcal{D}^{\alpha}_{t}\bigl{[}\Phi_{p}(x^{*})\bigr{]}\text{
weakly in }D^{\prime}(\mathcal{I}_{T}),$
$\displaystyle{}^{\mathcal{ABC}}\mathcal{D}^{\alpha}_{t}\bigl{[}\Phi_{p}(y^{n})\bigr{]}\rightharpoonup{}^{\mathcal{ABC}}\mathcal{D}^{\alpha}_{t}\bigl{[}\Phi_{p}(y^{*})\bigr{]}\text{
weakly in }D^{\prime}(\mathcal{I}_{T}).$

Expressing $x^{n}y^{n}-x^{*}y^{*}=(x^{n}-x^{*})y^{n}+x^{*}(y^{n}-y^{*})$,
based on (28) and the boundedness of $(x^{n})$ and $(y^{n})$, it follows that
$x^{n}y^{n}\to x^{*}y^{*}$ in $L^{2}(\mathcal{I}_{T})$. In addition, we
observe $u^{n}\to u^{*}$ in $L^{2}(\mathcal{I}_{T})$ along a subsequence of
$(u^{n})$, denoted again by $(u^{n})$. Using the closeness and convexity of
$U_{ad}$, it can be deduced that $U_{ad}$ is weakly closed. Consequently,
$u^{*}\in U_{ad}$, and similarly, as described earlier, $u^{n}x^{n}\to
u^{*}x^{*}$ in $L^{2}(\mathcal{I}_{T})$.

Continuing, we can take the limit in the system satisfied by $\vartheta^{n}$
as $n\to\infty$, which leads to the conclusion that the optimal solution for
(10)–(11) is given by $(\vartheta^{*},u^{*})$. ∎

###  6.2 Optimality conditions

Let
$\vartheta^{\varepsilon}=(x^{\varepsilon},y^{\varepsilon})=(x,y)(u^{\varepsilon})$
and $\vartheta^{*}=(x^{*},y^{*})=(x,y)(u^{*})$ be the approximate and optimal
solutions of (10)–(11) and (26) respectively, where
$u^{\varepsilon}=u^{*}+\varepsilon u\in U_{ad},\ \forall u\in U_{ad}$. We
subtract the system associated with $\vartheta^{*}$ from the one corresponding
to
$\vartheta^{\varepsilon}=\vartheta^{*}+\varepsilon\mathscr{Z}^{\varepsilon}$
(with
$\mathscr{Z}^{\varepsilon}=(\mathscr{Z}_{1}^{\varepsilon},\mathscr{Z}_{2}^{\varepsilon})$),
we get

$\begin{cases}{}^{\mathcal{ABC}}\mathcal{D}^{\alpha}_{t}\bigl{[}\Phi_{p}(\vartheta^{\varepsilon})-\Phi_{p}(\vartheta^{*})\bigr{]}=\Psi(\vartheta^{\varepsilon})-\Psi(\vartheta^{*}),&\text{
in }\mathcal{I}_{T},\\\ \mathscr{Z}^{\varepsilon}(0)=0.\end{cases}$ (29)

Since $\Phi_{p}(\vartheta^{\varepsilon}(0))-\Phi_{p}(\vartheta^{*}(0))=0$,
then applying (5) to the first equation of (29), we have

$\Phi_{p}(\vartheta^{\varepsilon})-\Phi_{p}(\vartheta^{*})={}^{\mathcal{AB}}I^{\alpha}\bigl{(}\Psi(\vartheta^{\varepsilon})-\Psi(\vartheta^{*})\bigr{)}.$

Remember that $\Phi_{p}$ is a regular function. Using the mean value theorem,
we find that there exists a constant $\zeta$ such that

$\displaystyle\Phi_{p}(\vartheta^{\varepsilon})-\Phi_{p}(\vartheta^{*})$
$\displaystyle=\Phi_{p}^{\prime}(\zeta)(\vartheta^{\varepsilon}-\vartheta^{*})=\varepsilon\Phi_{p}^{\prime}(\zeta)\mathscr{Z}^{\varepsilon}$
$\displaystyle=\varepsilon\Phi_{p}^{\prime}(\zeta)\mathscr{Z}^{\varepsilon}(t)-\varepsilon\Phi_{p}^{\prime}(\zeta)\mathscr{Z}^{\varepsilon}(0)$
$\displaystyle={}^{\mathcal{AB}}I^{\alpha}\bigl{(}\varepsilon\Phi_{p}^{\prime}(\zeta){}^{\mathcal{ABC}}\mathcal{D}^{\alpha}_{t}[\mathscr{Z}^{\varepsilon}(t)]\bigr{)}.$

Thus,

$\Phi_{p}^{\prime}(\zeta){}^{\mathcal{ABC}}\mathcal{D}^{\alpha}_{t}\bigl{[}\mathscr{Z}^{\varepsilon}(t)\bigr{]}=\frac{\Psi(\vartheta^{\varepsilon})-\Psi(\vartheta^{*})}{\varepsilon}.$

On the other side

$\frac{\Psi(\vartheta^{\varepsilon})-\Psi(\vartheta^{*})}{\varepsilon}=\mathscr{N}^{\varepsilon}\mathscr{Z}^{\varepsilon}+\mathscr{Y}u,$

with

$\mathscr{N}^{\varepsilon}=\begin{pmatrix}\mu-\beta
y^{\varepsilon}-u^{\varepsilon}-\eta&\mu-\beta x^{*}\\\ \beta
y^{\varepsilon}+u^{\varepsilon}&\beta x^{*}-\eta\end{pmatrix}\ \text{ and }\
\mathscr{Y}=\begin{pmatrix}-x^{*}\\\ x^{*}\end{pmatrix}.$

The components of the matrix $\mathscr{N}^{\varepsilon}$ are uniformly bounded
with respect to $\varepsilon$. By using Lemma 5 and a similar approach
described in Subsection 6.1, as $\varepsilon$ tends to 0 in (29), we obtain

$\begin{cases}\Phi_{p}^{\prime}(\vartheta^{*}){}^{\mathcal{ABC}}\mathcal{D}^{\alpha}_{t}\bigl{[}\mathscr{Z}(t)\bigr{]}=\mathscr{N}\mathscr{Z}+\mathscr{Y}u,&\text{
in }\mathcal{I}_{T},\\\ \mathscr{Z}(0)=0,\end{cases}$ (30)

where

$\mathscr{N}=\begin{pmatrix}\mu-\beta y^{*}-u^{*}-\eta&\mu-\beta x^{*}\\\
\beta y^{*}+u^{*}&\beta x^{*}-\eta\end{pmatrix}.$

For $p=2$, according to Theorem 4, we can state that the problem (30) has a
unique solution. We now introduce
$\mathscr{P}=(\mathscr{P}_{1},\mathscr{P}_{2})$ in such a way that

$\Phi_{p}^{\prime}(\vartheta^{*})\int_{0}^{T}\bigl{(}{}^{\mathcal{ABC}}\mathcal{D}^{\alpha}_{t}[\mathscr{Z}(t)]\bigr{)}\mathscr{P}\,dt=\int_{0}^{T}(\mathscr{N}\mathscr{Z}+\mathscr{Y}u)\mathscr{P}\,dt.$

By Lemma 5, we have

$\int_{0}^{T}({}^{\mathcal{ABC}}\mathcal{D}^{\alpha}_{t}\mathscr{Z})\
\mathscr{P}\,dt=-\int_{0}^{T}\mathscr{Z}\
({}^{\mathcal{ABC}}_{T}\mathcal{D}^{\alpha}_{t}\mathscr{P})\,dt+\frac{B(\alpha)}{1-\alpha}\mathscr{P}(T)\int_{0}^{T}\mathscr{Z}(t)E_{\alpha,\alpha}[-\gamma(T-t)^{\alpha}]\,dt.$

Consequently, the dual system corresponding to (10)–(11) can be formulated as

$\begin{cases}-\Phi_{p}^{\prime}(\vartheta^{*}){}^{\mathcal{ABC}}_{T}\mathcal{D}^{\alpha}_{t}\mathscr{P}-\mathscr{N}\mathscr{P}=\mathscr{W}^{*}\mathscr{W}\vartheta^{*},&\text{
in }\mathcal{I}_{T},\\\
\mathscr{P}(T)=\mathscr{W}^{*}\mathscr{W}\vartheta^{*}(T),\end{cases}$ (31)

with

$\mathscr{W}=\begin{pmatrix}0&0\\\ 0&1\end{pmatrix}.$

We can show that the problem (31) has a unique solution (Theorem 4 with
$p=2$).

######  Theorem 7.

Let $\mathscr{P}$ be the solution of (31). Then,

$u^{*}=\max\biggl{\\{}0,\min\Bigl{(}\frac{\mathscr{P}_{1}^{*}-\mathscr{P}_{2}^{*}}{\varrho}S^{*},1\Bigr{)}\biggr{\\}}.$
(32)

###### Proof.

Let $(S,I)$ be the solution of (10)–(11). The Hamiltonian associated with (26)
is defined by

$\begin{split}\mathscr{H}(S,I,\mathscr{P}_{1},\mathscr{P}_{1},u)&=\
\mathscr{G}I^{2}(t)+\frac{\varrho}{2}u^{2}(t)\\\ &\qquad+\mathscr{P}_{1}(\mu
N-\beta SI-uS-\eta S)+\mathscr{P}_{2}(\beta SI+uS-\eta I),\end{split}$ (33)

where $(\mathscr{P}_{1},\mathscr{P}_{2})$ is the adjoint variable. Following
[24], we obtain the necessary optimality conditions for (10)–(11) and (26)
such that

$\begin{split}&{}^{\mathcal{ABC}}\mathcal{D}^{\alpha}_{t}\bigl{[}\Phi_{p}(S(t))\bigr{]}=\frac{\partial\mathscr{H}}{\partial\mathscr{P}_{1}}(S,I,\mathscr{P}_{1},\mathscr{P}_{1},u),\\\
&{}^{\mathcal{ABC}}\mathcal{D}^{\alpha}_{t}\bigl{[}\Phi_{p}(I(t))\bigr{]}=\frac{\partial\mathscr{H}}{\partial\mathscr{P}_{2}}(S,I,\mathscr{P}_{1},\mathscr{P}_{1},u),\\\
&\frac{\partial\mathscr{H}}{\partial
u}(S^{*},I^{*},\mathscr{P}_{1}^{*},\mathscr{P}_{1}^{*},u^{*})=0.\end{split}$
(34)

Since $\frac{\partial\mathscr{H}}{\partial
u}(S^{*},I^{*},\mathscr{P}_{1}^{*},\mathscr{P}_{1}^{*},u^{*})=0$, then

$u^{*}=\max\biggl{\\{}0,\min\Bigl{(}\frac{\mathscr{P}_{1}^{*}-\mathscr{P}_{2}^{*}}{\varrho}S^{*},1\Bigr{)}\biggr{\\}}.$

∎

##  7 Numerical simulations

In this section, we present a numerical approach to solving the optimality
system over the entire time interval $\mathcal{I}_{T}$. Based on [14, 26, 35],
we provide specific (positive) parameters and initial conditions, as detailed
in Table 2, assuming the boundedness constraint (16).

Table 2: Initial conditions and parameters values. Symbol | Description | Value
---|---|---
$S_{0}$ | Initial susceptible individuals | 10
$I_{0}$ | Initial infected individuals | 4
$\mu$ | Birth rate | 0.04
$\beta$ | Effective contact rate | 0.09
$\eta$ | Natural mortality rate | 0.04
$T$ | Final time | 200

###  7.1 Forward-Backward Sweep Method

We show numerical simulations related to the above optimal control problem. A
MatLab code has been developed and several simulations have been performed
using different data. The optimality system is solved based on an iterative
discrete scheme that converges after an appropriate test similar to the one
related to the Forward-Backward Sweep Method (FBSM). Due to transversality
conditions, the state model is solved forward in time, followed by the
backward solution of the adjoint system. Then, the optimal control values are
updated using the state and adjoint variables obtained in the previous steps.
This iterative process is repeated until a tolerance criterion is satisfied.
The algorithm is summarized in Algorithm 1.

######  Property 1 ([11]).

Let $\mathscr{X}\colon\mathcal{I}_{T}\to\mathbb{R}$. Then,

${}^{\mathcal{ABC}}_{T}\mathcal{D}^{\alpha}_{t}\mathscr{X}(t)={}^{\mathcal{ABC}}\mathcal{D}^{\alpha}_{t}\mathscr{X}(T-t).$
(35)

Algorithm 1 FBSM.

0: Let $\mathcal{S}=\\{t_{i}=1+i\delta_{t}/i=0,\cdots\mathcal{N}-1\\}$ be a
uniform subdivision on $\mathcal{I}_{T}$, where $\delta_{t}$ is the step size
of this subdivision and $\mathcal{N}$ is the number of sub-intervals.
$\delta=10^{-3}$ is tolerance and $Rel\\_Err=-1$ is the relative error.

1: Variables initialization: For the state problem $(S_{\rm old},I_{\rm
old})$, the adjoint system $(\mathscr{P}_{rmold_{1}},\mathscr{P}_{\rm
old_{2}})$, and the control $u_{\rm old}$.

2: while $Rel\\_Err<0$ do

3: Solve (10) with initial condition (11) given in Table 2.

4: Under (35), solve the (31) for $(\mathscr{P}_{\rm old_{1}},\mathscr{P}_{\rm
old_{2}})$ using the transversality conditions $(\mathscr{P}_{\rm
old_{1}}(T),\mathscr{P}_{\rm old_{2}}(T))$ and $(S,I)$.

5: Use (32) to update the control $u$.

6: Calculate relative errors:

  * $\psi_{1}=\delta\|S\|-\|S-S_{\rm old}\|$, $\psi_{2}=\delta\|I\|-\|I-I_{\rm old}\|$,

  * $\psi_{3}=\delta\|\mathscr{P}_{1}\|-\|\mathscr{P}_{1}-\mathscr{P}_{\rm old_{1}}\|$, $\psi_{4}=\delta\|\mathscr{P}_{2}\|-\|\mathscr{P}_{2}-\mathscr{P}_{\rm old_{1}}\|$,

  * $\psi_{5}=\delta\|u\|-\|u-u_{\rm old}\|$.

7: $Rel\\_Err=\min\\{\psi_{1},\psi_{2},\psi_{3},\psi_{4},\psi_{5}\\}$.

8: end while

###  7.2 Numerical approximations in the absence of vaccination

The Figures 2, 3 and 4 show numerical approximations in the absence of
vaccination (NAAV), considering different values of $\alpha$ and $p$.
Throughout the observation period, the number of susceptible individuals
decreases as the number of infected individuals rapidly increases. In each of
these three cases ($\alpha=1$, $\alpha=0.9$, and $\alpha=0.8$), we can observe
that the delayed spread of the disease is directly related to the values of
$p$ and $\alpha$. Thus, whenever $p$ is large, it results in the slow spread
of the disease. In addition, the values of $\alpha$ are also related to how
quickly the disease spreads, because natural derivative values result in the
rapid spread of the disease, while fractional values result in the slow spread
of the disease.

Figure 2: Numerical simulations of the SI model (10)–(11) for $\alpha=1$ without control. Figure 3: Numerical simulations of the SI model (10)–(11) for $\alpha=0.9$ without control. Figure 4: Numerical simulations of the SI model (10)–(11) for $\alpha=0.8$ without control. Table 3: The cost of $\mathscr{T}$ without control. | $\alpha=0.8$ | $\alpha=0.9$ | $\alpha=1$
---|---|---|---
$p=2$ | $1.2641e^{+06}$ | $1.2579e^{+06}$ | $1.2672e^{+06}$
$p=3$ | $1.0124e^{+06}$ | $1.0021e^{+06}$ | $1.0174e^{+06}$
$p=4$ | $8.6621e^{+05}$ | $8.3237e^{+05}$ | $8.8288e^{+05}$
$p=5$ | $7.8946e^{+05}$ | $7.2172e^{+05}$ | $8.2216e^{+05}$

###  7.3 Numerical approximations in the presence of vaccination

The vaccination strategy is currently being implemented, with the assumption
that vaccination will begin on the first day. From Figures 5, 6 and 7, we can
see that the vaccination strategy has spectacularly fought the spread of the
epidemic, where the proportion of infected individuals has been reduced and,
conversely, the number of susceptible individuals is constantly increasing. We
can say that our proposed strategy is effective in controlling the spread of
the epidemic.

Figure 5: Numerical simulations of the SI model (10)–(11) for $\alpha=1$ with control. Figure 6: Numerical simulations of the SI model (10)–(11) for $\alpha=0.9$ with control. Figure 7: Numerical simulations of the SI model (10)–(11) for $\alpha=0.8$ with control. Table 4: The cost of $\mathscr{T}$ with control. | $\alpha=0.8$ | $\alpha=0.9$ | $\alpha=1$
---|---|---|---
$p=2$ | $2.9641e^{+04}$ | $2.3167e^{+04}$ | $4.1000e^{+04}$
$p=3$ | $2.4244e^{+04}$ | $2.0200e^{+04}$ | $3.1549e^{+04}$
$p=4$ | $1.9380e^{+04}$ | $1.8588e^{+04}$ | $2.2178e^{+04}$
$p=5$ | $1.9117e^{+04}$ | $1.7859e^{+04}$ | $2.0106e^{+04}$

Note that, depending on the tables 3 and 4, we observe that $\mathscr{T}$
decreases under the effect of vaccination for different values of $p$ and
$\alpha$. The value of $\mathscr{T}$ is optimal for $p=5$ and $\alpha=0.9$.

##  8 Conclusion

In this manuscript, we have presented an innovative application of
epidemiological models, where the interactions between susceptible and
infected individuals are captured by a system of FDEs using the
$\mathcal{ABC}$ fractional derivative and the $p$ Laplacian operator. Our
investigation aims to contribute to more realistic models of disease spread in
certain scenarios. We established the existence of a unique bounded non-
negative solution to our biological model. We also proved the Hyers-Ulam
stability and optimality conditions. We validated our theoretical findings
through numerical simulations, where the results obtained show that the rapid
spread of the disease is associated with $\alpha$ taking natural values and
$p$ taking small values. Moreover, the slow spread of the disease results from
the fractional and large values of both $\alpha$ and $p$. In conclusion, the
spread of the disease has been successfully controlled.

## Declarations

### CRediT author statement

A. Zinihi: Conceptualization, Methodology, Software, Formal analysis,
Investigation, Writing – Original Draft, Writing – Review & Editing,
Visualization.

M. R. Sidi Ammi: Conceptualization, Methodology, Validation, Formal analysis,
Investigation, Writing – Original Draft, Writing – Review & Editing,
Supervision.

M. Ehrhardt: Validation, Formal analysis, Investigation, Writing – Original
Draft, Writing – Review & Editing, Project administration.

### Data availability

All information analyzed or generated, which would support the results of this
work are available in this article. No data was used for the research
described in the article.

### Conflict of interest

The authors declare that there are no problems or conflicts of interest
between them that may affect the study in this paper.

## References

  * [1] Abdeljawad, T., and Baleanu, D. Integration by parts and its applications of a new nonlocal fractional derivative with Mittag-Leffler nonsingular kernel.  J. Nonlin. Sci. Appl. 10, 03 (2017), 1098–1107.
  * [2] Agarwal, R. P., O’Regan, D., and Staněk, S. Positive solutions for Dirichlet problems of singular nonlinear fractional differential equations.  J. Math. Anal. Appl. 371, 1 (2010), 57–68.
  * [3] Altaf Khan, M., Ismail, M., Ullah, S., and Farhan, M. Fractional order SIR model with generalized incidence rate.  AIMS Mathematics 5, 3 (2020), 1856–1880.
  * [4] Anderson, R. M., and May, R. M. Infectious diseases of humans: Dynamics and control.  Oxford University Press, 1991.
  * [5] Atangana, A., and Baleanu, D. New fractional derivatives with nonlocal and non-singular kernel: Theory and application to heat transfer model.  Therm. Sci. 20, 2 (2016), 763–769.
  * [6] Bai, C. Z., and Fang, J. X. The existence of a positive solution for a singular coupled system of nonlinear fractional differential equations.  Appl. Math. Comput. 150, 3 (2004), 611–621.
  * [7] Bai, Z., and Qiu, T. Existence of positive solution for singular fractional differential equation.  Appl. Math. Comput. 215, 7 (2009), 2761–2767.
  * [8] Chai, G. Positive solutions for boundary value problem of fractional differential equation with $p$-Laplacian operator.  Boundary Value Problems 2012, 1 (2012).
  * [9] Chen, T., Liu, W., and Hu, Z. A boundary value problem for fractional differential equation with $p$-Laplacian operator at resonance.  Nonlin. Anal.: Real World Appl. 75, 6 (2012), 3210–3217.
  * [10] Din, A., Li, Y., Khan, F. M., Khan, Z. U., and Liu, P. On analysis of fractional order mathematical model of Hepatitis B using Atangana–Baleanu Caputo (ABC) derivative.  Fractals 30, 01 (2021).
  * [11] Djida, J. D., Mophou, G., and Area, I. Optimal control of diffusion equation with fractional time derivative with nonlocal and nonsingular Mittag-Leffler kernel.  J. Optim. Theory Appl. 182, 2 (2018), 540–557.
  * [12] Fernandez, A., and Baleanu, D. The mean value theorem and Taylor’s theorem for fractional derivatives with Mittag-Leffler kernel.  Adv. Diff. Eqs. 2018, 1 (2018).
  * [13] Forti, G. L. Hyers-Ulam stability of functional equations in several variables.  Aequat. Math. 50, 1–2 (1995), 143–190.
  * [14] Gatto, N. M., and Schellhorn, H. Optimal control of the SIR model in the presence of transmission and treatment uncertainty.  Math. Biosci. 333 (2021), 108539.
  * [15] Ghoshal, S. S., Junca, S., and Parmar, A. Fractional regularity for conservation laws with discontinuous flux.  Nonlin. Anal.: Real World Appl. 75 (2024), 103960.
  * [16] Han, J., and Liu, C. Global weak solution for a chemotaxis Navier–Stokes system with $p$-Laplacian diffusion and singular sensitivity.  Nonlin. Anal.: Real World Appl. 73 (2023), 103898.
  * [17] Hilfer, R. Applications of fractional calculus in physics.  World Scientific, 2000.
  * [18] Hoang, M. T., and Ehrhardt, M. Differential equation models for infectious diseases: Mathematical modeling, qualitative analysis, numerical methods and applications.  IMACM preprint 24/02, University of Wuppertal, Germany, 2024.
  * [19] Hosseininia, M., Heydari, M. H., and Razzaghi, M. Meshless local Petrov-Galerkin method for 2D fractional Fokker-Planck equation involved with the ABC fractional derivative.  Comput. Math. Appl. 125 (2022), 176–192.
  * [20] Huo, J., Zhao, H., and Zhu, L. The effect of vaccines on backward bifurcation in a fractional order HIV model.  Nonlin. Anal.: Real World Appl. 26 (2015), 289–305.
  * [21] Jena, R. M., Chakraverty, S., and Baleanu, D. SIR epidemic model of childhood diseases through fractional operators with Mittag-Leffler and exponential kernels.  Math. Comput. Simul. 182 (2021), 514–534.
  * [22] Keeling, M. J., and Rohani, P. Modeling infectious diseases in humans and animals.  Princeton University Press, 2007.
  * [23] Kermack, W. O., and McKendrick, A. G. A contribution to the mathematical theory of epidemics.  Proc. Royal Soc. A 115, 772 (1927), 700–721.
  * [24] Khan, A., Zarin, R., Humphries, U. W., Akgül, A., Saeed, A., and Gul, T. Fractional optimal control of COVID-19 pandemic model with generalized Mittag-Leffler function.  Adv. Diff. Eqs. 2021, 1 (2021).
  * [25] Khan, H., Chen, W., and Sun, H. Analysis of positive solution and Hyers‐Ulam stability for a class of singular fractional differential equations with $p$‐Laplacian in Banach space.  Math. Meth. Appl. Sci. 41, 9 (2018), 3430–3440.
  * [26] Lin, F., Muthuraman, K., and Lawley, M. An optimal control theory approach to non-pharmaceutical interventions.  BMC Infect. Dis. 10, 1 (2010).
  * [27] Maamar, M., Ehrhardt, M., and Tabharit, L. A Nonstandard Finite Difference Scheme for a Time-Fractional Model of Zika Virus Transmission.  Math. Biosci. Engrg. 21, 01 (2024), 924–962.
  * [28] Mittag-Leffler, G. M. Une généralisation de l’intégrale de Laplace-Abel.  CR Acad. Sci. Paris (Ser. II) 137 (1903), 537–539.
  * [29] Niimi, K., Harding, D. R. K., Holmes, A. R., Lamping, E., Niimi, M., Tyndall, J. D. A., Cannon, R. D., and Monk, B. C. Specific interactions between the Candida albicans ABC transporter Cdr1p ectodomain and ad‐octapeptide derivative inhibitor.  Molecul. Microbiol. 85, 4 (2012), 747–767.
  * [30] Okazawa, N., and Yokota, T. Global existence and smoothing effect for the complex Ginzburg–Landau equation with p-Laplacian.  J. Diff. Eqs. 182, 2 (2002), 541–576.
  * [31] Phukan, J., and Dutta, H. Dynamic analysis of a fractional order SIR model with specific functional response and Holling type II treatment rate.  Chaos Solitons & Fractals 175 (2023), 114005.
  * [32] Podlubny, I. Fractional differential equations: An introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications.  Elsevier, 1998.
  * [33] Pu, H., Li, S., Dong, M., Jiao, S., and Shang, Y. Numerical method for coupled thermal analysis of the regenerative cooling structure.  J. Thermophys. Heat Transf. 32, 2 (2018), 326–336.
  * [34] Samko, S. G., Kilbas, A., and Marichev, O. Fractional Integrals and Derivatives: Theory and Applications.  Gordon and Breach, 1993.
  * [35] Sidi Ammi, M. R., Zinihi, A., Raezah, A. A., and Sabbar, Y. Optimal control of a spatiotemporal SIR model with reaction–diffusion involving $p$-Laplacian operator.  Res. Phys. 52 (2023), 106895.
  * [36] Ulam, S. M. A collection of mathematical problems.  Interscience, New York, 1964.
  * [37] Vong, S. Positive solutions of singular fractional differential equations with integral boundary conditions.  Math. Comput. Modell. 57, 5–6 (2013), 1053–1059.
  * [38] Wang, L., Wang, J., Zhong, Q., and Cheng, K. Concentration of solutions for a fractional relativistic Schrödinger-Choquard equation with critical growth.  Nonlin. Anal. 230 (2023), 113233.
  * [39] Zhao, K. Generalized UH-stability of a nonlinear fractional coupling $({p}_{1},{p}_{2})$-Laplacian system concerned with nonsingular Atangana-Baleanu fractional calculus.  J. Inequal. Appl. 2023, 1 (2023).