Effects of Fractional Derivatives with Different Orders in SIS Epidemic Models

Abstract

We study epidemic Susceptible–Infected–Susceptible (SIS) 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 the Caputo derivative, for which we establish existence results of the solutions. Furthermore, we investigate a model based on the Caputo–Fabrizio operator, for which we provide existence of solutions and a study of the equilibria. Both models can be framed in the context of SIS models with time-varying total population, in which the competition between birth and death rates is macroscopically described by the fractional orders of the derivatives. Numerical simulations for both models and a direct numerical comparison are also provided.

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 [1]. Their idea was to divide the population into groups that distinguish 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 [2], describing infections that do not confer immunity to recovery from illnesses such as influenza and common cold. This theory describes compartmental models, where the population is divided into groups depending on the state of individuals, that is, with respect to disease, distinguishing between infected and susceptible. The use of mathematical models for epidemiology is useful to predict the behavior of an infection and make 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 [3] for a recent review of fractional epidemic models.

In this paper, we consider two fractional SIS models. One model is based on the Caputo derivative, for which we establish the existence of results of the solutions and provide numerical simulations. The novelty is to allow the susceptible and infected population to evolve according to different fractional orders: the mathematical setting, in which this type of models has to be considered, needs to be deeply investigated in order to select problems that enjoy existence, uniqueness and closed-form solutions. We refer to [4] for a stability analysis of a Lotka–Volterra predator-prey fractional model with Caputo derivative and different orders of derivative, and to [5] for the case with the same order of derivative. in the framework of epidemic models, we refer to the paper [6], in which an HIV model is analyzed with different orders of fractional derivatives in the Caputo sense.

The other model is based on the Caputo–Fabrizio operator. Here we let the susceptible population evolve according to the Caputo–Fabrizio fractional operator, whereas the infected population dynamics are based on an ordinary differential equation. in this case, we rewrite the system as a system of ordinary differential equations, study the equilibria and present some simulations. We remark that the SIS-type model contains many difficulties that need to be analyzed before generalization to more complicated models: for instance, the closed solutions can be provided in only a few cases. One of them is the class of problems with Caputo–Fabrizio presented in Equation (2) below. Thanks to the explicit form of the solution, we are able to study, without any approximation, related properties of solutions, just looking at their behavior from a mathematical point of view.

More precisely, let $\alpha ,{\alpha }_1,{\alpha }_2\in (0,1]$. We consider the initial-value problem for Caputo derivatives with different orders

(1)$\left\{\begin{array}{c}{{}^{C}D}_t^{{\alpha }_1}S\left(t\right)=-\beta \frac{S\left(t\right)}{S\left(t\right)+I\left(t\right)}I\left(t\right)+\gamma I\left(t\right)\hfill \\ {{}^{C}D}_t^{{\alpha }_2}I\left(t\right)=\beta \frac{S\left(t\right)}{S\left(t\right)+I\left(t\right)}I\left(t\right)-\gamma I\left(t\right)\hfill \\ S\left(0\right)=S_0\hfill \\ I\left(0\right)=I_0\hfill \end{array}\right.$

and the initial-value problem for Caputo–Fabrizio operator

(2)$\left\{\begin{array}{c}{{}^{CF}D}_t^{\alpha }S\left(t\right)=-\left({\displaystyle \frac{\beta }{S\left(t\right)+I\left(t\right)}}S\left(t\right)-\gamma \right)I\left(t\right)\hfill \\ I^{\prime }\left(t\right)=\left({\displaystyle \frac{\beta }{S\left(t\right)+I\left(t\right)}}S\left(t\right)-\gamma \right)I\left(t\right)\hfill \\ S\left(0\right)=S_0\hfill \\ I\left(0\right)=I_0\hfill \end{array}\right.$

where $\beta ,\gamma >0$, $S_0,I_0\ge 0$. Note that results similar to those presented in this paper can be derived also in the case in which the derivatives of S and I in (2) are reversed. The analysis of general, different fractional orders appears to be much more involved than the particular case being examined, and it will be the object of future investigations.

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

$\begin{array}{c}\hfill {{}^{C}D}_t^{{\alpha }_i}u\left(t\right)=\frac{1}{\Gamma (1-{\alpha }_i)}{\int }_0^t\frac{du}{ds}\left(s\right){(t-s)}^{-{\alpha }_i}ds,i=1,2\end{array}$

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

$\begin{array}{c}\hfill \Gamma \left(z\right)={\int }_0^{\infty }{e}^{-s}{s}^{z-1}ds,z>0.\end{array}$

The Caputo derivative is well-defined for a function

$\begin{array}{c}\hfill u\in C(0,T),u^{\prime }\left(s\right){(t-s)}^{-{\alpha }_i}\in L^1\left((0,t)\right),\forall t\in (0,T),i=1,2.\end{array}$

The Caputo–Fabrizio operator in (2) is defined for $u\in H^1\left((a,b)\right)$, $a<b$, by 

$\begin{array}{c}\hfill {{}^{CF}D}_t^{\alpha }u\left(t\right)=\frac{M\left(\alpha \right)}{1-\alpha }{\int }_a^t{u}^{\prime }\left(\tau \right)e^{-\frac{\alpha }{1-\alpha }(t-\tau )}d\tau \end{array}$

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

The case ${\alpha }_1={\alpha }_2=1$ in (1) has been extensively studied in the literature. Beginning with the models introduced in [8,9,10] in which the logistic equation is used to model population dynamics, numerous researchers are working in the field. The case ${\alpha }_1={\alpha }_2\in (0,1)$ was earlier discussed, for instance, in [11,12,13] and the references therein. in this case, thanks to the property of the Caputo derivative

(3)$\begin{array}{c}\hfill {{}^{C}D}_t^{\alpha }u\left(t\right)=0\mathrm{if} \mathrm{and} \mathrm{only} \mathrm{if}u\left(t\right)=const\forall t\ge 0,\alpha \in (0,1)\end{array}$

the total population $N\left(t\right)=S\left(t\right)+I\left(t\right)$ is constant. 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 [14] 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 property (3) of the Caputo derivative. The present approach, based on choosing ${\alpha }_1$ possibly different from ${\alpha }_2$, allows us to study a time-varying total population just adapting the fractional order of derivation. in particular, the offset between ${\alpha }_1$ and ${\alpha }_2$ may provide a macroscopic description of the competition between the birth and the death rates, which characterizes the classical, ordinary SIS models with time-varying population size; see, for instance, [15], Chapter 2.2. When ${\alpha }_1\ne {\alpha }_2$, we deduce from (1) that for ${\alpha }_1,{\alpha }_2\in (0,1)$,

${{}^{C}D}_t^{{\alpha }_1}S\left(t\right)+{{}^{C}D}_t^{{\alpha }_1}I\left(t\right)=0$

and, by the definition of Caputo derivative, that

$\frac{d}{dt}\left(\frac{1}{\Gamma (1-{\alpha }_1)}\left(S\ast {\kappa }_1-\frac{S\left(0\right)}{1-{\alpha }_1}{t}^{1-{\alpha }_1}\right)+\frac{1}{\Gamma (1-{\alpha }_2)}\left(I\ast {\kappa }_2-\frac{I\left(0\right)}{(1-{\alpha }_2)}{t}^{1-{\alpha }_2}\right)\right)=0$

where ${\kappa }_i\left(t\right):=t^{-{\alpha }_i}$, $i=1,2$. Although the above identity has a theoretical interest, from a technical point view, its investigation appears to be hard. When considering a non-singular exponential kernel as the Caputo–Fabrizio operator in the system (2), we may look as well for constant quantities associated to the system. As a result, we find the relation

$\frac{d}{dt}\left(M\left(\alpha \right)S\left(t\right)+\alpha I\left(t\right)+(1-\alpha )\left(\frac{\beta }{S\left(t\right)+I\left(t\right)}S\left(t\right)-\gamma \right)I\left(t\right)\right)=0.$

This identity allows a qualitative and asymptotic analysis of the system (2).

Indeed, our effort is to study the fractional indices of derivation in which ${\alpha }_1$ may differ from ${\alpha }_2$. We mainly study the existence of a solution, and we propose numerical tests based on the method presented in [16], which 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) for future works. In particular, the numerical tests show that, for $\beta >\gamma$, if ${\alpha }_1<{\alpha }_2$, then the total population decreases for large times, and this behavior may take into account the death of part of the infected population. Conversely, if ${\alpha }_1>{\alpha }_2$, then the total population increases, and this can be explained by the birth of susceptible individuals. The results are symmetrically opposite for $\beta <\gamma$.

Recently, the Caputo–Fabrizio operator [17] has attracted many researchers; we refer to [18] (and the references therein) for a complete study of the inclusion of memory in epidemic models by using the Caputo–Fabrizio operator (of the same order). We refer to [19] for a SIR-type model using the Caputo–Fabrizio fractional operator with assumptions about births and deaths: the paper also contains numerical simulations for a selected set of initial conditions. The peculiarity of this fractional operator is the presence of a non-singular kernel in contrast with the Caputo derivative, in which a singular kernel appears in the definition.

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 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 ([20], Chapter 2). In Section 2.2, we complete the study with some numerical simulations. In Section 3, we address a SIS model with the fractional Caputo–Fabrizio operator. In Section 3.1, we establish existence of solutions 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), in particular in terms of well-posedness of the equations. We establish existence results for the solutions of (1), and we present some related numerical simulations. To this end, we introduce the notation

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

