# Effects of fractional derivatives in epidemic models

Caterina Balzotti1, Mirko D’Ovidio2, Anna Chiara Lai2, Paola Loreti2
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract.

We study epidemic Susceptible-Infected-Susceptible models in the fractional
setting. The novelty is to consider models in which the susceptible and
infected populations evolve according to different fractional orders. We study
a model based on Caputo derivative, for which we establish existence results
of the solutions. Also, we investigate a model based on Caputo-Fabrizio
operator, for which we provide existence of solutions and a study of the
equilibria. Numerical simulations for both models and a direct numerical
comparison are also provided.

###### Key words and phrases:

Caputo derivative; Caputo-Fabrizio operator; SIS models

###### 2020 Mathematics Subject Classification:

26A33; 34A34

1 Istituto per le Applicazioni del Calcolo,
Consiglio Nazionale delle Ricerche, Rome, Italy;
2 Dipartimento di Scienze di Base e Applicate per l’Ingegneria,
Sapienza Università di Roma, Rome, Italy

##  1\. Introduction

The interest of the scientific community in mathematical modeling for
epidemiology has grown considerably in recent years. The study of epidemic
models began in the early 1900s with the pioneering work of Kermack and
McKendrick [KM27]. Their idea was to divide the population into groups which
distinguish the individuals based on their status with respect to the
infection, giving rise to the compartmental modeling for epidemics. The
evolution in time of the disease is then described by a system of ordinary
differential equations for each considered class.

In this work we focus on the SIS (Susceptible-Infected-Susceptible) model
[Het89], describing infections which do not confer immunity to recovery from
illness, such as influenza and common cold. Such theory describes
compartmental models, where the population is divided into groups depending on
the state of individuals, that is with respect to disease, distinguishing
infected, susceptible. The use of mathematical models for epidemiology is
useful to predict the behavior of an infection and take strategic decisions in
emergency situations to limit the spread of the disease which is
microscopically modeled by the fractional order of the derivative.

In recent years, the use of fractional derivatives for epidemic models has
grown widely. The main advantage of fractional calculus is that it can
incorporate memory effects into the model. Moreover, fractional models have an
extra degree of freedom compared to classical models, which is particularly
useful for fitting real data when available. We refer to [CLYL21] for a recent
review of fractional epidemic models.

In this paper we consider two fractional SIS models. One model is based on
Caputo derivative, for which we establish existence results of the solutions
and provide numerical simulations. The novelty is to allow the susceptible and
infected population evolve according to different fractional orders. The other
model is based on Caputo-Fabrizio operator. Here we let the susceptible
population evolve according to the Caputo-Fabrizio fractional operator,
whereas the infected population dynamics is based on ordinary differential
equation. In this case, we rewrite the system as a system of ordinary
differential equations, we study the equilibria and present some simulations.

More precisely, let $\alpha,\alpha_{1},\alpha_{2}\in(0,1)$. We consider the
initial-value problem for Caputo derivative with different orders

(1)
$\begin{dcases}{{}^{C}D}^{\alpha_{1}}_{t}S(t)=-\beta\frac{S(t)}{S(t)+I(t)}I(t)+\gamma
I(t)\\\ {{}^{C}D}^{\alpha_{2}}_{t}I(t)=\beta\frac{S(t)}{S(t)+I(t)}I(t)-\gamma
I(t)\\\ S(0)=S_{0}\\\ I(0)=I_{0}.\end{dcases}$

and the initial-value problem for Caputo-Fabrizio operator

(2)
$\begin{cases}{{}^{CF}D}^{\alpha}_{t}S=-\left(\frac{\beta}{S+I}S-\gamma\right)I\\\
I^{\prime}=\left(\frac{\beta}{S+I}S-\gamma\right)I\\\ S(0)=S_{0}\\\
I(0)=I_{0}.\end{cases}$

where $\beta,\gamma>0$, $S_{0},I_{0}\geq 0$.

Here, in the formula (1) above we denote the Caputo derivative [Cap08] by

$\displaystyle{{}^{C}D}^{\alpha_{i}}_{t}u(t)=\frac{1}{\Gamma(1-\alpha_{i})}\int_{0}^{t}\frac{du}{ds}(s)\,(t-s)^{-\alpha_{i}}\,ds,\quad
i=1,2$

where $\Gamma(\cdot)$ is the gamma function defined as

$\displaystyle\Gamma(z)=\int_{0}^{\infty}e^{-s}s^{z-1}\,ds,\quad z>0.$

The Caputo derivative is well-defined for a function

$\displaystyle u\in C(0,T),\,u^{\prime}(s)(t-s)^{-\alpha_{i}}\in
L^{1}((0,t)),\quad\forall t\in(0,T),\quad i=1,2.$

The Caputo-Fabrizio operator in (2) is defined by

$\displaystyle{{}^{CF}D}^{\alpha}_{t}u(t)=\frac{M(\alpha)}{1-\alpha}\int_{a}^{t}u^{\prime}(\tau)e^{-\frac{\alpha}{1-\alpha}(t-\tau)}d\tau$

where $M(\alpha)$ is a non-negative scaling factor satisfying $M(0)=M(1)=1$

The case $\alpha_{1}=\alpha_{2}=1$ in (1) has been extensively studied in the
literature. Beginning with the models introduced in [Ver38, Ver45, Ver47] in
which the logistic equation is used to model population dynamics, a numerous
researchers work in the field.

The problem to find a solution to the fractional logistic equation attracted
many authors and many works on this topic have been written only providing
approximations for that solution. A contribution to this discussion is given
in [DL18] based on a series representation of the solution which involves
Euler’s numbers. However, this approach is not applicable in this new context
since it is based on the following property of the Caputo derivative

$\displaystyle{{}^{C}D}^{\alpha}_{t}u(t)=0\quad\textrm{if and only if}\quad\
u(t)=const\quad\forall t\geq 0,\quad\alpha\in(0,1).$

Indeed, our effort is to study the fractional indices of derivation in which
$\alpha_{1}$ may differ from $\alpha_{2}$. We mainly study existence of
solution and we propose numerical test based on the method presented in
[GLM20] that we illustrate by several pictures. The numerical tests are in
agreement with the theory developed in the present paper, and they offer new
perspectives (e.g. symmetries emerging in the evolution of the total
population $N(t)=I(t)+S(t)$) for future works.

This is not only a challenging problem from the mathematical point of view
since it is of interest in many applications in which we may act in different
ways in order to slowdown the process.

Recently the Caputo-Fabrizio operator [CF15] has attracted many researchers.
The peculiarity of this fractional operator is the presence of non-singular
kernel in contrast with the Caputo derivative in which singular kernel appears
in the definition. However, the Caputo-Fabrizio operator can be used to model
processes with memory.

Our work is to consider (2) in which we are able to reduce the problem to an
ordinary differential equation. Here we obtain the solution of the
differential equation and their equilibria. The analysis is concluded by
numerical simulations. Finally, in order to point out the effects of the
Caputo and the Caputo-Fabrizio differential operators, we show a direct
comparison between the evolutions of the numerical solutions to the system (1)
and the system (2).

We conclude with the plan of the paper. In Section 2 we introduce the SIS
model with fractional Caputo derivatives with different orders. In Section
2.1, we investigate the existence of solutions. To this end, we adopt a
constructive approach whose ideas are borrowed from the Carathéodory existence
theorem for ordinary differential equations, see [CL55, Chapter 2]. In Section
2.2 we complete the study with some numerical simulation. In Section 3 we
address a SIS model with fractional Caputo-Fabrizio operator. In Section 3.1
we establish existence of solution by rewriting the system as a system of
ordinary differential equations, and we characterize the associated
equilibria. Section 3.2 is devoted to numerical simulations, also including
numerical tests directly comparing the proposed Caputo and Caputo-Fabrizio
models. Finally in Section 4 we draw our conclusions.

##  2\. Caputo fractional epidemic models

In this section we investigate the properties of the Caputo fractional SIS
model with different fractional orders (1), which we rewrite here for reader’s
convenience.

(3)
$\begin{dcases}{{}^{C}D}^{\alpha_{1}}_{t}S(t)=-\beta\frac{S(t)}{S(t)+I(t)}I(t)+\gamma
I(t)\\\ {{}^{C}D}^{\alpha_{2}}_{t}I(t)=\beta\frac{S(t)}{S(t)+I(t)}I(t)-\gamma
I(t)\\\ S(0)=S_{0}\\\ I(0)=I_{0}.\end{dcases}$

