A fractional-order compartmental model for the spread of the COVID-19 pandemic.

We propose a time-fractional compartmental model (SEI A I S HRD) comprising of the susceptible, exposed, infected (asymptomatic and symptomatic), hospitalized, recovered and dead population for the COVID-19 pandemic. We study the properties and dynamics of the proposed model. The conditions under which the disease-free and endemic equilibrium points are asymptotically stable are discussed. Furthermore, we study the sensitivity of the parameters and use the data from Tennessee state (as a case study) to discuss identifiability of the parameters of the model. The non-negative parameters in the model are obtained by solving inverse problems with empirical data from California, Florida, Georgia, Maryland, Tennessee, Texas, Washington and Wisconsin. The basic reproduction number is seen to be slightly above the critical value of one suggesting that stricter measures such as the use of face-masks, social distancing, contact tracing, and even longer stay-at-home orders need to be enforced in order to mitigate the spread of the virus. As stay-at-home orders are rescinded in some of these states, we see that the number of cases began to increase almost immediately and may continue to rise until the end of the year 2020 unless stricter measures are taken.

Introduction

Fractional differential equations (FDEs) are used to model complex phenomena such as such as the modeling of memory-dependent phenomena (Di Giuseppe et al. [1], Baleanu et al. [2], Podlubny [3]), mechanical properties of materials (Caputo and Mainardi [4]), anomalous diffusion in porous media (Fomin et al. [5], Metzler and Klafter [6]), groundwater flow problems (Cloot and Botha [7], Iaffaldano et al. [8]), and control theory (Podlubny [9]), among others. They serve as a generalization of the integer-order differential equations and give more degree of freedom for modeling of biological and physical processes. FDEs have been applied in biological tissues [10], DNA sequencing [11], Pine Wilt disease [12], lung tissue mechanics and models [13], harmonic oscillators [14], Dengue fever [15], measles [16], human liver [17], diffusion processes [18], SEIR models [19]. Infectious disease outbreaks are one of the main causes of deaths in human. Their dynamics and spread are modeled and studied before the introduction of vaccines. COVID-19 began in December 2019 in China. Over the last few months, it has spread rapidly leading to over 400,000 deaths across the globe. The first occurrence in the United States was seen around mid January in Washington [20] and has spread across America with over 100,000 deaths and 1.5 million infected. The pandemic has disrupted the day-to-day activities of the human life with over six million jobs lost in the United States. Several actions and measures have been taken by the federal, state and local governments to mitigate the spread of the pandemic. The most prominent measures taken include social distancing, testing, use of face-masks and contact tracing. It is important to model this pandemic in order to better understand the spread and dynamics as well as address the challenges of the pandemic. In short, mathematical models are important to guide the decisions of health and government officials.

The goal of this study is to examine and analyze the spread of the pandemic using a modification of the Susceptible-Exposed-Infected-Recovered (SEIR) model with a time-fractional derivative. The use of fractional derivatives in the model stems from the fact that the spread of infectious diseases depends not only on the current state but also on its past states (history or memory dependency). Additionally, time-fractional order models reduce errors resulting from neglect of parameters in models. We shall focus on some selected states in the US. We note that models that consider the US as a whole may be misleading and have limited applicability as different states have different economical and political perspectives which determines the different control strategies used for each state. For example, while some states such as Maryland, New Jersey, New York, Connecticut, among others, enforced the use of masks in public places and longer stay-at-home order [21], other states do not enforce these measures thereby allowing for a possibility of increase of infected individuals in such states. There have been several models for the study of the pandemic. Lu et al. [22] considered a fractional-order SEIHDR model which incorporates intercity movements. Liu et al. [23] studied the dynamics of the pandemic by considering asymptomatic and symptomatic infected populations separately. Giordano et al. [24] studied the COVID-19 epidemic with intervention strategies in Italy. They proposed a model consisting of different stages: susceptible, infected, diagnosed, ailing, recognized, threatened, healed and extinct. They further discuss the long time behavior of the populations in which the susceptible, healed and extinct population remains. Stella et al. [25] studied the role of asymptomatic individuals via complex networks. In particular, they formulated a model which aims at studying the interactions in the population by means of complex networks. They further extended the model to a structured nonhomogenous version by means of the Watts-Strogatz complex network. Wu et al. [26] studied domestic and international spread of the pandemic by using different data sets. Zhao and Chen [27] discussed the dynamics of the pandemic by considering the Susceptible, unquarantined infected, quarantined infected and Confirmed infected (SUQC) model and parametrize the intervention effect of control measures. Zhang et al. [28] considered a fractional SEIR model with different order of the time-fractional derivative for each of the different population being studied. Tuan et al. [29] proposed a fractional-order model using the Caputo derivative for studying the transmission of COVID-19. They discussed the existence and uniqueness of solutions to the proposed model. They further used the generalized Adams–Bashforth–Moulton method for simulating the model. Bahloul et al. [30] proposed a fractional-order Susceptible-Exposed-Infected-Quarantined- Recovered-Death-Insusceptible (SEIQRDP) model for predicting the spread of COVID-19.