so that the above system (1) can be rewritten as

$\left\{\begin{array}{c}{{}^{C}D}_t^{{\alpha }_1}S\left(t\right)=f(S\left(t\right),I\left(t\right))\hfill \\ {{}^{C}D}_t^{{\alpha }_2}I\left(t\right)=-f(S\left(t\right),I\left(t\right))\hfill \\ S\left(0\right)=S_0\hfill \\ I\left(0\right)=I_0.\hfill \end{array}\right.$

2.1. Solutions

The existence of a solution to (1) in the time interval $[0,T)$ (with $T\le +\infty$) with (positive) initial data $(S_0,I_0)$ is implied by the existence of a couple of absolutely continuous functions $\left(S\right(t),I(t\left)\right)$ satisfying the Volterra-type fractional integral equation

(4)$\begin{array}{c}\hfill S\left(t\right)=S_0+\frac{1}{\Gamma \left({\alpha }_1\right)}{\int }_0^{t}f(S\left(s\right),I\left(s\right)){(t-s)}^{{\alpha }_1-1}ds\mathrm{and}\\ \hfill I\left(t\right)=I_0-\frac{1}{\Gamma \left({\alpha }_2\right)}{\int }_0^{t}f(S\left(s\right),I\left(s\right)){(t-s)}^{{\alpha }_2-1}ds\end{array}$

for all $t\in (0,T)$. Fix ${\alpha }_1,{\alpha }_2\in [0,1]$. in the following, we make use of the 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)\ge \frac{t^{{\alpha }_i}}{\Gamma ({\alpha }_i+1)}={\int }_0^t\frac{{(t-s)}^{{\alpha }_i-1}}{\Gamma \left({\alpha }_i\right)}ds,\forall t>0,i=1,2.$

In the next two results we establish some invariance properties for the field f in the cases $\gamma \ge \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\left(t\right),y\left(t\right)$ 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 (1), $(x_0,y_0)$ plays the role of initial data $(S_0,I_0)$, whereas $x\left(t\right)$ and $y\left(t\right)$ represent the (approximated) evolution of the susceptible population $S\left(t\right)$ and infected population $I\left(t\right)$, respectively.

Lemma 1. 

Let ${\alpha }_1,{\alpha }_2\in (0,1]$ and let $0<\beta \le \gamma$ and fix $\epsilon \in (0,1)$. Let $T=T_{\epsilon ,\gamma }$ be the positive real number such that

$G(T;{\alpha }_1,{\alpha }_2)(1+\epsilon )\gamma \le \epsilon .$

Fix $x_0,y_0>0$ and define

$Q_{x_0,y_0}:=[x_0,x_0+\epsilon y_0]\times [(1-\epsilon )y_0,(1+\epsilon )y_0].$

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

$\tilde{x}\left(t\right):=x_0+{\int }_0^{t}f(x\left(s\right),y\left(s\right))\frac{{(t-s)}^{{\alpha }_1-1}}{\Gamma \left({\alpha }_1\right)}ds$

$\tilde{y}\left(t\right):=y_0-{\int }_0^{t}f(x\left(s\right),y\left(s\right))\frac{{(t-s)}^{{\alpha }_2-1}}{\Gamma \left({\alpha }_2\right)}ds$

also satisfy $(\tilde{x}\left(t\right),\tilde{y}\left(t\right))\in Q_{x_0,y_0}$ for all $t\in [0,T]$.

Proof. 

