Fractional model for the spread of COVID-19 subject to government intervention and public perception

COVID-19 pandemic has impacted people all across the world. As a result, there has been a collective effort to monitor, predict, and control the spread of this disease. Among this effort is the development of mathematical models that could capture accurately the available data and simulate closely the futuristic scenarios. In this paper, a fractional-order memory-dependent model for simulating the spread of COVID-19 is proposed. In this model, the impact of governmental interventions and public perception are incorporated as part of the nonlinear time-varying transmission rate. In addition, an algorithm for approximating the optimal values of the fractional order and strength of governmental interventions is provided. This approach makes our model suitable for capturing the given data set and consequently reliable for future predictions. The model simulation is performed using the two-step generalized exponential time-differencing method and tested for data from Mainland China, Italy, Saudi Arabia and Brazil. The simulation results demonstrate that the fractional order model calibrates to the data better than its integer order counterpart. This observation is further endorsed by the calculated error metrics.


Introduction
The novel coronavirus pandemic, COVID-19, has impacted all aspects of our life and the disease continues to threaten global economy, public health, human interactions, and safety. Since the outbreak of the disease, there has been a collective effort to monitor, predict, and control its spread. For this purpose, different biological and non-biological techniques, tools, interventions, analysis, and experiments have been employed to understand the dynamics of disease spread.
Epidemiological mathematical models continue to play a pivotal role in understanding, capturing and predicting the spread of infectious diseases. In particular, first-order models such as SIR, SEIR, and many others, are intensively explored to mimic the collected data and make effective predictions. In the recent months, several such models have been developed to understand and predict the transmission of COVID-19 pandemic. See for example [1,2] .
Fractional-order models have been recognized as a powerful mathematical tool to study anomalous behaviors observed in many physical processes with prominent memory and hereditary properties. In general, such properties might not be satisfactorily captured by the integer-order counterpart models. This motivated researchers to consider fractional-order compartmental models for better understanding of dynamics of infectious diseases.
Since the outbreak of COVID-19 pandemic, a variety of fractional-order models for the spread of this disease are formulated [11][12][13][14] These models differ mainly in the considered compartments, mitigation effects, multi-fractional orders, and parameter forms.
Fractional-order epidemiological models are often obtained from their integer-order counterparts by replacing the firstorder derivative with a fractional one. However, this leads to a mismatch between the dimensions if the same parameters are used. This balance of dimensionality has been highlighted in some models but without much epidemiological interpretation of the new parameters. A common interpretation is that the parameters are observed in inhomogeneous time scale [15] corresponding to the fractional power of time unit. These parameters in the fractional model are commonly approximated by taking the fractional power of the classical parameters or by introducing an auxiliary parameter [16][17][18][19][20] .
In this article, we propose and simulate a system of fractional order differential equations for the dynamics of COVID-19 pandemic. The proposed model uses the nonlinear time varying transmission rate given in [1,21] which incorporates the effects of governmental interventions and public perception. In addition, an algorithm is developed so that for a given data set, the optimal fractional order of the model together with empirical values imitating the strength of governmental interventions are approximated. Consequently, reliable predictions for the future scenarios could be made.
The proposed model and algorithms are applied to the data from Mainland China, Saudi Arabia, Italy and Brazil. The simulation results suggest that the fractional order model has better agreement with the data than its integer order counterpart. This outcome is also supported by the error metrics considered in this paper. In addition, prediction trajectories for Saudi Arabia (a case with several changes in governmental interventions) and Brazil (a case with no governmental interventions) are included.
Since the focus is to provide adequate study and simulation of the fractional order model, we adopt the two-step generalized time differencing GETD2 scheme developed by Garrappa and Popolizio [22] to obtain the numerical results. However, to achieve a cost-effective implementation of this scheme, we use the partial fraction decomposition of the global Padé approximants of the two-parameter Mittag-Leffler function [23] .
The rest of this paper is arranged as follows. Next, the memory-dependent SEIR model is introduced. Section 3 is devoted to the discussion of reproduction number of the disease and local stability analysis. The computational algorithm is described in Section 4 . The simulation results are presented and discussed in Section 5 .