The remainder of the paper is organized as follows: In Section 2, we give some preliminary results and definitions from fractional Calculus. Furthermore, we discuss the formulation of the model by considering seven compartments: Susceptible-Exposed-Asymptomatic infected-Symptomatic infected, Hospitalized, Recovered and Dead populations (SEIIHRD). Section 3 details the properties and theoretical analysis of the model. Section 4 discusses the parameter sensitivity and identifiability analysis. In Section 5, we applied the model to observed data for some selected states in the US. In particular, this section details solving several inverse problems for parameter estimation and computation of the basic reproduction number. Finally, we give some concluding remarks in Section 6.

Model setup

Preliminaries

In this section, we present some preliminary results and definitions in fractional calculus.

Podlubny [3] The gamma function is defined by the integralwhich converges in the right half of the complex plane .

Caputo [31] For any the Caputo-fractional derivative of order of a function is defined as

Podlubny [3] The Mittag–Leffler function which generalizes the exponential function for fractional calculus is defined as

More generally, the two parameter Mittag–Leffler function is defined asIt has the following properties:

.

.

Li et al. [32] A point is said to be an equilibrium point of the system if and only if .

Lu et al. [22] An equilibrium point of the system is said to be asymptotically stable if all the eigenvalues of the Jacobian matrix evaluated at the equilibrium point, satisfies where are the eigenvalues of .

Model formulation

The model discussed in this work is a modification of the SEIR model having three additional compartments. We consider a SEIIHRD compartmental model which comprises of the susceptible, exposed, infected (asymptomatic and symptomatic), hospitalized, recovered and dead population. We assume that the natural death and birth rates are the same. We further assume that deaths in the S and R compartments are due to natural deaths and deaths in the other compartments are as a result of the pandemic. The schema below shows the transmission flow of the model (Fig. 1 ). Thus, the model consists of the following system of ODEsThe system has the associated initial dataThe total population is N which is further divided into and is the natural birth rate, and are the transmission rates to the susceptible population from the asymptomatic and symptomatic populations, respectively. is the incubation period of an exposed individual, denotes the fraction of the exposed population that becomes asymptomatic after the incubation period and the remaining of the population are symptomatic. and are the infectious rates of an asymptomatic and a symptomatic individual, respectively. is the recovery rate through hospitalization and, and are the mortality rates of the exposed, asymptomatic, symptomatic and hospitalized populations, respectively. We note that the death population in the model comprises of deaths during exposure infectious period ( and ) and hospitalization and assume that deaths due to other natural occurrences are negligible for this populations.

Fig. 1

Schematic diagram of the proposed SEIIHRD model.

We note that the parameters of the model are non-negative and have dimensions given by 1/time. This observation was originally noted in Diethelm [15]. To alleviate this difference in dimensions, we replace the parameters (except ) with a power of new parameters to obtain the new system of equations

Model analysis

In this section, we discuss the properties of the model beginning with the existence, uniqueness, non-negativity and boundedness of solutions of the model (1). For simplicity in analysis, we reduce the system (1) tosince is a linear combination of populations in some of the other compartments.

There exist a unique solution to the system (2) and the solution is non-negative and bounded for any given initial data .

See Appendix A. □

Stability analysis

Computation of the basic reproduction number

