A Mathematical Study of a Coronavirus Model with the Caputo Fractional-Order Derivative

: In this work, we introduce a minimal model for the current pandemic. It incorporates the basic compartments: exposed, and both symptomatic and asymptomatic infected. The dynamical system is formulated by means of fractional operators. The model equilibria are analyzed. The theoretical results indicate that their stability behavior is the same as for the corresponding system formulated via standard derivatives. This suggests that, at least in this case for the model presented here, the memory effects contained in the fractional operators apparently do not seem to play a relevant role. The numerical simulations instead reveal that the order of the fractional derivative has a deﬁnite inﬂuence on both the equilibrium population levels and the speed at which they are attained.


Introduction
Among the new tools that have recently come into the limelight in scientific research, fractional derivatives play an important role, in view of the fact that, in a sense, they embed memory effects in dynamical systems [1][2][3][4]. This concept has been successfully applied in mathematical physics [5], but more recently it has also been used in other contexts. In the biological framework, ecological models incorporating this type of operator have been introduced in [6], and for epidemic situations, in [7]. In view of the COVID-19 pandemic that has been affecting humanity for almost two years, in the last months several papers have been published on this topic. Some of them attempt to set up and analyze models for this epidemic based on the concept of fractional derivatives. The simplest epidemiological model, SI, consists of two classes, namely susceptibles S and infectives I. If possible disease recovery is allowed, the class of healed individuals is introduced as R. These individuals are immune from the disease, at least temporarily. The model thus becomes type SIR. If, additionally, the recovered individuals lose disease immunity, they again become prone to being infected and thus they return to the class of susceptibles. In such cases, therefore, disease relapses are allowed, and the SIRS model is obtained. A further step forward consists of introducing another compartment, hosting the individuals that have been exposed to the disease. It is indeed well-known that, in the case of measles, individuals that have contracted the disease do not show symptoms until about at least a week after being infected. Meanwhile, they are very much able to transmit the disease, especially in the very first few days after infection. The exposed class therefore has paramount importance in transmissible disease modeling. The simplest nontrivial compartment model that shows this feature is the SEIR model. It has been taken into consideration, together with a fractional derivative approach, in [8,9]. Reference [10] also considers this formulation, but in addition time delays are explicitly built into the model. Another extension of the SEIR system is considered in [11], where the class P of "deprived and marginalised" individuals is also considered. The latter are individuals that are less able to cope with strict distancing measures and therefore contribute to the further spread of the disease. Apparently, they are "recruited" from the infected by a simple linear transition. In other words, the surviving infected exit this class either by recovering or by becoming part of P. In [12], instead, three additional classes are introduced, thereby forming the SEIQRDP model. Apart from the class D of deceased individuals, quarantined individuals Q are introduced, as well as the "insusceptible", P. In the model, disease mortality is allowed only among quarantined, thereby most likely implying that infected that die of the disease are accounted for as having been recognized as disease carriers, and therefore put in isolation before dying. The Atagana-Baleanu fractional derivative is employed instead of the Caputo fractional derivative in [13]. Here, the infected may migrate to the quarantined Q class. Individuals in isolation may possibly become confirmed cases and therefore are lumped into the new class C. Notably, individuals belonging to both Q and C classes are assumed not to be able to disseminate the disease further, which may be a strong assumption. In Italy, at least in the early phases of the pandemic, most of the contagions involved doctors and nurses working in hospitals. This indicates that virus leakage was nevertheless possible in supposedly confined environments. After collecting data for various compartments from two countries, in [14] parameters are estimated, on the basis of which, initially a standard dynamical system and later on some fractional models, are built. In particular, the Caputo, Caputo-Fabrizio, Atangana-Baleanu and the fractal-fractional nonlocal operators are employed. Five classes of infected are considered, of which the first one essentially represents the exposed in our above terminology. The remaining ones consider symptomatic individuals that have or have not been tested and found to be positive or not. In [15], again the same previous five classes of infected individuals are considered, those having been tested or not, or those under treatment. The model also introduces a time-dependent transmission rate that tries to incorporate the social distancing measures, including the wearing of masks, and social contacts restrictions, that is, shop closures, curfew and so forth. This attempts to model the ways of fighting the epidemic diffusion that have been enforced by many governments worldwide, although the implementation times differed from country to country. Specifically, in [15], a model is proposed that extends the possible disease transmission mechanism of [14] to two other compartments which were not able to spread it. In addition to the equilibria analysis, a numerical scheme for the integration of the fractional system is developed. In both [14,15], the major novelty at the modeling level appears to be the introduction of the compartment of the vaccinated, V. In [16], an early model for treatment is also formulated using standard derivatives.
In all the papers considered above, however, no real asymptomatic class appears in the model. By this terminology, we mean people that are unrecognized carriers, in that they do not show symptoms. This class has enormous importance, given the fact that individuals in it believe that are not infectious and therefore think that they may mingle without restraint with susceptibles. However, not being recognized as disease-carriers, they contribute to the fast diffusion of the epidemic. A model with standard derivatives taking asymptomatics into account has been formulated in [17]. Symptomatic and asymptomatic infectious are partitioned into two different classes. The model is formulated via standard derivatives. Instead, in [18] a fractional operator model is presented, but the main aim of the paper consists of devising a suitable stable numerical scheme, which essentially relies on the use of the Laplace transform. The basic reproduction number is evaluated and, if smaller than one, it is shown that the epidemic is eradicated. Parameters for the simulations are taken from the literature or estimated by data fitting. The asymptomatic class is also considered in [19], where the main aim is the study of a numerical algorithm based on the Riesz wavelets is proposed.
In [20], quarantined individuals-coming from the exposed class, the symptomatic and asymptomatic infected, the isolated and treated-are considered in addition to the standard classes of the SEIR model. Notably, the model does not seem to account for disease mortality. For infected I the parameter ρ, according to the parameter list, should be interpreted to be the same as γ. It describes the migration from class I to treatment. At least both the infected and treated individuals should experience disease mortality. Indeed, if by the term "treated" the authors mean hospitalized, not including disease mortality does not appear to match observations. As a matter of fact, in the case of covid, hospitalized mortality has had a very big impact on societies worldwide. A similar remark may perhaps be made on the class P of "isolated" individuals into which the "quarantined" (being so after exposure) and infected migrate. After all, these people are still affected by the disease. In the model, the isolated and quarantined are assumed to be able to propagate the disease. However, treated individuals are not. In such cases, the same remark on the Italian hospitals' situation mentioned above may apply. The paper investigates the existence of the solutions and the feasibility and stability of the system equilibria.
We have therefore chosen a simple model that is formulated via fractional derivatives in the line of the above research lines. It includes only the basic epidemiological classes, because the aim of modeling is to elucidate, if possible, how very complex behaviors are generated, using the minimal feasible assumptions. These must contain the basic mechanisms underlying the variables' relationships, but should not remove the system's basic features (Occam's razor).
The model we propose differs in some more important or smaller aspects from each one of the systems presented in the above papers. In some cases it is simpler, which may be a shortcoming if real data need to be analyzed. On the other hand, however, it may represent an advantage in that it is well-known that minimal models sometimes help in highlighting system features whose real origin may be obscured in more complicated ones. This paper is organized as follows. The next Section 2 contains the model description. The basic analytical properties of the system are presented in Section 3. A numerical simulation is presented in Section 4 together with a brief final discussion.

Material and Methods
We consider the following fractional-order epidemiological model incorporating a spread of corona virus. We consider the following classes of individuals: susceptibles S, namely the individuals who have not yet been exposed to the virus; exposed E, people who have been infected by the virus but are in the incubation period, in which they cannot yet spread the disease; symptomatic infected I, individuals that manifest symptoms and can communicate the disease; asymptomatic infected A, those persons that can communicate the disease even without having explicit symptoms; the removed class R, which includes the people that recovered from the disease. We further let N(t) denote the total population. Thus, N(t) = S(t) + E(t) + I(t) + A(t) + R(t). The model reads as follows: where D q is the standard Caputo differentiation with q ∈ (0, 1). For the benefit of the reader, we recall the definition of the Caputo fractional derivative of order q [1,2]: In the above model, all individuals experience natural mortality at rate d p ; susceptibles are recruited at constant rate Λ, these being the only two demographic features of the model. The epidemic part of the model considers disease transmission among susceptibles and both types of infected, at rate β I , through mass action. Note that the factor k modulates the possible difference in the ability to spread the disease from asymptomatic people versus symptomatic ones. After contagion, a susceptible migrates into the class of the exposed, from which, after activation of the disease, it moves on to the infected compartments, at possibly different rates ω p and ω p ; α instead tells the proportion of exposed that will never show symptoms, that is, they move to class A. From this class, at rate ξ, individuals may still be able to become symptomatic and migrate to the I class. Both types of infected individuals may recover at rates γ p and γ p and are subject to disease-related mortality, for which we assume two possibly different rates, µ for symptomatic and ν for asymptomatic people. All the parameters are non negative and their meaning is summarized in Table 1. The basic mechanisms underlying the model are shown in Figure 1.
transmissibility ratio between asymptomatics and symptomatics, ν disease-related mortality for asymptomatics, µ disease-related mortality for symptomatic individuals, ω p progression rate from exposed to symptomatic, ω p progression rate from exposed to asymptomatic, α fraction of exposed that turn asymptomatic, ξ progression rate from asymptomatic to symptomatic, γ p recovery rate from symptomatic infection, γ p recovery rate from asymptomatic infection.
The system (1) is completed by the following initial conditions

Results
This section studies the existence, uniqueness, non-negativity and boundedness of the solutions of a fractional order epidemiological model (1). In addition, the stability analysis of fractional order epidemiological model (1) is also performed.

Existence and Uniqueness
The sufficient condition for the existence and uniqueness of the solution of a fractional order system (1) are as follows: Theorem 1. For each non-negative initial conditions, there exists a unique solution of fractional order system (1).
Proof. We seek a sufficient condition for the existence and uniqueness of the solutions of fractional order system (1) in the region Ω × (0, T] where The approach used in [6] is adopted. Consider a mapping For any X, X, it follows from (3) that Thus, G(X) satisfies the Lipschitz condition. Consequently, the existence and uniqueness of the solution of the fractional order system (1) follows.

Non-Negativity and Boundedness
The solutions of the fractional order system (1) represent the densities of the interacting populations and must therefore be non-negative and bounded. These features are ensured by the following results.

Theorem 2.
All the solutions of fractional order system (1), which start in R 5 + are uniformly bounded and non-negative.
Proof. We follow the approach used by [6]. Define the function By using the standard comparison theorem for fractional order [21], where E q is the Mittag-Leffler function. According to Lemma 5 and Corollary 6 in [21], Therefore, all the solutions of fractional order system (1) starting in R 5 + are confined to the region Σ, where Next, we show that the solutions of the fractional order system (1) are non-negative. From the first equation of system (1) we have According to the standard comparison theorem for fractional order [21], and the positivity of the Mittag-Leffler function E q,1 (t) > 0, for any q ∈ (0, 1) [22], From the second equation of system (1) we find From the third equation of system (1) we get From the fourth equation of system (1) we obtain From the last equation of system (1) we have Therefore, Thus, it has been proved that the solutions of system (1) are non-negative.

Equilibrium Points and Local Stability
In this section, we investigate the local stability of equilibrium points. The fractionalorder system (1) has the following equilibrium points: It is always feasible.

The Basic Reproduction Number
The basic reproduction number R 0 for our fractional model is found using the next generation matrix method [23].
The reduced system may be written in compact form as: The Jacobian matrices of F(X) and V(X) at the disease-free equilibrium point P 0 are We find that The next generation matrix is Next, we will discuss the stability of the equilibrium points of system (1). At the point P(S, E, I, A, R), the Jacobian matrix of system (1) is given by Using the Jacobian matrix (12) and the Matignon condition [4,7], the local stability of the equilibrium points of the fractional-order system (1) is investigated. We have the following results.
Theorem 3. The coronavirus-free equilibrium P 0 of the fractional-order system (1) is locally asymptotically stable if and unstable if Here, recall that R 0 denotes the basic reproduction number.

Proof.
Letting the Jacobian matrix (12) around the coronavirus-free equilibrium P 0 is The eigenvalues of the Jacobian matrix J(P 0 ) around the coronavirus-free equilibrium P 0 are λ 1 = −d p of multiplicity order two and the roots of the characteristic polynomial of the minor matrix of J(P 0 ) given by where a i , i = 0, . . . , 2 are given in (15). It is evident that a 2 > 0. From condition (13), the following equations are also satisfied: Thus, a i > 0, i = 0, . . . , 2 and a 2 a 1 > a 0 . From the Routh-Hurwitz criterion, all the roots λ i of the characteristic Equation (16) have negative real parts. By using Matignon's condition [4,7], it can be observed that |arg(λ i )| > q π 2 for all 0 < q < 1. Therefore, the coronavirus-free equilibrium point P 0 is locally asymptotically stable if condition (13) is satisfied. Under condition (14), we have a 0 < 0 and lim λ→+∞ C(λ) = +∞. Then a positive real root λ * > 0 of the characteristic Equation (16) exists; from Matignon's condition [4,7], it can be observed that |arg(λ * )| = 0 < q π 2 for all 0 < q < 1. Thus, the coronavirus-free equilibrium point P 0 is unstable. The coronavirus-free equilibrium point P 0 is locally asymptotically stable when the coronavirus-symptomatic-free equilibrium P I and the coronavirus endemic equilibrium P * do not exist.
Since we can deduce the stability of coronavirus symptomatic infected-free equilibrium P I from the stability of coronavirus endemic equilibrium P * by taking α = 1 and ξ = 0, the stability of the coronavirus endemic equilibrium P * is now discussed.

Theorem 4. Let
The coronavirus endemic equilibrium P * of the fractional-order system (1) is locally asymptotically stable if Proof. At the coronavirus endemic equilibrium P * , the Jacobian matrix (12) is given by The eigenvalues of the Jacobian matrix J(P * ) at the equilibrium P * are −d p < 0 and the roots of the characteristic polynomial arising from the remaining minor of order four: where c i , i = 0, . . . , 3 are given in (17). It is evident that c 3 > 0. From condition (18), the following inequalities on the coefficients are also satisfied.
From Theorem 4, we obtain the following result.

Remark 1. Letting
The coronavirus symptomatic-infected-free equilibrium P I of the system (1) is locally asymptotically stable if Proof. The result can easily be deduced from Theorem 4 by taking α = 1 and ξ = 0.

Global Stability
The global asymptotic stability of the equilibria P 0 , P I and P * of the fractional-order system (1) is now investigated by using a suitable constructed Lyapunov function, Lemma 4.6 in [24] and Lemma 3.1 in [25].
then the coronavirus-free equilibrium P 0 of the fractional-order system (1) is globally asymptotically stable in R 5 + .
Proof. First, the four equations in (1) are independent of R, therefore the last equation in (1) can be omitted without loss of generality. Hence, let us consider the following function: It is easily seen that the above function is non-negative and also V 0 = 0 if and only if S = S 0 , E = 0, I = 0 and A = 0. Further, calculating the q-order derivative of V 0 along the positive solutions of (1), we find: From condition (22), we can show that the coefficient of the term I + kA in the last equality are negative. Thus, we have D q V 0 ≤ 0 for all (S, E, I, A) ∈ R 4 + and D q V 0 = 0 if and only if (S, E, I, A) = (S 0 , 0, 0, 0). Thus, the only invariant set contained in R 4 + is {(S 0 , 0, 0, 0)}. Hence, by Lemma 4.6 of [24], it is enough to prove the convergence of the solutions (S, E, I, A) to (S 0 , 0, 0, 0). From the last equation of (1) we can easily show that R converges also to 0. Therefore, P 0 is globally asymptotically stable in R 5 + if R 0 < 1.
then the coronavirus-symptomatic-infected-free equilibrium P I of the fractional-order system (1) is globally asymptotically stable in R 5 Proof. Consider the function This function is positive and V 1 (S, E, I, A) = 0 if and only if (S, E, I, A) = (S I , E I , 0, A I ).
By calculating the q-order derivative of V 1 along the solution trajectories of system (1) and using Lemma 3.1 in [25], we obtain  (25), we obtain From the fourth equation of (1), we have Multiplying the Equation (25) by where F 1 (x, y, z, u) is a function to be determined later, yields Then, we obtain The function F 1 (x, y, z, u) is chosen such that the coefficient of E I is equal to zero. In this case, we obtain F 1 (x, y, z, u) = y u − 1.

Thus, inequality (27) is equivalent to
By the arithmetic mean-geometric mean inequality we have 3 − 1 x − xu y − y u ≤ 0 for all x ≥ 0, y ≥ 0 and u ≥ 0. Hence D q V 1 (S, E, I, A) ≤ 0, and D q V 1 (S, E, I, A) = 0 if and only if S = S I , I = 0 and u = y (i.e., E E I = A A I ). Since S must remain constant at S I , D q S is zero. This implies that A = A I and E = E I . Thus, by Lemma 4.6 in [24], it follows that the coronavirus-symptomatic-infected-free equilibrium P I is globally asymptotically stable in , (i.e., R 0 > 1), and either α = 1 or ξ = 0, (28) the fully endemic equilibrium P * of the fractional-order system (1) is globally asymptotically stable in R 5

Proof. Consider the function
This function is positive and V 2 (S, E, I, A) = 0 if and only if (S, E, I, A) = (S * , E * , I * , A * ). By calculating the q-order derivative of V 2 along the solution of system (1) and using Lemma 3.1 in [25], we obtain From the third and fourth equations of (1), we have Multiplying the Equations (30) and (31) respectively by where F 2 (x, y, z, u) and F 3 (x, y, z, u) are functions to be determined later, yields Then, we have The functions F 2 (x, y, z, u) and F 3 (x, y, z, u) are chosen such that the coefficients of ω p E * and I * are equal to zero. In this case, we obtain and F 3 (x, y, z, u) = y u − 1.
Thus, inequality (35) is equivalent to By the arithmetic mean-geometric mean inequality we have 3 − 1 ≤ 0 for all x ≥ 0, y ≥ 0, z ≥ 0 and u ≥ 0. Hence, D q V 1 (S, E, I, A) ≤ 0, and D q V 1 (S, E, I, A) = 0 if and only if S = S * and u = z = y (i.e., E E * = I I * = A A * ). Since S must remain constant at S * , D q S is zero. This implies that A = A * , I = I * and E = E * . Thus, by Lemma 4.6 in [24], it follows that the fully endemic equilibrium P * is globally asymptotically stable in R 5

Discussion
The proposed epidemic model for COVID-19 containing a fractional derivative formulation has been analyzed. The results show that the model is sound, because trajectories are confined to a compact set and the subpopulations corresponding to the various compartments into which the whole population is partitioned remain bounded. Unboundedness would indeed be biologically impossible. The model exhibits the endemic equilibrium P * at which the disease persists. It also contains the disease-free point with only susceptible individuals, P 0 , for which the pandemic is finally eradicated. This equilibrium should be the goal of the policies attempting epidemic diffusion confinement. There is also an "intermediate equilibrium", P I , which does not contain symptomatic infected, but it still harbors the disease among the asymptomatic individuals, so that it is in reality still endemic. The susceptible-only point P 0 is incompatible with both P I and P * , the latter being unfeasible when the former is stably attained. Note also that P I is a kind of particular case of coexistence P * , which arises if asymptomatics will never develop symptoms and all the exposed migrate into the asymptomatic class. In Section 3.5, global stability for each equilibrium in suitable conditions is shown. This result is relevant but also to be expected, in view of the incompatibility of multiple equilibria discussed above. Figures 2-4 show these three theoretical types of possible model behaviors. The comparison of these results with publicly available datasets on COVID-19 is not really our goal. In part, this is also due to the discussion that follows below, indicating that the population equilibrium levels change with the changes of the fractional derivative order.
Note that, of the two endemic points P * and P I , at equilibrium P I the disease remains hidden in the population. However, it still has its pernicious effects, since people will still "silently" die of it, at rate ν. In a sense, this equilibrium is worse than the explicit endemic one P * , because the latter keeps on showing to the people that the disease is not eradicated and therefore warns the population to maintain suitable safety measures. Finally, we observe that the stability results obtained analytically for the fractional operator model in this case do coincide with those of the model formulated using standard derivatives. Therefore, theoretically, the memory effects on this dynamical system are essentially scantly relevant, and the use of a classical standard formulation appears to be adequate. However, to further analyse in practice the effect of the order of the fractional derivatives, we have investigated how it may possibly affect the equilibrium levels. Simulations of the population values attained at equilibrium are also reported as a function of the fractional derivative order, in Figures 5-7, which correspond to the time series of the Figures 2-4 mentioned above. It is clearly seen that the order q of the fractional derivative does influence the final stage of the compartment levels. In Figure 5, for the endemic equilibrium P * , all compartments show a population decrease when q falls below 0.5, with the exception of S, which slightly increases. The same trend is also observed for the other endemic point P I , in Figure 7, although in this case the compartments are arranged differently in the two frames, to better show the details. For the disease-free point P 0 , in Figure 6, a low value of q decrements the susceptibles and increases the other compartments. These are very small, but in fact should vanish. However, the simulations have all been carried out with the same final time T = 10,000, so that the result indicates that a low q entails a longer time for disease eradication. (Right) infected and asymptomatic individuals are reduced to the single digits, essentially vanishing, recovered and exposed individuals remain confined respectively around the levels 100 and 50. In this situation, the system approaches the disease-free equilibrium P 0 . The parameter values are Λ = 500, β I = 0.0005, k = 0.1, d p = 0.082, α = 0.15, ω p = 0.1, ω p = 0.1, ξ = 0.1, µ = 0.001, γ p = 1.764, γ p = 0.6024, ν = 0.0005. A possible line of research to be pursued in the future is to compare the real-world scenario data with the simulation outcomes to assess the optimal fractional order to be used in the model, or if instead the standard derivatives are sufficient.  (Right) exposed, infected, asymptomatics and recovered. In this situation the system approaches the disease-free equilibrium P 0 . All the simulations have been stopped at the same final time, T = 10,000. For low values of the fractional derivative order, the populations E, I, A and R are not zero, but this feature is intentionally left in the figure. Thus, it denotes the different speeds at which the populations vanish as a function of the fractional derivative order. The parameter values are the same as those in Figure 3.