Caputo SIR model for COVID-19 under optimized fractional order

Everyone is talking about coronavirus from the last couple of months due to its exponential spread throughout the globe. Lives have become paralyzed, and as many as 180 countries have been so far affected with 928,287 (14 September 2020) deaths within a couple of months. Ironically, 29,185,779 are still active cases. Having seen such a drastic situation, a relatively simple epidemiological SIR model with Caputo derivative is suggested unlike more sophisticated models being proposed nowadays in the current literature. The major aim of the present research study is to look for possibilities and extents to which the SIR model fits the real data for the cases chosen from 1 April to 15 March 2020, Pakistan. To further analyze qualitative behavior of the Caputo SIR model, uniqueness conditions under the Banach contraction principle are discussed and stability analysis with basic reproduction number is investigated using Ulam–Hyers and its generalized version. The best parameters have been obtained via the nonlinear least-squares curve fitting technique. The infectious compartment of the Caputo SIR model fits the real data better than the classical version of the SIR model (Brauer et al. in Mathematical Models in Epidemiology 2019). Average absolute relative error under the Caputo operator is about 48% smaller than the one obtained in the classical case (ν=1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\nu =1$\end{document}). Time series and 3D contour plots offer social distancing to be the most effective measure to control the epidemic.


Introduction
Literature includes mathematical models for pathways for the transmission of infectious ailments. These models play an important role in quantifying and assessing the efficient control and preventive measures of infectious ailments [2][3][4]. It has been proved in multifarious ways that mathematical modeling is a very flexible and efficient way of researching the dynamics of transmission of infectious ailments. Mathematical analysis and numerical simulations can be used to create and evaluate control measures that are convincing. As for the model of compartmental ailments, there are several infectious ailments, beginning with the very classic SIR model to the more complex ones [1].
In December 2019, there was an emergence of a novel ailment in China which was later declared as pandemic named COVID-19 by the World Health Organization [5]. It is well known that the COVID-19 pandemic has given rise to fearful living and coexistence around the globe, with millions of people worldwide infected. The mortality-recovery ratio appears to be in a positive proportion. Nevertheless, due to the sensitivity of polymerase chain reaction, the absence or presence of the previously infected host is observed and the recovery rate appears to be promising in the absence of any curative vaccine. The challenge faced by health care professionals, the World Health Organization, and the Center for Disease Control in each quarter was whether reinfection could occur after a COVID-19 patient had been clinically treated. The subtle nature of the disease has led many scientists and medical practitioners to massively embark on multiple studies to combat and avoid the spread of the disease [6][7][8][9][10][11].
Out of the seven known human COVID-19, four are the known human influenza pathogens. SARS-CoV, MERS-CoV, and 2019-nCoV are responsible for severe respiratory diseases [12]. While COVID-19 has long been developed and studied by researchers and medical professionals, many people are lacking knowledge of the disease, and vaccines and antiviral drugs are still not available to directly prevent or treat the infection. The last major outbreak of SARS-CoV in China was in 2003. It was an acute respiratory infectious ailment with a large risk of demise. China managed the SARS outbreak through numerous inspections and effective preventive measures. The time of incubation for COVID-19 is substantial and very long compared to SARS. Many studies computed different periods of incubation for the COVID-19, as an example, 5.2 days [13], 3.0 days [14], and 4.75 days [15]. In [14], an incubation period of up to 24 days is reported, this extends to 38 days as reported in Enshi Tujia and Miao autonomous prefecture in Hubei province of China. Notice that people afflicted with asymptoms are very numerous [16], and in comparison to SARS-CoV and MERS-CoV the death rate is much inferior [14]. Genetic virus studies show that SARS-CoV and 2019-nCoV are 85% homologous [17], but 2019-nCoV binds ACE2 to an affinity higher than SARS-CoVS [18]. In the positive reported cases induced at COVID-19, SARS will be surpassed by the end of 29 January 2020. The variability of the asymptomatic cases, the incubation time and the super transmissibility of the virus bring considerable difficulties in controlling the epidemics. There is also very recent and novel research that discussed the dynamic and transmission of COVID-19 [19,20].
Over the last few decades, many scientists have shown that the fractional models can more accurately explain natural phenomena than the differential equations of the integer order. Because of this advantage, the fractional calculus has taken on the importance and popularity of modeling realistic cases, especially those with memory effects [21][22][23]. Furthermore, applications of the fractional calculus and mathematical modeling are found in many fields of social sciences, engineering, and mathematical biology [24][25][26][27][28][29][30][31][32]. Given this importance, we feel inspired to investigate and examine a new fractional version of the model involving Caputo operator. To the best of our knowledge, this is the first time Caputo operator has been employed for the model being considered.