We shall use the next generation matrix originally proposed by Diekmann et al. [33] and further elaborated on by van den Driesche and Watmough [34] and Diekmann et al. [35] to determine . The disease-free equilibrium point (DFE) is . Now, consider the four compartments containing the infected individuals and let be the DFE point. Since the DFE point exists and is stable (shown in the next section) in the absence of any disease, then the linearized equation at the DFE iswhere is the rate of appearance of new infections in compartment and is the rate of transfer of infections to and from compartment . We further defineThen is the basic reproduction number where is the spectral radius of and is the next generation matrix. For the system (2),and the basic reproduction number is given as

The fractional system (2) has at most two equilibrium points

a disease-free equilibrium point DFE =

an endemic equilibrium point EE =

where

The DFE point is locally asymptotically stable if and where

See Appendix B. □

The EE point is locally asymptotically stable if and where

See Appendix C. □

Parameter sensitivity and identifiability analysis

We discuss the sensitivity and identifiability of the parameters with respect to the proposed model.

Sensitivity analysis

The sensitivity analysis (SA) deals with the significance or importance of the parameters in the model. In particular, it finds the most influential parameters that drives the dynamics of the model. It also describes the extent to which parameter changes affects the result of the methods or models with the goal of identifying the best set of parameters that describes the process or phenomena in question.

There are several SA methods which are broadly classified as local and global methods. In this work, we shall focus on the Morris screening method (local method) and Sobol analysis method (global method).

Morris screening method

The Morris screening method is a local sensitivity measure that makes use of the first order derivative of an output function with respect to the input parameter . It measures the effect of the output when the input variable is perturbed one at a time around a nominal value. It serves as a first check, in most analysis, in screening parameters for identifiability. The method evaluates elementary effects [36], [37], [38] with the th parameter through the forward perturbationMorris [39] proposed two sensitivity measures, the mean () and the standard deviation () of the elementary effects. For non-monotonic models, may lead to a very small value due to cancellation effects. For this reason, Campolongo et al. [40] proposed the use of absolute values for evaluating the mean. In order to obtain a dimension-free sensitivity, we prefer the use of the sensitivity measure given in Brun et al. [41] aswhere N is the number of sample points andA common practice in the literature [38], [42], [43] is to plot the indices against the standard deviation. We see from Fig. 2 (a) that the fractional order has the highest influence on the model output over time. The transmission rates for the symptomatic and asymptomatic populations and, the fractional parameter are the next most influential parameters in the model. This is further corroborated by Fig. 2(b). The four parameters ( and ) denoted by red squares are the most important parameters, that is they have the largest and values. The parameters ( and ) represented by the blue squares have the least influence on the model output and can be considered unimportant. The other parameters represented by the green squares have more influence than the parameters represented by the blue squares.

Fig. 2

Morris screening test.

One major setback of the Morris screening test for sensitivity analysis is the consideration of each parameter individually and independently of the other parameters. In real applications, this is not true as parameters have collinearity and dependencies on one another.

Sobol analysis

The Sobol method is a variance-based sensitivity analysis method which unlike the Morris screening method takes into account the effect of the relationship between each parameters of the model. It uses the decomposition of variance to calculate Sobol’s sensitivity indices: first and total order sensitivity measures. The basic idea of the Sobol’s method is the decomposition of the model output function into summands of increasing dimensionality, that iswhere is the partial variance of the contribution of the parameter and is the partial variances caused by the interaction of the parameters for .

The first order sensitivity index measures the main effect of parameter on the model output; that is the partial contribution of to the variance .The index [44], [45] is defined asThe larger this index, the more sensitive the parameter is to the model output [44], [45]. Using the law of total variances [38], [45], the index can also be expressed asandwhere is the partial variance caused by and is the mean of the model output calculated by using all the values of the other parameters (except ) and is the total variance.

The total sensitivity indices [46] measures the effects of parameter and the interaction with the other parameters. It is defined asThe total variance, for this index is given asandThe mean and the variance can be evaluated using quasirandom sampling method [38], [47] and are given asandwhere and are two independent parameter sample matrices of dimensions . We shall use the python SALib package [48] to compute the first and total order variance indices. Fig. 3 shows that the fractional order has the highest interaction with the other parameters. These results are consistent with the results in the Morris screening test as important parameters of the test show high interaction with the other parameters.

Fig. 3

Sobol sensitivity indices.

Parameter identifiability