Memory-dependent COVID-19 Model
Following [1] , we assume the total population N at any given time to be distributed among the four compartments: susceptible, exposed, infectious, and recovered, labeled by S, E, I and R, respectively. That is, The assumed dynamics among these populations is depicted in Fig. 1 and can be describes as follows: • Susceptible individuals enter the exposed compartment at the rate βI/N, where β is a nonlinear transmission rate. • Exposed individuals progress to the infectious compartment at the rate σ, while the infectious individuals recover from the disease at a rate γ . • The population is assumed to have a natural death rate μ, which is given by μ = 1 ×365 , where denotes the lifeexpectancy.
• The birthrate of the population is assumed to be equal to the death rate.
Furthermore, we consider a variable P (t) to represent the public perception of the risk of the pandemic. This quantity increases with the number of infections at a rate φγ , and diminishes over time at a rate ω. The parameter φ is the proportion of severe cases.
The values σ −1 , γ −1 , and ω −1 can be interpreted as the mean latent period, mean infectious period, and mean duration of public reaction, respectively. These quantities together with μ, and β all have units of (time) −1 .
Based on the above assumptions, we propose the following model for the dynamics of COVID-19 pandemic: subject to the non-negative initial conditions The operator c D q denotes the Caputo fractional derivative of order q ∈ (0 , 1) defined by This derivative has been used successfully in the literature for compartmental models of different inf ectious diseases [3,4,6] .
In addition, the Caputo derivative allows to specify the initial data in a natural way.
In the special case q = 1 , the model in (1) reduces to a modification of COVID-19 model introduced in [1] . As observed in the right-hand side of (1) , the q th power of the parameters is used as suggested by Diethelm [3] and Angstmann et al. [18] . This is necessary to ensure the consistency in units since the left-hand side has the units of (time) −q .
For the transmission rate β in (1) , as in [1] , we use the formula In the integer order SEIR models, q = 1 , the parameters γ , μ, σ and β 0 describe the transfer rates between the different compartments and are measured per unit time, i.e. (time ) −1 . These parameters are the intuitive standard epidemiological rates such as, recovery, death, incubation, and transmission rates, which can be identified empirically or statistically. As for the fractional order model, q ∈ (0 , 1) , these transfer rates retain similar epidemiological meanings. However, they are to be measured on time scale (time ) −q which might not be feasible. Raising these parameters to the power q is an approach to approximate the corresponding transfer rates measured on time scale (time ) −q . An advantage of this approach is that the model reduces to the classical counterpart when q = 1 .

Basic reproduction number and equilibria stability
The basic reproduction number, R 0 , is the expected number of secondary infections that a single infectious individual will generate on average within a susceptible population. This number is arguably the most important quantity in determining the nature of an epidemic. In practice, a disease spreads whenever the reproduction number exceeds 1 ( R 0 > 1 ). On the other hand, an outbreak with a reproduction number below 1 will gradually disappear.
An equilibrium point is defined as follows.
Definition 3.1. [24] A point u * is an equilibrium point of the Caputo fractional dynamical system c D α Using this definition, we see that is an equilibrium point of (1) . This point is called a disease-free equilibrium. By applying the next-generation operator approach (see [25,26] ), the reproduction number R 0 is obtained as the spectral radius of the next generation matrix , which is the transmission rate β evaluated at u * . Thus, R 0 of (1) is given by Next, we analyze the local asymptotic stability of the disease-free equilibrium. For this purpose, we state the following result which follows by linearizing the system about the equilibrium point (as in Ahmed et al. [27] ) and then applying Theorem 2 of Matignon [28] .

Lemma 3.2.
Consider the following fractional-order system: then u * is locally asymptotically stable.
Based on Lemma 3.2 , we have the following stability result for system (1) .
Proof. The Jacobian of the right hand side of (1) at u * is given by The eigenvalues of the matrix J u * are given by and thus λ 6 is also negative. Consequently, if R 0 < 1 , then condition (3) is satisfied and Lemma 3.2 implies that u * is locally asymptotically stable.

Computational algorithm
Here, we describe briefly the numerical algorithm used to simulate model (1) . For this purpose, we consider the matrix form of the system, Table 1 Model parameters with possible range of values together with references.

days [1]
Proportion of severe cases, φ 0.2 [1] Mean duration of public reaction, ω −1 11.2 days [1,21] To solve this system numerically we use the GETD2 scheme developed in [22] . It has been shown that this scheme when applied to (4) , has a (1 + q ) −rate of convergence. In addition, it has elegant stability properties that makes it robust for stiff systems without need to impose severe restrictions on the step size.
where g j = g(v j ) . To initialize the scheme we take v 0 = u 0 and compute v 1 as For efficient implementation of the GETD2 scheme, we use a technique based on the partial fraction decomposition of the global Padé approximant of type (7,2) developed in [23] for the Mittag-Leffler function.