Let $(x\left(t\right),y\left(t\right))\in Q_{x_0,y_0}$ for all $t\in [0,T]$. in particular, one has $x\left(t\right)>0$ and $y\left(t\right)\in ((1-\epsilon )y_0,(1+\epsilon )y_0)$ for all $t\in [0,T]$. Then

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

and this implies

$\begin{array}{cc}\hfill \tilde{x}\left(t\right)& =x_0+{\int }_0^{t}f(x\left(s\right),y\left(s\right))\frac{{(t-s)}^{{\alpha }_1-1}}{\Gamma \left({\alpha }_1\right)}ds\le x_0+{\int }_0^t\gamma y\left(s\right)\frac{{(t-s)}^{{\alpha }_1-1}}{\Gamma \left({\alpha }_1\right)}ds\hfill \\ \hfill & \le x_0+\gamma (1+\epsilon )y_0{\int }_0^t\frac{{(t-s)}^{{\alpha }_1-1}}{\Gamma \left({\alpha }_1\right)}ds\le x_0+G(t;{\alpha }_1,{\alpha }_2)(1+\epsilon )\gamma y_0\le x_0+\epsilon y_0.\hfill \end{array}$

Since $\gamma \ge \beta$ then $f(x,y)>0$ for all $x,y>0$; therefore $\tilde{x}\left(t\right)>x_0>0$ for all $t\in [0,T]$ and this proves that $\tilde{x}\left(t\right)\in (0,x_0+\epsilon y_0)$ for all $t\in [0,T]$. On the other hand,

$\begin{array}{cc}\hfill |\tilde{y}\left(t\right)-y_0|& \le {\int }_0^t\left|f(x\left(s\right),y\left(s\right))\right|\frac{{(t-s)}^{{\alpha }_2-1}}{\Gamma \left({\alpha }_2\right)}ds<{\int }_0^t\gamma y\left(s\right)\frac{{(t-s)}^{{\alpha }_2-1}}{\Gamma \left({\alpha }_2\right)}ds\hfill \\ \hfill & \le \gamma (1+\epsilon )y_0{\int }_0^t\frac{{(t-s)}^{{\alpha }_2-1}}{\Gamma \left({\alpha }_2\right)}ds\le G(t;{\alpha }_1,{\alpha }_2)(1+\epsilon )\gamma y_0\le \epsilon y_0\hfill \end{array}$

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 \le \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 $\epsilon =\epsilon (x_0,y_0):=min\{\frac{1}{2},\frac{1}{2}\frac{x_0}{y_0}\}$ so that

$\frac{1}{2}{x}_0\ge \epsilon y_0.$

Let $T=T_{\beta ,\gamma ,x_0,y_0}$ be such that

(5) $G(T;{\alpha }_1,{\alpha }_2)(\beta +\gamma )(1+\epsilon )\le \epsilon$

and define

$Q_{x_0,y_0}:=[x_0-\epsilon y_0,x_0+\epsilon ]\times [(1-\epsilon )y_0,(1+\epsilon )y_0].$

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

$\tilde{x}\left(t\right):=x_0+{\int }_0^{t}f(x\left(s\right),y\left(s\right))\frac{{(t-s)}^{{\alpha }_1-1}}{\Gamma \left({\alpha }_1\right)}ds$

$\tilde{y}\left(t\right):=y_0-{\int }_0^{t}f(x\left(s\right),y\left(s\right))\frac{{(t-s)}^{{\alpha }_2-1}}{\Gamma \left({\alpha }_2\right)}ds$

also satisfy $(\tilde{x}\left(t\right),\tilde{y}\left(t\right))\in Q_{x_0,y_0}$ for all $t\in [0,T]$.

Proof. 

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