The concept of identifiability is dependent on sensitivity. It entails the selection of the subset of parameters of a model having little or no collinearity and uncertainty, and which can be identified uniquely from a given set of observed data or measurements. In other words, it answers the question ”Can the available data be described by the model and the selected subset of parameters?”. There are several techniques or tests for parameter identifiability. Most of the tests are based on the Fisher information matrix (FIM) where for a model output function . Cobelli and Di Stefano [49] showed that a sufficient condition for identifiability is the non-singularity of FIM. Burth et al. [50] proposed an iterative estimation process which implements a reduced-order estimation by finding parameters whose axis lie closest to the direction of FIM. The associated parameter values are then fixed at prior estimates during the iterated process. Brun et al. [41] studied parameter identifiability using two indices; a parameter importance ranking index and a collinearity index which depends on the smallest eigenvalues of submatrices of corresponding to the parameter subset . Cintrón-Arias et al. [51] explained the need for a good parameter subset for identifiability to satisfy the full rank test. They further introduced two indices; the selection score and the condition number of . The smaller these indices the lesser the collinearity and uncertainty in the parameter values of the subset. Finally, they used the coefficient of variation index to examine the effect of parameters in the parameter subset. In this work, we shall use the test proposed by Cintrón-Arias et al. in identifying the parameters. The algorithm can be summarized in the following steps

Algorithm for parameter subset selection [51].

To discuss the results in this section, we shall use the state of Tennessee as a case study to understand parameter identifiability. Furthermore, we used the following values obtained using a random search algorithm as the nominal parameter set for the model: and the nominal error variance . We further divide the parameters into three groups according to their importance rankings discussed in the previous section:where and are the most and least influential parameter sets, respectively, while contains more influential parameters than . We display some selections of the parameter subsets of size in Table 1 where we have chosen the subsets with the smallest score values. The entries in Table 1 are ordered with respect to the selection score for each subset of same cardinality. A high selection score and condition number for a parameter subset indicates substantial collinearity and linear dependence, and thus is poorly identifiable even if the parameter subsets contains that is contains the set of most influential parameters. We observe that most of the selections in Table 1 contains at least one element in each of the groups listed above. This shows that while parameter importance ranking is crucial in identifying parameters that drives the dynamics of a model, it does not have substantial effect in identifiability. Identifiability depends on proper selection of subsets including parameters in each of the three groups above that describes the measurement or data.

Table 1

Selection scores and condition numbers for some selected parameter subsets.

p	Parameter Subsets	κ(χ(θ))	ϑ(θ)
13	(b_1,β_A,β_S,σ,γ_1,γ_2,ρ,μ_1,μ_2,μ_3,μ_4,η,α)	9.530e + 03	8.481e + 01
12	(b_1,β_A,β_S,σ,γ_1,γ_2,ρ,μ_1,μ_3,μ_4,η,α)	9.513e + 03	6.147e + 00
11	(b_1,β_S,σ,γ_1,γ_2,ρ,μ_1,μ_3,μ_4,η,α)	4.814e + 03	2.479e−04
10	(b_1,β_S,σ,γ_1,γ_2,ρ,μ_1,μ_4,η,α)	3.696e + 03	3.852e−07
	(b_1,β_S,σ,γ_1,γ_2,ρ,μ_3,μ_4,η,α)	3.685e + 03	4.143e−06
9	(β_S,σ,γ_1,γ_2,ρ,μ_1,μ_4,η,α)	3.411e + 03	2.766e−07
	(b_1,β_S,σ,γ_1,γ_2,ρ,μ_4,η,α)	3.311e + 03	2.824e−07
	(β_S,γ_1,γ_2,ρ,μ_4,η,α)	2.718e + 03	5.344e−08
7	(β_S,γ_1,γ_2,ρ,μ_1,η,α)	2.947e + 03	5.388e−08
5	(γ_1,γ_2,ρ,μ_4,α)	6.171e + 01	2.006e−09
	(β_S,γ_1,ρ,μ_4,α)	6.389e + 01	2.006e−09
	(γ_1,ρ,μ_4,η,α)	6.989e + 01	2.010e−09
4	(γ_1,ρ,η,α)	6.258e + 01	3.458e−10
	(γ_1,γ_2,ρ,α)	5.289e + 01	3.491e−10
	(β_S,γ_1,ρ,α)	5.349e + 01	3.496e−10