with $\alpha_{1},\alpha_{2}\in[0,1]$. The particular case of fractional SIS
models of Caputo type with $\alpha_{1}=\alpha_{2}$ has been studied in several
works, see e.g. [BDL20, ES13, HOEK18] and references therein. The novelty of
our work is the use of different fractional indices. This approach has been
already proposed in [LWLT19], where the authors study an inverse problem to
calibrate the parameters for the dengue fever. The results obtained fit well
the real data, suggesting that mixed order fractional epidemic models are
needed in applications.

Here we establish existence results for the solutions of (3) and we present
some related numerical simulations.

We introduce the notation

$f(x,y):=\left(\gamma-\beta\frac{x}{x+y}\right)y$

so that above system (3) rewrites

$\begin{dcases}{{}^{C}D}^{\alpha_{1}}_{t}S(t)=f(S(t),I(t))\\\
{{}^{C}D}^{\alpha_{2}}_{t}I(t)=-f(S(t),I(t))\\\ S(0)=S_{0}\\\
I(0)=I_{0}.\end{dcases}$

###  2.1. Solutions

The existence of a solution to (3) in the time interval $[0,T)$ (with
$T\leq+\infty$) with (positive) initial data $(S_{0},I_{0})$ is implied by the
existence of a couple of absolutely continuous functions $(S(t),I(t))$
satisfying the Volterra-type fractional integral equation

(4)
$\begin{split}S(t)=S_{0}+\frac{1}{\Gamma(\alpha_{1})}\int_{0}^{t}f(S(s),I(s))(t-s)^{\alpha_{1}-1}ds\quad\text{and}\quad\\\
I(t)=I_{0}-\frac{1}{\Gamma(\alpha_{2})}\int_{0}^{t}f(S(s),I(s))(t-s)^{\alpha_{2}-1}ds.\end{split}$

for all $t\in(0,T)$. Fix $\alpha_{1},\alpha_{2}\in[0,1]$. In the following we
make use of the following function

$G(t;\alpha_{1},\alpha_{2}):=\max\left\\{\frac{t^{\alpha_{i}}}{\Gamma(\alpha_{i}+1)}|i=1,2\right\\}$

which satisfies

$G(t;\alpha_{1},\alpha_{2})\geq\frac{t^{\alpha_{i}}}{\Gamma(\alpha_{i}+1)}=\int_{0}^{t}\frac{(t-s)^{\alpha_{i}-1}}{\Gamma(\alpha_{i})}ds,\qquad\forall
t>0,\,i=1,2.$

In the next two results we establish some invariance properties for the field
$f$ in the cases $\gamma\geq\beta$ and $\gamma<\beta$, respectively.

Note that below we often make use of the symbols $(x_{0},y_{0})$ to denote
couples of positive real numbers and $x(t),y(t)$ to denote continuous real
valued functions. This choice is meant to lighten the notation and we point
out that, in the search of solutions of (3), $(x_{0},y_{0})$ plays the role of
initial data $(S_{0},I_{0})$ whereas $x(t)$ and $y(t)$ represent the
(approximated) evolution of the susceptible population $S(t)$ and infected
population $I(t)$, respectively.

######  Lemma 1.

Let $\alpha_{1},\alpha_{2}\in[0,1]$ and let $0<\beta\leq\gamma$ and fix
$\varepsilon\in(0,1)$. Let $T=T_{\varepsilon,\gamma}$ be the positive real
number such that

$G(T;\alpha_{1},\alpha_{2})(1+\varepsilon)\gamma\leq\varepsilon.$

Fix $x_{0},y_{0}>0$ and define

$Q_{x_{0},y_{0}}:=[x_{0},x_{0}+\varepsilon
y_{0}]\times[(1-\varepsilon)y_{0},(1+\varepsilon)y_{0}].$

Then for any couple of absolutely continuous functions $x,y$ such that
$(x(t),y(t))\in Q_{x_{0},y_{0}}$ for all $t\in[0,T]$, the associated functions

$\tilde{x}(t):=x_{0}+\int_{0}^{t}f(x(s),y(s))\frac{(t-s)^{\alpha_{1}-1}}{\Gamma(\alpha_{1})}ds$
$\tilde{y}(t):=y_{0}-\int_{0}^{t}f(x(s),y(s))\frac{(t-s)^{\alpha_{2}-1}}{\Gamma(\alpha_{2})}ds.$

also satisfy $(\tilde{x}(t),\tilde{y}(t))\in Q_{x_{0},y_{0}}$ for all
$t\in[0,T]$.

###### Proof.

Let $(x(t),y(t))\in Q_{x_{0},y_{0}}$ for all $t\in[0,T]$, In particular one
has $x(t)>0$ and $y(t)\in((1-\varepsilon)y_{0},(1+\varepsilon)y_{0})$ for all
$t\in[0,T]$. Then

$f(x(s),y(s))=\left(\gamma-\beta\frac{x(s)}{x(s)+y(s)}\right)y(s)<\gamma
y(s)\qquad\forall s\in[0,T]$

and this implies

$\displaystyle\tilde{x}(t)$
$\displaystyle=x_{0}+\int_{0}^{t}f(x(s),y(s))\frac{(t-s)^{\alpha_{1}-1}}{\Gamma(\alpha_{1})}ds\leq
x_{0}+\int_{0}^{t}\gamma
y(s)\frac{(t-s)^{\alpha_{1}-1}}{\Gamma(\alpha_{1})}ds$ $\displaystyle\leq
x_{0}+\gamma(1+\varepsilon)y_{0}\int_{0}^{t}\frac{(t-s)^{\alpha_{1}-1}}{\Gamma(\alpha_{1})}ds\leq
x_{0}+G(t)(1+\varepsilon)\gamma y_{0}\leq x_{0}+\varepsilon y_{0}.$

Since $\gamma\geq\beta$ then $f(x,y)>0$ for all $x,y>0$, therefore
$\tilde{x}(t)>x_{0}>0$ for all $t\in[0,T]$ and this proves that
$\tilde{x}(t)\in(0,x_{0}+\varepsilon y_{0})$ for all $t\in[0,T]$. On the other
hand,

$\displaystyle|\tilde{y}(t)-y_{0}|$
$\displaystyle\leq\int_{0}^{t}|f(x(s),y(s))|\frac{(t-s)^{\alpha_{2}-1}}{\Gamma(\alpha_{2})}ds<\int_{0}^{t}\gamma
y(s)\frac{(t-s)^{\alpha_{2}-1}}{\Gamma(\alpha_{2})}ds$
$\displaystyle\leq\gamma(1+\varepsilon)y_{0}\int_{0}^{t}\frac{(t-s)^{\alpha_{2}}}{\Gamma(\alpha_{2})}ds\leq
G(t)(1+\varepsilon)\gamma y_{0}\leq\varepsilon y_{0}$

and this concludes the proof. ∎

The next result deals with the case $\beta>\gamma$ and it posits some
invariance properties of $f$, similar to those in Lemma 1, in a time interval
$[0,T]$. The main difference with Lemma 1 is that in this case the time $T$
depends not only on system parameters ($\beta\leq\gamma,\alpha_{1}$ and
$\alpha_{2}$) but also on initial data.

######  Lemma 2.

Let $\alpha_{1},\alpha_{2}\in[0,1]$ and $\beta>\gamma>0$. Fix $x_{0},y_{0}>0$
and let
$\varepsilon=\varepsilon(x_{0},y_{0}):=\min\\{\frac{1}{2},\frac{1}{2}\frac{x_{0}}{y_{0}}\\}$
so that $\varepsilon\in(0,\frac{1}{2}]$ and

$x_{0}>\frac{1}{2}x_{0}\geq\varepsilon y_{0}.$

Let $T=T_{\beta,\gamma,x_{0},y_{0}}$ be such that

(5) $G(T;\alpha_{1},\alpha_{2})(\beta+\gamma)(1+\varepsilon)\leq\varepsilon$

and define

$Q_{x_{0},y_{0}}:=[x_{0}-\varepsilon
y_{0},x_{0}+\varepsilon]\times[(1-\varepsilon)y_{0},(1+\varepsilon)y_{0}].$