$\left|f\right(x\left(s\right),y\left(s\right)|\le (\beta +\gamma \left)\right|y\left(s\right)|\forall s\in [0,T],$

then one has

$\begin{array}{cc}\hfill |\tilde{x}\left(t\right)-x_0|& =|{\int }_0^{t}f(x\left(s\right),y\left(s\right))\frac{{(t-s)}^{{\alpha }_1-1}}{\Gamma \left({\alpha }_1\right)}ds|\le {\int }_0^t(\gamma +\beta )\left|y\left(s\right)\right|\frac{{(t-s)}^{{\alpha }_1-1}}{\Gamma \left({\alpha }_1\right)}ds|\hfill \\ \hfill & \le G(T;{\alpha }_1,{\alpha }_2)(\beta +\gamma )(1+\epsilon )y_0\le \epsilon y_0.\hfill \end{array}$

This implies $\tilde{x}\left(t\right)\in (x_0-\epsilon y_0,x_0+\epsilon y_0)$ for all $t\in [0,T]$. Arguing as above, one also deduces

$\begin{array}{cc}\hfill |\tilde{y}\left(t\right)-y_0|& =|{\int }_0^{t}f(x\left(s\right),y\left(s\right))\frac{{(t-s)}^{{\alpha }_2-1}}{\Gamma \left({\alpha }_2\right)}ds|\le \epsilon y_0\hfill \end{array}$

and this implies $\tilde{y}\left(t\right)\in [(1-\epsilon )y_0,(1+\epsilon )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 (1) in the case $\gamma \ge \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 \le \gamma$. Then for every initial data $S_0,I_0>0$ the system (1) admits a positive solution in $(0,+\infty )$.

Otherwise, if $0<\gamma <\beta$ then for every initial data $S_0,I_0>0$ the system (1) admits a positive solution in $(0,T]$ where $T=T_{x_0,y_0,\beta ,\gamma }$ is any positive number satisfying

(6) $G(T;{\alpha }_1,{\alpha }_2)\le \left\{\begin{array}{cc}\frac{1}{3(\beta +\gamma )}\hfill & \mathrm{if}{S}_0\ge I_0;\hfill \\ \frac{S_0}{S_0+2I_0}\frac{1}{(\beta +\gamma )}\hfill & \mathrm{if}{S}_0<I_0.\hfill \end{array}\right.$

Proof. 

In order to have a more light notation, let $x_0:=S_0>0$ and $y_0:=I_0>0$. Furthermore, fix $\epsilon :=1/2$. Define

$T:=\left\{\begin{array}{cc}{T}_{\epsilon ,\gamma }\hfill & \mathrm{as}\mathrm{in}\mathrm{Lemma}1\mathrm{if}\gamma \ge \beta \hfill \\ T_{x_0,y_0,\beta ,\gamma }\hfill & \mathrm{as}\mathrm{in}\left(6\right)\mathrm{if}\gamma <\beta .\hfill \end{array}\right.$

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 \ge \beta$ and $Q_{x_0,y_0}$ as in Lemma 2 if $\gamma <\beta$. We first prove that (1) admits a positive solution in $[0,T]$ – note that T is independent from $x_0,y_0$ when $\gamma \ge \beta$. To this end, consider the sequence of functions for $n\ge 1$

(7)${\tilde{x}}_n\left(t\right):=\left\{\begin{array}{cc}{x}_0\hfill & \mathrm{if}t\in [0,T/n]\hfill \\ x_0+{\int }_0^{t-T/n}f({\tilde{x}}_n\left(s\right),{\tilde{y}}_n\left(s\right))\frac{{(t-s)}^{{\alpha }_1-1}}{\Gamma \left({\alpha }_1\right)}ds\hfill & \mathrm{if}t\in (T/n,T];\hfill \end{array}\right.$

(8)${\tilde{y}}_n\left(t\right):=\left\{\begin{array}{cc}{y}_0\hfill & \mathrm{if}t\in [0,T/n]\hfill \\ y_0-{\int }_0^{t-T/n}f({\tilde{x}}_n\left(s\right),{\tilde{y}}_n\left(s\right))\frac{{(t-s)}^{{\alpha }_2-1}}{\Gamma \left({\alpha }_2\right)}ds\hfill & \mathrm{if}t\in (T/n,T].\hfill \end{array}\right.$

Set for brevity ${\tilde{z}}_n\left(t\right):=({\tilde{x}}_n\left(t\right),{\tilde{y}}_n\left(t\right))$.

Claim 1. The sequence of functions ${\left\{{\tilde{z}}_n\left(t\right)\right\}}_{n\ge 1}$ is well defined, equicontinuous and equibounded in $[0,T]$. in particular ${\tilde{z}}_n\left(t\right)\in Q_{x_0,y_0}$ for all $t\in [0,T]$ and for all $n\ge 1$.

We prove the claim by showing by induction that for all n and for all $k=1,\cdots ,n$, ${\tilde{z}}_n\left(t\right)$ is well defined in $t\in [0,kT/n]$, ${\tilde{z}}_n\left(t\right)\in Q_{x_0,y_0}$ for all $t\in [0,kT/n]$, and that

(9)$|{\tilde{z}}_n\left(t_1\right)-{\tilde{z}}_n\left(t_2\right){|\le 2\Vert f\Vert }_{L^{\infty }\left(Q_{x_0,y_0}\right)}\omega \left(\right|t_1-t_2\left|\right)\forall t_1,t_2\in [0,kT/n]$

where $\omega$ is the modulus of continuity of the function $G(·;{\alpha }_1,{\alpha }_2)$. If $k=1$, then by the definitions in (7) and (8) it follows that $z_n\left(t\right)=(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\left(t\right)$ is well defined on $[0,kT/n]$, ${\tilde{z}}_n\left(t\right)\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\left(t\right)$ continuously extends to $[0,(k+1)T/n]$ by setting for $t\in (kT/n,(k+1)T/n]$

${\tilde{z}}_n\left(t\right):=\left(x_0+{\int }_0^{t-T/n}f\left({\tilde{z}}_n\left(s\right)\right)\frac{{(t-s)}^{{\alpha }_1-1}}{\Gamma \left({\alpha }_1\right)}ds,y_0-{\int }_0^{t-T/n}f\left({\tilde{z}}_n\left(s\right)\right)\frac{{(t-s)}^{{\alpha }_2-1}}{\Gamma \left({\alpha }_2\right)}ds\right).$

By applying Lemma 1 or Lemma 2 (according to the cases $\gamma \ge \beta$ and $\gamma <\beta$, respectively) to $(x_n\left(s\right),y_n\left(s\right))=z_n\left(s\right):={\tilde{z}}_n(s+(k-1)T/n)$ (which satisfies $(x_n\left(s\right),y_n\left(s\right))\in Q_{x_0,y_0}$ for all $s\in [0,T/n]$ by inductive hypothesis) we have that ${\tilde{z}}_n\left(t\right)\in Q_{x_0,y_0}$ for all $t\in [kT,(k+1)T/n]$. It is left to show that

(10)$|{\tilde{z}}_n\left(t_1\right)-{\tilde{z}}_n\left(t_2\right){|\le 2\Vert f\Vert }_{L^{\infty }\left(Q_{x_0,y_0}\right)}\omega \left(\right|t_1-t_2\left|\right)\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\left(s\right)\in Q_{x_0,y_0}$ for all $s\in [0,(k+1)T/n]$, one has

$\begin{array}{cc}\hfill |{\tilde{z}}_n\left(t_2\right)-{\tilde{z}}_n\left(t_1\right)|& \le {\int }_{t_1-T/n}^{t_2-T/n}|f\left({\tilde{z}}_n\left(s\right)\right)|\frac{{(t-s)}^{{\alpha }_1-1}}{\Gamma \left({\alpha }_1\right)}ds+{\int }_{t_1-T/n}^{t_2-T/n}\left|\tilde{f}\left(z_n\left(s\right)\right)\right|\frac{{(t-s)}^{{\alpha }_2-1}}{\Gamma \left({\alpha }_2\right)}ds\hfill \\ \hfill & \le {\Vert f\Vert }_{L^{\infty }\left(Q_{x_0,y_0}\right)}{\int }_{t_1-T/n}^{t_2-T/n}\frac{{(t-s)}^{{\alpha }_1-1}}{\Gamma \left({\alpha }_1\right)}+\frac{{(t-s)}^{{\alpha }_2-1}}{\Gamma \left({\alpha }_2\right)}ds\hfill \\ \hfill & \le {2\Vert f\Vert }_{L^{\infty }\left(Q_{x_0,y_0}\right)}|G(t_2-T/n;{\alpha }_1,{\alpha }_2)-G(t_1-T/n;{\alpha }_1,{\alpha }_2)|\hfill \\ \hfill & \le {2\Vert f\Vert }_{L^{\infty }\left(Q_{x_0,y_0}\right)}\omega \left(\right|t_2-t_1\left|\right).\hfill \end{array}$

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\ge t_2-T/n=|t_2-T/n|$. Then arguing as above one gets

$\begin{array}{cc}\hfill |{\tilde{z}}_n\left(t_2\right)-{\tilde{z}}_n\left(t_1\right)|& =|{\tilde{z}}_n\left(t_2\right)-(x_0,y_0)|\le {\int }_0^{t_2-T/n}\left|f\left({\tilde{z}}_n\left(s\right)\right)\right|\frac{{(t-s)}^{{\alpha }_1-1}}{\Gamma \left({\alpha }_1\right)}ds\hfill \\ \hfill & \le {\int }_0^{t_2-T/n}\left|f\left({\tilde{z}}_n\left(s\right)\right)\right|\frac{{(t-s)}^{{\alpha }_1-1}}{\Gamma \left({\alpha }_1\right)}ds\hfill \\ \hfill & \le {2\Vert f\Vert }_{L^{\infty }\left(Q_{x_0,y_0}\right)}\omega \left(\right|t_2-T/n\left|\right)\hfill \\ \hfill & \le {2\Vert f\Vert }_{L^{\infty }\left(Q_{x_0,y_0}\right)}\omega \left(\right|t_2-t_1\left|\right).\hfill \end{array}$

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 $\left\{{\tilde{z}}_{n_k}\left(t\right)\right\}$ converging uniformly in $[0,T]$ to a continuous limit function $z\left(t\right)=\left(S\right(t),I(t\left)\right)$ satisfying $z\left(t\right)\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\left(t\right)>0$ and $I\left(t\right)>0$ for all $t\in [0,T]$. Since f is continuous, then

$f\left({\tilde{z}}_{n_k}\left(t\right)\right)\to f\left(z\left(t\right)\right)=f(S\left(t\right),I\left(t\right))\mathrm{as}k\to +\infty ,\forall t\in [0,T]$

moreover,

$|f\left({\tilde{z}}_{n_k}\left(t\right)\right){|\le \Vert f\Vert }_{L^{\infty }\left(Q_{x_0,y_0}\right)}\forall t\in [0,T].$

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

$\underset{k\to \infty }{lim}{\int }_0^{t}f\left({\tilde{z}}_{n_k}\left(s\right)\right)\frac{{(t-s)}^{{\alpha }_i-1}}{\Gamma \left({\alpha }_i\right)}ds={\int }_0^{t}f\left(z\left(s\right)\right)\frac{{(t-s)}^{{\alpha }_i-1}}{\Gamma \left({\alpha }_i\right)}ds={\int }_0^{t}f(S\left(s\right),I\left(s\right))\frac{{(t-s)}^{{\alpha }_i-1}}{\Gamma \left({\alpha }_i\right)}ds.$

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

$\begin{array}{cc}\hfill \underset{k\to +\infty }{lim}{\tilde{x}}_{n_k}& =\underset{k\to +\infty }{lim}{x}_0+{\int }_0^{t}f\left({\tilde{z}}_{n_k}\left(s\right)\right)\frac{{(t-s)}^{{\alpha }_1}}{\Gamma ({\alpha }_1+1)}ds-{\int }_{t-T/n_k}^{t}f\left({\tilde{z}}_{n_k}\left(s\right)\right)\frac{{(t-s)}^{{\alpha }_1}}{\Gamma ({\alpha }_1+1)}ds\hfill \\ & =x_0+{\int }_0^{t}f(S\left(s\right),I\left(s\right))\frac{{(t-s)}^{{\alpha }_1-1}}{\Gamma \left({\alpha }_1\right)}ds\hfill \end{array}$

and

$\begin{array}{cc}\hfill \underset{k\to +\infty }{lim}{\tilde{y}}_{n_k}& =\underset{k\to +\infty }{lim}{y}_0-{\int }_0^{t}f\left({\tilde{z}}_{n_k}\left(s\right)\right)\frac{{(t-s)}^{{\alpha }_2-1}}{\Gamma \left({\alpha }_2\right)}ds-{\int }_{t-T/n_k}^{t}f\left({\tilde{z}}_{n_k}\left(s\right)\right)\frac{{(t-s)}^{{\alpha }_2-1}}{\Gamma \left({\alpha }_2\right)}ds\hfill \\ & =y_0-{\int }_0^{t}f(S\left(s\right),I\left(s\right))\frac{{(t-s)}^{{\alpha }_2-1}}{\Gamma \left({\alpha }_2\right)}ds.\hfill \end{array}$

On the other hand $({\tilde{x}}_{n_k}\left(t\right),{\tilde{y}}_{n_k}\left(t\right))\to (S\left(t\right),I\left(t\right))$ 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\left(t\right)=S_0+{\int }_0^{t}f(S\left(s\right),I\left(s\right))\frac{{(t-s)}^{{\alpha }_1-1}}{\Gamma \left({\alpha }_1\right)}ds$

and

$I\left(t\right)=I_0-{\int }_0^{t}f(S\left(s\right),I\left(s\right))\frac{{(t-s)}^{{\alpha }_2-1}}{\Gamma \left({\alpha }_2\right)}ds.$

Then $\left(S\right(t),I(t\left)\right)$ is the required solution of (1) in $(0,T]$. If $\beta >\gamma$ then we are done, because we proved the local existence of a solution of (1). If $\beta \le \gamma$, then we can iteratively extend $\left(S\right(t),I(t\left)\right)$ 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 (1) follows [16]. Let us consider a general equation

(11)${{}^{C}D}_t^{\alpha }u\left(t\right)=f\left(u\left(t\right)\right)$

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\left(t^n\right)$, with $t^n=n\Delta t$. Let $\alpha \in (0,1)$, then the Caputo derivative can be approximated as

${{}^{C}D}_t^{\alpha }{u}^n=\frac{1}{\Gamma (2-\alpha )\Delta t^{\alpha }}\left(u^n-\sum _{j=0}^{n-1}{C}_{n,j}{u}^j\right),$

with

$\begin{array}{c}\hfill C_{n,0}=g\left(n\right),C_{n,j}=g(n-j)-g(n-(j-1))\mathrm{for}j=1,\cdots ,n-1\end{array}$

and $g\left(n\right):=n^{1-\alpha }-{(n-1)}^{1-\alpha }\mathrm{for}n\ge 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\left(u^n\right).$

Let us denote by $S^n=S\left(t^n\right)$ and $I^n=I\left(t^n\right)$. By applying this discretization to system (1) we obtain

$\begin{array}{cc}\hfill S^{n+1}& =\sum _{j=0}^{n-1}{C}_{n,j}{S}^j+\Gamma (2-{\alpha }_1)\Delta t^{{\alpha }_1}{f}_S(S^n,I^n)\hfill \\ \hfill I^{n+1}& =\sum _{j=0}^{n-1}{C}_{n,j}{I}^j+\Gamma (2-{\alpha }_2)\Delta t^{{\alpha }_2}(-f_S(S^n,I^n)),\hfill \end{array}$

with

$\begin{array}{c}\hfill f_S(S,I)=-\beta \frac{S}{S+I}I+\gamma I.\end{array}$

We show the evolution in time of $S\left(t\right)$, $I\left(t\right)$ and their sum $S\left(t\right)+I\left(t\right)$ as ${\alpha }_1,{\alpha }_2$ changes: in Figure 1 we considered a case with $\beta >\gamma$ and in Figure 2 with $\gamma <\beta$. Note in plots (a) and (b) that the steepness of the solutions $S\left(t\right)$ and $I\left(t\right)$ in the long run are related to the size of ${\alpha }_1$ and ${\alpha }_2$, respectively: this can be interpreted as a memory effect of the Caputo fractional operator [21]. 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 $S\left(t\right)+I\left(t\right)$ is in general not monotone. For instance, in Figure 1, the total population $S\left(t\right)+I\left(t\right)$ 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 [11,13]. A rigorous study of the symmetries emerging in the simulations and, in particular the intersections of all the functions $N\left(t\right)$ at the same time, is still under investigation.

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 [22,23,24,25] for some examples of epidemic models based on the Caputo–Fabrizio fractional operator.

More precisely, here we study the fractionary SIS system (2). The identity

(12)$\frac{d}{dt}{}^{\mathrm{CF}}{D}_t^{\alpha }u\left(t\right)=\frac{M\left(\alpha \right)}{1-\alpha }{u}^{\prime }\left(t\right)-\frac{\alpha }{1-\alpha }{}^{\mathrm{CF}}{D}_t^{\alpha }u\left(t\right)\forall t\in (a,b),\alpha \in (0,1)$

is crucial to reconduct (2) 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 to the system parameters:

(13)$\begin{array}{c}\hfill \gamma <\frac{\alpha }{1-\alpha },\end{array}$

(14)$\begin{array}{c}\hfill M\left(\alpha \right)\ge \alpha \forall \alpha \in (0,1].\end{array}$

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

3.1. Solutions and Equilibria

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

Theorem 2. 

Let $\alpha \in (0,1]$, $\beta ,\gamma >0$, $S_0,I_0\ge 0$, $I_0+S_0>0$ and assume the conditions (13) and (14). Let

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

and

$\begin{array}{cc}\hfill P_{\alpha }(S_0,I_0)& :=M\left(\alpha \right)S_0+\alpha I_0+(1-\alpha )\left(\beta -\gamma -\frac{\beta }{S_0+I_0}{I}_0\right)I_0.\hfill \end{array}$

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

$g_{\alpha }\left(x\right):=\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\left(\alpha \right)P_{\alpha }(S_0,I_0)x)}}{M\left(\alpha \right)}.$

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

(15) $\begin{array}{cc}\hfill S\left(t\right)& =g_{\alpha }\left(I\left(t\right)\right)\forall t\ge 0\hfill \end{array}$

and $I\left(t\right)$ solves

(16) $\left\{\begin{array}{c}{I}^{\prime }=\left(\beta -\gamma -{\displaystyle \frac{\beta I}{g_{\alpha }\left(I\right)+I}}\right)I\hfill \\ I\left(0\right)=I_0.\hfill \end{array}\right.$

Proof. 

Preliminarily remark that, by (12) and by the first equation of (2)

$\begin{array}{cc}\hfill \frac{d}{dt}{D}_{\alpha }^{\mathrm{CF}}S& =\frac{M\left(\alpha \right)}{1-\alpha }{S}^{\prime }-\frac{\alpha }{1-\alpha }{D}_{\alpha }^{CF}S\hfill \\ \hfill & =\frac{M\left(\alpha \right)}{1-\alpha }{S}^{\prime }+\frac{\alpha }{1-\alpha }\left(\frac{\beta }{S+I}S-\gamma \right)I.\hfill \end{array}$

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

(17)$\left\{\begin{array}{c}{\displaystyle \frac{M\left(\alpha \right)}{1-\alpha }}{S}^{\prime }+{\displaystyle \frac{\alpha }{1-\alpha }}\left({\displaystyle \frac{\beta }{S+I}}S-\gamma \right)I=-{\left(\left({\displaystyle \frac{\beta }{S+I}}S-\gamma \right)I\right)}^{\prime }\hfill \\ I^{\prime }=\left({\displaystyle \frac{\beta }{S+I}}S-\gamma \right)I\hfill \\ S\left(0\right)=S_0\hfill \\ I\left(0\right)=I_0.\hfill \end{array}\right.$

Now define

$\begin{array}{cc}\hfill P\left(t\right):=& M\left(\alpha \right)S\left(t\right)+\alpha I\left(t\right)+(1-\alpha )\left(\frac{\beta }{S\left(t\right)+I\left(t\right)}S\left(t\right)-\gamma \right)I\left(t\right)\hfill \\ \hfill =& M\left(\alpha \right)S\left(t\right)+\alpha I\left(t\right)+(1-\alpha )\left(\beta -\gamma -\frac{\beta }{S\left(t\right)+I\left(t\right)}I\left(t\right)\right)I\left(t\right)\forall t\ge 0\hfill \end{array}$

and remark that (17) implies $P^{\prime }\left(t\right)=0$ and consequently $P\left(t\right)=P\left(0\right)=P_{\alpha }(S_0,I_0)$ for all $t\ge 0$. Then, if  $I\left(t\right)$ and $S\left(t\right)$ are solutions of (17) (or equivalently, of (2)) then

(18)$M\left(\alpha \right)S\left(t\right)+\alpha I\left(t\right)+(1-\alpha )\left(\beta -\gamma -\frac{\beta }{S\left(t\right)+I\left(t\right)}I\left(t\right)\right)I\left(t\right)=P_{\alpha }(S_0,I_0).$

By solving above equation with respect to $S\left(t\right)$ and selecting the only possibly non-negative solution, we obtain for all $t\ge 0$

$\begin{array}{cc}\hfill S\left(t\right)=& \frac{1}{M\left(\alpha \right)}\left(-B_{\alpha }I\left(t\right)+P_{\alpha }(S_0,I_0)/2\right.\hfill \\ \hfill & \left.+\sqrt{{(B_{\alpha }I\left(t\right)-P_{\alpha }(S_0,I_0)/2)}^2-C_{\alpha }{I}^2\left(t\right)+M\left(\alpha \right)P_{\alpha }(S_0,I_0)I\left(t\right)}\right)\hfill \\ \hfill =& g_{\alpha }\left(I\left(t\right)\right).\hfill \end{array}$

Incidentally notice that if $\alpha =1$ then $P_1\left(I_0\right)=I_0+S_0=:N$, $B_1=C_1=1$ we recover from above relation the classical identity $S\left(t\right)=N-I\left(t\right)$. We check that $S\left(t\right)$, as a function $g_{\alpha }$ of $I\left(t\right)$, is well defined. To this end, note that $g_{\alpha }\left(x\right)$ is defined in $\mathbb{R}$, because 

$\begin{array}{cc}& {(B_{\alpha }x-P_{\alpha }(S_0,I_0)/2)}^2-C_{\alpha }{x}^2+M\left(\alpha \right)P_{\alpha }(S_0,I_0)x=\hfill \\ \hfill & (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\ge 0\hfill \end{array}$

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\left(\alpha \right)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 }\left(x\right)=N-x$.

By plugging the identity $S\left(t\right)=g_{\alpha }\left(I\left(t\right)\right)$ in the second equation of (2) we obtain (17).

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

$\begin{array}{cc}\hfill P_{\alpha }(S_0,I_0)& \ge \alpha I_0+(1-\alpha )(\beta -\gamma -\frac{\beta }{S_0+I_0}{I}_0)I_0\hfill \\ \hfill & \ge (\alpha -(1-\alpha )\gamma )I_0.\hfill \end{array}$

Moreover if $I\left(t\right)=P_{\alpha }(S_0,I_0)/(\alpha -(1-\alpha )\gamma )$ then $I^{\prime }\left(t\right)=-\gamma I\left(t\right)<0$ and, consequently $I\left(t\right)\le P_{\alpha }(S_0,I_0)/(\alpha -(1-\alpha )\gamma )$. Finally, remark that 0 is an equilibrium for (17) and that the velocity field $f_{\alpha }\left(I\right)=\left(\beta -\gamma -{\displaystyle \frac{\beta I}{g_{\alpha }\left(I\right)+I}}\right)I$ is locally Lipschitz continuous. Therefore, by the local uniqueness of the solutions of (17) if $I_0>0$ then $I\left(t\right)>0$ for all t.    □

Equilibria

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

Proposition 1. 

Let $\alpha \in (0,1]$ and assume conditions (13) and (14) and let $S_0,I_0\ge 0$. The equilibria of the system (16) are 0 and

$E_{\alpha }(S_0,I_0):=\frac{(\beta -\gamma )P_{\alpha }(S_0,I_0)}{\alpha (\beta -\gamma )+M\left(\alpha \right)\gamma }.$

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

$I\left(t\right)\to E_{\alpha }(S_0,I_0),S\left(t\right)\to \frac{\gamma }{\beta -\gamma }{E}_{\alpha }(S_0,I_0)=\frac{\gamma P_{\alpha }(S_0,I_0)}{\alpha (\beta -\gamma )+M\left(\alpha \right)\gamma }\mathrm{as}t\to +\infty .$

Finally, if $E_{\alpha }(S_0,I_0)\le 0$, that is if $R\le 1$, then

$I\left(t\right)\to 0,S\left(t\right)\to P_{\alpha }(S_0,I_0)\mathrm{as}t\to +\infty .$

Proof. 

The first part of the claim follows by a direct computation, using in particular (14) 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\le 0$. We proved in Theorem 2 that $I\left(t\right)$ is non-negative, and consequently $[0,+\infty )$ is an invariant set for the dynamics (17). Moreover, by a direct computation one can check that the function $V\left(x\right):={(x-x_R)}^2$ is a Lyapunov function for (17) in $[0,+\infty )$ and, consequently, $x_R$ is a globally asymptotically stable equilibrium for (17) 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 seen 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 behavior of the function $N\left(t\right):=S\left(t\right)+I\left(t\right)$, which is not constant unless we are in the ordinary case $\alpha =1$. As we show below, the asymptotic behavior of $N\left(t\right)$ can be used for inverse problems, i.e., to reconstruct the fractional order $\alpha$ of (2) from the observed data in 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\left(t\right):=I\left(t\right)+S\left(t\right)$ 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\left(\alpha \right)\gamma }$

as $t\to +\infty$. Otherwise, if either $R>1$ and $I_0>E_{\alpha }(S_0,I_0)$ or $R\le 1$, then $N\left(t\right)$ is strictly decreasing and tends to $P_{\alpha }(S_0,I_0)/M\left(\alpha \right)$ as $t\to +\infty$.

Proof. 

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

(19)$g_{\alpha }^{″}\left(x\right)=\frac{P_{\alpha }^2(S_0,I_0)(1-\alpha )\beta }{4M\left(\alpha \right)\sqrt{{(B_{\alpha }x-P_{\alpha }(S_0,I_0)/2)}^2-C_{\alpha }{x}^2+M\left(\alpha \right)P_{\alpha }(S_0,I_0)}}\ge 0\forall x\in \mathbb{R}.$

Moreover, by a direct computation,

$1+g_{\alpha }^{\prime }\left(I_0\right)=1+\frac{2(1-\alpha )\beta M\left(\alpha \right){(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

(20)$1+g_{\alpha }^{\prime }\left(x\right)>0\forall x\ge I_0.$

Furthermore, remark that, using the assumption (14), i.e., $M\left(\alpha \right)\ge \alpha$, we obtain

$\begin{array}{cc}\hfill 1+g_{\alpha }^{\prime }\left(E_{\alpha }(S_0,I_0)\right)=& 1+\frac{M\left(\alpha \right)\left((1-\alpha )(\beta -\gamma )\gamma -\alpha \beta \right)}{(1-\alpha ){(\beta -\gamma )}^2+\beta M\left(\alpha \right)}\hfill \\ \hfill \ge & 1+\frac{M\left(\alpha \right)\left((1-\alpha )(\beta -\gamma )\gamma -\beta \right)}{(1-\alpha ){(\beta -\gamma )}^2+\beta M\left(\alpha \right)}\hfill \\ \hfill =& \frac{(1-\alpha )(\beta -\gamma )\left(M\right(\alpha )\gamma +\beta -\gamma )}{(1-\alpha ){(\beta -\gamma )}^2+\beta M\left(\alpha \right)}\hfill \end{array}$

therefore

(21)$\beta >\gamma \Rightarrow 1+g_{\alpha }^{\prime }\left(x\right)>0\forall x\ge E_0(S_0,I_0).$

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

$1+g_{\alpha }^{\prime }\left(0\right)=(1-\alpha )(1-\beta +\gamma )$

therefore

(22)$\beta \le \gamma \Rightarrow 1+g_{\alpha }^{\prime }\left(x\right)>0\forall x\ge 0.$

Now, by Theorem 2

(23)$N^{\prime }\left(t\right)=I^{\prime }\left(t\right)+S^{\prime }\left(t\right)=I^{\prime }\left(t\right)(1+g_{\alpha }^{\prime }\left(I\left(t\right)\right).$

Assume $R>1\mathrm{and}{I}_0\in (0,E(S_0,I_0))$. Then

(24)$I^{\prime }\left(t\right)=\left(\beta -\gamma -\frac{\beta I\left(t\right)}{I\left(t\right)+g_{\alpha }\left(I\left(t\right)\right)}\right)I\left(t\right)>0\forall t\ge 0.$

In particular, $I\left(t\right)\ge I_0$ for all $t\ge 0$, and, in view of (20) and (23), we deduce $N^{\prime }\left(t\right)>0$ for all $t>0$; hence, $N\left(t\right)$ is strictly increasing.

Assume now that either $R\le 1\mathrm{or}{I}_0\ge 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 }\left(t\right)<0$, and, since $E_{\alpha }(S_0,I_0)$ is a (globally asymptotically stable) equilibrium by Proposition 1, then $I\left(t\right)>E_{\alpha }(S_0,I_0)$ for all $t\ge 0$. Then, in view of (21) and (23), we deduce $N^{\prime }\left(t\right)<0$ for all $t\ge 0$.

Finally, if $R\le 1$ then $I^{\prime }\left(t\right)<0$. Since $I\left(t\right)>0$ for all $t>0$ then, in view of (22) and (23), we deduce $N^{\prime }\left(t\right)<0$ for all $t\ge 0$, and this concludes the proof.    □

From the above result, we can estimate the fractional order $\alpha$ from the asymptotic behavior of the total population $N\left(t\right)$. Indeed, if $R>1$ and $N\left(t\right)\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\left(\alpha \right)\gamma }=N_{\infty }.$

By assuming, as in [26], $M\left(\alpha \right)\equiv 1$, and setting $N_0:=I_0+S_0$, the 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\left(t\right)$ is constant if and only if $\alpha =1$.

3.2. Numerical Simulations

The numerical discretization of (2) 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 (16) by a proper ODE solver and then we discretize (15). Specifically, we use the MATLAB tool ode23t to compute $I^{n+1}$ and then, following (15), we obtain $S^{n+1}=g_{\alpha }\left(I^{n+1}\right)$.

In Figure 3 and Figure 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\left(\alpha \right)$ and $\gamma$ are the same in all simulations, fixed as $S_0=6$, $I_0=4$, $M\left(\alpha \right)\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\left(t\right)=S\left(t\right)+I\left(t\right)$ is constant, since (2) coincides with the classical SIS model, while if $\alpha <1$, then the monotonicity of $N\left(t\right)$ is in agreement with Theorem 2.

In Figure 5 and Figure 6, we fix $\alpha =0.5$ and we set $S_0=10-I_0$ and let the initial data $I_0$ vary 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.

Comparison between Caputo SIS Model and Caputo–Fabrizio SIS Model

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 (1) and Caputo–Fabrizio fractional SIS model (2). To this end, we fix the initial data $S_0=6$ and $I_0=4$. We consider the set the fractional orders ${\alpha }_1$ and ${\alpha }_2$ in (1), the fractional order $\alpha$ in (2), 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 with $\gamma <\beta$, in particular $\gamma =0.2$ and $\beta =0.7$, whereas in Figure 8, we test a case with $\beta <\gamma$, in particular $\gamma =0.2$ and $\beta =0.1$. Note that the evolution of susceptible population $S\left(t\right)$ 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 (1) 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 and the setting in a mathematical framework of these generalized models. Not only is our work based on numerical simulations, but we also tried to give some mathematical answers, and we used the numerical simulations to complete the investigation. 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 take some further steps 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 the effects of each particular differential operator on the system emerge. We remark that the numerical simulations of the Caputo-type model showed a non-monotone nature of the solutions for small values of one of the two fractional orders. The sum of infected and susceptible populations is not constant, suggesting that the vital dynamics are hidden in the mixed fractional orders, and this fact deserves further investigations. The tests related to the Caputo–Fabrizio model confirm the theoretical results developed in Section 3. The two types of fractional derivative are finally compared in section, where we observed again the lack of monotonicity of Caputo solutions, in contrast to the monotonicity shown by Caputo–Fabrizio solutions.

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.

Author Contributions

Conceptualization; methodology; software; validation; formal analysis; investigation; writing—original draft preparation; writing—review and editing: C.B., M.D., A.C.L. and P.L. All authors have read and agreed to the published version of the manuscript.

Funding

This research was partially supported by Progetto Ateneo Sapienza Università di Roma "Fractional derivatives in Science and Engineering”, grant number RM118164367E9B16.

Data Availability Statement

Data is contained within this article.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures

Figure 1. Numerical solutions to (1) with β=0.7, γ=0.2, I0=2 and S0=8.

Figure 2. Numerical solutions to (1) with β=0.1, γ=0.2, I0=8 and S0=2.

Figure 3. Numerical solutions to (2) for β=0.7, γ=0.2, S0=6, I0=4, and M(α)≡1 as α changes.

Figure 4. Numerical solutions to (2) with β=0.1, γ=0.2, S0=6, I0=4, and M(α)≡1 as α changes.

Figure 5. Numerical solutions to (2) for α=0.5, β=0.7, γ=0.2 and M(α)≡1 as I0 and S0 change.

Figure 6. Numerical solutions to (2) for α=0.5, β=0.1, γ=0.2 and M(α)≡1 as I0 and S0 change.

Figure 7. Numerical solutions to (1) and (2) for α=α1=0.2 (top), α=α1=0.5 (center), and α=α1=0.8 (bottom) with α2=1, β=0.7, γ=0.2, S0=6, I0=4 and M(α)≡1.

Figure 8. Numerical solutions to (1) and (2) for α1=0.2 (top), α1=0.5 (center) and α1=0.8 (bottom) with α2=1, β=0.1, γ=0.2, S0=6, I0=4 and M(α)≡1.