3	(γ_2,ρ,α)	5.140e + 01	5.431e−11
	(β_S,ρ,α)	5.203e + 01	5.850e−11
	(ρ,η,α)	6.224e + 01	9.109e−11

To have an idea of the variations of the condition number and the selection score, we give a plot of these values for in Fig. 4 (with logarithmic scales). Good parameter combination in Fig. 4 corresponds to values in the lower left corner of the figure where the values, and are relatively small.

Fig. 4

The condition number against the parameter selection scores of the sensitivity matrices for all parameter subsets with . Logarithmic scales are used on both axis.

To further analyze the parameter identifiability of the model, we consider the parameters subsets:such that . The choice of these parameter subsets are due to their relative small condition numbers and selection scores. In other to create synthetic data, we assume the nominal parameter subsets and error variance (given at the beginning of this section) to be the true parameter vectors and true variance. Furthermore, we add random noise to the model output as follows:We solve seven inverse problems for each of the parameter subsets . We analyze the result using the coefficient of variation and standard error [51] given asandwhere and .

From Table 6 , it is seen that the parameters have standard errors that is approximately 20 times their estimates. This shows substantial uncertainty in these parameter values and any parameter subset containing these parameters may result in illogical parameter estimation from observations. The standard errors of in show improvements and implies lower linear dependence and collinearity for parameters in than those in . Thus, a substantial improvement in uncertainty quantification is seen from to . Further improvements are observed for each of the other parameter subsets as more parameters are removed. For instance, with the removal of and in it seen that the standard error for and dropped from 4% and 1% to approximately 0.07% and 0.003%, respectively, of their estimates. Other improvements in include where the standard error is reduced to at least one-tenth of its standard error in . We note that there is no substantial gain in the removal of from and in as seen in Table 6.

Table 6

Parameter estimates for solving seven inverse problems from a synthetic data generated using the given nominal parameters and variance. For each parameter subset, we display the estimate the standard error, E and the coefficient of variation, .

	b_1	β_A	β_S	σ	γ_1	γ_2	ρ	μ_1	μ_2	μ_3	μ_4	η	α
θ˜	1.13e−4	2.60e−4	1.23e + 0	9.96e−1	2.69e−2	1.00e + 0	6.18e−3	3.32e−3	1.59e−5	3.05e−4	2.13e−5	1.57e−1	1.00e + 0
E	3.16e−6	1.42e−4	5.74e−3	2.29e−3	7.39e−4	4.56e−3	1.01e−5	4.94e−3	7.37e−4	6.48e−3	5.52e−5	9.47e−4	9.35e−4
V	2.81e−2	5.49e−1	4.66e−3	2.30e−3	2.75e−2	4.56e−3	1.64e−3	1.49e + 0	4.64e + 1	2.12e + 1	2.59e + 0	6.04e−3	9.35e−4
θ˜	3.52e−4		1.23e + 0	9.96e−1	2.74e−2	1.00e + 0	6.17e−3	3.31e−3		2.26e−4	2.42e−5	1.57e−1	1.00e + 0
E	3.09e−6		1.46e−3	1.82e−3	1.37e−4	1.19e−3	4.11e−6	1.43e−3		1.75e−3	7.48e−6	4.59e−4	3.53e−4
V	8.78e−3		1.19e−3	1.83e−3	4.99e−3	1.19e−3	6.67e−4	4.31e−1		7.75e + 0	3.09e−1	2.94e−3	3.53e−4
θ˜			1.23e + 0	9.96e−1	2.73e−2	9.97e−1	6.16e−3	3.27e−3			2.15e−5	1.59e−1	9.96e−1
E			1.41e−3	1.72e−3	1.30e−4	1.11e−3	3.94e−6	1.76e−4			6.92e−6	4.13e−4	3.32e−4
V			1.15e−3	1.73e−3	4.73e−3	1.12e−3	6.39e−4	5.37e−2			3.21e−1	2.60e−3	3.33e−4
θ˜			1.23e + 0		2.73e−2	9.97e−1	6.16e−3				1.85e−5	1.59e−1	9.96e−1
E			1.39e−3		1.13e−4	1.06e−3	1.85e−6				1.23e−6	3.71e−4	1.04e−4
V			1.13e−3		4.14e−3	1.06e−3	3.00e−4				6.64e−2	2.33e−3	1.04e−4
θ˜					2.68e−2	9.98e−1	6.15e−3				1.98e−05		9.97e−1
E					1.94e−5	2.28e−5	1.35e−6				1.24e−06		4.57e−5
V					7.25e−4	2.28e−5	2.20e−4				6.24e−02		4.59e−5
θ˜					2.68e−2	9.99e−1	6.15e−3						9.98e−1
E					1.87e−5	2.27e−5	1.33e−6						4.55e−5
V					7.01e−4	2.27e−5	2.16e−4						4.57e−5
θ˜						9.98e−1	6.16e−3						9.97e−1
E						2.26e−5	1.33e−6						4.51e−5
V						2.26e−5	2.15e−4						4.53e−4