Formulation of the Caputo SIR model
Although there have been various complex models proposed recently for transmission dynamics of COVID-19 as discussed in Sect. 1, we use a relatively simple and most popular model used in the field of mathematical epidemiology. The SIR model [1] is the stepping stone of all the models proposed after it. This is the reason to choose the SIR model to investigate dynamics of the COVID-19 pandemic while using the Caputo differential operator. The SIR model consists of three mutually distinct categories. First are those people who could potentially catch the disease S(t), second are those who currently have the disease and can infect others I(t) (infectives), and then R(t) stands for the removed, and this is the group of people who have already caught the disease and have now either recovered from the disease or have died with the disease. With all mathematical models, we have to make various assumptions to simplify the real world phenomena because things are just too complicated to express everything in a set of simple equations, so the first assumption that we make here is that the epidemic is sufficiently short, so it does not last that long so that we can assume that the total population remains constant. The second assumption in the model relates to the way in which the disease is transmitted, and we assume that the rate of increase in the infectives is proportional to the contact between the susceptibles and the infectives, and we assume that this occurs at a constant rate. Finally, our third assumption relates to the removal rate and this is the category R(t), so there is a constant rate which could be a death rate or a recovery rate. Having made the assumptions, a set of equations (1)-(3) that are going to govern the model can be written down: where N = S + I + R and D stands for the classical integer-order differential operator. We now have three differential equations for three categories of people within the population.
So, the number of susceptibles is going to decrease according to the number of contacts (contact rate = β) between the infectives and the susceptibles. Similarly, the number of infectives will increase due to contact between people and decrease because of people either recovering or dying as a result of the disease. Finally, the removed category includes people that no longer can catch the disease either because they have recovered or died, and this is going to increase at the constant rate (recovery rate = γ ) depending on how many infectives there are. Next, we need some initial data before we can solve the system of differential equations. Therefore, There are many instances wherein classical epidemiological models have been reinvestigated using operators of nonlocal nature. It is because of their memory property which makes them the most suitable tool to capture dynamics for the spread of a disease since epidemics are known to have memory retaining characteristics.
Having been inspired by a plethora of research works carried out in fractional mathematical epidemiology (see, for example, [33][34][35] and most of the references cited therein), we have introduced the following Caputo-type SIR model wherein its distinguished feature is that the dimensional inconsistency of the model has been removed by carrying the fractional order ν in the power of biological parameters β and γ : where C D ν 0,t stands for the Caputo (noninteger-order) differential operator with order of system (5) to be ν ∈ (0, 1).

Existence and uniqueness results
Among the vital areas in the concept of noninteger-order differential equations is the concept of existence and uniqueness of solutions in a dynamical system. Recently, the theory has attracted many researchers' attention [36]. By means of a fixed point theorem, we report the existence and uniqueness of (5). Take into account the proposed model (5) in the form as comes next: with Therefore, (5) can accommodate the form with T representing the transpose. Now, we can write (8) as Consider the Banach space on [0, b], with [0, b] denoting a set of continuous functions from on R with the associated norm given by ϒ = sup t∈J |ϒ| be E = C([0, b]; R). Also, take into account the following theorem.
which is an equivalence of (5) accommodates a unique solution only if L R < 1, and .
Proof Take into account P : E → E expressed as (11) shows that the unique solution for (5) represents the fixed point of P. Moreover, take sup t∈J R(t, 0) = M 1 and κ ≥ ϒ 0 + M 1 . Therefore, it will suffice to verify PH κ ⊂ H κ , and the set given by H κ = {ϒ ∈ E : ϒ ≤ κ} is convex and closed. Now, for any ϒ ∈ H κ , we have which affirms the result. Further, for given ϒ 1 , ϒ 2 ∈ E, one reaches indicating that (Pϒ 1 ) -(Pϒ 2 ) ≤ L R ϒ 1 -ϒ 2 . Thus, as a consequence of the Banach contraction rule, (5) possesses a unique solution on J . Now, we want show the existence of solutions of (5) by means of the Schauder fixed point principle.

Lemma 3.2 Take M = ∅ as a bounded, convex, and closed subset of a Banach space E.
Consider P 1 , P 2 as operators in obedience of the following: • P 1 is continuous and compact; • P 2 is a contraction mapping.
Then ∃u ∈ M so that u = P 1 u + P 2 u.

Then (5) possesses at least one solution on J only if
Take into account P 1 , P 2 on B ζ given by Therefore, ∀ϒ 1 , ϒ 2 ∈ B ζ , we have Thus, Now, the contraction of P 2 will be proved.
Given any t ∈ J and ϒ 1 , ϒ 2 ∈ B ζ , it yields Having R as a continuous function, P 1 is continuous. Furthermore, ∀t ∈ J and ϒ 1 ∈ B ζ , P 1 ϒ ≤ K < +∞, P 1 is uniformly bounded. At last, we prove that P 1 is compact. Start with sup (t,ϒ)∈J ×B ζ |R(κ, ϒ(κ))| = R * , which yields Hence, P 1 is equicontinuous and so is relatively compact on B ζ . In accordance with the Arzelá-Ascoli principle, P 1 is compact on B ζ . Having satisfied all the claims of Lemma 3.2,