Model simulation
In this section, we perform numerical simulations of the fractional order COVID-19 model (1) using the GETD2 scheme (5) . We consider different case studies in which we calibrate the model to the data on confirmed COVID-19 infections in Mainland China, Italy, Saudi Arabia and Brazil. In addition, we present some predictions regarding possible spike and/or decay in the number of infections based on the intensity of behavioral response of population. For these purposes, we begin with a discussion about the parameters and their values.
The range of values of the parameters of the COVID-19 model (1) are listed in Table 1 . For our simulations, we use the following values.
On the other hand, as outlined in Algorithm 1 , the optimal values of the fractional order q and governmental intervention strength ν will be identified using the given data sets.
Based on the timeline of governmental interventions, the numerical values assumed by the parameter ν may vary on different time intervals. For the implementation of Algorithm 1 , feasible values of this parameter should reflect the intensity of the interventions. For example, periods with strict regulations should have no smaller values of ν than those with relaxed interventions. As such, at the initial phase of the pandemic and before any governmental intervention, we set ν to 0. On the other hand, based on the studies in [29,31] , substantial values of ν can be considered when accounting for the effect of mandatory use of face mask and physical distancing in public places.
The intensity κ was estimated in the study of Influenza pandemic by He et al. [21] as κ = 1117 . 3 . However, by studying and matching different data sets of COVID-19 infections, we observed that to better reflect the human behavioral response, piecewise adjustments to this value are to be considered to reflect the interventions timeline. In particular, we use κ = 0 at the early stage of the outbreak since it is assumed that adequate awareness and precautions are not yet in place. For approximating q and ν in Algorithm 1 , we use the value κ = 1117 . 3 for later time intervals. Nonetheless, when most of the governmental lockdown measures are being lifted, it is noted that the intensity of human behavioral response plays a major role in the decline of the number of infections. Thus, once the optimal values of ν and q are approximated, we fine-tune κ in order to achieve better matching.
Algorithm 1 Approximating the optimal q and ν pair Step 1 Create a list { q j } l j=1 ∈ [ a, 1] , a > 0 , for values of q .
Step 3 On each sub-interval S i , create a list V i of feasible values of ν for time-interval S i .
Step 4 Construct a collection C of all ordered combinations c Step 5 Define an empty list M . For q j , j = 1 , 2 , . . . , l Define an empty list E.
For c ∈ C Solve system (4) Use an error metric (say mean square error) to quantify the error betweenthe computed infection numbers and the reported ones. Store the triple ( q j , c, error) in E. Select from E the triple with minimum error value and store it in M . Select from M the triple with minimum error value and return the corresponding c and q j .
In all our simulations, the values of the RMSE between the approximated infected population and the actual data is provided for optimal q and q = 1 . Moreover, to give a perspective of the corresponding relative error, we also present the values of the mean absolute percentage error (MAPE).

Case study: outbreak in Mainland China
In Table 2 , we consider the reported data of cumulative COVID-19 infections in Mainland China from January 19, 2020 to March 6, 2020, as provided in [32] . To simulate the outbreak of the pandemic, we use the β 0 = 0 . 52 , μ = 1 / (76 . 91 × 365) , together with those given in (7) . Furthermore, the initial populations are As per the lockdown records of Chinese Administrative Division, we consider the governmental intervention in two phases. The First phase corresponds to the period Jan 20 -Feb 02 and the second one to the period Feb 03 -March 06.
Applying Algorithm 1 , we obtain q = 0 . 93 as an approximation of the optimal fractional order. The corresponding values of governmental intervention strength ν are given in Table 3 . For comparison purposes, we also provide the optimal values of ν for the integer order problem, q = 1 . Furthermore, the values of RMSE and MAPE in approximating the infected population are also appended. The error metrics in Table 3 suggest that the fractional order model performs better in terms of the RMSE. However, the MAPE seems to support the integer order model. This blended outcome can be understood in view of Fig. 2 . As can be observed, the absolute errors for the fractional order model are smaller throughout most of the period and thus smaller RMSE. On the other hand, although the relative errors are smaller in general for the fractional model, they attained some large values towards the end of the period, which in turn impacted the MAPE.   5 provides an insight on the effect of governmental interventions. As observed, the change in governmental intervention strength ν on Feb 3rd yields an instantaneous change in the rate of change of exposed population. This phenomenon is expected, because governmental interventions such as inter-city lockdown, curfew, school closure among others, are expected to mitigate the contacts between susceptible and infectious individuals. However, some increase in the number of infections is still observed. This accounts for the individuals who were incubating the virus before implementation of the governmental interventions.
Our numerical experiments also demonstrate the robustness of the fractional-order model as a generalization of the integer-order model. It is observed that the fractional order model gives extra degree of freedom, via the order of differentiation q, for data fitting.