Parameter identifiability might be misleading without the investigation of the residual of the model [51]. The Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC) indices make use of residuals to determine the quality of models in the presence of a given set of data. Table 2 shows the AIC and BIC estimates for each parameter set . It is seen that the best improvements occur from to for . Thus, the best case scenario of uncertainty quantification obtained for this analysis is that of .

Table 2

AIC and BIC metrics to estimate the quality of the model with different parameter sets.

θ˜	AIC	BIC
θ_1	164.80	198.67
θ_2	163.39	192.05
θ_3	150.84	174.28
θ_4	144.85	163.09
θ_5	142.63	155.68
θ_6	140.80	151.22
θ_7	141.68	149.50

Finally, we present results in Table 3 obtained from solving the inverse problem for using the data given in TDH [52]. The remaining parameters are fixed using the nominal parameters above and is the number of days used in the simulation

Table 3

Final parameter estimates from the data given in TDH [52] with .

Parameters	Estimates	Standard errors	Coefficients of variation
β_S	9.9792e−01	1.8476e−01	1.8514e−01
γ_1	7.2723e−02	6.2483e−03	8.5918e−02
γ_2	2.9536e−01	2.2578e−02	7.6443e−02
ρ	1.9201e−02	4.9837e−04	2.5956e−02
μ_4	2.5413e−04	2.7412e−04	1.0786e + 00
η	6.7561e−01	4.2262e−02	6.2554e−02
α	9.8887e−01	2.9311e−02	2.9641e−02

Results and discussions

Seven inverse problems for California, Florida, Georgia, Maryland, Texas, Washington and Wisconsin were solved to estimate the full parameter sets of the model. As seen in Fig. 5 , the fits are reasonably good even for states like Tennessee and Wisconsin whose current infected population begins to flatten. Table 7 (see Appendix) shows the fit parameter sets for each of the states with being the number of days used in the simulation. We see that the contact rates and for the asymptomatic and symptomatic population lies within 0.5–1.5 days a range suggested by Li et al. [53], Read et al. [54], Shen et al. [55], Eikenberry et al. [56]. The incubation period is seen to be relatively small, around 1–2 days. The average length of active infection and of the asymptomatic and symptomatic population, respectively, is seen to be around 1–27 days which is in line with suggested range of days in the literature [57], [58]. The recovery rate via hospitalization is seen to be very small for California, Florida, Georgia and Washington where the data records little or no recovery at all for infected patients. The death rates from the exposed, asymptomatic infected, symptomatic infected and hospitalized compartments are seen to be relatively small ranging from 0 to 0.01 day. With a value of around 0, Florida and Georgia has no death from the hospitalized compartments. This is reasonable since is also practically 0 indicating no recovery for this states as seen in the data (not shown here).

Fig. 5

Model fits of the compartmental model to infected, recovered and deaths.

Table 7

Model fit parameters for some selected states in the US, shows the number of days used for each state.

