We propose a mathematical model of the coronavirus disease 2019 (COVID-19) to investigate the transmission and control mechanism of the disease in the community of Nigeria. Using stability theory of differential equations, the qualitative behavior of model is studied. The pandemic indicator represented by basic reproductive number R0 is obtained from the largest eigenvalue of the next-generation matrix. Local as well as global asymptotic stability conditions for the disease-free and pandemic equilibrium are obtained which determines the conditions to stabilize the exponential spread of the disease. Further, we examined this model by using Atangana–Baleanu fractional derivative operator and existence criteria of solution for the operator is established. We consider the data of reported infection cases from April 1, 2020, till April 30, 2020, and parameterized the model. We have used one of the reliable and efficient method known as iterative Laplace transform to obtain numerical simulations. The impacts of various biological parameters on transmission dynamics of COVID-19 is examined. These results are based on different values of the fractional parameter and serve as a control parameter to identify the significant strategies for the control of the disease. In the end, the obtained results are demonstrated graphically to justify our theoretical findings.

Mathematical model COVID-19 Atangana-Baleanu fractional operator existence of solutions stability analysis numerical simulation
Introduction

The ongoing ravaging COVID-19 is a contagious disease instigated by SARS-CoV-2. The first case of the disease was reported in December 2019 in Wuhan, China, and has, within few weeks, spread across the globe, leading to the present 2020 COVID-19 pandemic . The coronavirus disease 2019 has been regarded as the largest global health crisis in human history as a result of the magnitude of confirmed cases, accompanied by the degree of fatalities across the continents . Reliable data had it that by April 2020, COVID-19 pandemic had led to over 3 million confirmed cases with 230,000 deaths and the disease has spread to over 210 nations globally . The symptoms and signs of COVID-19 develop within 2 to 14 days . When the disease is fully incubated, the infected individuals may exhibit fever, fatigue, cough and breathing disorder that is similar to those infections instigated by SARS-CoV and MERS-CoV . However, many COVID-19 acute cases and fatalities come from the elderly people (from the age of 65 upward) and individuals with severe health challenges (such as people with kidney disease, hypertension, diabetes, obesity and other health issues that deteriorate the immune system) .

The first confirmed case in Nigeria was reported on 27 February 2020, when a citizen from Italy is tested positive for the virus . The disease transmission raised gradually over the month of April 2020, after a substantial number of cases were noted for in the country. Since then, Nigeria is focused on spotting and referring identified infectious patients for treatment to devoted COVID-19 centers. As of 18 June 2020, Nigeria had reported 17,735 confirmed cases out of which 11,299 are active cases, 5,967 recovered peoples and 469 deaths due to COVID-19 infection .

The global scourge of COVID-19 pandemic has elicited the attention of scholars in different disciplines, prompting several proposals to examine and envisage the development of the pandemic . Ndairov et al.  propose a model for the transmissibility of COVID-19 in the presence of super-spreaders individuals. They perform the stability and sensitivity analyses of the model and discovered that daily reduction in the number of confirmed cases of COVID-19 is a function of the number of hospitalizations. Yang et al.  proposed a model to study the transmission pathways of COVID-19 in terms of human-to-human and environment-to-human spread. Their analysis confirms the tendency of COVID-19 to remain pandemic even with prevention and intervention measures.

A model for the dynamics of COVID-19 with parameter estimations, sensitivity analysis and data fitting is investigated in  while a model for COVID-19 infection that describes the impact of slow diagnosis on the dynamics of COVID-19 is also studied in . In , the researchers employ a statistical study of coronavirus disease data to calculate time-regulated risk for fatality from the COVID-19 in Wuhan. Their results indicate that movement restrictions and adequate social distancing procedures are capable of reducing the spread of the disease. Furthermore, a data-oriented model that includes behavioral impacts of humans and governmental efforts on the dynamics of COVID-19 in Wuhan is proposed in . A good number of mathematics and non-mathematics studies have also been conducted on COVID-19 .

In recent times, the integer order differential systems are generalized and improved in order to formulate several mathematical models using fractional differential operators. Since, demonstration of some real-world phenomena with the help of fractional derivative operator is more appropriate and useful for improving performance of numerous engineering and applied sciences systems . In this paper, initially we formulate the mathematical model in terms of integer order derivative and then apply the Atangana–Baleanu fractional derivative operator. The motivation behind utilizing the Atangana–Baleanu operator is that it has nonlocal and nonsingular kernel in the form of Mittag-Leffler function. Moreover, the complex behavior in the model must be ideal portrayed utilizing this operator. Literature pertaining to Atangana–Baleanu derivative and their applications to several systems arising in the field of applied sciences and engineering can be cited in .

Formulation of the Model

The population of human under consideration is divided into six compartments which are susceptible S(t), since the incubation period of COVID-19 is between two to fourteen days there are those who are infected without exhibiting any sign of symptoms and are undetected E(t). Individuals who are infected or suspected case of COVID-19 need to go through an incubation period before the suspected symptom can noticeable these categories are quarantined Q(t), there are those certain proportion of the population have been infected with sign and symptoms of COVID-19 and is highly infectious but not yet quarantined or isolated I(t). C(t) represent confirmed case of COVID-19 from Quarantine category. R(t) represent recovery after treatment. The susceptible population is increased by immigration or by birth at the rate θ , in each of the class, individuals can die a natural death and the rate θ , there is a force of infection between the susceptible population and exposed population this is represented by β1 , ε represent the progression from exposed class to highly infected class, the disease induced death rate for highly infected class, quarantine class and confirmed case of COVID-19 class is represented by δ , proportion of people identified as suspected case of COVID-19 are represented by β2 , after medical diagnosis, some of the suspected cases were confirmed, others that are not detected can return back to the susceptible population at the rate α . In the meantime, some highly infectious individuals will be moved to quarantine class at the rate γ . The progression rate from quarantine to confirm case after diagnosis is denoted by τ . The pictorial diagram illustrating the model is shown in Fig. 1 while the system of equations governing the model is given as:

The model’s flow diagram

dSdt=θβ1SIμS+αQ,

dEdt=β1SI(μ+β2+ε)E,

dIdt=εE(γ+μ+δ)I,

dQdt=γI+β2E(α+μ+δ+τ)Q,

dCdt=τQ(μ+δ+φ)C,

dRdt=φCμR.

The rest of the sections are organized as follows: Some valuable preliminaries dependent on the Atangana–Baleanu fractional operator is given in Section 3. In Section 4, we present stability analysis of the equilibria (drug-free equilibrium state and endemic equilibrium state). Analysis of fractional coronavirus model using the Atangana–Baleanu operator is given in Section 5. The Approximation technique and Numerical Simulation are given to reveal the behavior of dynamics components is accounted for in Section 6. The conclusion is finally drawn in the last Section 7.

Preliminaries

This section of the paper will convert some basic definitions and properties related to Fractional calculus. During the paper process, we are going to refer to the following given specific definitions and properties of the Atangana–Baleanu fractional derivatives of Caputo type  that are peculiar to our study.

Definition 3.1 The Caputo fractional derivative for order 0$]]>κ>0 is defined as acDtκf(t)=1Γ(nκ)at(tξ)nκ1f(n)(ξ)dξ, where n1<κn,nN,fCn1[0,t]. Definition 3.2 The Atangana–Baleanu fractional derivative for a given function for order κ in Caputo sense are defined as aABCDtκf(t)=B(κ)1κatdf(ξ)dξEκ[κ1κ(tξ)κ]dξ, where B(κ)=(1κ)+κΓ(κ) is a normalization function and Eα() is the Mittag-Leffler function. Definition 3.3 Atangana–Baleanu fractional integral order κ is defined as 0ABItκ(f(t))=1κB(κ)f(t)+κB(κ)Γ(κ)atf(ξ)(tξ)κ1dξ if f(t) is a constant, integral will be resulted with zero. Definition 3.4 The Laplace transforms for the Atangana–Baleanu fractional operator of order κ, where 0 <κ1 is given as L{aABCDκf(t)}(s)=B(κ)1κsκL{f(t}(s)sκ1f(a)sκ+κ1κ Theorem 3.1. The following time fractional ordinary differential equation 0ABCDκf(t)=ϑ(t) has a unique solution considering the inverse Laplace transform and the convolution property below f(t)=1κB(κ)ϑ(t)+κB(κ)Γ(κ)0tϑ(ξ)(tξ)κ1dξ. Basic Properties of the Model The Invariant Region The invariant region sets out the domain where the model’s solutions are both biologically and mathematically meaningful. Since the model deals with human population, all of the model’s variables and parameters are assumed to be positive. To achieve this, we consider first the total human population Nh, where Nh = S + E + I + Q + C + R. By differentiating with respect to t both side of the total population N dNhdt=dEdt+dIdt+dQdt+dCdt+dRdt. dNhdt=θμ(S+E+I+Q+C+R)δIδCδQθμNhδIδCδQ. In the absence of the disease induced death due to COVID-19 (δ=0) , Eq. (2) becomes dNhdt=θμNh. Integrating on both side dNhθμNhdt1μln(θμNh)t+kθμNhAeμt, with the initial condition Nh(0)=Nh0 , where A is constant. Applying the initial condition in Eq. (4), we get Nhθμ[θμNhμ]eμt. As t in Eq. (5), the total human population reduces to Nhθμ . In this regard, all the feasible solution sets for human population in Eq. (1) enters and remains in the region Z={(S,E,I,Q,C,R)+6:Nhθμ}. We therefore conclude that the proposed model is well posed and are both biologically and mathematically meaningful in the domain Z . Positivity of Solution For the COVID-19 model dynamics in Eq. (1) to be epidemiologically meaningful it is important to prove that all its state variables are positive for all time. Theorem 4.1 Let the Z={(S,E,I,Q,C,R)+6:S00,E0,I0,Q0,C0,R0} . Then the solutions {S,E,I,Q,C,R} are non-negative for t0 . Proof: First, we consider the susceptible compartment in Eq. (1) give as dSdt=θβ1SEμS+αQ. dS(t)dt(β1E+μ)S(t)dS(t)S(t)(β1E+μ)dtdS(t)S(t)(β1E+μ)dt. By applying the initial condition S0 and solving the above Eq. (7), we get S(t)S0e(β1E+μ)t0. By repeating the same procedure for E,I,Q,C,R respectively in the system Eq. (1), we obtained the following results E(t)S0e(μ+β2+ε)t0,I(t)I0e(γ+δ+μ)t0,Q(t)Q0e(α+μ+δ+τ)t0, C(t)C0e(μ+δ+φ)t0,R(t)R0eμt0. This shows that the solutions of the model are positive. Hence the proof. Disease-Free Equilibrium State (DFE) The COVID-19 model in Eq. (1) has a disease-free equilibrium DFE obtain by setting the right-hand side of Eq. (1) to zero. Therefore, ΩDFE=(S,E,I,Q,C,R)=(θμ,0,0,0,0,0) Existence of Endemic Equilibrium Point (EE) We present the existence of the COVID-19 endemic equilibrium states. It is a positive equilibrium state where the COVID-19 disease is persisting in the population. Theorem 4.2. Let there be a unique endemic equilibrium state when the basic reproduction number R0 > 1 in the COVID-19 periodically forced model in Eq. (1). Proof. Suppose ΩEE=(S,E,I,Q,C,R) is a nontrivial equilibrium state of system Eq. (1) which then connote that all the compartment of ΩEE are non-negative. By equating the left-hand side of Eq. (1) to zero we get the following endemic equilibrium states S=k1k2β1ε, E=k2k3(β1εθk1k2μ)εβ1(αβ2k2+αεγk1k2k3), I=k3(β1εθk1k2μ)β1(αβ2k2+αεγk1k2k3), Q=(β2k2+εγ)(β1εθk1k2μ)β1(αβ2k2+αεγk1k2k3), C=τφ(β2k2+εγ)(β1εθk1k2μ)β1k4(αβ2k2+αεγk1k2k3)k4, R=τφ(β2k2+εγ)(β1θk1μ)β1μεk4(αβ2k2+αεγk1k2k3). where k1=(μ+β2+ε),k2=(γ+μ+δ),k3=(α+μ+β),k4=(μ+δ+φ). Basic Reproduction Number <italic>R<sub>0</sub></italic> The basic reproductive ratio is a threshold quantity that represents the overall number of secondary diseases caused by a single infected individual created into a fully susceptible population throughout its infectious period. F and V are the matrices for the new infections generated and the terms of transition. Following the same approach as , we have f=(β1SIεEγI+β1EτQ),v=((μ+β2+ε)E(γ+μ+δ)I(α+μ+δ+τ)Q(μ+δ+φ)C). The Jacobian matrix of f and v computed at the disease-free equilibrium is given as F and V such that F=(0β1S00ε000β2γ0000τ0),V=((μ+β2+ε)0000(γ+μ+δ)0000(α+μ+δ+τ)0000(μ+δ+φ)), V1=(1k100001k200001k3000τ1k4).Therefore,FV1=[β1S000εk1000β2k1γk20000τk30] By substituting k1 and k2 from Eq. (11), we therefore obtain the basic reproduction number which is the spectral radius of the matrix FV1 as R0=θεβ1μ(μ+β2+ε)(γ+μ+δ). Global Stability of the Disease-Free Equilibrium Theorem 4.3. If R01 , then the disease-free equilibrium ΩDFE=(θμ,0,0,0,0,0) is globally asymptotically stable otherwise it is unstable. Proof. We consider the Lyapunov function of the type L=x1E+x2I and L=x1E+x2I where, x1=εk1k2,x2=1k2 L=εk1k2(β1SIk1E)+1k2(εEk2I)=εβ1SIk1k2εEk2+εEk2I=εβ1SIk1k2I=(εβ1SIk1k21)(εβ1θμk1k21)I L=I(R01) From the result in Eq. (12), we can conclude that L0 provided that R01 . In addition, L=0 provided that R0=1 or I=0 . Analysis of Fractional Coronavirus Model Using the Atangana-Baleanu Operator Let us consider the mathematical model given by an ordinary differential equation system Eq. (1) using Atangana–Baleanu fractional derivative operator as below 0ABCDκS(t)=Φ1(t,S(t))=θβ1SIμS+αQ 0ABCDκE(t)=Φ2(t,E(t))=β1SI(μ+β2+ε)E 0ABCDκI(t)=Φ3(t,I(t))=εE(γ+μ+δ)I 0ABCDκQ(t)=Φ4(t,Q(t))=γI+β2E(α+μ+δ+τ)Q 0ABCDκC(t)=Φ5(t,C(t))=τQ(μ+δ+φ)C 0ABCDκR(t)=Φ6(t,R(t))=φCμR where 0ABCDκ represents the fractional operator of type Atangana–Baleanu–Caputo (ABC) having fractional order κ, where 0 <κ1 , subject to initial conditions S0(t)=S(0),E0(t)=E(0),I0(t)=I(0),Q0(t)=Q(0), C0(t)=C(0),R0(t)=R(0). The system in Eq. (14) can be converted to the Volterra-type integral equation by using the ABC fractional integral. The model is written by referring Theorem 3.1 as below: S(t)S(0)=1κB(κ){θβ1S(t)I(t)μS(t)+αQ(t)}+κB(κ)Γ(κ)0t{θβ1S(ξ)I(ξ)μS(ξ)+αQ(ξ)}(tξ)κ1dξ, E(t)E(0)=1κB(κ){β1S(t)I(t)(μ+β2+ε)E(t)}+κB(κ)Γ(κ)0t{β1S(ξ)I(ξ)(μ+β2+ε)E(ξ)}(tξ)κ1dξ, I(t)I(0)=1κB(κ){εE(t)(γ+μ+δ)I(t)}+κB(κ)Γ(κ)0t{εE(ξ)(γ+μ+δ)I(ξ)}(tξ)κ1dξ, Q(t)Q(0)=1κB(κ){γI(t)+β2E(t)(α+μ+δ+τ)Q(t)}+κB(κ)Γ(κ)0t{γI(ξ)+β2E(ξ)(α+μ+δ+τ)Q(ξ)}(tξ)κ1dξ, C(t)C(0)=1κB(κ){τQ(t)(μ+δ+φ)C(t)}+κB(κ)Γ(κ)0t{τQ(ξ)(μ+δ+φ)C(ξ)}(tξ)κ1dξ, R(t)R(0)=1κB(κ){φC(t)μR(t)}+κB(κ)Γ(κ)0t{φC(ξ)μR(ξ)}(tξ)κ1dξ. Theorem 5.1 The kernels Φ1,Φ2,Φ3,Φ4,Φ5 and Φ6 given in Eq. (14) satisfy the Lipschitz condition and contraction if the following inequality holds: 0π1,π2,π3,π4,π5,π6<1. Proof: Let the kernel Φ1(t,S(t))=θβ1SIμS+αQ. Let S1 and S2 be two functions; then we obtain the following: Φ1(t,S1(t))Φ1(t,S2(t))=(θβ1S1IμS1+αQ)(θβ1S2IμS2+αQ) =(β1I+μ)(S1(t)S2(t))(β1I(t)+μ)S1(t)S2(t) (β1M1+μ)S1(t)S2(t) (β1M1+μ)S1(t)S2(t)π1S1(t)S2(t), where π1=β1M1+μ,M1=maxtJI(t) . Similarly, we get, Φ2(t,E1(t))Φ2(t,E2(t))π2E1(t)E2(t), Φ3(t,I1(t))Φ3(t,I2(t))π3I1(t)I2(t), Φ4(t,Q1(t))Φ4(t,Q2(t))π4Q1(t)Q2(t), Φ5(t,C1(t))Φ5(t,C2(t))π5C1(t)C2(t), Φ6(t,R1(t))Φ6(t,R2(t))π6R1(t)R2(t), where, π2=(μ+β2+ε),π3=γ+μ+δ,π4=α+μ+δ+τ,π5=μ+δ+φ , and π6=μ . Considering the kernels of the model, Eq. (16) can be rewritten as S(t)=S(0)+1κB(κ){Φ1(t,S(t))}+κB(κ)Γ(κ)0t{Φ1(ξ,S(ξ))}(tξ)κ1dξ, E(t)=E(0)+1κB(κ){Φ2(t,E(t))}+κB(κ)Γ(κ)0t{Φ2(ξ,E(ξ))}(tξ)κ1dξ, I(t)=I(0)+1κB(κ){Φ3(t,I(t))}+κB(κ)Γ(κ)0t{Φ3(ξ,I(ξ))}(tξ)κ1dξ, Q(t)=Q(0)+1κB(κ){Φ4(t,Q(t))}+κB(κ)Γ(κ)0t{Φ4(ξ,Q(ξ))}(tξ)κ1dξ, C(t)=C(0)+1κB(κ){Φ5(t,C(t))}+κB(κ)Γ(κ)0t{Φ5(ξ,C(ξ))}(tξ)κ1dξ, R(t)=R(0)+1κB(κ){Φ6(t,R(t))}+κB(κ)Γ(κ)0t{Φ6(ξ,R(ξ))}(tξ)κ1dξ. Therefore, we get the following recursive formula. Sn(t)=S(0)+1κB(κ){Φ1(t,Sn1(t))}+κB(κ)Γ(κ)0t{Φ1(ξ,Sn1(ξ))}(tξ)κ1dξ, En(t)=E(0)+1κB(κ){Φ2(t,En1(t))}+κB(κ)Γ(κ)0t{Φ2(ξ,En1(ξ))}(tξ)κ1dξ, In(t)=I(0)+1κB(κ){Φ3(t,In1(t))}+κB(κ)Γ(κ)0t{Φ3(ξ,In1(ξ))}(tξ)κ1dξ, Qn(t)=Q(0)+1κB(κ){Φ4(t,Qn1(t))}+κB(κ)Γ(κ)0t{Φ4(ξ,Qn1(ξ))}(tξ)κ1dξ, Cn(t)=C(0)+1κB(κ){Φ5(t,Cn1(t))}+κB(κ)Γ(κ)0t{Φ5(ξ,Cn1(ξ))}(tξ)κ1dξ, Rn(t)=R(0)+1κB(κ){Φ6(t,Rn1(t))}+κB(κ)Γ(κ)0t{Φ6(ξ,Rn1(ξ))}(tξ)κ1dξ. We next get the difference between the iterative terms in the expression Θ1n(t)=Sn(t)S(n1)(t)=1κB(κ){Φ1(t,Sn1(t))Φ1(t,Sn2(t))} +κB(κ)Γ(κ)0t{Φ1(ξ,Sn1(ξ))Φ1(ξ,Sn2(ξ))}(tξ)κ1dξ, Θ2n(t)=En(t)E(n1)(t)=1κB(κ){Φ2(t,En1(t))Φ2(t,En2(t))} +κB(κ)Γ(κ)0t{Φ2(ξ,En1(ξ))Φ2(ξ,En2(ξ))}(tξ)κ1dξ, Θ3n(t)=In(t)I(n1)(t)=1κB(κ){Φ3(t,In1(t))Φ3(t,In2(t))} +κB(κ)Γ(κ)0t{Φ3(ξ,In1(ξ))Φ3(ξ,In2(ξ))}(tξ)κ1dξ, Θ4n(t)=Qn(t)Q(n1)(t)=1κB(κ){Φ4(t,Qn1(t))Φ4(t,Qn2(t))} +κB(κ)Γ(κ)0t{Φ4(ξ,Qn1(ξ))Φ4(ξ,Qn2(ξ))}(tξ)κ1dξ, Θ5n(t)=Cn(t)C(n1)(t)=1κB(κ){Φ5(t,Cn1(t))Φ5(t,Cn2(t))} +κB(κ)Γ(κ)0t{Φ5(ξ,Cn1(ξ))Φ5(ξ,Cn2(ξ))}(tξ)κ1dξ, Θ6n(t)=Rn(t)R(n1)(t)=1κB(κ){Φ6(t,Rn1(t))Φ6(t,Rn2(t))} where Sn=m=0Θ1m,En=m=0Θ2m,In=m=0Θ3m,Qn=m=0Θ4m,Cn=m=0Θ5m,Rn=m=0Θ6m. Applying the norm of both sides and considering triangular inequality, the Eq. (21) becomes Θ1n(t)=Sn(t)S(n1)(t)1κB(κ)Φ1(t,Sn1(t))Φ1(t,Sn2(t)) +κB(κ)Γ(κ)0t{Φ1(ξ,Sn1(ξ))Φ1(ξ,Sn2(ξ))}(tξ)κ1dξ, Θ2n(t)=En(t)E(n1)(t)1κB(κ)Φ2(t,En1(t))Φ2(t,En2(t)) +κB(κ)Γ(κ)0t{Φ2(ξ,En1(ξ))Φ2(ξ,En2(ξ))}(tξ)κ1dξ, Θ3n(t)=In(t)I(n1)(t)1κB(κ)Φ3(t,In1(t))Φ3(t,In2(t)) +κB(κ)Γ(κ)0t{Φ3(ξ,In1(ξ))Φ3(ξ,In2(ξ))}(tξ)κ1dξ, Θ4n(t)=Qn(t)Q(n1)(t)1κB(κ)Φ4(t,Qn1(t))Φ4(t,Qn2(t)) +κB(κ)Γ(κ)0t{Φ4(ξ,Qn1(ξ))Φ4(ξ,Qn2(ξ))}(tξ)κ1dξ, Θ5n(t)=Cn(t)C(n1)(t)1κB(κ)Φ5(t,Cn1(t))Φ5(t,Cn2(t)) +κB(κ)Γ(κ)0t{Φ5(ξ,Cn1(ξ))Φ5(ξ,Cn2(ξ))}(tξ)κ1dξ, Θ6n(t)=Rn(t)R(n1)(t)1κB(κ)Φ6(t,Rn1(t))Φ6(t,Rn2(t)) +κB(κ)Γ(κ)0t{Φ6(ξ,Rn1(ξ))Φ6(ξ,Rn2(ξ))}(tξ)κ1dξ. Since the kernels satisfy the Lipschitz condition, we get the following Θ1n(t)1κB(κ)π1Θ1(n1)(t)+κB(κ)Γ(κ)π10t(tξ)κ1Θ1(n1)(ξ)dξ, Θ2n(t)1κB(κ)π2Θ2(n1)(t)+κB(κ)Γ(κ)π20t(tξ)κ1Θ2(n1)(ξ)dξ, Θ3n(t)1κB(κ)π3Θ3(n1)(t)+κB(κ)Γ(κ)π30t(tξ)κ1Θ3(n1)(ξ)dξ, Θ4n(t)1κB(κ)π4Θ4(n1)(t)+κB(κ)Γ(κ)π40t(tξ)κ1Θ4(n1)(ξ)dξ, Θ5n(t)1κB(κ)π5Θ5(n1)(t)+κB(κ)Γ(κ)π50t(tξ)κ1Θ5(n1)(ξ)dξ, Θ6n(t)1κB(κ)π6Θ6(n1)(t)+κB(κ)Γ(κ)π60t(tξ)κ1Θ6(n1)(ξ)dξ. This completes the proof of the theorem. Theorem 5.2 (Existence of the Solution). The system given by Eq. (14) has a solution under the conditions that we can find tmax satisfying πjB(κ)(1κ+tmaxκΓ(κ))<1,forj{1,2,...6} . Proof: Let the functions S(t), E(t), I(t), Q(t), C(t) and R(t) are bounded. From Eq. (23) we get the following relations. Θ1n(t)S(0){π1B(κ)(1κ+tmaxκΓ(κ))}n, Θ2n(t)E(0){π2B(κ)(1κ+tmaxκΓ(κ))}n, Θ3n(t)I(0){π3B(κ)(1κ+tmaxκΓ(κ))}n, Θ4n(t)Q(0){π4B(κ)(1κ+tmaxκΓ(κ))}n, Θ5n(t)C(0){π5B(κ)(1κ+tmaxκΓ(κ))}n, Θ6n(t)R(0){π6B(κ)(1κ+tmaxκΓ(κ))}n. Hence, the functions Θ1n(t),Θ2n(t),Θ3n(t),Θ4n(t),Θ5n(t) and Θ6n(t) given in Eq. (24) exist and are smooth. Moreover, to show that the functions in Eq. (24) are the solutions of Eq. (14), we assume that S(t)S(0)=Sn(t)Δ1(n)(t), E(t)E(0)=En(t)Δ2(n)(t), I(t)I(0)=In(t)Δ3(n)(t), Q(t)Q(0)=Qn(t)Δ4(n)(t), C(t)C(0)=Cn(t)Δ5(n)(t), R(t)R(0)=Rn(t)Δ6(n)(t). where Δ1(n)(t),Δ2(n)(t), Δ3(n)(t),Δ4(n)(t),Δ5(n)(t), and Δ6(n)(t), are reminder terms of series solution. Then, we must show that these terms approach to zero at infinity, that is, Δ1()(t)0 , Δ2()(t)0,Δ3()(t)0,Δ4()(t)0,Δ5()(t)0 and Δ6()(t)0. Thus, for the term Δ1(n)(t) +κB(κ)Γ(κ)0t(tξ)κ1Φ1(ξ,S(ξ))Φ1(ξ,Sn1(ξ))dξ π1B(κ)(1κ+tmaxκΓ(κ))S(t)Sn1(t). Continuing this way recursively, we get {1B(κ)(1κ+tmaxκΓ(κ))}n+1 where B=S(t)Sn1(t) . When we take the limit of both sides as n tends to infinity, we get Δj()(t)0,j{1,2,3,4,5,6}. The Approximation Technique and Numerical Simulation Approximation Technique Consider the coronavirus model Eq. (14) along with initial conditions Eq. (15). The terms SI in this model is nonlinear. Implementing the Laplace transform on both sides of Eq. (14), we obtain, B(κ)1κsκL{S(t)}sκ1S(0)sκ+κ1κ=L{θβ1SIμS+αQ}, B(κ)1κsκL{E(t)}sκ1E(0)sκ+κ1κ=L{β1SI(μ+β2+ε)E}, B(κ)1κsκL{I(t)}sκ1I(0)sκ+κ1κ=L{εE(γ+μ+δ)I}, B(κ)1κsκL{Q(t)}sκ1Q(0)sκ+κ1κ=L{γI+β2E(α+μ+δ+τ)Q}, B(κ)1κsκL{C(t)}sκ1C(0)sκ+κ1κ=L{τQ(μ+δ+φ)C}, B(κ)1κsκL{R(t)}sκ1R(0)sκ+κ1κ=L{φCμR}. Rearranging, we get L{S(t)}=S(0)s+sκ(1κ)+κsκ(B(κ))L{θβ1SIμS+αQ}, L{E(t)}=E(0)s+sκ(1κ)+κsκ(B(κ))L{β1SI(μ+β2+ε)E}, L{I(t)}=I(0)s+sκ(1κ)+κsκ(B(κ))L{εE(γ+μ+δ)I}, L{Q(t)}=Q(0)s+sκ(1κ)+κsκ(B(κ))L{γI+β2E(α+μ+δ+τ)Q}, L{C(t)}=C(0)s+sκ(1κ)+κsκ(B(κ))L{τQ(μ+δ+φ)C}, L{R(t)}=R(0)s+sκ(1κ)+κsκ(B(κ))L{φCμR}. Further, the inverse Laplace transform on Eq. (29), yields S(t)=S(0)+L1[(sκ(1κ)+κsκ(B(κ)))L{θβ1SIμS+αQ}], E(t)=E(0)+L1[(sκ(1κ)+κsκ(B(κ)))L{β1SI(μ+β2+ε)E}], I(t)=I(0)+L1[(sκ(1κ)+κsκ(B(κ)))L{εE(γ+μ+δ)I}], Q(t)=Q(0)+L1[(sκ(1κ)+κsκ(B(κ)))L{γI+β2E(α+μ+δ+τ)Q}], C(t)=C(0)+L1[(sκ(1κ)+κsκ(B(κ)))L{τQ(μ+δ+φ)C}], R(t)=R(0)+L1[(sκ(1κ)+κsκ(B(κ)))L{φCμR}]. The series solutions achieved by the method are given by, S=n=0Sn,E=n=0En,I=n=0In,Q=n=0Qn,C=n=0Cn,R=n=0Rn. The nonlinearity of SI is written as SI=n=0Gn, whereas Gn is further decomposed as follows  Gn=j=0nSjj=0nIjj=0n1Sjj=0n1Ij. Using initial conditions, we get the recursive formula given by Sn+1(t)=Sn(0)+L1[(sκ(1κ)+κsκ(B(κ)))L{θβ1SnInμSn+αQn}], En+1(t)=En(0)+L1[(sκ(1κ)+κsκ(B(κ)))L{β1SnIn(μ+β2+ε)En}], In+1(t)=In(0)+L1[(sκ(1κ)+κsκ(B(κ)))L{{εEn(γ+μ+δ)In}}], Qn+1(t)=Qn(0)+L1[(sκ(1κ)+κsκ(B(κ)))L{γIn+β2En(α+μ+δ+τ)Qn}], Cn+1(t)=Cn(0)+L1[(sκ(1κ)+κsκ(B(κ)))L{τQn(μ+δ+φ)Cn}], Rn+1(t)=Rn(0)+L1[(sκ(1κ)+κsκ(B(κ)))L{φCnμRn}]. where S0(t)=S(0),E0(t)=E(0),I0(t)=I(0),Q0(t)=Q(0),C0(t)=C(0), R0(t)=R(0). The approximate solution is assumed to obtain as a limit when n tends to infinity. S(t)=limnSn(t),E(t)=limnEn(t),I(t)=limnIn(t),Q(t)=limnQn(t), C(t)=limnCn(t),R(t)=limnRn(t). Numerical Simulations In this section, we have presented data fitting, numerical simulations and graphical demonstration of the Atangana Baleanu COVID-19 model Eq. (14) for the population of Nigeria. We consider the available cumulative infection cases for April 1, 2020, till April 30, 2020 and parameterized the model . The parameters were estimated based on the some assumptions and facts given in Tab. 1 which plays significant role in estimating R0 (basic reproductive number). R0 is expected number of cases directly generated by one individual in a population. When R0 > 1 the infection will be able to start spreading in a population, but if R0 < 1 the disease will die out. Details Defination of Variables and Parameters Parameter Description Value Source θ recruitment rate into susceptible class 1.3  μ natural death rate 4.317 × 10−5  δ covid-19 infection death rate 0.00618 Fitted β1 force of infection 2.51×107 Fitted β2 proportion of people identified as suspected cases 0.04 Fitted ε progression rate from exposed class to highly infected class 0.00916 Fitted α Rate of returning to susceptible class after diagnosis 0.068  γ Progression rate from highly infectious class to quarantine class 0.001 Fitted τ progression rate from quarantine to confirm case after diagnosis 0.002  φ rate of recovery from infection 0.0714 Fitted Next, we evaluate and present the number of cumulative infectious cases in different compartments with respect to time in days using various plots for population of Nigeria. We begin estimation of our model from initial time (t = t0 ) as per data reported on April 1, 2020 . Hence, the required initial values are S(0) = 205773342 ; E(0) = 15,000; I(0) = 100; Q(0) = 100; C(0) =175; R(0) = 31. By applying iterative Laplace transform using Eq. (31) successively up to four terms we get series form approximate solution of the fractional COVID-19 model Eq. (14) as given below S(t)=2057733420.0000382421(1κ+κGamma[κ])2+0.0000764842κ(1κ+κGamma[κ])20.0000382421κ2(1κ+κGamma[κ])2+2.61κ+κGamma[κ]2.6κ1κ+κGamma[κ]+0.000012912814395020359tκ(1κ+κGamma[κ])3(κ+(1+κ)Gamma[κ])+ E(t)=15000+0.000032630000000000004(1κ+κGamma[κ])20.00006526κ(1κ+κGamma[κ])2+0.000032630000000000004κ2(1κ+κGamma[κ])20.000012912814395020362tκ(1κ+κGamma[κ])3(κ+(1+κ)Gamma[κ])+0.00003873844318506108tκκ(1κ+κGamma[κ])3(κ+(1+κ)Gamma[κ])..... I(t)=100+0.28430815806309984tκ(1κ+κGamma[κ])2(κ+(1+κ)Gamma[κ])0.5686163161261999tκκ(1κ+κGamma[κ])2(κ+(1+κ)Gamma[κ])+0.28430815806309984tκκ2(1κ+κGamma[κ])2(κ+(1+κ)Gamma[κ])39.573442828747645tκ(1κ+κGamma[κ])(κ+(1+κ)Gamma[κ])+.... Q(t)=100+10.024125186027561tκ(1κ+κGamma[κ])2(κ+(1+κ)Gamma[κ])20.04825037205514tκκ(1κ+κGamma[κ])2(κ+(1+κ)Gamma[κ])+10.024125186027561tκκ2(1κ+κGamma[κ])2(κ+(1+κ)Gamma[κ])132.0967231202757tκ(1κ+κGamma[κ])(κ+(1+κ)Gamma[κ])+ C(t)=1752.222828365962386tκ(1κ+κGamma[κ])(κ+(1+κ)Gamma[κ])+2.2228283659623855tκκ(1κ+κGamma[κ])(κ+(1+κ)Gamma[κ])2.222828365962386Gamma[κ](1κ+κGamma[κ])(κ+(1+κ)Gamma[κ])+4.445656731924771κGamma[κ](1κ+κGamma[κ])(κ+(1+κ)Gamma[κ]).... R(t)=31+0.9551899812522688tκ(1κ+κGamma[κ])(κ+(1+κ)Gamma[κ])0.9551899812522688tκκ(1κ+κGamma[κ])(κ+(1+κ)Gamma[κ])+0.9551899812522688Gamma[κ](1κ+κGamma[κ])(κ+(1+κ)Gamma[κ])1.9103799625045377κGamma[κ](1κ+κGamma[κ])(κ+(1+κ)Gamma[κ])+.... In Tab. 1, we have estimated some needed biological parameter values related with basic reproduction number R0 corresponds to model (1) like, covid-19 infection death rate ( δ ), force of infection ( β1 ), proportion of people identified as suspected cases ( β2 ), progression rate from exposed class to highly infected class ( ε ) and other parameters are fitted from previous literature. Estimating these parameters to suitable values decreases rate of infection meaningfully. It is examined and noted that if R0 is near to 2 the number of cumulative infected population rises very quickly. Also reducing the value of R0 near to one lowers the number of infected cases rapidly. Figs. 2a7b shows the dynamical behaviour of various classes of mathematical model like susceptible population S(t), symptomatic and undetected population E(t), highly infectious but not yet quarantined or isolated I(t), individuals who are infected or suspected and quarantined Q(t), confirmed and quarantine population C(t), recovered population R(t) with R0 = 1.96 and R0 = 1.16 respectively for various values of κ = 1, 0.9, 0.8, 0.7 verses time in days. It is observed that as the value of κ decreases from 1 the effect of the fractional derivative order becomes prominent. We clearly observe the significant variance in both type of plots (a) and (b) in each figure when value of fractional parameter κ changes. From Figs. 3a and 3b, it is noticed that as the time increase with respect to time the symptomatic and undetected population E(t) also increases due to spread of infection in the population and their interaction with infected people. It is clear that reducing the value of fractional parameter κ significantly affect COVID-19 dynamics and greatly decrease the fraction of asymptomatic carriers with level of infection in the population. Hence, the number of cumulative cases of infections in every class is continuously depends upon the values of fractional order. Plots of susceptible population for various values of <inline-formula id="ieqn-104"> <alternatives><inline-graphic xlink:href="ieqn-104.png"/><tex-math id="tex-ieqn-104"><![CDATA[$\bi{\kappa }$]]></tex-math><mml:math id="mml-ieqn-104"><mml:mi mathvariant="bold-italic">κ</mml:mi></mml:math> </alternatives></inline-formula> with respect to time t in days and different R<sub>0</sub> [Fig. 2(a): R<sub>0</sub> = 1.96, Fig. 2(b): R<sub>0</sub> = 1.16] Plots of asymptotic and undected population for various values of <inline-formula id="ieqn-105"> <alternatives><inline-graphic xlink:href="ieqn-105.png"/><tex-math id="tex-ieqn-105"><![CDATA[$\bi {\kappa }$]]></tex-math><mml:math id="mml-ieqn-105"><mml:mi mathvariant="bold-italic">κ</mml:mi></mml:math> </alternatives></inline-formula> with respect to time t in days and different R<sub>0</sub> [Fig. 3(a): R<sub>0</sub> = 1.96, Fig. 3(b): R<sub>0</sub> = 1.16] Plots of highly infectitious population for various values of <inline-formula id="ieqn-106"> <alternatives><inline-graphic xlink:href="ieqn-106.png"/><tex-math id="tex-ieqn-106"><![CDATA[$\bi {\kappa }$]]></tex-math><mml:math id="mml-ieqn-106"><mml:mi mathvariant="bold-italic">κ</mml:mi></mml:math> </alternatives></inline-formula> with respect to time t in days and different R<sub>0</sub> [Fig. 4(a): R<sub>0</sub> = 1.96, Fig. 4(b): R<sub>0</sub> = 1.16] Plots of symptomatics population for various values of <inline-formula id="ieqn-107"> <alternatives><inline-graphic xlink:href="ieqn-107.png"/><tex-math id="tex-ieqn-107"><![CDATA[$\bi {\kappa }$]]></tex-math><mml:math id="mml-ieqn-107"><mml:mi mathvariant="bold-italic">κ</mml:mi></mml:math> </alternatives></inline-formula> with respect to time t in days and different R<sub>0</sub> [Fig. 5(a): R<sub>0</sub> = 1.96, Fig. 5(b): R<sub>0</sub> = 1.16] Plots of confirmed infected population for various values of <inline-formula id="ieqn-108"> <alternatives><inline-graphic xlink:href="ieqn-108.png"/><tex-math id="tex-ieqn-108"><![CDATA[$\bi {\kappa }$]]></tex-math><mml:math id="mml-ieqn-108"><mml:mi mathvariant="bold-italic">κ</mml:mi></mml:math> </alternatives></inline-formula> with respect to time t in days and different R<sub>0</sub> [Fig. 6(a): R<sub>0</sub> = 1.96, Fig. 6(b): R<sub>0</sub> = 1.16] Plots of recovered population for various values of <inline-formula id="ieqn-109"> <alternatives><inline-graphic xlink:href="ieqn-109.png"/><tex-math id="tex-ieqn-109"><![CDATA[$\bi {\kappa }\$]]></tex-math><mml:math id="mml-ieqn-109"><mml:mi mathvariant="bold-italic">κ</mml:mi></mml:math> </alternatives></inline-formula> with respect to time t in days and different R<sub>0</sub> [Fig. 7(a): R<sub>0</sub> = 1.96, Fig. 7(b): R<sub>0</sub> = 1.16]

Also, Fig. 8a shows the comparison between estimated and actual number of cumulative cases of confirmed and quarantine population C(t) for data of available infection cases from April 1, 2020, till April 30, 2020 of Nigeria. This clearly shows that estimated and actual cases of infections are very near to each other. Fig. 8b shows the plots of total infected population in various compartments of model Eq. (1) versus time (days).

(a) Cumulative number of cases of infected population C(t) with respect to time(days) (b) Plot of total population in various classes of model versus time (days)
Conclusions

In this manuscript, we have analysed and examined transmission dynamics of COVID-19 infection formulated in terms of mathematical model based on fractional differential system. We have used Atangana–Baleanu fractional derivative operator to obtain existence criteria of solution of mathematical model for the operator. The numerical simulations are carried out using iterative Laplace transform method. The essential axioms of the proposed model have been studied to observe the biological and mathematical feasibility. Further, we have examined the sensitivity analysis by finding the basic reproductive number R0 which explains the significant of every biological parameter involved in the proposed model. Moreover, the local and global asymptotic stability conditions for the disease-free and endemic equilibrium are obtained which determines the conditions to stabilize the exponential spread of the disease. It is noted from this analysis, that parameters β1 and ε strengthen the outbreak of the infection at large extent and needs notable consideration to implement some control strategies to keep this under control. Towards the end all the hypothetical results are supported with the assistance of graphical portrayal by numerical investigation which would be beneficial for researchers to contemplate the dynamics of the COVID-19.

Funding Statement: The authors received no specific funding for this study.

Conflicts of Interest: “The authors declare that they have no conflicts of interest to report regarding the present study.”