Case study: outbreak in Saudi Arabia
Since the outbreak of COVID-19 in Saudi Arabia, the Kingdom has witnessed several phases of governmental interventions aimed at mitigating spread of the diseases. Due to seasonality, special breaks and festivals, there have been a number of relaxations and tightening of these interventions. To simulate the outbreak of the Pandemic within the Kingdom, we only adopt those interventions that were implemented at a national scale, as provided in https://www.garda.com/crisis24/news-alerts/334921/saudi-arabia-curfewhours-modified-for-ramadan-april-21-update-29 . We consider the parameters in (7) together with β 0 = 0 . 54 , μ = 1 / (75 . 13 × 365) .

The initial conditions are
Based on the data provided in https://ourworldindata.org/coronavirus-source-data , our simulation is performed for the period of March 03 to July 17 excluding the days 4th, 5th and 7th of March, where there was no reported data.
Using Algorithm 1 , we obtained an approximation of the optimal q and ν pair, with q = 0 . 87 . For comparison purposes, we present in Table 4 the ν values corresponding to both q = 0 . 87 and q = 1 . In addition, the values of the RMSE and MAPE in approximating the infected population for these values of q are presented. Brief remarks on the governmental actions and the values of intensity κ are also provided.
We observe that the values of ν presented in Table 4 reflect the intensity of the control measures in the sense that periods with more restrictive interventions have higher values of ν. From an epidemiology perspective, these interventions such as imposing a curfew, restricting the gatherings size, among others, reduce interactions and consequently mitigate the spread of the disease. As a result, the fitted values of ν tend to reflect the epidemiological reality since bigger values of ν will imply smaller values of the transmission rate β in (2) . Fig. 6 includes plots of the daily infections obtained from the fractional model (1) and plots of the parameters ν, κ and β. It is observed that the polynomial behavior of the data at the early stage of the outbreak is better captured by the power-law corresponding to q = 0 . 87 than by the exponential or almost exponential law, q = 1 or q = 0 . 98 . However, the final stage of simulation is better captured with q = 0 . 98 . This Further illustrates the effectiveness of fractional order models  for making predictions. The outcome of our simulation for the cumulative infections is presented in Fig. 7 . We note that q = 0 . 87 yields a satisfactory matching. Using the calibrated model (1) , predictions for the spread of the disease subject to different values of κ are shown in Fig. 8 . We assume that the governmental intervention continues to be the same as recent period. It is shown that, a drop in the intensity κ below the value between 21-June and 17-July will yield a spike in the infection numbers. On the other hand, maintaining same intensity will help obtain a slow decline of the infection numbers. However, a notable decay in the number of infections could be achieved if there is a significant increase in the intensity of the behavioral response. This corresponds to improved compliance with proper hand-washing, face masking, physical distancing, among others.

Case study: outbreak in Italy
In the early outbreak of COVID-19 in Italy, the country was hit badly and a number of governmental interventions were initiated. To simulate this outbreak using (1) , we consider   together with those provided in (7) . The initial conditions are The timeline of governmental interventions considered for this simulation are based on the report by Briscese et al. [33] , up to 3-April, together with other reopening information available on https://www.bbc.com/news/world-europe-52435273 . Our simulation covers the period of 31-January 2020 to 17-July 2020, based on the data provided in https://ourworldindata.org/ coronavirus-source-data . However for the period February 01-21, 2020, no new infections were reported. Consequently, we excluded these dates from our simulation.
Algorithm 1 applied to the fractional order model (1) yields an approximation for the optimal q and ν pair, with q = 0 . 88 . The corresponding values of ν are provided in Table 5 . In addition, brief descriptions of the governmental actions together with the adopted values of κ are also provided. Similar to Saudi Arabia case, we observe that the fitted values of ν reflect the intensity of governmental actions and are consistent from an epidemiological point of view.
Plots of the daily and cumulative infections obtained from the model (1) are presented in Figs. 9 and 10 , respectively. We observe from these plots that the fractional order model has a better matching. This observation is also supported by the values of the error metrics included in Table 5 .