θ	California	Florida	Georgia	Maryland	Texas	Washington	Wisconsin
b_1	1.0000e−15	1.0000e−15	4.6495e−02	1.0080e−05	1.4196e−05	2.0116e−02	1.0000e−15
β_A	8.8973e−05	9.8043e−01	1.0939e + 00	1.3360e + 00	5.6047e−01	9.5295e−01	6.3325e−01
β_S	1.1213e + 00	6.9601e−01	5.1042e−01	1.3320e + 00	1.1989e−02	1.0000e + 00	1.0773e + 00
σ	9.9999e−01	1.0000e + 00	7.5663e−01	1.0000e + 00	1.0000e + 00	7.3519e−01	7.1669e−01
γ_1	3.5769e−02	8.3577e−01	9.8459e−01	9.9686e−01	1.2284e−01	8.9911e−01	1.0000e + 00
γ_2	9.1747e−01	7.8446e−01	7.6914e−01	9.6675e−01	6.6645e−02	8.9098e−01	1.0000e + 00
ρ	1.3366e−07	1.0000e−15	1.0000e−15	4.5126e−05	5.4691e−02	2.9082e−06	1.0325e−03
μ_1	2.7409e−03	6.7600e−09	6.5427e−04	6.6585e−04	1.0000e−15	6.9886e−04	8.7380e−04
μ_2	1.1988e−05	8.1025e−04	3.4154e−04	1.0000e−15	1.2534e−02	3.1211e−06	8.6293e−02
μ_3	6.8398e−05	6.9384e−05	3.2020e−08	1.0000e−15	1.0000e−15	3.6933e−06	6.0737e−04
μ_4	1.0000e−15	1.0000e−15	1.0000e−15	.0000e−15	1.0000e−15	1.2343e−05	1.0000e−15
η	1.1097e−01	9.4378e−01	1.0000e + 00	9.9986e−01	3.0315e−01	3.9910e−01	4.6680e−05
α	9.9980e−01	8.7834e−01	9.9967e−01	7.7423e−01	9.8525e−01	1.0000e + 00	1.0000e + 00
T^*	150	110	110	110	130	150	140

Table 4 shows the basic reproduction number values computed for our model using (3) and the parameters listed in Table 7. The epidemic is expected to continue indefinitely if as predicted for all states. This suggest that stricter measures such as the use of masks in public places, social distancing, contact tracing and even longer stay-at-home orders need to be enforced in order to eradicate the epidemic. Fig. 5 shows plot of the infected (asymptomatic plus symptomatic), recovered and deaths. We see that shape of most curves are similar except for Wisconsin and Tennessee where the infected population begin to flatten over time. In particular, spikes of the curves in Fig. 5 begin around late March or early April 2020 and continues to rise.

Table 4

The basic reproduction number values for the selected states.

States	R_0
California	1.0837
Florida	1.1340
Georgia	1.1096
Maryland	1.2501
Texas	1.2367
Tennessee	1.0343
Washington	1.0964
Wisconsin	1.0754

Table 5 shows a computational comparison of the proposed model and the integer-order counterpart using the sum squared errors (SSE) obtained for the infected, recovered and dead populations. As seen in the table, the fractional-order model shows superiority over the integer-order model in fitting the COVID-19 pandemic. In Fig. 6 , we show the cumulative infected and hospitalized population. This plot also shows that while the infected population flatten over time for states like California, Florida, Maryland, Tennessee, Texas and Wisconsin, it increases for Georgia and Washington. If drastic measures are not taken, the model suggest that many will be infected, at some point, before the end of 2021 in these states.

Table 5

Comparison of sum squared errors (SSE) of the proposed model to the integer-order counterpart.

States	Integer-order model	Fractional-order model
California	8.323e + 09	2.898e + 09
Florida	3.688e + 09	2.544e + 09
Georgia	2.995e + 09	7.464e + 08
Maryland	6.277e + 08	6.551e + 08
Tennessee	3.308e + 08	1.054e + 08
Texas	5.473e + 09	2.712e + 09
Washington	8.780e + 08	3.103e + 08
Wisconsin	7.404e + 08	3.165e + 08

Fig. 6

Model prediction for cumulative infected and hospitalized populations.

Conclusions

A fractional-order compartmental model is proposed to study the spread and dynamics of the COVID-19 pandemic. We studied the properties of the model and discussed the parameters of the model in the context of sensitivity and identifiability. We solve several inverse problems to estimate the fit parameters of the model using data from John Hopkins University [59] for some selected states in the US. The basic reproduction number of the model was computed and shows that the states considered have values slightly greater than the critical value of one. The model suggests that stricter or aggressive measures need be enforced in order to slow the spread of the virus.

CRediT authorship contribution statement

T.A. Biala: Conceptualization, Methodology, Investigation, Data curation, Writing - original draft. A.Q.M. Khaliq: Supervision, Validation, Software, Writing - review & editing.

Declaration of Competing Interest

The authors declare that there is no conflict of interest regarding the preparation and submission of the manuscript.