Then for any couple of continuous functions $x(t),y(t)$ such that
$(x(t),y(t))\in Q_{x_{0},y_{0}}$ for all $t\in[0,T]$, the associated functions

$\tilde{x}(t):=x_{0}+\int_{0}^{t}f(x(s),y(s))\frac{(t-s)^{\alpha_{1}-1}}{\Gamma(\alpha_{1})}ds$
$\tilde{y}(t):=y_{0}-\int_{0}^{t}f(x(s),y(s))\frac{(t-s)^{\alpha_{2}-1}}{\Gamma(\alpha_{2})}ds.$

also satisfy $(\tilde{x}(t),\tilde{y}(t))\in Q_{x_{0},y_{0}}$ for all
$t\in[0,T]$.

###### Proof.

Fix $x_{0},y_{0}>0$. Let $x(t),y(t)$ be two continuous functions satisfying
$(x(t),y(t))\in Q_{x_{0},y_{0}}$ for all $t\in[0,T]$ with
$T=T_{\beta,\gamma,x_{0},y_{0}}$ satisfying (5). Since

$|f(x(s),y(s)|\leq(\beta+\gamma)|y(s)|\qquad\forall s\in[0,T],$

then one has

$\displaystyle|\tilde{x}(t)-x_{0}|$
$\displaystyle=|\int_{0}^{t}f(x(s),y(s))\frac{(t-s)^{\alpha_{1}-1}}{\Gamma(\alpha_{1})}ds|\leq\int_{0}^{t}(\gamma+\beta)|y(s)|\frac{(t-s)^{\alpha_{1}-1}}{\Gamma(\alpha_{1})}ds|$
$\displaystyle\leq
G(T;\alpha_{1},\alpha_{2})(\beta+\gamma)(1+\varepsilon)y_{0}\leq\varepsilon
y_{0}.$

This implies $\tilde{x}(t)\in(x_{0}-\varepsilon y_{0},x_{0}+\varepsilon
y_{0})$ for all $t\in[0,T]$. Arguing as above, one also deduces

$\displaystyle|\tilde{y}(t)-y_{0}|$
$\displaystyle=|\int_{0}^{t}f(x(s),y(s))\frac{(t-s)^{\alpha_{2}-1}}{\Gamma(\alpha_{2})}ds|\leq\varepsilon
y_{0}.$

and this implies $\tilde{y}(t)\in[(1-\varepsilon)y_{0},(1+\varepsilon)y_{0}]$
for all $t\in[0,T]$. ∎

We are finally in position to state the main result of the section, namely a
global existence for the solutions of (3) in the case $\gamma\geq\beta$ and a
local existence result in the case $\beta>\gamma$.

######  Theorem 1.

Let $\alpha_{1},\alpha_{2}\in[0,1]$ and let $0<\beta\leq\gamma$. Then for
every initial data $S_{0},I_{0}>0$ the system (3) admists a positive solution
in $(0,+\infty)$.

If otherwise $0<\gamma<\beta$ then for every initial data $S_{0},I_{0}>0$ the
system (3) admists a positive solution in $(0,T]$ where
$T=T_{x_{0},y_{0},\beta,\gamma}>0$ satisfies

(6)
$\max\left\\{\frac{T^{\alpha_{1}}}{\Gamma(\alpha_{1})},\frac{T^{\alpha_{2}}}{\Gamma(\alpha_{2})}\right\\}=G(T;\alpha_{1},\alpha_{2})=\begin{cases}\frac{1}{3(\beta+\gamma)}&\text{if
}S_{0}\geq I_{0};\\\
\frac{S_{0}}{S_{0}+2I_{0}}\frac{1}{(\beta+\gamma)}&\text{if
}S_{0}<I_{0}.\end{cases}$

###### Proof.

In order to have a more light notation, let $x_{0}:=S_{0}>0$ and
$y_{0}:=I_{0}>0$. Also fix $\varepsilon:=1/2$. Define

$T:=\begin{cases}T_{\varepsilon,\beta,\gamma}&\text{ as in Lemma \ref{l1} if
$\gamma\geq\beta$}\\\ T_{x_{0},y_{0},\beta,\gamma}&\text{ as in \eqref{Test}
if $\gamma<\beta$}.\end{cases}$

Note that, if $\gamma<\beta$ then $T$ meets the hypothesis of Lemma 2.

Consider $Q_{x_{0},y_{0}}$ as in Lemma 1 if $\gamma\geq\beta$ and
$Q_{x_{0},y_{0}}$ as in Lemma 2 if $\gamma<\beta$. We first prove that (3)
admits a positive solution in $[0,T]$ – note that $T$ is independent from
$x_{0},y_{0}$ when $\gamma\geq\beta$. To this end, consider the sequence of
functions for $n\geq 1$

(7) $\tilde{x}_{n}(t):=\begin{cases}x_{0}&\text{if }{t\in[0,T/n]}\\\
x_{0}+\int_{0}^{t-T/n}f(\tilde{x}(s),\tilde{y}(s))\frac{(t-s)^{\alpha_{1}-1}}{\Gamma(\alpha_{1})}ds&\text{if
}{t\in(T/n,T]};\end{cases}$ (8) $\tilde{y}_{n}(t):=\begin{cases}y_{0}&\text{if
}{t\in[0,T/n]}\\\
y_{0}-\int_{0}^{t-T/n}f(\tilde{x}(s),\tilde{y}(s))\frac{(t-s)^{\alpha_{2}-1}}{\Gamma(\alpha_{2})}ds&\text{if
}{t\in(T/n,T]}.\end{cases}$

Set for brevity $\tilde{z}_{n}(t):=(\tilde{x}_{n}(t),\tilde{y}_{n}(t))$.

_Claim 1._ The sequence of functions $\\{\tilde{z}_{n}(t)\\}_{n\geq 1}$ is
well defined, equicontinuous and equibounded in $[0,T]$. In particular
$\tilde{z}_{n}(t)\in Q_{x_{0},y_{0}}$ for all $t\in[0,T]$ and for all $n\geq
1$.
We prove the claim by showing by induction that for all $n$ and for all
$k=1,\dots,n$, $\tilde{z}_{n}(t)$ is well defined in $t\in[0,kT/n]$,
$\tilde{z}_{n}(t)\in Q_{x_{0},y_{0}}$ for all $t\in[0,kT/n]$, and that

(9) $|z_{n}(t_{1})-z_{n}(t_{2})|\leq 2\parallel
f\parallel_{L^{\infty}(Q_{x_{0},y_{0}})}\omega(|t_{1}-t_{2}|)\qquad\forall
t_{1},t_{2}\in[0,kT/n]$

where $\omega$ is the modulus of continuity of the function
$G(\cdot;\alpha_{1},\alpha_{2})$. If $k=1$, then by the definitions in (7) and
(8) it follows that $z_{n}(t)=(x_{0},y_{0})$ for all $t\in[0,T/n]$ and the
base of the induction readily follows. Fix $k$ such that $1<k<n$ and assume
that $\tilde{z}_{n}(t)$ is well defined on $[0,kT/n]$, $\tilde{z}_{n}(t)\in
Q_{x_{0},y_{0}}$ for every $t\in[0,kT/n]$, and that (9) holds. Then using the
second lines of (7) and (8) we have that the definition of $\tilde{z}_{n}(t)$
continuously extends to $[0,(k+1)T/n]$ by setting for $t\in(kT/n,(k+1)T/n]$

$\tilde{z}_{n}(t):=\left(x_{0}+\int_{0}^{t-T/n}f(\tilde{z}_{n}(s))\frac{(t-s)^{\alpha_{1}}}{\Gamma(\alpha_{1}+1)}ds,y_{0}-\int_{0}^{t-T/n}f(\tilde{z}_{n}(s))\frac{(t-s)^{\alpha_{2}}}{\Gamma(\alpha_{2}+1)}ds\right),$

By applying Lemma 1 or Lemma 2 (according to the cases $\gamma\geq\beta$ and
$\gamma<\beta$, respectively) to $(x(s),y(s))=z(s):=\tilde{z}_{n}(s+(k-1)T/n)$
(which satisfies $(x(s),y(s))\in Q_{x_{0},y_{0}}$ for all $s\in[0,T/n]$ by
inductive hypothesis) we have that $\tilde{z}_{n}(t)\in Q_{x_{0},y_{0}}$ for
all $t\in[kT,(k+1)T/n]$. It is left to show that

(10) $|z_{n}(t_{1})-z_{n}(t_{2})|\leq 2\parallel
f\parallel_{L^{\infty}(Q_{x_{0},y_{0}})}\omega(|t_{1}-t_{2}|)\qquad\forall
t_{1},t_{2}\in[0,(k+1)T/n]$

For all $t_{1},t_{2}\in[T/n,(k+1)T/n]$, since $\tilde{z}_{n}(s)\in
Q_{x_{0},y_{0}}$ for all $s\in[0,(k+1)T/n]$, one has

$\displaystyle|\tilde{z}_{n}(t_{2})-\tilde{z}_{n}(t_{1})|$
$\displaystyle\leq\int_{t_{1}-T/n}^{t_{2}-T/n}|f(\tilde{z}_{n}(s))|\frac{(t-s)^{\alpha_{1}}}{\Gamma(\alpha_{1}+1)}ds+\int_{t_{1}-T/n}^{t_{2}-T/n}|\tilde{f}(z_{n}(s))|\frac{(t-s)^{\alpha_{2}}}{\Gamma(\alpha_{2}+1)}ds$
$\displaystyle\leq\parallel
f\parallel_{L^{\infty}(Q_{x_{0},y_{0}})}\int_{t_{1}-T/n}^{t_{2}-T/n}\frac{(t-s)^{\alpha_{1}}}{\Gamma(\alpha_{1}+1)}+\frac{(t-s)^{\alpha_{2}}}{\Gamma(\alpha_{2}+1)}ds$
$\displaystyle\leq 2\parallel
f\parallel_{L^{\infty}(Q_{x_{0},y_{0}})}|G(t_{2}-T/n;\alpha_{1},\alpha_{2})-G(t_{1}-T/n;\alpha_{1},\alpha_{2})|$
$\displaystyle\leq 2\parallel
f\parallel_{L^{\infty}(Q_{x_{0},y_{0}})}\omega(|t_{2}-t_{1}|).$

The case in which $t_{1},t_{2}\in[0,T/n]$ is trivial, because $\tilde{z}_{n}$
is constant in that interval. It is left to discuss the case in which
$t_{1}\in[0,T/n]$ and $t_{2}\in[T/n,(k+1)T/n]$. In this case
$|t_{2}-t_{1}|=t_{2}-t_{1}\geq t_{2}-T/n=|t_{2}-T/n|$. Then arguing as above
one gets

$\displaystyle|\tilde{z}_{n}(t_{2})-\tilde{z}_{n}(t_{1})|$
$\displaystyle=|\tilde{z}_{n}(t_{2})-(x_{0},y_{0})|\leq\int_{0}^{t_{2}-T/n}|f(\tilde{z}_{n}(s))|\frac{(t-s)^{\alpha_{1}}}{\Gamma(\alpha_{1}+1)}ds$
$\displaystyle\leq\int_{0}^{t_{2}-T/n}|f(\tilde{z}_{n}(s))|\frac{(t-s)^{\alpha_{1}}}{\Gamma(\alpha_{1}+1)}ds$
$\displaystyle\leq 2\parallel
f\parallel_{L^{\infty}(Q_{x_{0},y_{0}})}\omega(|t_{2}-T/n|)$
$\displaystyle\leq 2\parallel
f\parallel_{L^{\infty}(Q_{x_{0},y_{0}})}\omega(|t_{2}-t_{1}|).$

This concludes the proof of the inductive step and, consequently, of the Claim
1.

Now, by Claim 1 and by Ascoli-Arzela’s theorem, there exists a subsequence
$\\{\tilde{z}_{n_{k}}(t)\\}$ converging uniformly in $[0,T]$ to a continuous
limit function $z(t)=(S(t),I(t))$ satisfying $z(t)\in{Q_{x_{0},y_{0}}}$ –
recall indeed that $Q_{x_{0},y_{0}}\subset(0,+\infty)\times(0,+\infty)$ is a
compact set. This implies in particular $S(t)>0$ and $I(t)>0$ for all
$t\in[0,T]$. Since $f$ is continuous, then

$f(\tilde{z}_{n_{k}}(t))\to f(z(t))=f(S(t),I(t))\qquad\text{as
$k\to+\infty$},\quad\forall t\in[0,T]$

moreover,

$|f(\tilde{z}_{n_{k}}(t))|\leq\parallel
f\parallel_{L^{\infty}({Q_{x_{0},y_{0}}})}\quad\forall t\in[0,T].$

Then, by Lebesgue’s dominated convergence theorem, for every fixed $t\in(0,T]$
and for $i=1,2$

$\lim_{k\to\infty}\int_{0}^{t}f(\tilde{z}_{n_{k}}(s))\frac{(t-s)^{\alpha_{i}-1}}{\Gamma(\alpha_{i})}ds=\int_{0}^{t}f(z(s))\frac{(t-s)^{\alpha_{i}-1}}{\Gamma(\alpha_{i})}ds=\int_{0}^{t}f(S(s),I(s))\frac{(t-s)^{\alpha_{i}-1}}{\Gamma(\alpha_{i})}ds$

Therefore, for all $t\in(0,T]$, choosing $k$ sufficiently large to have
$t>T/n_{k}$, one has
$\tilde{z}_{n_{k}}(t)=(\tilde{x}_{n_{k}}(t),\tilde{y}_{n_{k}}(t))$ where
$\tilde{x}_{n_{k}}$ and $\tilde{y}_{n_{k}}$ satisfy

$\displaystyle\lim_{k\to+\infty}\tilde{x}_{n_{k}}$
$\displaystyle=\lim_{k\to+\infty}x_{0}+\int_{0}^{t}f(\tilde{z}_{n_{k}}(s))\frac{(t-s)^{\alpha_{1}}}{\Gamma(\alpha_{1}+1)}ds-\int^{t}_{t-T/n_{k}}f(\tilde{z}_{n_{k}}(s))\frac{(t-s)^{\alpha_{1}}}{\Gamma(\alpha_{1}+1)}ds$
$\displaystyle=x_{0}+\int_{0}^{t}f(S(s),I(s))\frac{(t-s)^{\alpha_{1}-1}}{\Gamma(\alpha_{1})}ds$

and

$\displaystyle\lim_{k\to+\infty}\tilde{y}_{n_{k}}$
$\displaystyle=\lim_{k\to+\infty}y_{0}-\int_{0}^{t}f(\tilde{z}_{n_{k}}(s))\frac{(t-s)^{\alpha_{2}-1}}{\Gamma(\alpha_{2})}ds-\int^{t}_{t-T/n_{k}}f(\tilde{z}_{n_{k}}(s))\frac{(t-s)^{\alpha_{2}-1}}{\Gamma(\alpha_{2})}ds$
$\displaystyle=y_{0}-\int_{0}^{t}f(S(s),I(s))\frac{(t-s)^{\alpha_{2}-1}}{\Gamma(\alpha_{2})}ds.$

On the other hand $(\tilde{x}_{n_{k}}(t),\tilde{y}_{n_{k}}(t))\to(S(t),I(t))$
as $k\to\infty$ for all $t\in(0,T]$, and recalling $x_{0}=S_{0}$ and
$y_{0}=I_{0}$ one deduces that

$S(t)=S_{0}+\int_{0}^{t}f(S(s),I(s))\frac{(t-s)^{\alpha_{1}-1}}{\Gamma(\alpha_{1})}ds$

and

$I(t)=I_{0}-\int_{0}^{t}f(S(s),I(s))\frac{(t-s)^{\alpha_{2}-1}}{\Gamma(\alpha_{2})}ds.$

Then $(S(t),I(t))$ is the required solution of (3) in $(0,T]$. If
$\beta>\gamma$ then we are done, because we proved the local existence of a
solution of (3). If $\beta\leq\gamma$, then we can iteratively extend
$(S(t),I(t))$ to $(0,+\infty)$. For instance by applying the result to the
initial datum $(x_{0},y_{0}):=(S(T/2),I(T/2))$ we then obtain a solution
defined in $[0,T/2+T]$ and so on.

∎

###  2.2. Numerical simulations

The numerical discretization of system (3) follows [GLM20]. Let us consider a
general equation

(11) ${{}^{C}D}^{\alpha}_{t}u(t)=f(u(t)).$

and consider a numerical grid which uniformly divides the time interval
$[0,T]$ into $N_{t}$ steps of length $\Delta t$. We denote by
$u^{n}=u(t^{n})$, with $t^{n}=n\Delta t$. Let $\alpha\in(0,1)$, then the
Caputo derivative can be approximated as

${{}^{C}D}^{\alpha}_{t}u^{n}=\frac{1}{\Gamma(2-\alpha)\Delta
t^{\alpha}}\Big{(}u^{n}-\sum_{j=0}^{n-1}C_{n,j}u^{j}\Big{)},$

with

$\displaystyle C_{n,0}=g(n),\qquad C_{n,j}=g(n-j)-g(n-(j-1))\quad\text{ for
$j=1,\dots,n-1$}$

and $g(r)=r^{1-\alpha}-(r-1)^{1-\alpha}\text{ for $r\geq 1$}$. The numerical
scheme to solve (11) is then given by

$u^{n+1}=\sum_{j=0}^{n-1}C_{n,j}u^{j}+\Gamma(2-\alpha)\Delta
t^{\alpha}f(u^{n}).$

Let us denote by $S^{n}=S(t^{n})$ and $I^{n}=I(t^{n})$. By applying this
discretization to system (3) we obtain

$\displaystyle S^{n+1}$
$\displaystyle=\sum_{j=0}^{n-1}C_{n,j}S^{j}+\Gamma(2-\alpha_{1})\Delta
t^{\alpha_{1}}f_{S}(S^{n},I^{n})$ $\displaystyle I^{n+1}$
$\displaystyle=\sum_{j=0}^{n-1}C_{n,j}S^{j}+\Gamma(2-\alpha_{2})\Delta
t^{\alpha_{2}}f_{I}(S^{n},I^{n}),$

with

$\displaystyle f_{S}(S,I)=-\beta\frac{S}{S+I}I+\gamma I\qquad\text{and}\qquad
f_{I}(S,I)=\beta\frac{S}{S+I}I-\gamma I.$

We show the evolution in time of $S(t)$, $I(t)$ and their sum $S(t)+I(t)$ as
$\alpha_{1},\alpha_{2}$ changes: in Figures 1 we considered a case in which
$\beta>\gamma$ and in 2 $\gamma<\beta$. Note in plots (a) and (b) that the
steepness of the solutions $S(t)$ and $I(t)$ in the long run are related to
the size of $\alpha_{1}$ and $\alpha_{2}$, respectively: this can be
interpreted as a time delay effect of the Caputo fractional operator. An
interesting phenomenon is also the lack of monotonicity of the solutions, in
contrast with ordinary SIS models.

The sum of the two classes, see plots (c), shows that $N(t)=S(t)+I(t)$ is in
general not monotone. For instance, in Figure 1, $N(t)$ first decreases and
then increases when $\alpha_{1}>\alpha_{2}$, while it first increases and then
decreases when $\alpha_{1}<\alpha_{2}$ (a symmetrical behavior emerges in
Figure 2). In case $\alpha_{1}=\alpha_{2}$ then the sum is constant and we
recover the theory developed in [BDL20, HOEK18]. A rigorous study of the
symmetries emerging in the simulations and, in particular the intersections of
all the functions $N(t)$ at the same time, is still under investigation.

(a)  (b)  (c)

Figure 1. Numerical solutions to (3) with $\beta=0.7$, $\gamma=0.2$, $I_{0}=2$
and $S_{0}=8$.

(a)  (b)  (c)

Figure 2. Numerical solutions to (3) with $\beta=0.1$, $\gamma=0.2$, $I_{0}=8$
and $S_{0}=2$.

##  3\. Caputo-Fabrizio fractional epidemic models

In this section we are concerned with the SIS model using Caputo-Fabrizio type
fractional derivatives. We refer to [HA20, KSK+21, MSK19, UKF+20] for some
examples of epidemic models based on the Caputo-Fabrizio fractional operator.

More precisely, here we study the the fractionary SIS system (2), which we
recall to be

(12)
$\begin{cases}{}^{\text{CF}}D^{\alpha}_{t}S=-\left(\frac{\beta}{S+I}S-\gamma\right)I\\\
I^{\prime}=\left(\frac{\beta}{S+I}S-\gamma\right)I\\\ S(0)=S_{0}\\\
I(0)=I_{0}.\end{cases}$

where $\beta,\gamma>0$, $S_{0},I_{0}\geq 0$, $I_{0}+S_{0}>0$. Also recall that
the Caputo-Fabrizio operator of order $\alpha\in[0,1]$ for a function $u\in
H^{1}((a,b))$, $a<b$ is

${}^{\text{CF}}D^{\alpha}_{t}u(t):=\frac{M(\alpha)}{1-\alpha}\int_{a}^{t}u^{\prime}(\tau)e^{-\frac{\alpha}{1-\alpha}(t-\tau)}d\tau.$

where $M(\alpha)$ is a non-negative scaling factor satisfying $M(0)=M(1)=1$.
The indentity

(13)
$\frac{d}{dt}\,^{\text{CF}}D^{\alpha}_{t}u(t)=\frac{M(\alpha)}{1-\alpha}u^{\prime}(t)-\frac{\alpha}{1-\alpha}\,^{\text{CF}}D^{\alpha}_{t}u(t)\quad\forall
t\in(a,b),\alpha\in[0,1).$

is crucial to reconduct (12) to a system of ordinary differential equations,
which is the first step of our investigation. We make the following key
assumption, relating the order of derivation with the system parameters:

(14) $\displaystyle\gamma<\frac{\alpha}{1-\alpha},$ (15) $\displaystyle
M(\alpha)\geq\alpha\quad\forall\alpha\in[0,1].$

Note that, for any $\gamma>0$ the condition (14) is satisfied by choosing
$\alpha$ sufficiently close to $1$ whereas (15) trivially holds choosing
$M(\alpha)\equiv 1$, which is a setting earlier explored and motivated in
[LN15].

###  3.1. Solutions and equilibria

As anticipated above, we begin by rewriting (12) as a system of ordinary
equations for $\alpha\in[0,1)$.

######  Theorem 2.

Let $\beta,\gamma>0$, $S_{0},I_{0}\geq 0$, $I_{0}+S_{0}>0$ and assume the
conditions (14) and (15). Let

$B_{\alpha}:=\frac{1}{2}(\alpha+M(\alpha)+(1-\alpha)(\beta-\gamma)),\qquad
C_{\alpha}:=M(\alpha)(\alpha-(1-\alpha)\gamma)$

and

$\displaystyle P_{\alpha}(S_{0},I_{0})$ $\displaystyle:=M(\alpha)S_{0}+\alpha
I_{0}+(1-\alpha)\left(\beta-\gamma-\frac{\beta}{S_{0}+I_{0}}I_{0}\right)I_{0}.$

Also define the function $g_{\alpha}:\mathbb{R}\to\mathbb{R}$

$g_{\alpha}(x):=\frac{-B_{\alpha}x+P_{\alpha}(S_{0},I_{0})/2+\sqrt{(B_{\alpha}x-P_{\alpha}(S_{0},I_{0})/2)^{2}-C_{\alpha}x^{2}+M(\alpha)P_{\alpha}(S_{0},I_{0})x)}}{M(\alpha)}.$