Case study: outbreak in Brazil
Information from various sources indicate that there were no governmental control measures to mitigate the COVID-19 outbreak in Brazil. Some of the events in the early days of the outbreak can be found in an article by Reuters https://www.reuters.com/article/us-health-coronavirus-brazil-timeline/timeline-key-moments-inbolsonaros-handling-of-covid-19-crisis-idINKBN2482UL . Thus, for simulating the outbreak we set ν = 0 . The simulation covers the period from 26th February to 17th July using the data provided in https://ourworldindata.org/ coronavirus-source-data . We exclude the days February 27-29 and March 02-04, 08 and 10, since no infections were reported. The other parameters used in simulating (1) for the outbreak in Brazil are β 0 = 0 . 75 , μ = 1 / (75 . 88 × 365) , together with those provided in (7) . The initial conditions are With ν = 0 for the whole simulation time, implementation of Algorithm 1 reduces to approximating only the optimal fractional order q, which is approximated as q = 0 . 70 .
Our examination of the infection data in Fig. 11 shows that there is a significant change in the reported numbers around day sixty. Under no governmental interventions, the parameter in our model that can be used to account for this dramatic change is the human behavioral response κ. For this purpose, we split our time interval into [0, 60] and [61,134]. Then, we fit appropriate values of κ for both q = 0 . 70 (the optimal value of q ) and q = 1 . The calculated values are shown in Table 6 together with the RMSE and MAPE in approximating the daily number of infections.
Plots of the daily infections obtained from (1) are presented in Fig. 12 . As can be observed, the integer order model could not provide a reasonable matching of the data. On the other hand, the fractional order model with q = 0 . 70 is able to capture well the initial spread of the disease. Plots of the cumulative infections are also presented in 13 .
Overall, we remark that modeling the COVID-19 outbreak in Brazil is not a trivial task due to the erratic nature of the reported data. However, using the extra degree of freedom offered by the fractional order model through the fractional order q, we are able to find a possible fit for daily infection data.  Table 6 for q = 0 . 70 . Then we predict possible  infection numbers for the period from July 18 to August 16. The plots suggest that for a significant reduction in infections, a substantial increase in κ is required. This corresponds to improved compliance with proper hand-washing, face masking, physical distancing, among other public health measures.

Discussion
The proposed model aims for improved fitting of the real data and thus more realistic predictions. It provides extra degrees of freedom through the fractional order q and the transmission rate β(t) that includes the government intervention and intensity of human behavioral response.   From the outcomes of our numerical simulations, we note the following observations which can help in policy making for public health interventions.
• The data matching reveals that strict governmental lockdown measures result in large values of the parameter ν, and thus rapid decrease in the transmission rate. This is notably observed in the case of China.
• The intensity of human behavioral response, such as adherence to masking, hand washing, physical distancing, among others, compliments the governmental interventions in mitigating the transmission of the disease. This is demonstrated by the predictions in Fig. 8 , where we see that for fixed governmental interventions, the number of infections is inversely proportional to the behavioral response κ.
• Complete absence of governmental interventions results in growing transmission of the disease even with strong human behavioral response. This scenario is demonstrated in the Brazil model, where we set ν to zero since there is no record of governmental interventions. As can be observed, the human behavioral response alone is not adequate to mitigate the spread of the disease. Rather, it complements the governmental interventions. • As a consequence of the above observations, gradual lifting of governmental measures is advisable. For example, interventions such as mandatory physical distancing and masking in public places should be in place. These less strict interventions can still be well complemented by the behavioral responses and consequently help flatten the curve. The prediction for Saudi Arabia in Fig. 8 is a typical reference.

Conclusion
In this paper, we proposed a nonlinear fractional order model to study the dynamics of COVID-19 pandemic. Using the data of the outbreak in Mainland China, Saudi Arabia, Italy, and Brazil, the effectiveness of the proposed model in simulating the spread of the disease is demonstrated. Predictions based on varying intensity of human behavioral responses are provided. By examining the mean square error of the predicted number of infections and the reported data, we found that the proposed fractional model yields better matching than the integer order model. From the perspective of relative error, the values of mean absolute percentage error confirm this conclusion for the simulation of outbreak in Saudi Arabia, Italy and Brazil. On the other hand, for China, these two error metrics have mixed outcomes.