Fractional SIS Epidemic Models

Abstract

In this paper, we consider the fractional SIS (susceptible-infectious-susceptible) epidemic model (α-SIS model) in the case of constant population size. We provide a representation of the explicit solution to the fractional model and we illustrate the results by numerical schemes. A comparison with the limit case when the fractional order α converges to 1 (the SIS model) is also given. We analyze the effects of the fractional derivatives by comparing the SIS and the α-SIS models.

1. Introduction

The study of mathematical models for epidemiology has a long history, dating back to the early 1900s with the theory developed by Kermack and McKendrick [1]. Such theory describes compartmental models, where the population is divided into groups depending on the state of individuals with respect to disease, distinguishing between groups. The dynamic of the disease is then described by a system of ordinary differential equations for each class of individuals. The use of mathematical models for epidemiology is particularly useful to predict the progress of an infection and to take strategy to limit the spread of the disease. In this work, we focus on the α-SIS (susceptible-infectious-susceptible) epidemiological model. The SIS model has a long history, too [2]. It describes the spread of human viruses, such as influenza. The SIS model with constant population is particularly appropriate to describe some bacterial agent diseases, such as gonorrhea, meningitis, and streptococcal sore throat. SIS is a model without immunity, where the individual recovered from the infection comes back into the class of susceptibles.

1.1. Statement of the Problem

We propose an α-SIS model with constant population size. The novelty concerns the SIS equations with the time fractional Caputo derivative in place of time standard derivative and their explicit solutions in terms of Euler’s numbers and Euler’s Gamma functions.

Let us consider the Caputo fractional derivative introduced in (Section 2.1) below. We provide an explicit representation of the solution to