Then any couple non-negative absolutely functions continuous $(S(t),I(t))$ is
a solution of fractional system (12) if and only if

(16) $\begin{split}S(t)&=g_{\alpha}(I(t))\qquad\forall t\geq 0\end{split}$

and $I(t)$ solves

(17) $\begin{cases}I^{\prime}=\bigg{(}\beta-\gamma-\dfrac{\beta
I}{g_{\alpha}(I)+I}\bigg{)}I\\\ I(0)=I_{0}.\end{cases}$

###### Proof.

Preliminarly remark that, by (13) and by the first equation of (12)

$\displaystyle\frac{d}{dt}D^{\text{CF}}_{\alpha}S$
$\displaystyle=\frac{M(\alpha)}{1-\alpha}S^{\prime}-\frac{\alpha}{1-\alpha}D^{CF}_{\alpha}S$
$\displaystyle=\frac{M(\alpha)}{1-\alpha}S^{\prime}+\frac{\alpha}{1-\alpha}\left(\frac{\beta}{S+I}S-\gamma\right)I.$

By differentiating both sides of the first equation of (12), we deduce that
(12) is equivalent to the ordinary system

(18)
$\begin{cases}\frac{M(\alpha)}{1-\alpha}S^{\prime}+\frac{\alpha}{1-\alpha}\left(\frac{\beta}{S+I}S-\gamma\right)I=-\left(\left(\frac{\beta}{S+I}S-\gamma\right)I\right)^{\prime}\\\
I^{\prime}=\left(\frac{\beta}{S+I}S-\gamma\right)I\\\ S(0)=S_{0}\\\
I(0)=I_{0}.\end{cases}$

