A fractional-order compartmental model for predicting the spread of the Covid-19 pandemic
AA fractional-order compartmental model for predictingthe spread of the Covid-19 pandemic
T. A. Biala*, A. Q. M. Khaliq
Center for Computational Science and Department of Mathematical Sciences,Middle Tennessee State University,Murfreesboro, TN 37132-0001, USA.
Abstract
We propose a time-fractional compartmental model (SEI A I S HRD) comprisingof the susceptible, exposed, infected (asymptomatic and symptomatic), hospi-talized, recovered and dead population for the Covid-19 pandemic. We studythe properties and dynamics of the proposed model. The conditions under whichthe disease-free and endemic equilibrium points are asymptotically stable arediscussed. Furthermore, we study the sensitivity of the parameters and usethe data from Tennessee state (as a case study) to discuss identifiability of theparameters of the model. The non-negative parameters in the model are ob-tained by solving inverse problems with empirical data from California, Florida,Georgia, Maryland, Tennessee, Texas, Washington and Wisconsin. The basicreproduction number is seen to be slightly above the critical value of one sug-gesting 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 inorder to mitigate the spread of the virus. As stay-at-home orders are rescindedin some of these states, we see that the number of cases began to increase al-most immediately and may continue to rise until the end of the year 2020 unlessstricter measures are taken.
Keywords:
Time-fractional model, SEIR model, Covid-19, Sensitivityanalysis, Parameter estimation and identifiability.
1. Introduction
Fractional differential equations (FDEs) are used to model complex phenom-ena such as such as the modeling of memory-dependent phenomena (DiGuiseppe et al. [1], Baleanu et al. [2], Podlubny [3]), mechanical properties of mate-rials (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
Email address: *[email protected] (T. A. Biala*)
Preprint submitted to arXiv July 9, 2020 a r X i v : . [ q - b i o . P E ] J u l erve as a generalization of the integer-order differential equations and give moredegree of freedom for modeling of biological and physical processes. FDEs havebeen applied in biological tissues [10], DNA sequencing [11], Pine Wilt disease[12], lung tissue mechanics and models [13] , harmonic oscillators [14], Denguefever [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 ofvaccines. The novel coronavirus began in December 2019 in China. Over thelast few months, it has spread rapidly leading to over 400,000 deaths across theglobe. The first occurrence in the United States was seen around mid Januaryin Washington [20] and has spread across America with over 100,000 deaths and1.5 million infected. The pandemic has disrupted the day-to-day activities ofthe human life with over six million jobs lost in the United States. Several ac-tions and measures have been taken by the federal, state and local governmentsto mitigate the spread of the pandemic. The most prominent measures takeninclude social distancing, testing, use of face-masks and contact tracing. It isimportant to model this pandemic in order to better understand the spread anddynamics as well as address the challenges of the pandemic. In short, mathe-matical models are important to guide the decisions of health and governmentofficials.The goal of this study is to examine and analyze the spread of the pandemic us-ing a modification of the Susceptible-Exposed-Infected-Recovered (SEIR) modelwith a time-fractional derivative. The use of fractional derivatives in the modelstems from the fact that the spread of infectious diseases depends not only onthe current state but also on its past states (history or memory dependency).Additionally, time-fractional order models reduce errors resulting from neglectof parameters in models. We shall focus on some selected states in the US. Wenote that models that consider the US as a whole may be misleading and havelimited applicability as different states have different economical and politicalimpacts on the pandemic. For example, while some states such as Maryland,New Jersey, New York, Connecticut, among others, enforced the use of masksin public places and longer stay-at-home order [21], other states do not enforcethese measures thereby allowing for a possibility of increase of infected indi-viduals in such states. There have been several models for the study of thepandemic. Lu et al. [22] considered a fractional-order SEIHDR model whichincorporates intercity movements. Liu et al. [23] studied the dynamics of thepandemic by considering asymptomatic and symptomatic infected populationsseparately. Wu et al. [24] studied domestic and international spread of the pan-demic by using different data sets. Zhao and Chen [25] discussed the dynamicsof the pandemic by considering the Susceptible, unquarantined infected, quar-antined infected and Confirmed infected (SUQC) model and parametrize theintervention effect of control measures. Zhang et al. [26] considered a fractionalSEIR model with different order of the time-fractional derivative for each of thedifferent population being studied.The remainder of the paper is organized as follows: In section 2, we givesome preliminary results and definitions from fractional Calculus. Furthermore,2e discuss the formulation of the model by considering seven compartments:Susceptible-Exposed-Asymptomatic infected-Symptomatic infected, Hospital-ized, Recovered and Dead populations (SEI A I S HRD). Section 3 details the prop-erties and theoretical analysis of the model. Section 4 discusses the parametersensitivity and identifiability analysis. In section 5, we applied the model to ob-served data for some selected states in the US. In particular, this section detailssolving several inverse problems for parameter estimation and computation ofthe basic reproduction number. Finally, we give some concluding remarks inSection 6.
2. Model Setup
In this section, we present some preliminary results and definitions in frac-tional calculus.
Definition 2.1. [3] The gamma function Γ( α ) is defined by the integral Γ( α ) = (cid:90) ∞ e − x x α dx which converges in the right half of the complex plane R e ( z ) > . Definition 2.2. [27] For any t > , the Caputo-fractional derivative of order α, ( n < α ≤ n − of a function f ( t ) is defined as t D αt f ( t ) = 1Γ( n − α ) (cid:90) tt ( t − τ ) n − α − f ( n ) ( τ ) dτ. Definition 2.3. [3] The Mittag-Leffler function which generalizes the exponen-tial function for fractional calculus is defined as E α,β ( z ) = ∞ (cid:88) k =0 z k Γ( αk + 1) , α ∈ R + , z ∈ C . Remark 2.1.
More generally, the two parameter Mittag-Leffler function is defined as E α ( z ) = ∞ (cid:88) k =0 z k Γ( αk + β ) , α, β ∈ R + , z ∈ C . It has the following properties: E α,β ( z ) = zE α,α + β ( z ) + 1Γ( β ) . D αt e λt = t − α E , − α ( λt ) . D αt E α, ( λt α ) = λE α, ( λt α ) . efinition 2.4. [28] A point x ∗ is said to be an equilibrium point of the system t D αt = f ( t, x ( t )) , x ( t ) > if and only if f ( t, x ∗ ( t )) = 0 . Definition 2.5. [22] An equilibrium point x ∗ of the system t D αt x ( t ) = f ( t, x ( t )) , x ( t ) > is said to be asymptotically stable if all he eigenvalues of the Jacobian matrix J = ∂f /∂x , evaluated at the equilibrium point, satisfies | arg ( λ i ) | > απ , where λ i are the eigenvalues of J .2.2. Model Formulation The model discussed in this work is a modification of the SEIR model havingthree additional compartments. We consider a SEI A I S HRD compartmentalmodel which comprises of the susceptible, exposed, infected (asymptomatic andsymptomatic), hospitalized, recovered and dead population. We assume thatthe natural death and birth rates are the same. We further assume that deathsin the S and R compartments are due to natural deaths and deaths in the othercompartments are as a result of the pandemic. The schema below shows thetransmission flow of the model.
Figure 1: Schematic diagram of the proposed SEI A I S HRD model D αt S ( t ) = b N − S ( t ) N ( β A I A ( t ) + β S I S ( t )) − b S ( t ) , D αt E ( t ) = S ( t ) N ( β A I A ( t ) + β S I S ( t )) − ( σ + µ ) E ( t ) , D αt I A ( t ) = ησE ( t ) − ( γ + µ ) I A ( t ) , D αt I S ( t ) = (1 − η ) σE ( t ) − ( γ + µ ) I S ( t ) , D αt H ( t ) = γ I A ( t ) + γ I S ( t ) − ( ρ + µ ) H ( t ) , D αt R ( t ) = ρH ( t ) − b R ( t ) , D αt D ( t ) = µ E ( t ) + µ I A ( t ) + µ I S ( t ) + µ H ( t ) . The system has the associated initial data S (0) = S ≥ , E (0) = E ≥ , I A (0) = I A ≥ ,I S (0) = I S ≥ , H (0) = H ≥ , R (0) = R ≥ ,D (0) = D ≥ . The total population is N which is further divided into S ( t ) , E ( t ) , I A ( t ) , I S ( t ), H ( t ) and R ( t ), b is the natural birth rate, β A and β S are the transmission ratesto the susceptible population from the asymptomatic and symptomatic popu-lations, respectively. 1 /σ is the incubation period of an exposed individual, η denotes the fraction of the exposed population that becomes asymptomatic afterthe incubation period and the remaining (1 − η ) of the population are symp-tomatic. γ and γ are the infectious rates of an asymptomatic and a symp-tomatic individual, respectively. ρ is the recovery rate through hospitalizationand, µ , µ , µ and µ are the mortality rates of the exposed, asymptomatic,symptomatic and hospitalized populations, respectively. We note that the deathpopulation in the model comprises of deaths during exposure µ E ( t ), infectiousperiod ( µ I A and µ I S ) and hospitalization µ H , and assume that deaths dueto other natural occurrences are negligible for this populations.We note that the parameters of the model are non-negative and have dimen-sions 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 D αt S ( t ) = b α N − S ( t ) N ( β αA I A ( t ) + β αS I S ( t )) − b α S ( t ) , D αt E ( t ) = S ( t ) N ( β αA I A ( t ) + β αS I S ( t )) − ( σ α + µ α ) E ( t ) , D αt I A ( t ) = ησ α E ( t ) − ( γ α + µ α ) I A ( t ) , D αt I S ( t ) = (1 − η ) σ α E ( t ) − ( γ α + µ α ) I S ( t ) , D αt H ( t ) = γ α I A ( t ) + γ α I S ( t ) − ( ρ α + µ α ) H ( t ) , D αt R ( t ) = ρ α H ( t ) − b α R ( t ) , D αt D ( t ) = µ α E ( t ) + µ α I A ( t ) + µ α I S ( t ) + µ α H ( t ) . (1)5 . Model Analysis In this section, we discuss the properties of the model beginning with theexistence, uniqueness, non-negativity and boundedness of solutions of the model(1). For simplicity in analysis, we reduce the system (1) to D αt S ( t ) = b α N − S ( t ) N ( β αA I A ( t ) + β αS I S ( t )) − b α S ( t ) , D αt E ( t ) = S ( t ) N ( β αA I A ( t ) + β αS I S ( t )) − ( σ α + µ α ) E ( t ) , D αt I A ( t ) = ησ α E ( t ) − ( γ α + µ α ) I A ( t ) , D αt I S ( t ) = (1 − η ) σ α E ( t ) − ( γ α + µ α ) I S ( t ) , D αt H ( t ) = γ α I A ( t ) + γ α I S ( t ) − ( ρ α + µ α ) H ( t ) , D αt R ( t ) = ρ α H ( t ) − b α R ( t ) , (2)since D is a linear combination of populations in some of the other compart-ments. Theorem 3.1.
There exist a unique solution to the system (2) and the solutionis non-negative and bounded for any given initial data ( S , E , I A , I S , H , R ) ≥ ∈ R .Proof. See Appendix A. R We shall use the next generation matrix originally proposed by Diekmann et al. [29] and further elaborated on by van den Driesche and Watmough [30]and Diekmann et al. [31] to determine R . According to system (2), the disease-free equilibrium point (DFE) is ( N, , , , , Y = ( Y , Y , Y , Y ) = ( E, I A , I S , H ) containing the infected individualsand let Y ∗ be the DFE point. Since the DFE point exists and is stable (shownin the next section) in the absence of any disease, then the linearized equationat the DFE is D αt Y i = F i ( Y ) − V i ( Y ) , i = 1(1)4 , where F i ( Y ) is the rate of appearance of new infections in compartment i and V i ( Y ) is the rate of transfer of infections to and from compartment i . We furtherdefine F = ∂ F ( Y ) ∂Y j (cid:12)(cid:12) Y = Y ∗ and V = ∂ V i ( Y ) ∂Y j (cid:12)(cid:12) Y = Y ∗ , i, j = 1(1)4 . Then ρ ( F V − ) is the basic reproduction number R , where ρ ( x ) is the spectralradius of x and F V − is the next generation matrix. For the system (2), F = β αA β αS
00 0 0 00 0 0 00 0 0 0 , V = σ α + µ α − ησ α γ α + µ α − (1 − η ) σ α γ α + µ α − γ α − γ α ( ρ α + µ α ) R = σ α [ ηβ αA ( γ α + µ α ) + (1 − η ) β αS ( γ α + µ α )]( σ α + µ α )( γ α + µ α )( γ α + µ α ) . (3) Lemma 3.1.
The fractional system (2) has at most two equilibrium points a disease-free equilibrium point DFE = ( N, , , , , an endemic equilibrium point EE = ( S ∗ , E ∗ , I ∗ A , I ∗ S , H ∗ , R ∗ ) ,where S ∗ = S R ,E ∗ = b α S R ( σ α + µ α ) (cid:18) − R (cid:19) ,I ∗ A = ηb α σ α S R ( σ α + µ α )( γ α + µ α ) (cid:18) − R (cid:19) ,I ∗ S = (1 − η ) b α σ α S R ( σ α + µ α )( γ α + µ α ) (cid:18) − R (cid:19) ,H ∗ = b α σ α S R ( σ α + µ α )( ρ α + µ α ) (cid:18) − R (cid:19) (cid:18) ηγ α γ α + (1 − η ) γ α γ α + µ α (cid:19) ,R ∗ = ρ α b α σ α S R ( σ α + µ α )( ρ α + µ α ) (cid:18) − R (cid:19) (cid:18) ηγ α γ α + (1 − η ) γ α γ α + µ α (cid:19) . Theorem 3.2.
The DFE point is locally asymptotically stable if R < and K < , where K = σ α [ ηβ αA ( σ α + γ α + µ α + µ α ) + (1 − η ) β αS ( σ α + γ α + µ α + µ α )]( σ α + γ α + µ α + µ α ) ( σ α + γ α + µ α + µ α ) ( γ α + γ α + µ α + µ α ) . Proof.
See Appendix B.
Theorem 3.3.
The EE point is locally asymptotically stable if R > and ( A A − A ) A ≥ A A , where A = b α (2 R −
1) + R ( σ α + γ α + γ α + µ α + µ α + µ α ) ,A = − σ α R ( ηβ αA + (1 − η ) β αS ) + b α R (2 R − σ α + γ α + γ α + µ α + µ α + µ α )+ R [( γ α + µ α )( σ α + µ α ) + ( γ α + µ )( γ α + µ α ) + ( σ α + µ α )( γ α + µ α )] ,A = b α R (2 R − γ α + µ α )( σ α + µ α ) + ( γ α + µ )( γ α + µ α ) + ( σ α + µ α )( γ α + µ α )]+ R ( σ α + µ α )( γ α + µ α )( γ α + µ ) − σ α ((1 − η )( γ α + µ α ) β αS + η ( γ α + µ α ) β αA )) ,A = b α R (2 R − σ α + µ α )( γ α + µ α )( γ α + µ α ) − σ α ((1 − η ) β αS ( γ α + µ α )+ ηβ αA ( γ α + µ α ))] Proof.
See Appendix C. 7 . Parameter Sensitivity and Identifiability Analysis
We discuss the sensitivity and identifiability of the parameters with respectto the proposed model.
The sensitivity analysis (SA) deals with the significance or importance of theparameters in the model. In particular, it finds the most influential parametersthat drives the dynamics of the model. It also describes the extent to whichparameter changes affects the result of the methods or models with the goal ofidentifying the best set of parameters that describes the process or phenomenain question.There are several SA methods which are broadly classified as local and globalmethods. In this work, we shall focus on the Morris screening method (localmethod) and Sobol analysis method (global method).
The Morris screening method is a local sensitivity measure that makes useof the first order derivative of an output function y = f ( θ ) = f ( θ , · · · , θ p ) withrespect to the input parameter θ . It measures the effect of the output when theinput variable is perturbed one at a time around a nominal value. It serves asa first check, in most analysis, in screening parameters for identifiability. Themethod evaluates elementary effects [32, 33, 34] with the i th parameter throughthe forward perturbation g i ( θ ) = f ( θ , θ , · · · , θ i + ∆ θ i , · · · , θ p ) − f ( θ i , · · · , θ p )∆ θ i , i = 1(1) p Morris [35] proposed two sensitivity measures, the mean ( µ ) and the standarddeviation (˜ σ ) of the elementary effects. For non-monotonic models, µ may leadto a very small value due to cancellation effects. For this reason, Campolongo et al. [36] proposed the use of absolute values for evaluating the mean. In orderto obtain a dimension-free sensitivity, we prefer the use of the sensitivity measure δ given in Brun et al. [37] as δ i = (cid:118)(cid:117)(cid:117)(cid:116) N N (cid:88) j =1 ˜ g ij , i = 1(1) p, and j = 1(1) N, where N is the number of sample points and˜ g i ( θ ) = f ( θ , θ , · · · , θ i + ∆ θ i , · · · , θ p ) − f ( θ , · · · , θ p )∆ θ i θ i f ( θ , · · · , θ p ) . A common practice in the literature [38, 39, 34] is to plot the indices δ against ˜ σ ,the standard deviation. We see from Fig. 2(a) that the fractional order α has thehighest influence on the model output over time. The transmission rates for the8ymptomatic and asymptomatic populations and, the fractional parameter η arethe next most influential parameters in the model. This is further corroboratedby Fig. 2(b). The four parameters ( α β S , β A and η ) denoted by red squares arethe most important parameters, that is they have the largest δ and ˜ σ values.The parameters ( b , ρ, µ , µ , µ and µ ) represented by the blue squares havethe least influence on the model output and can be considered unimportant.The other parameters represented by the green squares have more influencethan the parameters represented by the blue squares.One major setback of the Morris screening test for sensitivity analysis is theconsideration of each parameter individually and independently of the otherparameters. In real applications, this is not true as parameters have collinearityand dependencies on one another. (a) Sensitivity of the parameters over time(b) Parameter Importance Figure 2: Morris screening test
The Sobol method is a variance-based sensitivity analysis method which un-like the Morris screening method takes into account the effect of the relationshipbetween each parameters of the model. It uses the decomposition of variance to9alculate Sobol’s sensitivity indices: first and total order sensitivity measures.The basic idea of the Sobol’s method is the decomposition of the model outputfunction y = f ( θ , · · · , θ p ) into summands of increasing dimensionality, that is V ( y ) = V , ··· ,p + p (cid:88) i =1 V i + p (cid:88) i =1 p (cid:88) j> V i,j + · · · where V i is the partial variance of the contribution of the parameter θ i and V i, ··· ,s is the partial variances caused by the interaction of the parameters ( θ , · · · , θ s )for s ≤ p .The first order sensitivity index measures the main effect of parameter θ i on themodel output; that is the partial contribution of θ i to the variance V ( y ).Theindex [40, 41] is defined as S i = V i V ( y ) . The larger this index, the more sensitive the parameter is to the model output[40, 41]. Using the law of total variances [41, 34], the index can also be expressedas V ( y ) = V θ i ( E θ ∼ i ( y | θ i )) + E θ i ( V θ ∼ i ( y | θ i ))and S i = V θ i ( E θ ∼ i ( y | θ i )) V ( y )where V θ i ( E θ ∼ i ( y | θ i )) is the partial variance caused by θ i and E θ ∼ i ( y | θ i ) isthe mean of the model output calculated by using all the values of the otherparameters θ ∼ i (except θ i ) and V ( y ) is the total variance.The total sensitivity indices [42] measures the effects of parameter θ i and theinteraction with the other parameters. It is defined as S T i = V i + V i,j + · · · + V i,j, ··· ,p V ( y ) . The total variance, V ( y ), for this index is given as V ( y ) = V θ ∼ i ( E θ i ( y | θ ∼ i )) + E θ ∼ i ( V θ i ( y | θ ∼ i ))and S T i = E θ ∼ i ( V θ i ( y | θ ∼ i )) V ( y )= V ( y ) − V θ ∼ i ( E θ i ( y | θ ∼ i )) V ( y ) . The mean and the variance can be evaluated using quasirandom samplingmethod [43, 34] and are given as V θ i ( E θ ∼ i ( y | θ i )) = 1 N N (cid:88) j =1 f ( B j ) (cid:0) f ( A i B ,j ) − f ( A j ) (cid:1) , E θ ∼ i ( V θ i ( y | θ ∼ i )) = 12 N N (cid:88) j =1 (cid:0) f ( A j ) − f ( A i B ,j ) (cid:1) , where A and B are two independent parameter sample matrices of dimensions N × p . We shall use the python SALib package [44] to compute the first andtotal order variance indices. Fig. (3) shows that the fractional order β S has thehighest interaction with the other parameters. These results are consistent withthe results in the Morris screening test as important parameters of the test showhigh interaction with the other parameters. Figure 3: Sobol Sensitivity Indices
The concept of identifiability is dependent on sensitivity. It entails the se-lection of the subset of parameters of a model having little or no collinearityand uncertainty, and which can be identified uniquely from a given set of ob-served data or measurements. In other words, it answers the question ”Can theavailable data be described by the model and the selected subset of parame-ters?”. There are several techniques or tests for parameter identifiability. Mostof the tests are based on the Fisher information matrix (FIM) F = χ T χ where χ = ∂y/∂θ for a model output function y . Cobelli and Di Stefano [45] showedthat a sufficient condition for identifiability is the non-singularity of FIM. Burth et al. [46] proposed an iterative estimation process which implements a reduced-order estimation by finding parameters whose axis lie closest to the directionof FIM. The associated parameter values are then fixed at prior estimates dur-ing the iterated process. Brun et al. [37] studied parameter identifiability usingtwo indices; a parameter importance ranking index δ and a collinearity index γ K χ T χ correspondingto the parameter subset K . Cintr´on-Arias et al. [47] explained the need fora good parameter subset for identifiability to satisfy the full rank test. Theyfurther introduced two indices; the selection score and the condition number of χ T χ . The smaller these indices the lesser the collinearity and uncertainty inthe parameter values of the subset. Finally, they used the coefficient of varia-tion index to examine the effect of parameters in the parameter subset. In thiswork, we shall use the test proposed by Cintr´on-Arias et al. in identifying theparameters. The algorithm can be summarized in the following steps Algorithm 1
Algorithm for Parameter subset Selection [47] Perform a combinatorial search for all possible parameter subsets. Let S p = { θ = ( λ , λ , · · · , λ p ) ∈ R p (cid:12)(cid:12) λ k ∈ I and λ k (cid:54) = λ m ∀ k, m = 1 , · · · , p } , where I = { b , β A , β S , σ, γ , γ , ρ, µ , µ , µ , µ , η, α } . Select parameter subsets that pass the full rank test; that isΘ p = { θ (cid:12)(cid:12) θ ∈ S p ⊂ R p , Rank( χ ( θ )) = p } . For each θ ∈ Θ p , calculate the parameter selection score ζ ( θ ) = | ϑ ( θ ) | where ϑ = (cid:112) Σ( θ ) ii θ i , i = 1 , · · · , p, and Σ( θ ) = σ (cid:2) χ T ( θ ) χ ( θ ) (cid:3) − ∈ R p . Calculate the condition number κ ( χ ( θ )) for each parameter subset θ ∈ Θ p .The smaller the values of κ ( χ ( θ )) and ϑ ( θ ), the lower the uncertainty pos-sibilities in the estimate.To discuss the results in this section, we shall use the state of Tennesseeas a case study to understand parameter identifiability. Furthermore, we usedthe following values obtained using a random search algorithm as the nominalparameter set θ for the model: b = 8 . , β A = 1 . , β S = 1 . , σ = 9 . ,γ = 2 . , γ = 1 . , ρ = 6 . , µ = 3 . ,µ = 1 . , µ = 3 . , µ = 2 . , η = 1 . ,α = 1 . σ = 10. We further divide the parameters intothree groups according to their importance rankings discussed in the previoussection: S = ( β A , β S , η, α ) ,S = ( σ, γ , γ ) ,S = ( b , ρ, µ , µ , µ , µ ) , Figure 4: The condition number κ ( χ ( θ )) against the parameter selection scores ϑ ( θ ) of the N × p sensitivity matrices for all parameter subsets θ = Θ p with p = 2. Logarithmic scalesare used on both axis where S and S are the most and least influential parameter sets, respectively,while S contains more influential parameters than S . We display some selec-tions of the parameter subsets of size p in Table 1 where we have chosen thesubsets with the smallest score values. The entries in Table 1 are ordered withrespect to the selection score ϑ ( θ ) for each subset of same cardinality. A highselection score and condition number for a parameter subset indicates substan-tial collinearity and linear dependence, and thus is poorly identifiable even ifthe parameter subsets contains S , that is contains the set of most influentialparameters. We observe that most of the selections in Table 1 contains at leastone element in each of the groups listed above. This shows that while parameterimportance ranking is crucial in identifying parameters that drives the dynamicsof a model, it does not have substantial effect in identifiability. Identifiabilitydepends on proper selection of subsets including parameters in each of the threegroups above that describes the measurement or data.To have an idea of the variations of the condition number and the selectionscore, we give a plot of these values for p = 2 in fig. 4 (with logarithmic scales).Good parameter combination in fig. 4 corresponds to values in the lower leftcorner of the figure where the values, ϑ ( θ ) and κ ( χ ( θ )), are relatively small.13 Parameter Subsets κ ( χ ( θ )) ϑ ( θ )13 ( b , β A , β S , σ, γ , γ , ρ, µ , µ , µ , µ , η, α ) 9.530e+03 8.481e+0112 ( b , β A , β S , σ, γ , γ , ρ, µ , µ , µ , η, α ) 9.513e+03 6.147e+0011 ( b , β S , σ, γ , γ , ρ, µ , µ , µ , η, α ) 4.814e+03 2.479e-0410 ( b , β S , σ, γ , γ , ρ, µ , µ , η, α ) 3.696e+03 3.852e-07( b , β S , σ, γ , γ , ρ, µ , µ , η, α ) 3.685e+03 4.143e-069 ( β S , σ, γ , γ , ρ, µ , µ , η, α ) 3.411e+03 2.766e-07( b , β S , σ, γ , γ , ρ, µ , η, α ) 3.311e+03 2.824e-077 ( β S , γ , γ , ρ, µ , η, α ) 2.718e+03 5.344e-08( β S , γ , γ , ρ, µ , η, α ) 2.947e+03 5.388e-085 ( γ , γ , ρ, µ , α ) 6.171e+01 2.006e-09( β S , γ , ρ, µ , α ) 6.389e+01 2.006e-09( γ , ρ, µ , η, α ) 6.989e+01 2.010e-094 ( γ , ρ, η, α ) 6.258e+01 3.458e-10( γ , γ , ρ, α ) 5.289e+01 3.491e-10( β S , γ , ρ, α ) 5.349e+01 3.496e-103 ( γ , ρ, α ) 5.140e+01 5.431e-11( β S , ρ, α ) 5.203e+01 5.850e-11( ρ, η, α ) 6.224e+01 9.109e-11 Table 1: Selection scores and condition numbers for some selected parameter subsets
To further analyze the parameter identifiability of the model, we considerthe parameters subsets: θ = ( b , β A , β S , σ, γ , γ , ρ, µ , µ , µ , µ , η, α ) ,θ = ( b , β S , σ, γ , γ , ρ, µ , µ , µ , η, α ) ,θ = ( β S , σ, γ , γ , ρ, µ , µ , η, α ) ,θ = ( β S , γ , γ , ρ, µ , η, α ) ,θ = ( γ , γ , ρ, µ , α ) ,θ = ( γ , γ , ρ, α ) ,θ = ( γ , ρ, α )such that θ i +1 ⊂ θ i , i = 1 , · · · ,
6. The choice of these parameter subsets aredue to their relative small condition numbers and selection scores. In otherto create synthetic data, we assume the nominal parameter subsets and errorvariance (given at the beginning of this section) to be the true parameter vectorsand true variance. Furthermore, we add random noise to the model output asfollows: Y j = z ( t j , θ ) + σ N (0 , , j = 1 , · · · , N. We solve seven inverse problems for each of the parameter subsets θ i , i =1 , · · · ,
7. We analyze the result using the coefficient of variation and standarderror [47] given as SE j (˜ θ ) == (cid:113) ˜Σ j,j , j = 1 . · · · , p v j (˜ θ ) = SE j (˜ θ ) θ j , j = 1 . · · · , p, where ˜Σ j,j = ˜ σ (cid:104) χ (˜ θ ) T χ (˜ θ ) (cid:105) − and ˜ σ = 1 n − p | Y − z (˜ θ ) | .15 β A β S σ γ γ ρ µ µ µ µ η α ˜ θ . e - . e - . e + . e - . e - . e + . e - . e - . e - . e - . e - . e - . e + E . e - . e - . e - . e - . e - . e - . e - . e - . e - . e - . e - . e - . e - V . e - . e - . e - . e - . e - . e - . e - . e + . e + . e + . e + . e - . e - θ . e - . e + . e - . e - . e + . e - . e - . e - . e - . e - . e + E . e - . e - . e - . e - . e - . e - . e - . e - . e - . e - . e - V . e - . e - . e - . e - . e - . e - . e - . e + . e - . e - . e - θ . e + . e - . e - . e - . e - . e - . e - . e - . e - E . e - . e - . e - . e - . e - . e - . e - . e - . e - V . e - . e - . e - . e - . e - . e - . e - . e - . e - θ . e + . e - . e - . e - . e - . e - . e - E . e - . e - . e - . e - . e - . e - . e - V . e - . e - . e - . e - . e - . e - . e - θ . e - . e - . e - . e - . e - E . e - . e - . e - . e - . e - V . e - . e - . e - . e - . e - θ . e - . e - . e - . e - E . e - . e - . e - . e - V . e - . e - . e - . e - θ . e - . e - . e - E . e - . e - . e - V . e - . e - . e - T a b l e : P a r a m e t e r e s t i m a t e s f o r s o l v i n g s e v e n i n v e r s e p r o b l e m s f r o m a s y n t h e t i c d a t ag e n e r a t e du s i n g t h e g i v e nn o m i n a l p a r a m e t e r s a nd v a r i a n ce . F o r e a c hp a r a m e t e r s ub s e t , w e d i s p l a y t h ee s t i m a t e ( ˜ θ ) , t h e s t a nd a r d e rr o r , E a nd t h ec o e ffi c i e n t o f v a r i a t i o n , V = E / ˜ θ . θ AIC BIC θ θ θ θ θ θ θ Table 3: AIC and BIC metrics to estimate the quality of the model with different parametersets.
From Table 2, it is seen that the parameters ( µ , µ ) ⊂ θ have standarderrors that is approximately 20 times their estimates. This shows substantial un-certainty in these parameter values and any parameter subset containing theseparameters may result in illogical parameter estimation from observations. Thestandard errors of β S , γ , ρ, µ , µ , η, α in θ show improvements and implieslower 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 sub-sets as more parameters are removed. For instance, with the removal of β S and η in θ , it seen that the standard error for γ and γ dropped from 4% and1% to approximately 0.07% and 0.003%, respectively, of their estimates. Otherimprovements in θ include α where the standard error is reduced to at leastone-tenth of its standard error in θ . We note that there is no substantial gainin the removal of µ from θ and γ in θ as seen in Table 2.Parameter identifiability might be misleading without the investigation of theresidual of the model [47]. The Akaike Information Criterion (AIC) and BayesianInformation Criterion (BIC) indices make use of residuals to determine the qual-ity of models in the presence of a given set of data. Table 3 shows the AIC andBIC estimates for each parameter set θ i , i = 1 , · · · ,
7. It is seen that the bestimprovements occur from θ i to θ i +1 for i = 2 ,
3. Thus, the best case scenario ofuncertainty quantification obtained for this analysis is that of θ .Finally, we present results in Table 4 obtained from solving the inverse problemfor θ using the data given in [48]. The remaining parameters are fixed using thenominal parameters above and T ∗ is the number of days used in the simulation
5. Results and Discussions
Seven inverse problems for California, Florida, Georgia, Maryland, Texas,Washington and Wisconsin were solved to estimate the full parameter sets ofthe model. As seen in fig. 5, the fits are reasonably good even for states likeTennessee and Wisconsin whose current infected population begins to flatten.Table 5 shows the fit parameter sets for each of the states with T ∗ being the17arameters Estimates Standard Errors Coefficients of Variation β S γ γ ρ µ η α Table 4: Final parameter estimates from the data given in [48] with T ∗ = 100 number of days used in the simulation. We see that the contact rates β αA and β αS for the asymptomatic and symptomatic population lies within 0.5–1.5 days − , arange suggested by Li et al. [49], Read et al. [50], Shen et al. [51], Eikenberry et al. [52]. The incubation period 1 /σ α is seen to be relatively small, around1–2 days. The average length of active infection 1 /γ α and 1 /γ α of the asymp-tomatic and symptomatic population, respectively, is seen to be around 1–27days which is in line with suggested range of days in the literature [53, 54].The recovery rate via hospitalization is seen to be very small for California,Florida, Georgia and Washington where the data records little or no recoveryat all for infected patients. The death rates µ α , µ α , µ α , µ α from the exposed,asymptomatic infected, symptomatic infected and hospitalized compartmentsare seen to be relatively small ranging from 0–0.01 day − . With a µ α valueof around 0, Florida and Georgia has no death from the hospitalized compart-ments. This is reasonable since ρ α is also practically 0 indicating no recoveryfor this states as seen in the data (not shown here).Table 6 shows the basic reproduction number R values computed for our modelusing (3) and the parameters listed in Table 5. The epidemic is expected to con-tinue indefinitely if R > C a li f o r n i a F l o r i d a G e o r g i a M a r y l a nd T e x a s W a s h i n g t o n W i s c o n s i n b . e - . e - . e - . e - . e - . e - . e - β A . e - . e - . e + . e + . e - . e - . e - β S . e + . e - . e - . e + . e - . e + . e + σ . e - . e + . e - . e + . e + . e - . e - γ . e - . e - . e - . e - . e - . e - . e + γ . e - . e - . e - . e - . e - . e - . e + ρ . e - . e - . e - . e - . e - . e - . e - µ . e - . e - . e - . e - . e - . e - . e - µ . e - . e - . e - . e - . e - . e - . e - µ . e - . e - . e - . e - . e - . e - . e - µ . e - . e - . e - . e - . e - . e - . e - η . e - . e - . e + . e - . e - . e - . e - α . e - . e - . e - . e - . e - . e + . e + T ∗ T a b l e : M o d e l fi t p a r a m e t e r s f o r s o m e s e l ec t e d s t a t e s i n t h e U S , T ∗ s h o w s t h e nu m b e r o f d a y s u s e d f o r e a c h s t a t e . R California 1.0837Florida 1.1340Georgia 1.1096Maryland 1.2501Texas 1.2367Tennessee 1.0343Washington 1.0964Wisconsin 1.0754
Table 6: The basic reproduction number R values for the selected states.
6. Conclusions
A fractional-order compartmental model is proposed to study the spread anddynamics of the Covid-19 pandemic. We studied the properties of the modeland discussed the parameters of the model in the context of sensitivity andidentifiability. We solve several inverse problems to estimate the fit parametersof the model using data from John Hopkins University [55] for some selectedstates in the US. The basic reproduction number R of the model was computedand shows that the states considered have R values slightly greater than thecritical value of one. The model suggests that stricter or aggressive measuresneed be enforced in order to slow the spread of the virus.21 igure 5: Model fits of the compartmental model to infected, recovered and deaths. igure 6: Model prediction for selected US states until mid 2021 igure 7: Model prediction for cumulative infected and hospitalized populations eferences arXiv:2004.12308 .[23] Z. Liu, P. Magal, O. Seydi, G. Webb, Understanding Unreported Cases inthe COVID-19 Epidemic Outbreak in Wuhan, China, and the Importanceof Major Public Health Interventions, Biology 9 (50) (2020).[24] J. T. Wu, K. Leung, G. M. Leung, Nowcasting and forecasting the potentialdomestic and international spread of the 2019-nCoV outbreak originatingin Wuhan, China: a modelling study, The Lancet 395 (10225) (2020) 689– 697.[25] S. Zhao, H. Chen, Modeling the epidemic dynamics and control of COVID-19 outbreak in China, Quantitative Biology 8 (1) (2020) 11.[26] Y. Zhang, X. Yu, H. Sun, G. R. Tick, W. Wei, B. Jin, Applicability of timefractional derivative models for simulating the dynamics and mitigationscenarios of COVID-19, Chaos, Solitons & Fractals 138 (2020) 109959.[27] M. Caputo, Linear model of dissipation whose q is almost frequency inde-pendent. ii, Geophysical Journal International 13 (5) (1967) 529539.2628] Y. Li, Y. Chen, I. Podlubny, Mittagleffler stability of fractional order non-linear dynamic systems, Automatica 45 (8) (2009) 1965 – 1969.[29] O. Diekmann, J. Heesterbeek, J. Metz, On the definition and the compu-tation of the basic reproduction ratio R ppendicesAppendix A. Proof of Theorem 3.1. Proof.
By applying [56, Theorem 3.1], we obtain the existence of the solutions.To show the uniqueness and boundedness of solutions, it suffices to show, by[56, Remark 3.2] and Rademacher’s theorem, that F = ( f , f , f , f , f , f ) islocally Lipschitz continuous where f = b α N − SN ( β αA I A + β αS I S ) − b α S,f = SN ( β αA I A + β αS I S ) − ( σ α + µ α ) E,f = ησ α E − ( γ α + µ α ) I A ,f = (1 − η ) σ α E − ( γ α + µ α ) I S ,f = γ α I A + γ α I S − ( ρ α + µ α ) H,f = ρ α H − b α R. Let X = ( S, E, I A , I S , H, R ), ˜ X = ( ˜ S, ˜ E, ˜ I A , ˜ I S , ˜ H, ˜ R ) and || · || denote the L norm, then || F ( X ) − F ( ˜ X ) || ≤ || f ( X ) − f ( ˜ X ) || + || f ( X ) − f ( ˜ X ) || + || f ( X ) − f ( ˜ X ) || + || f ( X ) − f ( ˜ X ) || + || f ( X ) − f ( ˜ X ) || + || f ( X ) − f ( ˜ X ) ||≤ L || X − ˜ X || , where L = max ≤ i ≤ L i and L = b α ( N + 1) + β αA + β αS , L = β αA + β αS + σ α + µ α , L = ησ α + γ α + µ α , L = (1 − η ) σ α + γ α + µ α , L = γ α + γ α + ρ α + µ α and L = ρ α + b α . Thus, F satisfies the local Lipschitz conditions with respect to X which proves the uniqueness and boundedness of solution to (2). Next we showthe non-negativity of solutions. At first, we consider moving along the S -axis,that is E (0) = I A (0) = I S (0) = H (0) = R (0) = 0 and 0 < S (0) = S ≤ N , then D αt S ( t ) = b α N − b α S whose solution is given as S ( t ) = S E α, ( − b α t α ) + b α N t α E α,α +1 ( − b α t α ) > b α > t >
0. In a similar manner, moving along each of the otherrespective axis (that is all initial conditions are zeros except for the axis beingconsidered), it is easy to show that E ( t ) = E α, ( − ( σ α + µ α ) t α ) E ≥ I A ( t ) = E α, ( − ( γ α + µ α ) t α ) I A ≥ I S ( t ) = E α, ( − ( γ α + µ α ) t α ) I S ≥ H ( t ) = E α, ( − ( ρ α + µ α ) t α ) H ≥ R ( t ) = E α, ( − b α t α ) R ≥ . E − I A − I S − H − R plane, then let S ( t ∗ ) = 0, E ( t ∗ ) > I A ( t ∗ ) > I S ( t ∗ ) > H ( t ∗ ) > R ( t ∗ ) > t ∗ such that S ( t ) S ( t ) − S ( t ∗ ) = 1Γ( α ) D αt ( τ )( t − t ∗ ) α for some τ ∈ [ t ∗ , t ), we see that S ( t ) > S ( t ∗ ). This contradicts our previousstatement. Similar arguments can be used for each of the remaining populationvariables. Appendix B. Proof of Theorem 3.2
Proof.
The Jacobian matrix evaluated at the DFE point is given as J = − b α − β αA − β αS − ( σ α + µ α ) β αA β αS ησ α − ( γ α + µ α ) 0 0 00 (1 − η ) σ α − ( γ α + µ α ) 0 00 0 γ α γ α − ( ρ α + µ α ) 00 0 0 0 ρ α − b α . The eigenvalues of the Jacobian matrix are given as λ = λ = − b α , λ = − ( ρ α + µ α ) and the roots of the equation z + Az + Bz + C , where A = ( σ α + γ α + γ α + µ α + µ α + µ α ), B = − σ α ( ηβ αA + (1 − η ) β αS ) + ( σ α + µ α )( γ α + µ α ) + ( σ α + µ α )( γ α + µ α ) + ( γ α + µ α )( γ α + µ α ), C = − σ α ( η ( γ α + µ α ) β αA + (1 − η )( γ α + µ α ) β αS ) + ( σ α + µ α )( γ α + µ α )( γ α + µ α ).To show stability, we apply the Routh-Hurwitz criterion. Clearly A >
C > R <
1. It is easy to show that
AB > C if K <
Appendix C. Proof of Theorem 3.3
Proof.
The Jacobian matrix evaluated at the EE point is given as J = − b α − b α (cid:18) − R (cid:19) − β αA R − β αS R b α (cid:18) − R (cid:19) − ( σ α + µ α ) β αA R β αS R ησ α − ( γ α + µ α ) 0 0 00 (1 − η ) σ α − ( γ α + µ α ) 0 00 0 γ α γ α − ( ρ α + µ α ) 00 0 0 0 ρ α − b α . λ = − b α , λ = − ( ρ α + µ α ) and theroots of the equation z + A z + A z + A z + A . By Applying the Routh-Hurwitz criterion, the EE is stable if A > A > A A − A ) A − A A ≥
0. Clearly A > R >
1. Also, A > b α R (2 R − σ α + µ α )( γ α + µ α )( γ α + µ α ) > σ α ((1 − η ) β αS ( γ α + µ α )+ ηβ αA ( γ α + µ α ))] lpha ( γ α + µ α )) ⇒ b α R (2 R − > σ α b α ((1 − η ) β αS ( γ α + µ α ) + ηβ αA ( γ α + µ α ))( σ α + µ α )( γ α + µ α )( γ α + µ α ) ⇒ b α R (2 R − > R > Appendix D. Data Sets and Methods
Data for cumulative confirmed cases, recovered and deaths were obtainedfrom the Tennessee Department of Health (TDH) and Johns Hopkins University(JHU) Center for Systems Science and Engineering which are made available tothe public on TDH’s website [48] and Github [55], respectively. The raw filesare converted into Panda data frames and stored for ease of access. For the filesobtained from [55], the data are sorted according to counties, so we aggregatethe cases for each of the recorded counties to obtain the total number of casesin each day for the states. For the total number of population for each states,we used the data obtained from the US Census Bureau [57].All numerical simulations were done in python using our numerical scheme [58]from which we obtain the solution of the proposed model at each time step as1. Predictor: S p = S j + τ α Γ(1 + α ) F ( t j , S j , E j , I A,j , I
S,j , H j , R j , D j ) + ˜ H ,j E p = E j + τ α Γ(1 + α ) F ( t j , S j , E j , I A,j , I
S,j , H j , R j , D j ) + ˜ H ,j I A,p = I A,j + τ α Γ(1 + α ) F ( t j , S j , E j , I A,j , I
S,j , H j , R j , D j ) + ˜ H ,j I S,p = I S,j + τ α Γ(1 + α ) F ( t j , S j , E j , I A,j , I
S,j , H j , R j , D j ) + ˜ H ,j H p = H j + τ α Γ(1 + α ) F ( t j , S j , E j , I A,j , I
S,j , H j , R j , D j ) + ˜ H ,j R p = R j + τ α Γ(1 + α ) F ( t j , S j , E j , I A,j , I
S,j , H j , R j , D j ) + ˜ H ,j D p = D j + τ α Γ(1 + α ) F ( t j , S j , E j , I A,j , I
S,j , H j , R j , D j ) + ˜ H ,j
2. Corrector: 32 j +1 = S j + τ α Γ(2 + α ) (cid:16) α F ( t j , S j , E j , I A,j , I
S,j , H j , R j , D j )+ F ( t j +1 , S p , E p , I A,p , I
S,p , H p , R p , D p ) (cid:17) + ˜ H ,j ,E j +1 = E j + τ α Γ(2 + α ) (cid:16) α F ( t j , S j , E j , I A,j , I
S,j , H j , R j , D j )+ F ( t j +1 , S p , E p , I A,p , I
S,p , H p , R p , D p ) (cid:17) + ˜ H ,j ,I A,j +1 = I A,j + τ α Γ(2 + α ) (cid:16) α F ( t j , S j , E j , I A,j , I
S,j , H j , R j , D j )+ F ( t j +1 , S p , E p , I A,p , I
S,p , H p , R p , D p ) (cid:17) + ˜ H ,j ,I S,j +1 = I S,j + τ α Γ(2 + α ) (cid:16) α F ( t j , S j , E j , I A,j , I
S,j , H j , R j , D j )+ F ( t j +1 , S p , E p , I A,p , I
S,p , H p , R p , D p ) (cid:17) + ˜ H ,j ,H j +1 = H j + τ α Γ(2 + α ) (cid:16) α F ( t j , S j , E j , I A,j , I
S,j , H j , R j , D j )+ F ( t j +1 , S p , E p , I A,p , I
S,p , H p , R p , D p ) (cid:17) + ˜ H ,j ,R j +1 = R j + τ α Γ(2 + α ) (cid:16) α F ( t j , S j , E j , I A,j , I
S,j , H j , R j , D j )+ F ( t j +1 , S p , E p , I A,p , I
S,p , H p , R p , D p ) (cid:17) + ˜ H ,j ,D j +1 = D j + τ α Γ(2 + α ) (cid:16) α F ( t j , S j , E j , I A,j , I
S,j , H j , R j , D j )+ F ( t j +1 , S p , E p , I A,p , I
S,p , H p , R p , D p ) (cid:17) + ˜ H ,j , where F ( t j , S j , E j , I A,j , I
S,j , H j , R j , D j ) = b α N − S j N ( β αA I A,j + β αS I S,j ) − b α S j ,F ( t j , S j , E j , I A,j , I
S,j , H j , R j , D j ) = S j N ( β αA I A,j + β αS I S,j ) − ( σ α + µ α ) E j ,F ( t j , S j , E j , I A,j , I
S,j , H j , R j , D j ) = ησ α E j − ( γ α + µ α ) I A,j ,F ( t j , S j , E j , I A,j , I
S,j , H j , R j , D j ) = (1 − η ) σ α E j − ( γ α + µ α ) I S,j ,F ( t j , S j , E j , I A,j , I
S,j , H j , R j , D j ) = γ α I A,j + γ α I S,j − ( ρ α + µ α ) H j ,F ( t j , S j , E j , I A,j , I
S,j , H j , R j , D j ) = ρ α H j − b α R j ,F ( t j , S j , E j , I A,j , I
S,j , H j , R j , D j ) = µ α E j + µ α I A,j + µ α I S,j + µ α H j , H ,j = τ α Γ(2 + α ) j (cid:88) l =0 a l,j F ( t l , S l , E l , I A,l , I
S,l , H l , R l , D l ) , ˜ H ,j = τ α Γ(2 + α ) j (cid:88) l =0 a l,j F ( t l , S l , E l , I A,l , I
S,l , H l , R l , D l ) , ˜ H ,j = τ α Γ(2 + α ) j (cid:88) l =0 a l,j F ( t l , S l , E l , I A,l , I
S,l , H l , R l , D l ) , ˜ H ,j = τ α Γ(2 + α ) j (cid:88) l =0 a l,j F ( t l , S l , E l , I A,l , I
S,l , H l , R l , D l ) , ˜ H ,j = τ α Γ(2 + α ) j (cid:88) l =0 a l,j F ( t l , S l , E l , I A,l , I
S,l , H l , R l , D l ) , ˜ H ,j = τ α Γ(2 + α ) j (cid:88) l =0 a l,j F ( t l , S l , E l , I A,l , I
S,l , H l , R l , D l ) , ˜ H ,j = τ α Γ(2 + α ) j (cid:88) l =0 a l,j F ( t l , S l , E l , I A,l , I
S,l , H l , R l , D l )are the memory terms of the respective population variables and a l,j = τ γ Γ( γ + 2) − ( j − γ )( j + 1) γ + j γ (2 j − γ − − ( j − γ +1 , l = 0 , ( j − l + 2) γ +1 − j − l + 1) γ +1 + 3( j − l ) γ +1 − ( j − l − γ +1 , ≤ l ≤ j − , γ +1 − γ − , l = j. The fitted parameters were obtained by using python’s scipy.optimize.minimizescipy.optimize.minimize