Stability analysis
Here, the stability of (5) will scrutinize the aspect of Ulam-Hyers and generalized Ulam-Hyers [37,38]. It has been proved that stability analysis is important for an approximate solution. Assume ε > 0 and take into account the inequality and ε = max(ε j ) T , j = 1, . . . 5.
Definition 1 (8), which is the equivalence of (5), is Ulam-Hyers stable if ∃ X R > 0, so that ∀ ε > 0 and a solutionῩ ∈ E holds for (1), there is the unique solution ϒ ∈ E for (8), with Definition 2 (8), which is the equivalence of (5), is referred to as generalized Ulam-Hyers stable if ∃ a continuous function ϑ R : R + → R + , with ϑ R (0) = 0, so that ∀ solutionῩ ∈ E of (17), there is the unique solution ϒ ∈ E for (8), so that Remark 1 A functionῩ ∈ E satisfies (17) if and only if ∃ a function h ∈ E with the property as comes next: Lemma 4.1 Surmise thatῩ ∈ E holds for (17), thenῩ holds for the following: Proof By utilizing (ii) of 1,

,Ῡ(t) + h(t)
and Lemma 4.1, we have that Considering (i) of 1 yields Hence, we complete the result.

Estimation of parameters
A disease model like the one being investigated in the present research study is widely accepted and validated if its simulations for infectious compartment agree well with those of available real cases of the disease. In order to do so, one needs to obtain values of the parameters of the model in a best possible way so that they cause the model to be in good agreement with the available data. This brings into attention a few existing methods for getting best values of such parameters. The Bayesian technique, probability plotting, maximum likelihood estimation, and least-squares method are among the methods that exist for parameter estimation. In this study, we have used the least-squares method to compute the parameters (β = transmission rate and γ = recovery rate) of model (5).
Since there are only two non-demographical parameters for the proposed model, they are best fitted using the real cases of COVID-19 pandemic throughout Pakistan (source http://covid.gov.pk/stats/pakistan). Daily cases of the disease are taken from 1 April to 15 May 2020 when the present research study was being conceived and prepared.
In the least-squares approach, the objective function requires to be minimized. The objective is based upon adjusting the parameters for a model function to fit a data set in a best possible way. A simple data set consists of n points (data pairs) (x k , y k ), k = 1, . . . , n, where x k shows an independent variable and y k stands for a dependent variable whose value is computed using observation. The model function possesses the form g(x, p), where m adjustable parameters are kept in the vector p. The goal is to notice the parameters for the model that fit the data in a best possible way. Such a fit of a model to a data point is measured by its residual, defined as the difference between the real available value of the dependent variable and the value predicted by the model: The least-squares method obtains the optimal parameter values by minimizing the sum of squared residuals as shown below: For initial conditions, the overall population of Pakistan is found to be N(0) = 212.2M, the initial susceptible population is estimated to be S(0) = 212, 197, 711, the initial exposed population is estimated to be E(0) = 2, 289, and the initial recovered population is estimated to be R(0) = 0. There are only two biological parameters best fitted via the leastsquares fitting method thereby yielded best fit of the model's solution to the real cases chosen from Pakistan as depicted by Fig. 1 under the classical situation, that is, when ν = 1 wherein the best parameters are as follows: β = 3.0918 and γ = 3.0190 for the transmission and recovery rate, respectively. The average absolute relative error between the real cases and the model's simulations for the infectious compartment is decreased. Such a value for the error is approximately 3.3603e -02. Figure 1 shows the real cases by black circles, while the best fitted curve is shown by the blue solid line and the residuals are shown by side. The biological parameters included in the model are listed in Table 1 along with their best estimated values obtained via the least-squares technique. These parameters have finally produced the value of the basic reproduction number equal to R 0 = 1.0241  for the real COVID-19 cases under the classical situation in Pakistan from 1 April to 15 May 2020. Moreover, similar analysis for the Caputo fractional-order model is carried out wherein the best fitted parameters are found as β = 3.9405 and γ = 3.8010 for the transmission and recovery rate, respectively. One distinguished feature in the present research study is the computation of an optimum value of the fractional-order parameter ν which is 8.561591e-01. Hence, the basic reproduction number is obtained as R 0 = 1.0313 and the corresponding average absolute relative error is 1.7360e-02 for the Caputo model, whereas 3.3603e-02 is the average absolute relative error when ν = 1. This clearly reveals the superiority of the fractional operator over the classical one. Figure 2 shows the real cases by black circles, while the best fitted curve obtained under the Caputo fractional case is shown by the blue solid line and the residuals are shown by side.

Sensitivity analysis
In this portion, the concept of sensitivity analysis is employed to discover the robust significance of the generic parameters that are present in the basic reproduction number R 0 . Further, with the aid of parameter values from reliable assumptions, both analytic and numerical values of the parameters in R 0 are obtained. The analytic expressions obtained can be used to shed some light on how to control the onset of the model in variant localities if and only if the dynamics follow model (1). The threshold value R 0 is a quantity of which lowering the number to less than unity is considered as the major way of curtailing and aborting the spread of the ailment. The sensitivity index technique is used to measure the most sensitive parameters in the model, those with positive sign are considered as highly and proportionally sensitive for increasing the value of R 0 , while those with negative sign  are less sensitive for the decrease of R 0 and the other category are neutrally sensitive (with zero relative sensitivity). It is popularly known that the cause of raping transmission is associated directly with the basic reproduction number R 0 . The elasticity indices of R 0 to the associated parameters in the model are defined as follows [39]: where R 0 denotes the basic reproduction ratio and P i is as stated above. R 0 for model (5) under consideration is defined by the following expression: Thus, after some computation we reach The numerical values showing the relative significance of the R 0 parameters are given in Table 2. Positive relationship parameter is β, while that of negative relations is γ . A negative relationship indicates that an increase in this parameter's value would help to reduce the brutality of the disease. While a positive relationship indicates that an increase in the values of that parameter would have a substantial effect on the frequency of the spread of  Table 2.

Numerical simulations
This is the place where we get deep insights into the model's dynamical behavior. The present section offers the numerical simulations of Caputo model (5) while using the biological parameters as listed in Table 1. We have employed the FDE12 technique in order to get the solution of the nonlinear model shown in equations (5) and achieved the graphical results based upon parameters which are taken to be variables. Susceptible S(t), infectious I(t), and recovered R(t) populations are investigated with different values of the parameters. It may also be noted that the optimized value of the fractional order ν = 8.561591e-01 is used during each simulation. As shown in (a) plot of Fig. 5, the susceptible individuals decrease as the transmission rate increases; in contrast, the susceptible individuals decrease as the recovery rate decreases as shown by (b) plot of Fig. 5. Similarly, the slightly increased value of the transmission rate β is responsible for the spread of the virus (Fig. 6(a)), and the number of infectious individuals decreases when there is an improvement in their recovery which is observed in (b) plot of Fig. 6. Surprisingly, an increasing value of β causes an increase in the recovered population after a certain level of time interval as shown in (a) plot of Fig. 7, and similar surprising behavior is noted in (b) plot of Fig. 7. However, this should be obvious since with rapid transmission rate there will be more effective strategies used by public health sectors and concerned government to overcome the situation  Fig. 7.
We have also investigated the effects of fractional order ν on the dynamical behavior of Caputo model (5) for all three compartments. As can be seen in Fig. 8, the susceptible population decreases if ν increases, whereas Fig. 9 shows an increase in the infectious population for increasing values of ν, and this makes sense because the susceptible population, after being infected, is shifted to I compartment thereby leading to an increase therein. It is observed in Fig. 10 that the recovered population propels forward if ν gets large. Finally, we have also shown the dynamics of the basic reproduction number R 0 under different effects of the transmission rate β and the recovery rate γ in Fig. 11, wherein even a slightly larger value of β brings the reproduction number near to 1, which is clearly an alarming situation for policy makers to devise an effective approach to prevent R 0 to be greater than 1.

Conclusion
Using Caputo differential operator famous to have nonlocal nature, which is perfectly suitable to investigate transmission dynamics of a disease, we have fractionalized the SIR epidemic model having order ν without violating dimensional consistency between parameters and the operator itself. The proposed Caputo SIR model is used to comprehend the behavior of the devastating disease of COVID-19 that has recently emerged in the entire community of humankind. The model has shown promising results for it is proved to have a unique solution on the basis of the Banach contraction principle. Using Ulam-Hyers and its generalized version, the Caputo model is shown to be stable. Fitted parameters of the model, using real incidence cases of the virus from 1 April to 15 May 2020, are obtained under the least-squares approach wherein the fractional order ν is also optimized to be 8.561591e-01: one of the major contributions in the present work. Using these values, it has been proved that the Caputo model outperforms its classical version by 48% while the basic reproductive number is R 0 = 1.0313. Sensitivity of the parameters β and γ relating with R 0 is also investigated. Various numerical simulations carried out suggest that the epidemic can effectively be controlled only if its contact rate is lowered, and this is possible when people observe social distancing and wear protective masks. On the other hand, an underdeveloped country like Pakistan will be in huge trouble if such measures are not strictly followed upon. Our future study will explore the effects of nonsingular differential operators on the standard SIR model.