Now define

$\displaystyle P(t):=$ $\displaystyle M(\alpha)S(t)+\alpha
I(t)+(1-\alpha)\left(\frac{\beta}{S(t)+I(t)}S(t)-\gamma\right)I(t)$
$\displaystyle=$ $\displaystyle M(\alpha)S(t)+\alpha
I(t)+(1-\alpha)\left(\beta-\gamma-\frac{\beta}{S(t)+I(t)}I(t)\right)I(t)\quad\forall
t\geq 0$

and remark that (18) implies $P^{\prime}(t)=0$ and consequently
$P(t)=P(0)=P_{\alpha}(S_{0},I_{0})$ for all $t\geq 0$. Then, if $I(t)$ and
$S(t)$ are solutions of (18) (or equivalently, of (12)) then

(19) $M(\alpha)S(t)+\alpha
I(t)+(1-\alpha)\left(\beta-\gamma-\frac{\beta}{S(t)+I(t)}I(t)\right)I(t)=P_{\alpha}(S_{0},I_{0}).$

Solving above equation with respect to $S(t)$ and selecting the only possibly
non-negative solution, we obtain for all $t\geq 0$

$\begin{split}S(t)=&\frac{1}{M(\alpha)}\left(-B_{\alpha}I(t)+P_{\alpha}(S_{0},I_{0})/2\right.\\\
&\left.+\sqrt{(B_{\alpha}I(t)-P_{\alpha}(S_{0},I_{0})/2)^{2}-C_{\alpha}I^{2}(t)+M(\alpha)P_{\alpha}(S_{0},I_{0})I(t)}\right)\\\
=&g_{\alpha}(I(t)).\end{split}$