$$\begin{array}{cc}\hfill & \left\{\begin{array}{c}{\mathcal{D}}_t^{\alpha }S\left(t\right)=\mu -\beta S\left(t\right)I\left(t\right)+\gamma I\left(t\right)-\mu S\left(t\right)\hfill \\ {\mathcal{D}}_t^{\alpha }I\left(t\right)=\beta S\left(t\right)I\left(t\right)-\gamma I\left(t\right)-\mu I\left(t\right)\hfill \end{array}\right.\hfill \\ \hfill & \mathrm{with}S\left(t\right)+I\left(t\right)=N\left(t\right),S\left(0\right)=S_0\mathrm{and}I\left(0\right)=I_0\hfill \end{array}$$

for the constant population case $N\left(t\right)=1$, $\forall t$, where $\alpha \in (0,1)$ is the order of the Caputo fractional derivative, $\mu$ is the birth rate and the death removal rate, $\beta$ is the contact rate, and $\gamma$ is the recovery removal rate. The unknown functions $S\left(t\right)$ and $I\left(t\right)$ represent the percentage of susceptible and infected people at time $t>0$ with initial data $S_0$ and $I_0$. As far as we know, although, in the numerical literature, it is an unknown formula for the solution. By using a series representation for the solution to the fractional logistic equation we may give an explicit formula for the unknown functions S and I. From the numerical point of view, we validate the goodness of the theoretical formulas by applying two different numerical schemes. Then, we compare the fractional case results ($0<\alpha <1$) with the well-known standard case taking the limit $\alpha$ converges to 1, and we analyze the effects produced by the fractional derivatives.

1.2. Motivations

Let us consider an infective disease which does not confer immunity and which is transmitted through contact between people. We divide the population into two disjoint classes which evolve in time: the susceptibles and the infectives. The first class contains the individuals which are not yet infected but who can contract the disease; the second class contains the infected population which can transmit the disease. The SIS model [2] is a simple disease model without immunity, where the individuals recovered from the infection come back into the class of susceptibles. Such a model is used to describe the dynamic of infections which do not confer a long immunity, such as a cold or influenza. Fractional calculus is therefore considered in biological models to take into account macroscopic effect. The use of fractional derivatives in the model means that some global effect may produce slowdown in the process. This is verified and discussed in the validation of the model.

1.3. State of the Art

The logistic function was introduced by Pierre Francois Verhulst [3] to model the population growth. At the beginning of the process, the growth of the population is fast; then, as saturation process begins, the growth slows, and then growth is close to be flat. The problem to give a solution of the fractional logistic equation was unsolved and several attempts have been done (see, for instance, Reference [4,5,6,7]). Concerning the fractional SIS model, some works can be listed about numerical solutions obtained by considering different methods. The existence and uniqueness of the solution have been discussed in Reference [8] in case of constant population size. Moreover, in that work, the authors have studied a numerical solution by variational iteration method and Euler method. In the previous paper of Reference [9], the case of variable population size has been considered and the stability of the equilibrium point has been investigated together with the existence of the solution. From the technical point of view our result take advantage of the explicit representation by series of the solution of fractional logistic equation solved in the recent paper of Reference [10]. Thanks to a fruitful formulation of the SIS model, we are able to adapt the results obtained for fractional logistic equation in Reference [10] and to give the solution of fractional SIS model by series. For the physical derivation of the fractional model, we refer to the interesting paper of Reference [11], in which the authors have considered a probabilistic approach in terms of continuous-time random walk in order to obtain the fractional SIR model.

In recent years, the study of epidemiological models using fractional calculus has spread widely. In Reference [12], the authors prove via numerical simulations that the proposed fractional model gives better results than the classical theory, when compared to real data. Moreover, for some diseases, it is necessary to take into account the history of the system (see, for example, Reference [13]); thus, non-locality and memory become important to model real data. Indeed, fractional operators consider the entire history of the biological process, and we are able to model non-local effects often encountered in biological phenomena.

1.4. Main Results

We provide an explicit representation of the solution to

$$\begin{array}{cc}\hfill & \left\{\begin{array}{c}{\mathcal{D}}_t^{\alpha }S\left(t\right)=\mu -\beta S\left(t\right)I\left(t\right)+\gamma I\left(t\right)-\mu S\left(t\right)\hfill \\ {\mathcal{D}}_t^{\alpha }I\left(t\right)=\beta S\left(t\right)I\left(t\right)-\gamma I\left(t\right)-\mu I\left(t\right)\hfill \end{array}\right.\hfill \\ \hfill & \mathrm{with}S\left(t\right)+I\left(t\right)=1,S\left(0\right)=S_0\mathrm{and}I\left(0\right)=I_0\hfill \end{array}$$

(1)

in terms of uniformly convergent series on compact sets.

Let us introduce the basic reproduction number [14], i.e., the expected number of secondary infections produced during the period of infection, which is given by

$$\sigma =\frac{\beta }{\gamma +\mu },$$

(2)

where $\gamma +\mu$ is the infection period. Let

$${\displaystyle c=\frac{\sigma -1}{\sigma }}$$

(3)

be the so-called carrying capacity and define $b=\beta c$. The problem (1) can be solved by considering, for $c\ne 0$, the fractional logistic equation

$${\mathcal{D}}_t^{\alpha }I\left(t\right)=bI\left(t\right)\left(1-\frac{1}{c}I\left(t\right)\right).$$

(4)

In the following theorems, $B(x,y)$ denotes the Beta function, $\Gamma \left(x\right)$ denotes the Euler Gamma function, and $E_k^{\alpha }$ is the $\alpha$-Euler’s number introduced in Reference [10].

Theorem 1. 

Let $\alpha \in (0,1)$, $c\ne 0$ and $b^{1/\alpha }<1$. An explicit representation of the solution of the fractional $SIS$ model (1) with initial condition $I_0=c/2$ and $S_0=1-I_0$ is given by

$$\begin{array}{cc}\hfill I\left(t\right)& =c\sum _{k\ge 0}{E}_k^{\alpha }{b}^{\alpha k}\frac{t^{\alpha k}}{\Gamma (\alpha k+1)}\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill S\left(t\right)& =1-I\left(t\right),\hfill \end{array}$$

(6)

with

$$\begin{array}{c}\hfill E_0^{\alpha }=\frac{1}{2},E_1^{\alpha }=E_0^{\alpha }-{\left(E_0^{\alpha }\right)}^2\end{array}$$

and $\forall k\ge 1$

$$\begin{array}{c}\hfill {\displaystyle E_{2k}^{\alpha }=0,E_{2k+1}^{\alpha }=-\frac{1}{\alpha k+1}\sum _{\begin{array}{c}i,j\\ i+j=k\end{array}}\frac{E_i^{\alpha }{E}_j^{\alpha }}{B(\alpha i+1,\alpha j+1)}.}\end{array}$$

The series is uniformly convergent on any compact subset $K\subseteq (0,r_{\alpha })$, where

$$r_{\alpha }=\frac{1}{b^{1/\alpha }}{\left(\frac{\Gamma (\alpha +1)\Gamma (3\alpha +1)}{\Gamma (2\alpha +1)}\right)}^{\frac{1}{2\alpha }}.$$

(7)

Theorem 2. 

Let $\alpha \in (0,1)$, $c=0$. An explicit representation of the solution of the fractional SIS model (1) with initial condition $I_0=1/\left(2\beta \right)$ and $S_0=1-I_0$ is given by

$$\begin{array}{cc}\hfill I\left(t\right)& =\frac{1}{\beta }\sum _{k\ge 0}{A}_k^{\alpha }\frac{t^{\alpha k}}{\Gamma (\alpha k+1)},\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill S\left(t\right)& =1-I\left(t\right),\hfill \end{array}$$

(9)

with $A_0^{\alpha }=\frac{1}{2},A_1^{\alpha }=-{\left(A_0^{\alpha }\right)}^2$ and

$$A_{k+1}^{\alpha }=-\frac{1}{\alpha k+1}\sum _{\begin{array}{c}i,j\\ i+j=k\end{array}}\frac{A_i^{\alpha }{A}_j^{\alpha }}{B(\alpha i+1,\alpha j+1)}\forall k\ge 1.$$

The series converges uniformly in $K\subset (0,r_{\alpha })$ with $r_{\alpha }\le {(1/2)}^{1/\alpha }$.

1.5. Outline

The paper is organized as follows. In Section 2, we introduce the fractional α-SIS model with constant population size. In Section 3, we prove the main results of the paper. In Section 4, we validate the model using two numerical schemes, and we provide some numerical tests also comparing the α-SIS model with the SIS one.

2. The Settings

In this section, we briefly review the mathematical background on fractional derivatives and its connection with the SIS model.

2.1. The Fractional Derivatives

Fractional Calculus has a long history, starting from some works by Leibniz (1695) or Abel (1823), to where it has been developed up to today. The literature is vast and many definitions of fractional derivatives has been given. We recall the well-known derivatives of Caputo and Riemann-Liouville given by following the definitions we will deal with throughout. The Caputo Derivative of a function $u\left(t\right)$ is written as

$$\begin{array}{c}\hfill {\mathcal{D}}_t^{\alpha }u\left(t\right):=\frac{1}{\Gamma (1-\alpha )}{\int }_0^t\frac{u^{\prime }\left(s\right)}{{(t-s)}^{\alpha }}ds,t>0,\end{array}$$

(10)

whereas the Riemann-Liouville derivative of $u\left(t\right)$ is defined as follows:

$${\mathcal{D}}_t^{\alpha }u\left(t\right)=\frac{1}{\Gamma (1-\alpha )}\frac{d}{dt}{\int }_0^t\frac{u\left(s\right)}{{(t-s)}^{\alpha }}ds.$$

(11)

Notice that, for $a<b$, if $u\in L^1(a,b)$ such that $u^{\prime }\in L^1(a,b)$ and $|u^{\prime }\left(t\right)|\le t^{\gamma -1}$ a.e. with $\gamma >0$, then we have that for $t\in (a,b)$

$$\begin{array}{c}\hfill \left|{\mathcal{D}}_t^{\alpha }u\left(t\right)\right|\le \frac{1}{\Gamma (1-\alpha )}{\int }_0^t{s}^{\gamma -1}{(t-s)}^{1-\alpha -1}ds=\frac{B(\gamma ,1-\alpha )}{\Gamma (1-\alpha )},\end{array}$$

where

$$\begin{array}{c}\hfill B(\alpha ,\beta )=\frac{\Gamma \left(\alpha \right)\Gamma \left(\beta \right)}{\Gamma (\alpha +\beta )},\alpha >0,\beta >0\end{array}$$

is the Beta function and $\Gamma \left(\alpha \right)={\int }_0^{\infty }{e}^{-s}{s}^{\alpha -1}ds$, $\alpha >0$ is the Euler’s gamma function. The Caputo and the Riemann-Liouville fractional derivatives are linked by the following formula:

$$\begin{array}{c}\hfill {\mathcal{D}}_t^{\alpha }u\left(t\right)={\mathcal{D}}_t^{\alpha }u\left(t\right)-\frac{t^{-\alpha }}{\Gamma (1-\alpha )}u\left(0\right)={\mathcal{D}}_t^{\alpha }\left(u\left(t\right)-u\left(0\right)\right),\end{array}$$

(12)

which will be useful further on. We list some useful properties of the Caputo derivative:

(P1)

Let u be a constant function. Then ${\mathcal{D}}_t^{\alpha }u\left(t\right)=0$.

(P2)

Let $u:[a,b]\to \mathbb{R}$ such that $u\left(a\right)=0$ and ${\mathcal{D}}_t^{\alpha }u$, ${\mathcal{D}}_t^{\alpha }u$ exist almost everywhere. Then, ${\mathcal{D}}_t^{\alpha }u={\mathcal{D}}_t^{\alpha }u$.

(P3)

Let $u,v:[a,b]\to \mathbb{R}$ be such that ${\mathcal{D}}_t^{\alpha }u\left(t\right)$ and ${\mathcal{D}}_t^{\alpha }v\left(t\right)$ exist almost everywhere in $[a,b]$. Let $c,d\in \mathbb{R}$. Then, ${\mathcal{D}}_t^{\alpha }(cu\left(t\right)+dv\left(t\right))$ exists almost everywhere in $[a,b]$. In particular,

$$\begin{array}{c}\hfill {\mathcal{D}}_t^{\alpha }(cu\left(t\right)+dv\left(t\right))=c{\mathcal{D}}_t^{\alpha }u\left(t\right)+d{\mathcal{D}}_t^{\alpha }v\left(t\right).\end{array}$$

(P4)

Let $u\in C^1\left([a,b]\right)$. Then,

$${\mathcal{D}}_t^{\alpha }u\left(t\right)\to u^{\prime }\left(t\right),\mathrm{as}\alpha \to 1^{-}$$

pointwise in $(a,b]$.

(P1) and (P3) are immediate consequences of the definition of the Caputo derivative. (P2) can be obtained from (12). (P4) follows from the definition given for $\alpha \in (0,1)$. Our discussion here is based on the result in [15] (Theorem 2.20) for the Riemann-Liouville derivative and the definition (12) above of the Caputo derivative. The interested reader can also consult [16] (p. 20) in which the connection with the Marchaud derivative is considered.

Let us consider the equation ${\mathcal{D}}_t^{\alpha }u+au=0$ on $K=[0,\infty )$ with $u\left(0\right)=1$ where $a\in \mathbb{R}$. Then, u is the Mittag-Leffler function

$$u\left(t\right)=E_{\alpha }(-a{t}^{\alpha })=\sum _{k\ge 0}{(-a)}^k\frac{t^{\alpha k}}{\Gamma (\alpha k+1)},t\in K.$$

(13)

For the reader’s convenience, we write below the proof of this standard result. From the Laplace transform

$$\begin{array}{c}\hfill {\int }_0^{\infty }{e}^{-\lambda t}{\mathcal{D}}_t^{\alpha }u\left(t\right)dt={\lambda }^{\alpha }\tilde{u}\left(\lambda \right)-{\lambda }^{\alpha -1}u\left(0\right),\end{array}$$

where $\tilde{u}\left(\lambda \right)={\int }_0^{\infty }{e}^{-\lambda t}u\left(t\right)dt$, the equation takes the form ${\lambda }^{\alpha }\tilde{u}\left(\lambda \right)-{\lambda }^{\alpha -1}u\left(0\right)=a\tilde{u}\left(\lambda \right)$, that is

$$\begin{array}{c}\hfill \tilde{u}\left(\lambda \right)=u\left(0\right)\frac{{\lambda }^{\alpha -1}}{a+{\lambda }^{\alpha }}={\int }_0^{\infty }{e}^{-\lambda t}{E}_{\alpha }(-a{t}^{\alpha })dt,\lambda >0,\end{array}$$

since $u\left(0\right)=1$. From the Stirling’s formula for Gamma function, we have

$$\begin{array}{c}\hfill {\left(\frac{a^k}{\Gamma (\alpha k+1)}\right)}^{1/k}\sim a{\left(\frac{e}{\alpha k+1}\right)}^{\frac{\alpha +1}{k}}{\left(2\pi (\alpha k+1)\right)}^{-1/\left(2k\right)}(1+o\left(1\right)).\end{array}$$

Thus, we get that

$$\begin{array}{c}\hfill {\left(\frac{a^k}{\Gamma (\alpha k+1)}\right)}^{1/k}\to 0\mathrm{as}k\to \infty .\end{array}$$

Thus, by the root criterion, we get an infinite radius of convergence.

2.2. The Fractional SIS Model

In the discussion above the symbols $S\left(t\right)$ and $I\left(t\right)$ have been used denoting percentages. Indeed, $N\left(t\right)=1$ is a constant function for any t. Denoting by $\mathcal{S}\left(t\right)$ and $\mathcal{I}\left(t\right)$ the number of susceptibles and infectives, respectively, at time t, the fractional SIS model with non-constant population (see Reference [17,18] for $\alpha =1$, that is, the non-fractional case, but we say SIS model) is written as

$$\begin{array}{cc}\hfill & \left\{\begin{array}{c}{\displaystyle {\mathcal{D}}_t^{\alpha }\mathcal{S}\left(t\right)=\mathsf{\Lambda }\mathcal{N}\left(t\right)-\beta \frac{\mathcal{S}\left(t\right)\mathcal{I}\left(t\right)}{\mathcal{N}\left(t\right)}+\gamma \mathcal{I}\left(t\right)-\mu \mathcal{S}\left(t\right)}\hfill \\ {\displaystyle {\mathcal{D}}_t^{\alpha }\mathcal{I}\left(t\right)=\beta \frac{\mathcal{S}\left(t\right)\mathcal{I}\left(t\right)}{\mathcal{N}\left(t\right)}-\gamma \mathcal{I}\left(t\right)-\mu \mathcal{I}\left(t\right)}\hfill \end{array}\right.\hfill \\ \hfill & \mathrm{with}\mathcal{N}\left(t\right)=\mathcal{S}\left(t\right)+\mathcal{I}\left(t\right),S\left(0\right)={\mathcal{S}}_0\mathrm{and}I\left(0\right)={\mathcal{I}}_0,\hfill \end{array}$$

(14)

where $\mathsf{\Lambda }$ is the birth rate, $\mu$ is the death removal rate, $\beta$ is the contact rate, and $\gamma$ is the recovery removal rate. The sum of susceptibles and infectives is defined by $\mathcal{N}\left(t\right)$.

The problem to solve (14) is challenging for many reasons. To overcome such difficulties we introduce the difference between the susceptible and infective populations given by

$$\mathcal{Z}\left(t\right)=\mathcal{S}\left(t\right)-\mathcal{I}\left(t\right),$$

(15)

from which we are able to recover the functions $\mathcal{S}$ and $\mathcal{I}$ as follows:

$$\mathcal{S}\left(t\right)=\frac{\mathcal{N}\left(t\right)+\mathcal{Z}\left(t\right)}{2}\mathrm{and}\mathcal{I}\left(t\right)=\frac{\mathcal{N}\left(t\right)-\mathcal{Z}\left(t\right)}{2}.$$

By the linearity of the Caputo derivative, (see (P3)) the problem takes the form

$$\begin{array}{cc}\hfill {\mathcal{D}}_t^{\alpha }\mathcal{N}\left(t\right)& =(\mathsf{\Lambda }-\mu )\mathcal{N}\left(t\right)\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill {\mathcal{D}}_t^{\alpha }\mathcal{Z}\left(t\right)& =\left(\mathsf{\Lambda }-\frac{\beta }{2}+\gamma \right)\mathcal{N}\left(t\right)-(\gamma +\mu )\mathcal{Z}\left(t\right)\left(1-\frac{\beta }{2\mathcal{N}\left(t\right)(\gamma +\mu )}\mathcal{Z}\left(t\right)\right).\hfill \end{array}$$

(17)

In this new formulation, we are able to solve (16) by using standard results.

Proposition 1. 

The solution to (16) with initial datum ${\mathcal{N}}_0={\mathcal{S}}_0+{\mathcal{I}}_0$ is

$$\mathcal{N}\left(t\right)={\mathcal{N}}_0{E}_{\alpha }\left((\mathsf{\Lambda }-\mu )t^{\alpha }\right),$$

(18)

where $E_{\alpha }$ is the Mittag-Leffler function, defined in (13).

Notice that $\mathcal{N}\left(t\right)\ge 0$ is an increasing function as $\mathsf{\Lambda }-\mu >0$, whereas it exhibits a decreasing behavior for $\mathsf{\Lambda }-\mu <0$. Thus, we can write the non-obvious relation

$$\begin{array}{c}\hfill {\mathcal{D}}_t^{\alpha }\mathcal{N}\left(t\right)>0\mathrm{if}\mathsf{\Lambda }>\mu \mathrm{and}\mathcal{N}\left(\mathrm{t}\right)\mathrm{is}\mathrm{increasing},\\ \hfill {\mathcal{D}}_t^{\alpha }\mathcal{N}\left(t\right)<0\mathrm{if}\mathsf{\Lambda }<\mu \mathrm{and}\mathcal{N}\left(\mathrm{t}\right)\mathrm{is}\mathrm{decreasing}.\end{array}$$

We underline that the fractional derivative is a non-local operator, and we do not have a direct information about the behavior of the function under investigation.

The Equation (17) can be treated as a fractional logistic equation with a forcing term. We decided to focus on this equation in a different work. Although the problem can be studied from a numerical point of view, proceeding with a general approach seems to be hard.

Our results can be regarded as the special case $\mathsf{\Lambda }=\mu$, that is, constant population $\mathcal{N}\left(t\right)$, $t>0$. Indeed, for the Mittag-Leffler function, we have $E_{\alpha }\left(0\right)=1$, $\forall \alpha \in (0,1)$. Thus, we turn our problem in studying the fractional logistic equation. In particular, assuming $\mathsf{\Lambda }=\mu$, the problem reduces to

$$\begin{array}{cc}\hfill {\mathcal{D}}_t^{\alpha }\mathcal{N}\left(t\right)& =0\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill {\mathcal{D}}_t^{\alpha }\mathcal{Z}\left(t\right)& =\left(\mathsf{\Lambda }-\frac{\beta }{2}+\gamma \right)\mathcal{N}\left(t\right)-(\gamma +\mu )\mathcal{Z}\left(t\right)\left(1-\frac{\beta }{2\mathcal{N}\left(t\right)(\gamma +\mu )}\mathcal{Z}\left(t\right)\right),\hfill \end{array}$$

(20)

that is, $\mathcal{N}\left(t\right)$ is constant and satisfies (P1), as we can see from the first equation, and the second equation is the fractional logistic equation we are interested in with the suitable characterization of all parameters. Indeed, by considering $\mathcal{N}\left(t\right)=C$ with the corresponding compartmental $\mathsf{\Lambda }C$, $\beta C$, $\mu C$, $\gamma C$, the equations above take the form:

$$\begin{array}{cc}\hfill & \left\{\begin{array}{c}{\displaystyle C{\mathcal{D}}_t^{\alpha }\frac{\mathcal{S}\left(t\right)}{C}=\mathsf{\Lambda }C-\beta C\frac{\mathcal{S}\left(t\right)}{C}\frac{\mathcal{I}\left(t\right)}{C}+\gamma C\frac{\mathcal{I}\left(t\right)}{C}-\mu C\frac{\mathcal{S}\left(t\right)}{C}}\hfill \\ {\displaystyle C{\mathcal{D}}_t^{\alpha }\frac{\mathcal{I}\left(t\right)}{C}=\beta C\frac{\mathcal{S}\left(t\right)}{C}\frac{\mathcal{I}\left(t\right)}{C}-\gamma C\frac{\mathcal{I}\left(t\right)}{C}-\mu C\frac{\mathcal{I}\left(t\right)}{C}}\hfill \end{array}\right.\hfill \\ \hfill & \mathrm{with}C=\mathcal{S}\left(t\right)+\mathcal{I}\left(t\right),S\left(0\right)={\mathcal{S}}_0\mathrm{and}I\left(0\right)={\mathcal{I}}_0,\hfill \end{array}$$

(21)

and we get

$$\begin{array}{cc}\hfill & \left\{\begin{array}{c}{\displaystyle C{\mathcal{D}}_t^{\alpha }S\left(t\right)=\mathsf{\Lambda }C-\beta CS\left(t\right)I\left(t\right)+\gamma CI\left(t\right)-\mu CS\left(t\right)}\hfill \\ {\displaystyle C{\mathcal{D}}_t^{\alpha }I\left(t\right)=\beta CS\left(t\right)I\left(t\right)-\gamma CI\left(t\right)-\mu CI\left(t\right)}\hfill \end{array}\right.\hfill \\ \hfill & \mathrm{with}1=S\left(t\right)+I\left(t\right),S\left(0\right)=S_0\mathrm{and}I\left(0\right)=I_0,\hfill \end{array}$$

(22)

where $S_0={\mathcal{S}}_0/C$ and $I_0={\mathcal{I}}_0/C$. Remember that $I\left(t\right)=\mathcal{I}\left(t\right)/C$ is a percentage; by recalling that $\mathcal{Z}\left(t\right)=C-2\mathcal{I}\left(t\right)$ and $\mathsf{\Lambda }=\mu$, we obtain

$$\begin{array}{cc}\hfill {\displaystyle -2{\mathcal{D}}_t^{\alpha }\mathcal{I}\left(t\right)=}& \left(\mu +\gamma -\frac{\beta }{2}\right)C-(\gamma +\mu )(C-2\mathcal{I}\left(t\right))\left(1-\frac{\beta }{2C(\gamma +\mu )}(C-2\mathcal{I}\left(t\right))\right)\hfill \\ \hfill =& -\frac{\beta }{2}C+2(\gamma +\mu )\mathcal{I}\left(t\right)+\frac{\beta }{2C}{(C-2\mathcal{I}\left(t\right))}^2\hfill \\ \hfill =& 2(\gamma +\mu -\beta )\mathcal{I}\left(t\right)+2\frac{\beta }{C}{\mathcal{I}}^2\left(t\right),\hfill \end{array}$$

that is

$$\begin{array}{c}\hfill {\displaystyle -2C{\mathcal{D}}_t^{\alpha }I\left(t\right)=-2(\beta -(\gamma +\mu ))CI\left(t\right)+2\beta C{I}^2\left(t\right),}\end{array}$$

from which we recover

$$\begin{array}{c}\hfill {\displaystyle {\mathcal{D}}_t^{\alpha }\mathcal{I}\left(t\right)=\beta cI\left(t\right)-\beta I^2\left(t\right),}\end{array}$$

which is (4). We notice that in this characterization the carrying capacity c merits further investigations. Indeed, it must be $c\ne 1$. We are led to study both cases $c=0$ and $c\ne 0$. Since, in our formulation, $\mathcal{N}\left(t\right)=1$, we refer to $\mathcal{S}\left(t\right)$ and $\mathcal{I}\left(t\right)$ as percentages and use the symbol $S\left(t\right)$ and $I\left(t\right)$.

For $\alpha =1$, the Mittag-Leffler becomes the exponential $E_1\left((\mathsf{\Lambda }-\mu )t\right)=e^{(\mathsf{\Lambda }-\mu )t}$, whereas, for $\alpha \in (0,1)$, we have the following asymptotic behaviors for $\mathsf{\Lambda }\le \mu$,

$$\begin{array}{c}\hfill \frac{E_{\alpha }\left((\mathsf{\Lambda }-\mu )t^{\alpha }\right)}{e_0\left((\mathsf{\Lambda }-\mu )t^{\alpha }\right)}\to 1,\mathrm{as}t\to 0\mathrm{and}\frac{E_{\alpha }\left((\mathsf{\Lambda }-\mu )t^{\alpha }\right)}{e_{\infty }\left((\mathsf{\Lambda }-\mu )t^{\alpha }\right)}\to 1,\mathrm{as}t\to \infty ,\end{array}$$

where

$$\begin{array}{c}\hfill e_0\left((\mathsf{\Lambda }-\mu )t^{\alpha }\right)=exp\left(-|\mathsf{\Lambda }-\mu |\frac{t^{\alpha }}{\Gamma (1+\alpha )}\right),\mathrm{and}{e}_{\infty }\left((\mathsf{\Lambda }-\mu )t^{\alpha }\right)=\frac{1}{|\mathsf{\Lambda }-\mu |}\frac{t^{-\alpha }}{\Gamma (1-\alpha )}.\end{array}$$

For $\mathsf{\Lambda }>\mu$, the Mittag-Leffler (18) is an increasing function.

3. Proof of the Main Results

In this section, we collect the proof of the results presented in the work. From the theory of power series, we know that to each series representation with coefficients ${\left\{{\psi }_k\right\}}_k$ corresponds a radius of convergence $r_{\alpha }\in [0,\infty ]$ such that the series converges uniformly in $(0,r)$ for every $r<r_{\alpha }$. By the root test, we also have that

$$\begin{array}{c}\hfill r_{\alpha }={\left(\underset{k\to \infty }{lim}sup|\frac{{\psi }_k}{\Gamma (\alpha k+1)}{|}^{1/k}\right)}^{-1/\alpha },\end{array}$$

(23)

and the radius $r_{\alpha }$ obviously depends on the sequence ${\left\{{\psi }_k\right\}}_k$ and the order $\alpha \in (0,1)$ of the fractional derivative.

Proof of Theorem 1.

Similarly to the classical case, by the linearity (P3) of the Caputo derivative, we exploit $S\left(t\right)=1-I\left(t\right)$ to reduce problem (1) to

$${\mathcal{D}}_t^{\alpha }I\left(t\right)=\beta cI\left(t\right)\left(1-\frac{I\left(t\right)}{c}\right).$$

(24)

We rewrite (24) as

$${\mathcal{D}}_t^{\alpha }v\left(t\right)=\frac{1}{M^{\alpha }}v\left(t\right)(1-v\left(t\right)),$$

(25)

where $v\left(t\right)=I\left(t\right)/c$ and $M={\left(\beta c\right)}^{-1/\alpha }=b^{-1/\alpha }$. Equation (25) is the fractional logistic equation investigated in Reference [10], where the explicit solution is given for $M>1$ and $v\left(0\right)=1/2$ as

$$v\left(t\right)=\sum _{k\ge 0}\frac{E_k^{\alpha }}{M^{\alpha k}}\frac{t^{\alpha k}}{\Gamma (\alpha k+1)}.$$

(26)

In particular, the authors proved an estimate by below of the convergence ray $r_{\alpha }$. From (26), we recover $I\left(t\right)=cv\left(t\right)$, the solution of the α − SIS model.

Proof of Theorem 2.

By the linearity (P3) of the Caputo derivative and the fact that $S\left(t\right)=1-I\left(t\right)$, the problem (1) reduces to

$${\mathcal{D}}_t^{\alpha }I\left(t\right)=-\beta I^2\left(t\right).$$

(27)

Setting $u\left(t\right)=\beta I\left(t\right)$, we have that

$${\mathcal{D}}_t^{\alpha }u\left(t\right)=\beta {\mathcal{D}}_t^{\alpha }I\left(t\right)=-{\beta }^2{I}^2\left(t\right)=-u^2\left(t\right).$$

(28)

We prove that

$$u\left(t\right)=\sum _{k=0}^{\infty }{A}_k^{\alpha }\frac{t^{\alpha k}}{\Gamma (\alpha k+1)}$$

(29)

solves (28); hence, $I\left(t\right)=u\left(t\right)/\beta$ is the solution to (27).

To this end, we compute the Riemann-Liouville fractional derivative of $u\left(t\right)$ in (29), which is

$$\begin{array}{cc}\hfill {\mathcal{D}}_t^{\alpha }u\left(t\right)& =\sum _{k=0}^{\infty }{A}_k^{\alpha }\frac{t^{\alpha k-\alpha }}{\Gamma (\alpha k-\alpha +1)}\hfill \\ \hfill & =A_0^{\alpha }\frac{t^{-\alpha }}{\Gamma (1-\alpha )}+\sum _{k=0}^{\infty }{A}_{k+1}^{\alpha }\frac{t^{\alpha k}}{\Gamma (\alpha k+1)}\hfill \\ \hfill & =A_0^{\alpha }\frac{t^{-\alpha }}{\Gamma (1-\alpha )}+A_1^{\alpha }+A_2^{\alpha }\frac{t^{\alpha }}{\Gamma (\alpha +1)}+A_3^{\alpha }\frac{t^{2\alpha }}{\Gamma (2\alpha +1)}+A_4^{\alpha }\frac{t^{3\alpha }}{\Gamma (3\alpha +1)}+A_5^{\alpha }\frac{t^{4\alpha }}{\Gamma (4\alpha +1)}+\cdots .\hfill \end{array}$$

By (12), we have

$${\mathcal{D}}_t^{\alpha }u\left(t\right)=A_1^{\alpha }+A_2^{\alpha }\frac{t^{\alpha }}{\Gamma (\alpha +1)}+A_3^{\alpha }\frac{t^{2\alpha }}{\Gamma (2\alpha +1)}+A_4^{\alpha }\frac{t^{3\alpha }}{\Gamma (3\alpha +1)}+A_5^{\alpha }\frac{t^{4\alpha }}{\Gamma (4\alpha +1)}+\cdots .$$

(30)

Now, we compute $u^2\left(t\right)$:

$$\begin{array}{cc}\hfill u^2\left(t\right)& =\sum _{k=0}^{\infty }\sum _{s=0}^{\infty }{A}_k^{\alpha }{A}_s^{\alpha }\frac{t^{\alpha (k+s)}}{\Gamma (\alpha k+1)\Gamma (\alpha s+1)}\hfill \\ \hfill & =A_0^{\alpha }{A}_0^{\alpha }\hfill \\ \hfill & +\frac{2A_1^{\alpha }{A}_0^{\alpha }}{\Gamma (\alpha +1)}{t}^{\alpha }\hfill \\ \hfill & +\left(\frac{A_1^{\alpha }{A}_1^{\alpha }}{\Gamma (\alpha +1)\Gamma (\alpha +1)}+\frac{2A_0^{\alpha }{A}_2^{\alpha }}{\Gamma (2\alpha +1)}\right)t^{2\alpha }\hfill \\ \hfill & +\left(\frac{A_1^{\alpha }{A}_2^{\alpha }}{\Gamma (\alpha +1)\Gamma (2\alpha +1)}+\frac{2A_0^{\alpha }{A}_3^{\alpha }}{\Gamma (3\alpha +1)}\right)t^{3\alpha }\hfill \\ \hfill & +\left(\frac{A_2^{\alpha }{A}_2^{\alpha }}{\Gamma (2\alpha +1)\Gamma (2\alpha +1)}+\frac{2A_1^{\alpha }{A}_3^{\alpha }}{\Gamma (\alpha +1)\Gamma (3\alpha +1)}+\frac{2A_0^{\alpha }{A}_4^{\alpha }}{\Gamma (4\alpha +1)}\right)t^{4\alpha }+\cdots \hfill \end{array}$$

(31)

By (30) and (31), and by $A_0^{\alpha }=1/2$, we have

$$\begin{array}{cc}\hfill A_1^{\alpha }& =-A_0^{\alpha }{A}_0^{\alpha }=-1/4\hfill \\ \hfill A_2^{\alpha }& =-2A_1^{\alpha }{A}_0^{\alpha }\frac{\Gamma (\alpha +1)}{\Gamma (\alpha +1)}\hfill \\ \hfill A_3^{\alpha }& =A_1^{\alpha }{A}_1^{\alpha }\frac{\Gamma (2\alpha +1)}{\Gamma (\alpha +1)\Gamma (\alpha +1)}+2A_0^{\alpha }{A}_2^{\alpha }\frac{\Gamma (2\alpha +1)}{\Gamma (2\alpha +1)}\hfill \\ \hfill A_4^{\alpha }& =A_1^{\alpha }{A}_2^{\alpha }\frac{\Gamma (3\alpha +1)}{\Gamma (\alpha +1)\Gamma (2\alpha +1)}+2A_0^{\alpha }{A}_3^{\alpha }\frac{\Gamma (3\alpha +1)}{\Gamma (3\alpha +1)}\hfill \\ \hfill A_5^{\alpha }& =A_2^{\alpha }{A}_2^{\alpha }\frac{\Gamma (4\alpha +1)}{\Gamma (2\alpha +1)\Gamma (2\alpha +1)}+2A_1^{\alpha }{A}_3^{\alpha }\frac{\Gamma (4\alpha +1)}{\Gamma (\alpha +1)\Gamma (3\alpha +1)}+2A_0^{\alpha }{A}_4^{\alpha }\frac{\Gamma (4\alpha +1)}{\Gamma (4\alpha +1)};\hfill \end{array}$$

thus,

$$A_{k+1}^{\alpha }=-\sum _{j=0}^k\frac{\Gamma (k\alpha +1)}{\Gamma \left(\right(k-j)\alpha +1)\Gamma (j\alpha +1)}{A}_j^{\alpha }{A}_{k-j}^{\alpha }.$$

(32)

We use the fact that $\forall k\in \{0,1,\dots ,\}$,

$$\begin{array}{c}\hfill \frac{\Gamma (k\alpha +1)}{\Gamma \left(\right(k-j)\alpha +1)\Gamma (j\alpha +1)}=:R_k\le \Gamma (k\alpha +1).\end{array}$$

From the definition above of the coefficients ${\left\{A_k^{\alpha }\right\}}_k$, we get

$$\begin{array}{cc}\hfill \left|\frac{A_{k+1}^{\alpha }}{\Gamma \left(\right(k+1)\alpha +1)}\right|\le & \frac{1}{\Gamma \left(\right(k+1)\alpha +1)}\sum _{j=0}^k{R}_j\left|A_j^{\alpha }{A}_{k-j}^{\alpha }\right|\hfill \\ \hfill \le & \frac{\Gamma (k\alpha +1)}{\Gamma \left(\right(k+1)\alpha +1)}\sum _{j=0}^k\left|A_j^{\alpha }{A}_{k-j}^{\alpha }\right|.\hfill \end{array}$$

By iteration, we obtain that $A_k^{\alpha }\sim {\left|A_0^{\alpha }\right|}^k$. Since $(0,1)\ni A_0^{\alpha }\le 1/A_0^{\alpha }$, we write

$$\begin{array}{cc}\hfill \left|\frac{A_{k+1}^{\alpha }}{\Gamma \left(\right(k+1)\alpha +1)}\right|\le & \frac{\Gamma (k\alpha +1)}{\Gamma \left(\right(k+1)\alpha +1)}(k+1){\left(\frac{1}{A_0^{\alpha }}\right)}^k=:{\vartheta }_k,k\in {\mathbb{N}}_0.\hfill \end{array}$$

We now consider the fact that

$$\begin{array}{c}\hfill \frac{x^{x-\gamma }}{e^{x-1}}<\Gamma \left(x\right)<\frac{x^{x-1/2}}{e^{x-1}},x>1\end{array}$$

(where $\gamma \approx 0.5$ is the Mascheroni constant), and we get

$$\begin{array}{c}\hfill \sqrt[k]{|{\vartheta }_k|}\sim \frac{1}{|A_0^{\alpha }|}{\left((k+1)\frac{{(k\alpha +1)}^{k\alpha +1/2}}{{((k+1)\alpha +1)}^{(k+1)\alpha +1-\gamma }}\right)}^{1/k}.\end{array}$$

Since

$$\begin{array}{cc}\hfill {(k\alpha +1)}^{\frac{1}{k}(k\alpha +1/2)}\sim & exp\left(\left(\alpha +\frac{1}{2k}\right)ln(k\alpha +1)\right)\hfill \end{array}$$

and

$$\begin{array}{c}\hfill {((k+1)\alpha +1)}^{\frac{1}{k}((k+1)\alpha +1-\gamma )}\sim exp\left(\left(\alpha +\frac{1-\gamma }{k}\right)ln((k+1)\alpha +1)\right),\end{array}$$

we get that

$$\begin{array}{c}\hfill \sqrt[k]{|{\vartheta }_k|}\sim \frac{1}{|A_0^{\alpha }|}.\end{array}$$

Thus, we get the radius of convergence

$$\begin{array}{c}\hfill r_{\alpha }^{\vartheta }={\left(\underset{k\to \infty }{lim}|{\vartheta }_k{|}^{1/k}\right)}^{-1/\alpha }={\left(|A_0^{\alpha }|\right)}^{1/\alpha }\end{array}$$

for the series

$$\begin{array}{c}\hfill \sum _{k\ge 0}{\vartheta }_k.\end{array}$$

The convergence of the majorant series determines the uniform convergence in $(0,r_{\alpha })\subset (0,r_{\alpha }^{\vartheta })$ of the series we are interested in. This concludes the proof by considering $I=u/\beta$.

Remark 1. 

The solution in Theorem 1 has been given only for the initial datum $c/2$. This is because of the representation given in Reference [10] in terms of Euler polynomials. Let us focus on the solution in Theorem 2. Taking $A_0^{\alpha }\in (0,1)$, we see that, setting

$$\begin{array}{c}\hfill v\left(t\right)=u(t/2^q)=\sum _{n\ge 0}{A}_k^{\alpha }\frac{{(t/2^q)}^{n\alpha }}{\Gamma (n\alpha +1)},t\in K^q\subseteq (0,r_{\alpha }^q),\end{array}$$

where

$$\begin{array}{c}\hfill q=\left\{\begin{array}{cc}{\displaystyle \frac{1}{A_0^{\alpha }},}\hfill & A_0^{\alpha }<\frac{1}{2}\hfill \\ {\displaystyle 4+\frac{1}{2}\left(\frac{1}{A_0^{\alpha }}-4\right),}\hfill & A_0^{\alpha }\ge \frac{1}{2},\hfill \end{array}\right.\end{array}$$

we obtain $r_{\alpha }^q=2^q{\left(|A_0^{\alpha }|\right)}^{1/\alpha }$. This is the solution in $(0,r_{\alpha }^q)$ to

$$\begin{array}{c}\hfill {\mathcal{D}}_t^{\alpha }v=-\frac{1}{2^q}{v}^2,v\left(0\right)=A_0^{\alpha }\in (0,1)\end{array}$$

(see the proof of Theorem 3.1 in Reference [10]). In the special case $\alpha =1$, we know that

$$\begin{array}{c}\hfill w\left(t\right)={\left(\frac{1}{A_0}-t\right)}^{-1}=A_0\sum _{k\ge 0}{(-A_0)}^k{t}^{k}t\in (0,1/A_0)\end{array}$$

solves $w^{\prime }=-w^2$ with $w\left(0\right)=A_0\in (0,1)$. In particular, for $A_0^{\alpha }=A_0=1/2$, we obtain convergence in any compact sets $K\subset (0,2)$ for both solutions v and w. This underlines the fact that, by introducing non-locality, we may deal with solutions quite far from their non-linear analogues.

4. Numerical Comparison

In this section, we proceed with the validation of the previous results on the fractional SIS model by means of numerical approximations, and we analyze the effects of fractional derivatives by comparing the ordinary and fractional SIS model.

4.1. Numerical Approximation

The explicit solution (5)–(6) to the fractional SIS model (1) for $c\ne 0$ is defined for $b^{1/\alpha }<1$ and initial datum $I_0=c/2$. The explicit solution (8)–(9) to the fractional SIS model (1) for $c=0$ is defined for the initial datum $I_0=1/\left(2\beta \right)$. In order to compute the solution to the fractional SIS model for any set of parameters and any initial datum, we propose and compare two numerical schemes to approximate (1). To this end, let us consider the following problem:

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

(33)

on a time interval $[0,T]$ uniformly divided into $N+1$ time steps of length $\Delta t$. Our aim is to define the discrete solution $u_n=u\left(t_n\right)$ for $n=1,\cdots ,N$, where $t_n=n\Delta t$, and $u_0$ is known.

We refer to the following method as the Method 1. Following Reference [19], we observe that

$$\begin{array}{cc}\hfill {\mathcal{I}}^{1-\alpha }{u}^{\prime }& =f\left(u\right)\hfill \\ \hfill {\mathcal{I}}^{\alpha }{\mathcal{I}}^{1-\alpha }{u}^{\prime }& ={\mathcal{I}}^{\alpha }f\left(u\right)\hfill \\ \hfill {\mathcal{I}}^1{u}^{\prime }& ={\mathcal{I}}^{\alpha }f\left(u\right);\hfill \end{array}$$

thus, we rewrite (33) as

$$u\left(t\right)=u\left(0\right)+{\mathcal{I}}^{\alpha }f\left(u\right).$$

(34)

We introduce a Predictor-Evaluate-Corrector-Predictor (PECE) method [20]. Specifically, we use the implicit one-step Adams-Moulton method [21] (Chapter 11), i.e.,

$$u_{n+1}=u_0+\frac{1}{\Gamma \left(\alpha \right)}\left(\sum _{j=0}^n{a}_{j,n+1}f\left(u_j\right)+a_{n+1,n+1}f\left({\tilde{u}}_{n+1}\right)\right),$$

(35)

where the coefficients $a_{j,n+1}$ and ${\tilde{u}}_{n+1}$ are defined below.

First of all, we compute the term ${\tilde{u}}_{n+1}$ with the one-step Adams-Bashforth method. We introduce $g\left(s\right)=f\left(u\right(s\left)\right)$ and $g_{n+1}$ as a piece-wise linear function which interpolates g on the nodes $t_j$, $j=0,\cdots ,n+1$. We approximate the integral term of (34) with the product rectangle rule, i.e.,

$${\int }_{t_0}^{t_{n+1}}{(t_{n+1}-s)}^{\alpha -1}g\left(s\right)ds\approx \sum _{j=0}^n{b}_{j,n+1}g\left(t_j\right),$$

where

$$b_{j,n+1}={\int }_{t_j}^{t_{j+1}}{(t_{n+1}-s)}^{\alpha -1}ds=\frac{1}{\alpha }({(t_{n+1}-t_j)}^{\alpha }-{(t_{n+1}-t_{j+1})}^{\alpha }).$$

In particular, for our uniform discretization of the time interval $[0,T]$, we have

$$b_{j,n+1}=\frac{\Delta t^{\alpha }}{\alpha }({(n+1-j)}^{\alpha }-{(n-j)}^{\alpha }).$$

Therefore,

$${\tilde{u}}_{n+1}=u_0+\frac{1}{\Gamma \left(\alpha \right)}\sum _{j=0}^n{b}_{j,n+1}f\left(u_j\right).$$

(36)

Now, we compute the coefficients $a_{j,n+1}$; thus, we approximate ${\mathcal{I}}^{\alpha }g$ as

$${\int }_{t_0}^{t_{n+1}}{(t_{n+1}-s)}^{\alpha -1}g\left(s\right)ds\approx {\int }_{t_0}^{t_{n+1}}{(t_{n+1}-s)}^{\alpha -1}{g}_{n+1}\left(s\right)ds.$$

(37)

By using the product trapezoidal quadrature formula on the nodes $t_j$, Equation (37) becomes

$${\int }_{t_0}^{t_{n+1}}{(t_{n+1}-s)}^{\alpha -1}{g}_{n+1}\left(s\right)ds=\sum _{j=0}^{n+1}{a}_{j,n+1}g\left(t_j\right),$$

where $a_{j,n+1}$ are defined as

$$a_{j,n+1}={\int }_{t_{j-1}}^{t_j}\frac{s-t_{j-1}}{t_j-t_{j-1}}{(t_{n+1}-s)}^{\alpha -1}ds+{\int }_{t_j}^{t_{j+1}}\frac{t_{j+1}-s}{t_{j+1}-t_j}{(t_{n+1}-s)}^{\alpha -1}ds.$$

We observe that, from integration by parts, we have

$$\begin{array}{cc}\hfill {\int }_{t_{j-1}}^{t_j}\frac{s-t_{j-1}}{t_j-t_{j-1}}{(t_{n+1}-s)}^{\alpha -1}ds& =-\frac{{(t_{n+1}-t_j)}^{\alpha }}{\alpha }+{\int }_{t_{j-1}}^{t_j}\frac{{(t_{n+1}-s)}^{\alpha }}{\alpha (t_j-t_{j-1})}ds\hfill \\ \hfill {\int }_{t_j}^{t_{j+1}}\frac{t_{j+1}-s}{t_{j+1}-t_j}{(t_{n+1}-s)}^{\alpha -1}ds& =\frac{{(t_{n+1}-t_j)}^{\alpha }}{\alpha }-{\int }_{t_j}^{t_{j+1}}\frac{{(t_{n+1}-s)}^{\alpha }}{\alpha (t_{j+1}-t_j)}ds,\hfill \end{array}$$

and, therefore,

$$\begin{array}{cc}\hfill a_{0,n+1}& =\frac{{(t_{n+1}-t_0)}^{\alpha }}{\alpha }-{\int }_{t_0}^{t_1}\frac{{(t_{n+1}-s)}^{\alpha }}{\alpha (t_1-t_0)}ds\hfill \\ \hfill a_{n+1,n+1}& ={\int }_{t_n}^{t_{n+1}}\frac{{(t_{n+1}-s)}^{\alpha }}{\alpha (t_{n+1}-t_n)}ds\hfill \\ \hfill a_{j,n+1}& ={\int }_{t_{j-1}}^{t_j}\frac{{(t_{n+1}-s)}^{\alpha }}{\alpha (t_j-t_{j-1})}ds-{\int }_{t_j}^{t_{j+1}}\frac{{(t_{n+1}-s)}^{\alpha }}{\alpha (t_{j+1}-t_j)}ds\mathrm{for}j=1,\cdots ,n.\hfill \end{array}$$

Finally, in our uniform grid, the coefficients are

$$\begin{array}{cc}\hfill a_{0,n+1}& =\frac{\Delta t^{\alpha }}{\alpha (\alpha +1)}(n^{\alpha +1}-(n-\alpha ){(n+1)}^{\alpha })\hfill \end{array}$$

(38)

$$\begin{array}{cc}\hfill a_{n+1,n+1}& =\frac{\Delta t^{\alpha }}{\alpha (\alpha +1)}\hfill \end{array}$$

(39)

$$\begin{array}{cc}\hfill a_{j,n+1}& =\frac{\Delta t^{\alpha }}{\alpha (\alpha +1)}({(n-j+2)}^{\alpha +1}-2{(n-j+1)}^{\alpha +1}+{(n-j)}^{\alpha +1})\mathrm{for}j=1,\cdots ,n.\hfill \end{array}$$

(40)

Remark 2. 

The numerical scheme described above works for any $\alpha \in [0,1]$.

We now introduce a method to which we refer as Method 2. Let $\alpha \in (0,1)$. In Reference [22], the authors give the following approximation of the Caputo derivative

$${\mathcal{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),$$

(41)

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(r\right)=r^{1-\alpha }-{(r-1)}^{1-\alpha }\mathrm{for}r\ge 1$. The numerical scheme to solve (33) 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).$$

(42)

We refer to Reference [22] for further details on the properties of the scheme.

Remark 3. 

The numerical scheme above described works for $\alpha \in (0,1)$, with the extreme values excluded.

To summarize, in this section, we have introduced two numerical schemes, which we denote here by $M_1$ and $M_2$ for notational convenience. The solution to the fractional SIS model (1) with the first numerical scheme (that is Method 1) is

$$\begin{array}{cc}\hfill I\left(t_{n+1}\right)& =M_1\left(I\left(t_n\right)\right)\hfill \end{array}$$

(43)

$$\begin{array}{cc}\hfill S\left(t_{n+1}\right)& =1-I\left(t_{n+1}\right),\hfill \end{array}$$

(44)

where $M_1$ is defined in (35), and the solution with the second numerical scheme (that is Method 2) is

$$\begin{array}{cc}\hfill I\left(t_{n+1}\right)& =M_2\left(I\left(t_n\right)\right)\hfill \end{array}$$

(45)

$$\begin{array}{cc}\hfill S\left(t_{n+1}\right)& =1-I\left(t_{n+1}\right),\hfill \end{array}$$

(46)

where $M_2$ is defined in (42) and $n=1,\cdots ,N$. Note that the function $f\left(u\right)$ in (33), used for both the numerical schemes, is defined as $f\left(u\right)=\beta cu-\beta u^2$, while $u_0=I_0$.

Remark 4. 

The computational complexity is essentially due to the fact that fractional derivatives are non-local operators. In order to improve the speed of the iterative method the short memory principle of Podlubny [23] (p. 203) for Caputo derivative can be helpfully considered. However, the price for the reduced complexity is given by the loss in the order of accuracy (see, for example, Reference [24,25]).

4.2. Numerical Tests

In this section we compare the solutions to the fractional SIS model (1) computed with the explicit representation and the two numerical schemes, testing both the case $c\ne 0$ and $c=0$. In what follows, we denote by

• $I^C,S^C$ the solutions to the SIS model, our aim is to show the correspondence with the case $\alpha =1$,
• $I^F,S^F$ the solutions (5)–(6) or (8)–(9) to the fractional SIS model (1) defined by Theorems 1 or 2 respectively (depending on the carrying capacity c),
• $I_1^N,S_1^N$ the numerical solutions (43)–(44) computed with the methodology proposed as Method 1,
• $I_2^N,S_2^N$ the numerical solutions (45)–(46) computed with the methodology proposed as Method 2.

4.2.1. Test with $c\ne 0$

We start our numerical analysis with the case of carrying capacity $c\ne 0$. We fix this set of parameters: $\beta =0.7$, $\gamma =0.05$, $\mu =0.12$, $\sigma =4$, and $c=0.75$. The initial data are $I\left(0\right)=c/2$ and $S\left(0\right)=1-I\left(0\right)$, the final time is $T=5$, and the time step $\Delta t=0.05$.

First of all, we compare the exact fractional solutions (5)–(6) and the two numerical solutions (43)–(44) and (45)–(46) for $\alpha =0.99$, which approximately corresponds to the classical derivative. Note that we do not use $\alpha \equiv 1$ since the second numerical scheme works for $\alpha \in (0,1)$, as already observed in Remark 3. In Figure 1, we show the results. As expected, the exact fractional solution and the two numerical solutions to (1) overlap the solution for $\alpha =1$.

In Figure 2 and Figure 3, we show the results obtained with $\alpha =0.7$ and $\alpha =0.3$. In the first case, the two density curves are closer each other and the intersection point between them slightly moves to the right with respect to the solution shown in Figure 1. Such behavior is further emphasized by lower values of $\alpha$, as shown for example in Figure 3. Note that, in both cases, the three methodologies produces almost identical results.

To further investigate on the three methodologies, we compute the $L^{\infty }$-norm of the difference between the exact fractional solutions (5)–(6) and the two numerical solutions (43)–(44) and (45)–(46) and between the two numerical solutions each others, as shown in

Table 1. Comparison of the $L^{\infty }$-norm between the solutions computed with the three methodologies for different values of $\alpha$.

($\mathit{\alpha }$)	${∥{\mathit{I}}^{\mathit{F}}-{\mathit{I}}_1^{\mathit{N}}∥}_{\mathit{\infty }}$	${∥{\mathit{I}}^{\mathit{F}}-{\mathit{I}}_2^{\mathit{N}}∥}_{\mathit{\infty }}$	${∥{\mathit{I}}_1^{\mathit{N}}-{\mathit{I}}_2^{\mathit{N}}∥}_{\mathit{\infty }}$
0.99	1 × 10${}^{-5}$	9 × 10${}^{-4}$	9 × 10${}^{-4}$
0.7	1 × 10${}^{-5}$	2 × 10${}^{-3}$	2 × 10${}^{-4}$
0.3	3 × 10${}^{-5}$	8 × 10${}^{-3}$	8 × 10${}^{-3}$

. We observe that the errors range from orders of ${10}^{-5}$ to ${10}^{-3}$, increasing with respect to the decrease of $\alpha$. This fact further certifies the similarity between the three proposed methodologies.

4.2.2. Test with $c=0$

We focus now on the case of carrying capacity $c=0$. We fix this set of parameters: $\beta =0.7$, $\gamma =0.07$, $\mu =0.63$, $\sigma =1$ and $c=0$. Moreover, the initial data are $I\left(0\right)=1/\left(2\beta \right)$ and $S\left(0\right)=1-I\left(0\right)$, the final time is $T=1$, and the time step $\Delta t=0.01$.

In Figure 4, we compare the exact fractional solutions (8)–(9) and the two numerical solutions (43)–(44) and (45)–(46) for $\alpha =0.99$. Again, we observe that the fractional solutions, both explicit and numerical, perfectly overlap the solution to the SIS model. In Figure 5, we show the results obtained with $\alpha =0.7$. Analogously to the example with $c\ne 0$, the point of intersection between the two densities of population slightly moves to the right with respect to the solution shown in Figure 4. Moreover, the three different methodologies produce again almost identical results. Finally, in Figure 6, we show the results obtained with $\alpha =0.5$. In this case, the explicit fractional solutions (8)–(9) blow up in finite time, since the final time T is greater than the radius of convergence, while the two numerical solutions show that the intersection point between the two curves further moves to the right with respect to Figure 5.

5. Conclusions

In this work, we studied the fractional SIS model with constant population size. We proposed an explicit representation of the solution to the fractional model under particular assumptions on parameters and initial data. By considering the basic reproduction number, we rearrange the SIS model and obtain a logistic equation. In the new formulation of the problem, the carrying capacity has a new meaning based on the parameters of the SIS model. We exploit such a formulation in order to study the fractional SIS model and obtain a fruitful characterization of the problem, despite many difficulties introduced by non-locality. In our formulation, the carrying capacity can equal zero, and this brings our attention to a different non-linear problem which, in turn, is related to the underlined SIS model. We introduced two different numerical schemes to approximate the model and perform numerical simulations, with which we tested the proposed explicit solution.