Dynamic model of COVID-19 and citizens reaction using fractional derivative

This paper presents dynamic model of COVID-19 and citizens reaction from a fraction of the population in Nigeria using fractional derivative. We consider the reported cases from February to June 2020 using fractional derivative. The Stability analysis of the model was carried out and the basic reproduction number was calculated via the next generation matrix. The fractional derivative model was solved numerically and many graphs presented accordingly which could serve as a yard-stick of reducing this menacing virus and policy making.


Introduction
The present Corona Virus pandemic ravaging almost the entire nations of the earth probably caught the world off guard. The events as seen globally have again revealed how exponentially a new disease can 'sprout' and spread. Meanwhile, these events come with an associated explosion of clinical and epidemiological information worthy of research.
Disease outbreaks (epidemic) is not something uncommon and there is nothing to suggest there won't be more outbreaks in the foreseeable future. Severe acute respiratory syndrome corona virus (SARS-CoV) broke out in 2003 while the world again witnessed another form of corona virus epidemic known as the Middle Eastern Respiratory Syndrome corona virus (MERS-CoV) in 2012 with the first case reported in Saudi Arabia [13]. A total of 8,098 SARS cases were reported with a fatality rate of 11 percent in 17 countries, the largest fraction of case reported coming from southern mainland China and Hong Kong [34]. However, exactly 6 months since China reported the world's first case of the new disease, mortality rate of corona virus disease 2019 (COVID-19) has largely surpassed the two previous corona virus outbreaks within the last 2 decades with about 503,682 deaths globally as at the time (June 30, 2020) of writing this paper (see [33]). According to www.-covid19.ncdc.gov.ng [32], there are about 25,694 reported cases with 590 deaths reported in Nigeria as at June 30, 2020.
The emergence of COVID-19 coincided with the largest annual human migration in the world from Chinese new year festival according to Lin et al. [27]. It is stated that the rapid national and global spread of this deadly virus was as a result of large travel for spring festival season in China. The first case of corona virus was reported in the city of Wuham, China on 31st December, 2019 [16]. Since the first reported case, the virus has spread all over America, Europe, Africa and the world at large (CNN). Treatments or rather vaccines such as Hydroxychloroquine which is a class of drug for antimalarial, treatments for rheumatoid arthritis and systemic lupus erythematosus have been used but show no effect (WHO).
The first case of the corona virus disease was first confirmed in Nigeria on 28th February, 2020 [32]. It was stated that the symptoms include coughing incessantly, difficulties in breathing, high fever, pharyngitis and the appearance of the symptoms takes between 2 to 14 days. Since the first confirmed case in Nigeria, government has put in place strong measures to stop the spread of the virus among the population. These measures include ban on inter-state and international travels, restriction on mass gathering (social distancing), mandatory use of nose mask, total lockdown of markets, worship centres, schools and regular use of hand sanitizer after washing of hands.
These measures will substantially reduce the spread of this virus among the communities in this country as evidence from other countries like Italy, USA, China etc. have shown. However, the reality is that some of these measures have been observed to heighten the plight of approximately 48% of population ravaged by extreme poverty, marginalized by basic amenities, economy, health facilities etc. Most of these citizens are displaced because of continuous ethnic cleansing still taking place in their communities. Some also reside in a makeshift camp while some are even homeless. The indication shows that these group of population depends on daily hand-outs and small sales on daily basis to make ends meet. In fact, their survival solely depends on most of the things that the Nigerian government has placed restriction on. How on earth can this fraction of the population survive this present time if they abide by these measures? Our aim in this paper is to attempt to present a reasonable mathematical analysis of the spread of corona virus in Nigeria incorporated with the reaction from the fraction of the population.
In our study, we consider the dynamic model as a fractional derivative by applying Caputo fractional derivative of order q. The use of Caputo fractional derivative is because its data presentation has a well understood physical meaning compared to other fractional derivatives such as Riemann-Liouville. Fractional derivatives have played an important role in modelling physical problems. Though old and yet a novel topic, the first monograph is credited to [15]. They claimed that this branch of mathematics is applicable to transmission line theory, chemical analysis of aqueous solutions, design of heat-flux meters, rheology of soils, growth of intergranular grooves at metal surfaces, quantum mechanical calculations and dissemination of atmospheric pollutants. Other basic works in various aspect of fractional calculus can be found in [6,10,17,19,28,29,31,40], extensive approach by [39], fractional analysis by [4,8] and biological applications [20,[23][24][25] with references therein. For numerical analysis of fractional order in Caputo sense, see [9,11,14]. The convergence and accuracy of the numerical method were discussed by [12]. Mathematical modelling and analysis of infectious disease dynamics with more realistic assumptions were considered by [3,38].
Nowadays, many applications of fractional derivatives in epidemiology have exponentially increased. Atangana-Beleanu fractional derivative was used by [16] to model the dynamics of the novel coronavirus (2019-nCov). It was mentioned in their paper, that the use of this particular fractional derivatives is due to its unique properties such as their kernel is nonlocal and nonsingular and would best describe their model properly. Also Atangana and Owolabi [21] developed 'new numerical approach for fractional differential equations' for mathematical modelling of natural phenemenon. Ndolane [18] applied fractional derivative in the modelling of SIR epidemic with Mittag-Leffler fractional derivative where stability analysis of their model were discussed. Analysis of the financial chaotic model with the aid of fractional order was also produced by [22]. Mouaouine et al. [26] used Caputo fractional derivative to model SIR epidemic with nonlinear incidence rate. Also Caputo fractional derivative was used by [2] to analyze the SEIR model with fractional derivative using vertical transmission in a population with density-dependent death rate.
The rest of the paper is organized as follows: Formulation of the fractional dynamical model which includes preliminaries of fractional calculus is presented in Section Two. Section Three captures the systematic analysis of the dynamic model where the disease-free equilibrium state, basic reproduction number, existence of the endemic equilibrium state and stability analysis of the model based on its equilibria. In section Four, the fractional predictor-corrector algorithm which is a numerical method that will aid on the analysis of the fractional model is presented. Section Five consists of data fitting and numerical simulations of the model. Also, discussion and concluding remarks were presented.

Formulation of dynamic model
In this section, we seek the population size, N be divided into a set of distinct compartments or classes -Susceptible, Exposed, Infectious, Removed ðSEIRÞ framework with one extra compartment P which represents fraction of the population deprived and marginalized by wealth, basic amenities, security, health care etc. This particular population are expected by the government to adhere to the strict measures mentioned above in introduction irrespective of their struggles and difficult circumstances. The dire situation of this particular population shows that it will be extremely difficult for someone without proper financial backing, steady water supply, health care facilities etc. to follow this strict rules. Evidently, they will surely try to survive by any means possible and hence reacting negatively and breaking the necessary rules put in place to minimize the spread of the virus. Thus, it is also the intention of this paper to demonstrate that this particular population could rapidly spread this virus if something substantial is not done to address their needs.
We note that the compartment in which an individual resides at time t depends on that individual experience with respect to the virus. With respect to the SEIR framework and the extra compartment P, we define the state variables accordingly. SðtÞ = number of susceptible individuals at time t. EðtÞ = number of exposed individuals at time t. IðtÞ = number of infected individuals at time t. RðtÞ = number of post-infective individuals removed from the population (due to immunity, quarantine, or death) at time t. P ðtÞ = number of marginalized and deprived individuals at time t.
Birth rate, m > 0 and death rate, l > 0 occurs at equal rate and all the new born (immigration) are susceptible. Susceptible individuals are able to contact the virus since any encounter with the infected will result in transmission through the quantity HSI N (rate of new infection), where H is the force of infection (the infection rate). After the susceptible made contact with the infected, they move into the exposed compartment (those infected but not yet infectious) for a period of time. At the end of the incubation period, the transition rate from E compartment to I compartment is proportional to E, that is E where parameter is the progression rate to symptoms development. The parameterdenotes the death rate caused by the virus. At the end of the infectious period, the individuals from the infectious compartment I enters the R compartment at the removal rate denoted by u. The individuals on the removed compartment are made up of those who are unable to acquire or spread the virus since they are either immune, quarantined or dead. The compartment P mimics fraction of the population reacting negatively because of the deprivation and marginalization they suffer. This particular population is also susceptible to the virus but their chances of being infected are very high. The quantity wSP denotes the susceptible individuals infected after communicating with P compartment, where parameter w is the higher rate of virus transmission. The parameter ! denotes the proportion of the fraction of population with severe conditions while parameter / represents the duration of the reaction from this population. The quantity 1=/ represents the average length of the reaction.
The dynamic model (referred to as the SEIRP model) describing the temporal evolution of the sizes of the compartment is based on the following assumptions: 1. The population of the susceptible class will keep increasing if there is no contact with the infectious class. Also, the population of the susceptible class will decrease at a rate proportional to its current population if there is a contact with the infectious class. 2. The rate at which the susceptible becomes infected is proportional to the number of encounters between susceptible and infected individuals. This in turn is proportional to the product of the two population, HSI. Suppose we have larger values of H, it will correspond to the higher contact rates between infected and susceptible. In the same way, smaller values of H will correspond to lower contact rates between the two compartments. 3. The rate of transition from class I to class R is proportional to I, which is uI where u is that 1=u is the average length of the infectious period. 4. We also assume the population is completely homogeneous and there are deaths caused by the virus.
The simple truth is that COVID-19 is still rampaging globally and particularly in Nigeria where gradual increase of this menacing virus are recorded on daily basis. With the above assumptions, the fractional differential equations, equipped with Caputo fractional derivative that describes the number of individuals in five compartments including the compartment P representing the reaction of deprived and marginalized individuals are  (1) is subject to the following initial conditions respectively.
The papers of K. Diethelm and N.J Ford [8] have shown that our fractional dynamic model (1) of order 0 < q 1 with its respective initial conditions (3) is equivalent to the following Volterra-integral equations of second kind CðqÞ R t 0 ðt À sÞ qÀ1 F s; SðsÞ ð Þds; EðtÞ CðqÞ R t 0 ðt À sÞ qÀ1 G s; EðsÞ ð Þ ds; RðtÞ We assume that the functions in our fractional dynamic model (1) are continuous with respect to both of its argument in the domain of integration and it also satisfies the Lipschitz condition with respect to its second argument. However, the use of Caputo fractional derivative as stated in [7,8] is based on the fact that D q Ã has m-dimensional kernel and certainly there is need to specify the m initial conditions in order to obtain a unique solution of our fractional dynamic model (1). According to Kumar and Agrawal [7], Riemann-Liouville based fractional differential equations require fractional derivatives and integrals of the unknown function at t ¼ 0. In other words, the physical meaning of this type of conditions are unknown and its practical application appears impossible.

Preliminaries and fundamentals
We state the fundamental definitions of fractional calculus required for our model. We present the boundedness and the positivity of our solutions which shows that system (1) is both mathematically and biologically well-posed. Definition 2.1 [5]. The fractional integral of order q > 0 of a function f : R þ ! R of order q 2 R þ is defined by provided the integral part is integrable in R þ and the symbol C is the Gamma function defined by For simplicity, where fðtÞ ¼ 1, : where m is a positive integer uniquely defined by m À 1 < q 6 m. For a special case where 0 < q 6 1 which is a particular interest to this paper, the above Caputo fractional derivative of order q is reduced to Definition 2.3 [4]. The Mittag-Leffler of two-parameter type function is defined by the series expansion It follows from (6) that which is the classical exponential function.
Definition 2.4 [5]. The Caputo fractional derivative (5) treated by Laplace transform technique is given by where L is the Laplace transform operator. Note that when q ¼ 1, Eq. (7) reduces to the natural generalization of the corresponding classical result.
is positively invariant with respect to the governing model, 2. every solution of system (1) starting from the initial point S 0 ; E 0 ; I 0 ; R 0 and P 0 remain positive for all t P 0.
Proof. Adopting similar appproach from [1,18], let us consider the following expression which is the acclaimed human population at the time t. Applying Caputo fractional derivative into (9) and substituting the equivalent expressions from system (1) we get Applying Laplace transform in (10) Using partial fraction, we re-write (11) as follows 6 m l À m l À Nð0Þ E q Àlt q ð Þ: It follows that as t ! 1 which is the condition under which Eq. (1) except the last equation are bounded and shows that feasible region exist. h Positivity of the solution: Following similar approach from [1], we assume by contradiction that the second equation in fractional dynamic model is not true. Let t Ã ¼ minft : SðtÞEðtÞIðtÞPðtÞ ¼ 0g. Now if Eðt Ã Þ ¼ 0, it then implies that SðtÞ P 0; IðtÞ P 0; PðtÞ P 0 for all 0 6 t 6 t Ã .
Let us consider the following expression and assume it exists, It follows that D q Ã EðtÞ À B 1 EðtÞ P 0: Applying Laplace transform we get s qẼ ðsÞ À s qÀ1 Eð0Þ À B 1Ẽ ðsÞ P 0; so that EðsÞ P Eð0Þ s qÀ1 s q À B 1 ¼ Eð0Þ Applying inverse Laplace transform and expressed as a Mittag-Leffler function of one parameter, we get EðtÞ P Eð0Þ where the solution of (12) is satisfied by (13). Thus, the positivity of the solution E is represented by Similarly, suppose Sðt Ã Þ ¼ 0 which also implies that EðtÞ P 0; IðtÞ P 0; PðtÞ P 0 for all 0 6 t 6 t Ã . Let Applying Laplace transform we get s qS ðsÞ À s qÀ1 Sð0Þ À B 2S ðsÞ P 0; so that SðsÞ P Sð0Þ s qÀ1 s q À B 2 ¼ Sð0Þ Applying inverse Laplace transform and expressed as a Mittag-Leffler function of one parameter, we get SðtÞ P Sð0Þ where the solution of (14) is satisfied by (15). Thus, the positivity of the solution S is represented by which contradicts Sðt Ã Þ ¼ 0.
Following the same process, suppose Iðt Ã Þ ¼ 0 which also implies that EðtÞ P 0; SðtÞ P 0; PðtÞ P 0 for all 0 6 t 6 t Ã . Let such that D q Ã IðtÞ À B 3 IðtÞ P 0: Applying Laplace transform we get s qĨ ðsÞ À s qÀ1 Ið0Þ À B 3Ĩ ðsÞ P 0; so that IðsÞ P Sð0Þ s qÀ1 s q À B 3 ¼ Ið0Þ Applying inverse Laplace transform and expressed as a Mittag-Leffler function of one parameter, we get where the solution of (16) is satisfied by (17). Thus, the positivity of the solution I is represented by which contradicts Iðt Ã Þ ¼ 0. Also, suppose Pðt Ã Þ ¼ 0 which also implies that EðtÞ P 0; SðtÞ P 0; IðtÞ P 0 for all 0 6 t 6 t Ã . Let Applying Laplace transform we get s qP ðsÞ À s qÀ1 Ið0Þ À B 4P ðsÞ P 0; from which PðsÞ P Sð0Þ s qÀ1 s q À B 4 ¼ Pð0Þ Applying inverse Laplace transform and expressed as a Mittag-Leffler function of one parameter, we get PðtÞ P Pð0Þ where the solution of (18) is satisfied by (19). Thus, the positivity of the solution P is represented by which contradicts Pðt Ã Þ ¼ 0. Moreover, from the fourth equation in system (1), one can see that D q Ã RðtÞ P ÀlRðtÞ: ð20Þ Using similar approach from above, we have the following relation RðtÞ P Rð0ÞE q Àlt q ð Þ; from which the positivity of the solution R is represented by RðtÞ P Rð0ÞE q Àlt q ð ÞP 0: We therefore conclude that any solution of SðtÞ; EðtÞ; IðtÞ; RðtÞ; PðtÞ ð Þ is positive. Hence since the solution is bounded and positive, we can reasonably say that X is positively invariant with respect to system (1).

Systematic analysis of the dynamic model
Here, we give a vivid analysis of our model (1) with respect to existence of disease-free equilibrium state (W 0 ), basic reproduction number (R 0 ), which is the expected number of secondary cases generated by a single typical index case in a completely susceptible population. We also discuss existence of endemic equilibrium state and stability analysis of the model based on its equilibria.

Disease-free equilibrium state
Disease-free equilibrium state, W 0 simply means the point where there is total absence of COVID-19, which is represented by E ¼ I ¼ 0.
Lemma 1. The disease-free state of system (1) exists at the point Proof. The rate of change of each variable at equilibrium state is given by It then follows from system (1) that ð22Þ uI Ã À lR Ã ¼ 0; u!I Ã À /P Ã ¼ 0: is the disease-free equilibrium state of system (1). h

Basic reproduction number, R 0
Using the next generation matrix or the spectral radius method obtained from [37], the basic reproduction number R 0 ¼ qðFV À1 Þ, which is known as the spectral radius, is the largest eigenvalues of the matrix jFV À1 j. Here matrix F is the rate of new infections in different compartment differentiated with respect to E; I and P and evaluated at W 0 which is the disease free equilibrium (21). Matrix V defines the rate of transfer of infectives from one compartment to another. Therefore, the Jacobian matrix F evaluated at the disease equilibrium and the Jacobian matix V are ; where the inverse matrix of V is Thus we obtain k 1 ¼ k 2 ¼ 0 and k 3 ¼ ð/Hþwm!uÞ /ðlþÞðuþ-þlÞ , from which the basic reproduction number of our model (1) is

Existence of the endemic equilibrium state
Following similar approach from [18], we prove the existence and uniqueness of endemic equilibrium state. We will also apply Intermediate Value Theorem which plays a big role in many dynamic model analysis.
Proof. Let the rate of change in system (1) be such that W Ã ¼ ðS Ã ; E Ã ; I Ã ; R Ã ; P Ã Þ is the non-trivial solution of Eq. (24). Using third to last equations in system (1), we have where we have plugged (25) into (26) and (27) respectively. In another note For the existence of E Ã , we assume that g is a continuous function defined as g : R þ ! R, then after a simple algebraic manipulation, which is positive if R 0 > 1. Thus from the intermediate value theorem, there exists a unique number E Ã such that ð0 6 E Ã 6 m l þ and gðE Ã Þ ¼ 0: The above condition shows that solutions I Ã ; R Ã ; P Ã ; S Ã defined respectively by Eqs. (25)-(28) exists and unique. Hence W Ã ¼ ðS Ã ; E Ã ; I Ã ; R Ã ; P Ã Þ is the unique endemic equilibrium state of system (1). h

Stability analysis of the equilibria
Here, we consider the stabilities of the disease-free equilibrium W 0 and the endemic equilibrium state W Ã of our fractional dynamic model (1).
Theorem 2. The disease-free equilibrium W 0 is locally asymptotically stable if R 0 < 1 and unstable if R 0 > 1.
Proof. The Jacobian matrix (J) evaluated at the disease-free equilibrium state, W 0 is given by differentiating the independent variables in (2) with respect to S; E; I; R and P, so that Proof. To prove this theorem, we will exclude the fourth equation from system (1), since R has no impact in the rest of the fractional dynamic system and following the method from [36], we consider the following Lyapunov function Applying Caputo fractional derivatives of order q we have where b ¼ lþ . Now substituting the equivalent expressions of the fractional derivatives from (1) into the above equation and using the fact that m ¼ Let us note that This is due to the fact that Hence by Lassalles's invariance principle, our result follows. h

Numerical analysis
In this section, we consider the predictor-corrector method for the Caputo fractional derivative proposed by [12] and modified by [35] with the aim of reducing computational cost. The scheme which was derived from a well known Adams-Bashforth linear multi-step method for solving classical differential equations is very remarkable since it represents a good compromise between accuracy and ease of implementation. Extensive treatment of convergence and accuracy of the predictorcorrector method can be seen in [12].
Let 0 ¼ t 0 < t 1 < t 2 < . . . < t n < . . . ¼ t be the partition of t > 0. Let h be the step-size on the grid where the equispaced grid, t n ¼ t 0 þ nh. By Lagrange interpolation, the fractional predictor-corrector numerical algorithm for each of the fractional models, (1) and the weights are given as follows: where F j ¼ Fðt j ; S j Þ and S n is the numerical approximation to the exact solution Sðt n Þ.
where G j ¼ Gðt j ; S j Þ and E n is the numerical approximation to the exact solution Eðt n Þ.
where H j ¼ Hðt j ; S j Þ and I n is the numerical approximation to the exact solution Iðt n Þ.
where K j ¼ Kðt j ; S j Þ and R n is the numerical approximation to the exact solution Rðt n Þ.
where Q j ¼ Qðt j ; S j Þ and P n is the numerical approximation to the exact solution Pðt n Þ.
The fractional weights x n ; a n;0 and b n for each the numerical solutions to the model is given by ; n ¼ 1; 2; . . . ; 8 < : a n;0 ¼ ðnÀ1Þ qþ1 Àn q ðnÀqÀ1Þ

Data fitting and numerical simulations
In this section we present all numerical simulations produced using the fractional predictor-corrector method. According to [30], the current population of Nigeria as of June is 205,740,406 which is equivalent to 2:64% of the total world population. The total initial population N ¼ 106; 985; 011 which we consider to be 52% of the total population of Nigeria. This is because the 52% of the population is urban. The life expectancy for Nigeria in 2020 according to World Health Organisation (W.H.O.) is 54.81 years. It is quite reasonable to suggest that this is due to the sub standard of living and basic needs. Hence we estimate the natural death rate l ¼ 1 54:81 . The birth rate is estimated as m ¼ l Â N, which we consider to be limited population in the absence of any infection. It is also noted that there are approximately 10,578 confirmed cases of COVID-19, 3122 discharged and 299 deaths confirmed.
The parameters are listed as follows: (see Table 1). The average monthly cases from February-June, 2020 is shown in the Fig. 1. The total number of reported cases from February-June, 2020 is displayed in the Fig. 2. The plots for SEIRP compartments with different values of q and initial values were displayed in the Figs. 3-12.

Discussion of results and concluding remarks
The various and different simulations of each compartments produces distinctive and robust dynamics of the epidemic which is relatively in agreement with [16]. Application of non-integer order into the model has led to the manufacturing of visible effect on the dynamics of our model. The results generated via SðtÞ; EðtÞ; IðtÞ; RðtÞ and PðtÞ for different initial values were captured in the Figs. 3, 5, 7, 9 and 11, respectively. Figs. 4, 6, 8, 10 and 12 show the robust behaviour of each compartment made up of the susceptible population SðtÞ, the exposed population RðtÞ, the infective population IðtÞ, the removed population RðtÞ which consist of the recovered, dead and the fraction of the population reacting negatively PðtÞ against the time t in days ð0 6 t 6 10Þ for distinct values of q; ð0 < q 6 1Þ.
The rapid decrease in susceptible compartment which converges to zero as the values of q increase is evident as seen in Fig. 4. This is expected because of the interaction of the population in this compartment with the infected population. The exposed compartment which shows as there is decrease in the values of q, then the rate of rapid increase will decline as shown in Fig. 6. It is observed from Fig. 8 that the infective compart-    Dynamic model of COVID-19 and citizens reaction ment shows a slow growth initially but accelerated with time as the values of q increases. It is also observed from the Fig. 10 that, as the values of q decreases there is a sharp increase in the removed compartment. It is seen after reaching a certain      peak, it begins to flatten as time increases. The similarity strike between the susceptible compartment and the reaction compartment is not just a coincident but as expected. In fact this compartment shows a rapid decrease in the population and converges to zero with time as the values of q increases as shown in the Fig. 12.
As already witnessed, we have proposed a fractional dynamic model for the spread of COVID-19 and the reaction of some citizens. We have extensively discussed the positivity of the solution which depend on the fractional derivative order, stability of our model with respect to its equilibria etc. In no way do we claim that our model is perfect instead our model provides policy makers a yard-stick to a better solution to the war against COVID-19. Increase in contact tracing, strict lockdown and providing the necessary needs for the population reacting negatively due to their difficult circumstances can lead to decrease in the numbers of infection of COVID-19. It will also create a real avenue for the population to adhere to the governmental guidelines and thereby leading to exact estimates of the parameters and the progress of the virus.
The numerical simulation was performed to witness the interaction of the COVID-19 against the population (five different compartments) in Nigeria. The parameters applied were mostly either fitted, estimated or obtained from already existing literature. This is due to the limited access to real data. Since this is the case, the model is certain to fluctuate thereby creating uncertainties in policy making and estimates.