Incidentally notice that if $\alpha=1$ then $P_{1}(I_{0})=I_{0}+S_{0}=:N$,
$B_{1}=C_{1}=1$ we recover from above relation the classical identity
$S(t)=N-I(t)$. We check that $S(t)$, as a function $g_{\alpha}$ of $I(t)$, is
well defined. To this end set note that $g_{\alpha}(x)$ is defined in
$\mathbb{R}$, because

$\displaystyle(B_{\alpha}x-P_{\alpha}(S_{0},I_{0})/2)^{2}-C_{\alpha}x^{2}+M(\alpha)P_{\alpha}(S_{0},I_{0})x=$
$\displaystyle(B_{\alpha}^{2}-C_{\alpha})x^{2}-(B_{\alpha}-M_{\alpha})P_{\alpha}(S_{0},I_{0})x+P_{\alpha}^{2}(S_{0},I_{0})/4\geq
0$

for all $x\in\mathbb{R}$. More precisely, above inequality holds because the
discriminant of above polynomial reads

$(B_{\alpha}-M_{\alpha})^{2}P_{\alpha}^{2}(S_{0},I_{0})-(B_{\alpha}^{2}-C_{\alpha})P_{\alpha}^{2}(S_{0},I_{0})=-(1-\alpha)\beta
M(\alpha)P_{\alpha}^{2}(S_{0},I_{0})<0$

for all $\alpha\in[0,1)$, whereas we remarked above that if $\alpha=1$ then
$g_{\alpha}(x)=N-x$.

Plugging the identity $S(t)=g_{\alpha}(I(t))$ in the second equation of (12)
we obtain (18).

Finally we check that $S(t)$ and $I(t)$ are non-negative. By (16), $S(t)\geq
0$ if and only if
$I(t)\in[0,P_{\alpha}(S_{0},I_{0})/(\alpha-(1-\alpha)\gamma)]$ for all $t\geq
0$. Remark that this condition is satisfied for $t=0$, indeed we have

$\displaystyle P_{\alpha}(S_{0},I_{0})$ $\displaystyle\geq\alpha
I_{0}+(1-\alpha)(\beta-\gamma-\frac{\beta}{S_{0}+I_{0}}I_{0})I_{0}$
$\displaystyle\geq(\alpha-(1-\alpha)\gamma)I_{0}.$

Moreover if $I(t)=P_{\alpha}(S_{0},I_{0})/(\alpha-(1-\alpha)\gamma)$ then
$I^{\prime}(t)=-\gamma I(t)<0$ and, consequently $I(t)\leq
P_{\alpha}(S_{0},I_{0})/(\alpha-(1-\alpha)\gamma)$. Finally, remark that $0$
is an equilibrium for (18) and that the velocity field
$f_{\alpha}(I)=\bigg{(}\beta-\gamma-\dfrac{\beta I}{g_{\alpha}(I)+I}\bigg{)}I$
is locally Lipschitz continuous. Therefore, by the local uniqueness of the
solutions of (18) if $I_{0}>0$ then $I(t)>0$ for all $t$. ∎

####  3.1.1. Equilibria

We now characterize the equilibria and we study the asymptotic behavior of
(12).

######  Proposition 1.

Assume conditions (14) and (15) and let $S_{0},I_{0}\geq 0$. The equilibria of
the system (17) are $0$ and

$E_{\alpha}(S_{0},I_{0}):=\frac{(\beta-\gamma)P_{\alpha}(S_{0},I_{0})}{\alpha(\beta-\gamma)+M(\alpha)\gamma}.$

Setting $R:=\beta/\gamma$, $E_{\alpha}(I_{0})>0$ if and only if $R>1$. In this
case

$I(t)\to E_{\alpha}(S_{0},I_{0}),\quad
S(t)\to\frac{\gamma}{\beta-\gamma}E_{\alpha}(S_{0},I_{0})=\frac{\gamma
P_{\alpha}(S_{0},I_{0})}{\alpha(\beta-\gamma)+M(\alpha)\gamma}\qquad\text{as
}t\to+\infty.$

Finally, if $E_{\alpha}(S_{0},I_{0})\leq 0$, that is if $R\leq 1$, then

$I(t)\to 0,\quad S(t)\to P_{\alpha}(S_{0},I_{0})\qquad\text{as }t\to+\infty.$

###### Proof.

The first part of the claim follows by a direct computation, using in
particular (15) and the fact that $\beta,\gamma>0$ implies $R>0$. Let
$x_{R}=E_{\alpha}(S_{0},I_{0})$ if $R>1$ and $x_{R}=0$ if $R\leq 0$. We proved
that in Theorem 2 that $I(t)$ is non-negative, and consequently $[0,+\infty)$
is an invariant set for the dynamics (18). Moreover, by a direct computation
one can check that the function $V(x):=(x-x_{R})^{2}$ is a Lyapunov function
for (18) in $[0,+\infty)$ and, consequently, $x_{R}$ is a globally
asymptocally stable equilibrium for (18) and this concludes the proof. ∎

Note that, as in the classical case, the qualitative properties of the system
strongly depend on the reproduction number $R=\beta/\gamma$. The value
$E_{\alpha}(S_{0},I_{0})$ can be viewed as a fractional generalization of the
endemic equilibrium of classical SIS models, which is a stable equilibrium
when $R>1$.

The next result deals with the monotonicity and the asymptotic behaviour of
the function $N(t):=S(t)+I(t)$, which is not constant unless we are in the
ordinary case $\alpha=1$. As we show below, the asymptotic behaviour of $N(t)$
can be used for inverse problems, i.e., the goal to reconstruct the fractional
order $\alpha$ of (12) from the observed data on the long range.

######  Proposition 2.

Assume $S_{0},I_{0}>0$ and $\alpha\in[0,1)$. If $R>1$ and if
$I_{0}\in(0,E(S_{0},I_{0}))$ then $N(t):=I(t)+S(t)$ is a strictly increasing
function converging to

$N_{\alpha}(S_{0},I_{0}):=\frac{\beta}{\beta-\gamma}E_{\alpha}(S_{0},I_{0})=\frac{\beta
P_{\alpha}(S_{0},I_{0})}{\alpha(\beta-\gamma)+M(\alpha)\gamma}$

as $t\to+\infty$. If otherwise either $R>1$ and
$I_{0}>E_{\alpha}(S_{0},I_{0})$ or $R\leq 1$, then $N(t)$ is strictly
decreasing and tends to $P_{\alpha}(S_{0},I_{0})/M(\alpha)$ as $t\to+\infty$.

###### Proof.

We preliminary remark that $g_{\alpha}(x)$ is a convex function, indeed

(20)
$g_{\alpha}^{\prime\prime}(x)=\frac{P_{\alpha}^{2}(S_{0},I_{0})(1-\alpha)\beta}{4M(\alpha)\sqrt{(B_{\alpha}x-P_{\alpha}(S_{0},I_{0})/2)^{2}-C_{\alpha}x^{2}+M(\alpha)P_{\alpha}(S_{0},I_{0})}}\geq
0\quad\forall x\in\mathbb{R}.$

Moreover, by a direct computation,

$1+g^{\prime}_{\alpha}(I_{0})=1+\frac{2(1-\alpha)\beta
M(\alpha)(I_{0}+S_{0})^{3}P_{\alpha}(S_{0},I_{0})^{2}}{\left|I_{0}(I_{0}+S_{0})^{2}+(1-\alpha)I_{0}^{2}\beta\right|^{3}}>0$

and this, together with the convexity of $g_{\alpha}$, implies

(21) $1+g^{\prime}_{\alpha}(x)>0\quad\forall x\geq I_{0}.$

Also remark that, using the assumption (15), i.e., $M(\alpha)\geq\alpha$, we
obtain

$\displaystyle 1+g^{\prime}_{\alpha}(E_{\alpha}(S_{0},I_{0}))=$ $\displaystyle
1+\frac{M(\alpha)\left((1-\alpha)(\beta-\gamma)\gamma-\alpha\beta\right)}{(1-\alpha)(\beta-\gamma)^{2}+\beta
M(\alpha)}$ $\displaystyle\geq$ $\displaystyle
1+\frac{M(\alpha)\left((1-\alpha)(\beta-\gamma)\gamma-\beta\right)}{(1-\alpha)(\beta-\gamma)^{2}+\beta
M(\alpha)}$ $\displaystyle=$
$\displaystyle\frac{(1-\alpha)(\beta-\gamma)(M(\alpha)\gamma+\beta-\gamma)}{(1-\alpha)(\beta-\gamma)^{2}+\beta
M(\alpha)}$

therefore

(22) $\beta>\gamma\quad\Rightarrow\quad 1+g^{\prime}_{\alpha}(x)>0\quad\forall
x\geq E_{0}(S_{0},I_{0}).$

We conclude this preliminary study on $g_{\alpha}$ by noticing that

$1+g^{\prime}_{\alpha}(0)=(1-\alpha)(1-\beta+\gamma)$

therefore

(23) $\beta\leq\gamma\quad\Rightarrow\quad
1+g^{\prime}_{\alpha}(x)>0\quad\forall x\geq 0.$

Now, by Theorem 2

(24)
$N^{\prime}(t)=I^{\prime}(t)+S^{\prime}(t)=I^{\prime}(t)(1+g^{\prime}_{\alpha}(I(t)).$

Assume $R>1\text{ and }I_{0}\in(0,E(S_{0},I_{0}))$. Then

(25) $I^{\prime}(t)=\left(\beta-\gamma-\frac{\beta
I(t)}{I(t)+g_{\alpha}(I(t))}\right)I(t)>0\qquad\forall t\geq 0$

In particular $I(t)\geq I_{0}$ for all $t\geq 0$ and, in view of (21) and of
(24), we deduce $N^{\prime}(t)>0$ for all $t>0$, hence $N(t)$ is strictly
increasing.

Assume now that either $R\leq 1\text{ or }I_{0}\geq E(S_{0},I_{0})$.

First we discuss the case in which $R>1$ and $I_{0}>E(S_{0},I_{0})$. Then
$I^{\prime}(t)<0$ and, since $E_{\alpha}(S_{0},I_{0})$ is a (globally
asymptotically stable) equilibrium by Proposition 1, then
$I(t)>E_{\alpha}(S_{0},I_{0})$ for all $t\geq 0$. Then, in view of (22) and of
(24), we deduce $N^{\prime}(t)<0$ for all $t\geq 0$.

Finally, if $R\leq 1$ then $I^{\prime}(t)<0$. Since $I(t)>0$ for all $t>0$
then, in view of (23) and of (24), we deduce $N^{\prime}(t)<0$ for all $t\geq
0$ and this concludes the proof. ∎

From above result, we can estimate the fractional order $\alpha$ from the
asymptotic behaviour of the total population $N(t)$. Indeed, if $R>1$ and
$N(t)\to N_{\infty}$ then the corresponding fractional order $\alpha$ is the
solution of the equation

$\frac{\beta
P_{\alpha}(S_{0},I_{0})}{\alpha(\beta-\gamma)+M(\alpha)\gamma}=N_{\infty}.$

Assuming as in [LN15] $M(\alpha)\equiv 1$, and setting $N_{0}:=I_{0}+S_{0}$,
above equation reduces to the explicit formula

$\alpha=1-\frac{(N_{\infty}-N_{0})N_{0}\beta}{\beta I_{0}S_{0}-\beta
N_{0}(I_{0}+\gamma I_{0}-N_{\infty})-\gamma N_{0}N_{\infty}}.$

Note in particular that $\alpha=1$ if and only if $N_{\infty}=N_{0}$,
confirming the fact that, in the proposed mixed fractional model, $N(t)$ is
constant if and only if $\alpha=1$.

###  3.2. Numerical simulations

The numerical discretization of (12) is based on the results of Theorem 2. Let
us consider again the numerical grid introduced in Section 2.2. To compute the
discrete evolution in time of $I$ and $S$ we first solve system (17) by a
proper ODE solver and then we discretize (16). Specifically, we use the MATLAB
tool ode23t to compute $I^{n+1}$ and then, following (16), we obtain
$S^{n+1}=g_{\alpha}(I^{n+1})$.

(a)  (b)  (c)

Figure 3. Numerical solutions to (12) for $\beta=0.7$, $\gamma=0.2$,
$S_{0}=6$, $I_{0}=4$, and $M(\alpha)\equiv 1$ as $\alpha$ changes.

(a)  (b)  (c)

Figure 4. Numerical solutions to (12) with $\beta=0.1$, $\gamma=0.2$,
$S_{0}=6$, $I_{0}=4$, and $M(\alpha)\equiv 1$ as $\alpha$ changes.

In Figures 3 and 4 we show the results obtained as $\alpha$ changes in
$\\{0.2,0.4,0.6,0.8,1\\}$. The parameters $S_{0}$, $I_{0}$, $M(\alpha)$ and
$\gamma$ are the same in all simulations, fixed as $S_{0}=6$, $I_{0}=4$,
$M(\alpha)\equiv 1$ and $\gamma=0.2$. The parameter $\beta$, instead, is set
to $\beta=0.7$ in Figure 3, which implies $R>1$, while $\beta=0.1$ in Figure
4, and thus $R<1$. The dashed lines represent the equilibria estimated in
Proposition 1. Both $S$ and $I$ monotonically converge to their equilibria.
Note that the time delay induced by the fractional order emerges in the fact
that the smallest is $\alpha$, the slowest is the convergence of the related
solutions to equilibria. Finally remark that for $\alpha=1$ the sum
$N(t)=S(t)+I(t)$ is constant, since (12) coincides with the classical SIS
model, where if $\alpha<1$ then the monotonicity of $N(t)$ is in agreement
with Theorem 2.

(a)  (b)  (c)

Figure 5. Numerical solutions to (12) for $\alpha=0.5$, $\beta=0.7$,
$\gamma=0.2$ and $M(\alpha)\equiv 1$ as $I_{0}$ and $S_{0}$ change.

(a)  (b)  (c)

Figure 6. Numerical solutions to (12) for $\alpha=0.5$, $\beta=0.1$,
$\gamma=0.2$ and $M(\alpha)\equiv 1$ as $I_{0}$ and $S_{0}$ change.

In Figures 5 and 6 we fix $\alpha=0.5$ and we set $S_{0}=10-I_{0}$ and let
vary the initial data $I_{0}$ in $\\{2,4,6,8\\}$. We note that the qualitative
behavior of the solutions is in agreement with the theory developed here and
it is not much affected by initial data.

####  3.2.1. Comparison between Caputo SIS model and Caputo-Fabrizio SIS model

(a)  (b)  (c)

(d)  (e)  (f)

(g)  (h)  (i)

Figure 7. Numerical solutions to (3) and (12) for $\alpha=\alpha_{1}=0.2$
(top), $\alpha=\alpha_{1}=0.5$ (center) and $\alpha=\alpha_{1}=0.8$ (bottom)
with $\alpha_{2}=1$, $\beta=0.7$, $\gamma=0.2$, $S_{0}=6$, $I_{0}=4$ and
$M(\alpha)\equiv 1$.

(a)  (b)  (c)

(d)  (e)  (f)

(g)  (h)  (i)

Figure 8. Numerical solutions to (3) and (12) for $\alpha_{1}=0.2$ (top),
$\alpha_{1}=0.5$ (center) and $\alpha_{1}=0.8$ (bottom) with $\alpha_{2}=1$,
$\beta=0.1$, $\gamma=0.2$, $S_{0}=6$, $I_{0}=4$ and $M(\alpha)\equiv 1$.

We conclude this section with some tests pointing out the effect of the choice
of a particular fractional operator on SIS models. We directly compare the
numerical solution to the system with Caputo fractional SIS model (3) and
Caputo-Fabrizio fractional SIS model (12). To this end we fix the intial data
$S_{0}=6$ and $I_{0}=4$. We consider the set the fractional orders
$\alpha_{1}$ and $\alpha_{2}$ in (3), the fractional order $\alpha$ in (12)
and we set $\alpha_{1}=\alpha$, $\alpha_{2}=1$, letting $\alpha$ vary in the
set $\\{0.2,0.5,0.8\\}$. Figure 7 depicts the results in a case in which
$\gamma<\beta$, in particular $\gamma=0.2$ and $\beta=0.7$, whereas in Figure
8 we test a case in which $\beta<\gamma$, in particular $\gamma=0.2$ and
$\beta=0.1$. Note that the evolution of susceptible population $S(t)$ is
mostly affected by the change of differential operator and this effect is
augmented by choosing small fractional orders. Finally, note that the Caputo-
Fabrizio operator preserves the monotonicity properties of the ordinary SIS
case, while the Caputo operator used in (3) yields a more complex structure.

##  4\. Conclusions

We explored the effects of fractional differential operators on SIS epidemic
models. The main novelty of the present paper consists in letting the
susceptible and the infected population evolve with different orders of
fractional differential operators. We presented two fractional SIS models with
mixed fractional orders. One of them involves the Caputo fractional
derivative, characterized by a singular, power law kernel. The other proposed
model relies on the recently introduced Caputo-Fabrizio differential operator,
which has a non-singular, exponential kernel. For both models we established
existence results for the solutions and we conducted a qualitative analysis by
means of numerical simulations. In the case of the Caputo-Fabrizio operator,
under the assumption that the fractional behavior is restricted to the
susceptible population whereas the infected population evolves according to an
ordinary differential equation, we were able to move some further step forward
in the analysis of the system. Indeed we characterized the equilibria, noticed
their strong dependence on the reproduction number according to classical
theory, and proposed a method for inverse problems. Finally, we numerically,
directly compared the proposed Caputo and Caputo-Fabrizio SIS models, in order
to let emerge the effects of each particular differential operator on the
system.

We believe that the possibility of tuning the memory effects in a single
compartment of the population (susceptible/infected) by means of ad hoc
fractional orders, may provide finer and effective tools in data fitting and
mathematical modeling of epidemic dynamics.

Further possible extensions of the present work include the extension of the
proposed methods to other epidemic models, for instance the SIR model, and to
controlled, mixed fractional epidemic dynamics.

## References

  * [BDL20] Caterina Balzotti, Mirko D’Ovidio, and Paola Loreti.  Fractional sis epidemic models.  Fractal Fract., 4(3):44, 2020.
  * [Cap08] Michele Caputo.  Linear models of dissipation whose $Q$ is almost frequency independent. II.  Fract. Calc. Appl. Anal., 11(1):4–14, 2008.  Reprinted from Geophys. J. R. Astr. Soc. 13 (1967), no. 5, 529–539.
  * [CF15] Michele Caputo and Mauro Fabrizio.  A new definition of fractional derivative without singular kernel.  Progr. Fract. Differ. Appl, 1(2):1–13, 2015.
  * [CL55] Earl A Coddington and Norman Levinson.  Theory of ordinary differential equations.  Tata McGraw-Hill Education, 1955.
  * [CLYL21] Yuli Chen, Fawang Liu, Qiang Yu, and Tianzeng Li.  Review of fractional epidemic models.  Appl. Math. Model., 97:281–307, 2021.
  * [DL18] Mirko D’Ovidio and Paola Loreti.  Solutions of fractional logistic equations by Euler’s numbers.  Phys. A, 506:1081–1092, 2018.
  * [ES13] HAA El-Saka.  The fractional-order sir and sirs epidemic models with variable population size.  Math. Sci. Lett., 2(3):195, 2013.
  * [GLM20] Yoshikazu Giga, Qing Liu, and Hiroyoshi Mitake.  On a discrete scheme for time fractional fully nonlinear evolution equations.  Asymptotic Analysis, 120(1-2):151–162, 2020.
  * [HA20] M Higazy and Maryam Ahmed Alyami.  New caputo-fabrizio fractional order seiasqeqhr model for covid-19 epidemic transmission with genetic algorithm based control strategy.  Alexandria Engineering Journal, 59(6):4719–4736, 2020.
  * [Het89] Herbert W. Hethcote.  Three basic epidemiological models.  In Applied mathematical ecology (Trieste, 1986), volume 18 of Biomathematics, pages 119–144. Springer, Berlin, 1989.
  * [HOEK18] M. Hassouna, A. Ouhadan, and E. H. El Kinani.  On the solution of fractional order SIS epidemic model.  Chaos Solitons Fractals, 117:168–174, 2018.
  * [KM27] William Ogilvy Kermack and Anderson G McKendrick.  A contribution to the mathematical theory of epidemics.  P. R. Soc. Lond. A, 115(772):700–721, 1927.
  * [KSK+21] Sajjad Ali Khan, Kamal Shah, Poom Kumam, Aly Seadawy, Gul Zaman, and Zahir Shah.  Study of mathematical model of Hepatitis B under Caputo-Fabrizo derivative.  AIMS MATHEMATICS, 6(1):195–209, 2021.
  * [LN15] Jorge Losada and Juan J Nieto.  Properties of a new fractional derivative without singular kernel.  Progr. Fract. Differ. Appl, 1(2):87–92, 2015.
  * [LWLT19] T. Li, Y. Wang, F. Liu, and I. Turner.  Novel parameter estimation techniques for a multi-term fractional dynamical epidemic model of dengue fever.  Numer. Algorithms, 82(4):1467–1495, 2019.
  * [MSK19] Elvin J. Moore, Sekson Sirisubtawee, and Sanoe Koonprasert.  A Caputo-Fabrizio fractional differential equation model for HIV/AIDS with treatment compartment.  Adv. Difference Equ., pages Paper No. 200, 20, 2019.
  * [UKF+20] Saif Ullah, Muhammad Altaf Khan, Muhammad Farooq, Zakia Hammouch, and Dumitru Baleanu.  A fractional model for the dynamics of tuberculosis infection using Caputo-Fabrizio derivative.  Discrete Contin. Dyn. Syst. Ser. S, 13(3):975–993, 2020.
  * [Ver38] Pierre-François Verhulst.  Notice sur la loi que la population suit dans son accroissement.  Corresp. Math. Phys., 10:113–126, 1838.
  * [Ver45] Pierre-Françios Verhulst.  Recherches mathématiques sur la loi d’accroissement de la population.  Journal des économistes, 12:276, 1845.
  * [Ver47] Pierre-François Verhulst.  Deuxième mémoire sur la loi d’accroissement de la population.  Mémoires de l’académie royale des sciences, des lettres et des beaux-arts de Belgique, 20:1–32, 1847.