Fractal Fractional Order Operators in Computational Techniques for Mathematical Models in Epidemiology

1  Introduction

Coronavirus (COVID-19) is a new phenomenon in recent days, which affected the entire world while it was emerging. According to reference [1], it is reported that a mysterious outbreak of atypical pneumonia was traced to a seafood market in Wuhan, China. In December 2019 the first case of novel coronavirus was reported. The symptoms of coronavirus are dry cough, fever, fatigue and in severe cases, acute respiratory syndrome that appears in 2–10 days which further causes pneumonia, kidney failure and even death [2]. In between March and April, coronavirus became a global phenomenon with the whole world facing an emergency situation. Initial cases were reported in the wet seafood market of Wuhan, China [3]. That is why some researchers thought that it is transmitted by animals to humans. This virus is transmitted from one person to another through physical contact, droplets during sneezing and coughing [4]. Researchers in the field of epidemiology and other fields of biology are trying hard to develop the cure based on ongoing clinical trials, but different researching companies from different countries have developed the vaccine for COVID-19 in [5]. Developed countries like the USA, UK, Italy, Spain and many others are affected very badly, most of the global deaths are being reported from these countries [6]. Mathematical modeling is used to understand the dynamics and behavior of disease and then develop the procedures for the treatment of disease. For this purpose, many researchers developed the COVID-19 models (see [7–12]). Reproductive number has a notable role in the analysis of mathematical models. Reproductive number explains the behavior of the simulation of COVID-19. The fractional order mathematical models of a few more infectious diseases have recently been studied in [13].

In order to address problems in the real world, fractional calculus (FC) is essential. It is widely utilized in a variety of scientific, engineering, and financial sectors. The key characteristics of FC are fractional integrals and derivatives of fractional order. Researchers’ interest in fractional calculus and the numerous aspects of that study under inquiry has grown in recent years. This is due to the fact that genetic mutations are a crucial tool for characterizing the dynamic operation of diverse biological systems. These component operators’ non-local properties, which are absent from the integer separator operator, give them their power [14–16]. For the development of an artificial pancreas, Farman et al. [17] employed an Atangana Baleanu derivative to manage glucose levels in insulin treatments. Differential equations with various generalized fractional derivatives have been solved using a variety of numerical techniques [18–20]. In [21], a generalization of the squared remainder minimization method for resolving multi-term fractional differential equations was developed. The Caputo time-fractional derivative and redefined extended B-spline functions have been used for the time and spatial discretization, respectively in [22–25] and some details are also given in [26–28]. COVID-19 outbreaks have been well modeled in [29–31] for a variety of geographic locations. Additionally, some publications [32–34] explored the impact of quarantine and social isolation on the viral load in the environment. For the COVID-19 epidemic, some researchers suggested the best control approaches including cost-effectiveness assessments [35,36].

2  Basic Concepts of Fractional Operators

In this section, we consider some definition related to fractal fractional operator given in [31,34,37,38].

Definition 2.1: Let 0≤σ,σ1≤1, then with power law type kernel the fractal-fractional derivative is described by:

FFPJ0,tσ,σ1(g(t))=1Γ(m−σ)ddtσ1∫0t(t−s)m−σ−1g(s)ds,(1)

ddsσ1g(s)=limt→sg(t)−g(s)tσ1−sσ1.

Definition 2.2: Let 0≤σ,σ1≤1, then with the exponential-decay type kernel the fractal fractional derivative is described by:

FFPJ0,tσ,σ1(g(t))=M(σ)Γ(m−σ)ddtσ1∫0texp[−σ1−σ(t−s)ν−σ−1]g(s)ds,(2)

where σ>0, σ1≤ν∈N, and M(0)=M(1)=1.

Definition 2.3: Let 0≤σ,σ1≤1, then with the generalized Mittag-Leffler type kernel the fractal fractional derivative is described by:

FFMJ0,tσ,σ1(g(t))=AB(σ)1−σddtσ1∫0tEσ[−σ1−σ(t−s)σ]g(s)ds,(3)

where σ>0, σ1≤1, and AB(σ)=1−σ+σΓ(σ).

Definition 2.4: The function g(t) of order (σ,σ1), for fractal-fractional integral with power law type kernel is described by:

FFPI0,tσ,σ1(g(t))=1Γ(σ1)∫0ts1−σ1(t−s)σ−1g(s)ds.(4)

Definition 2.5: The function g(t) of order (σ,σ1), for fractal-fractional integral with exponential-decay type kernel is described by:

FFPI0,tσ,σ1(g(t))=σ1(1−σ)tσ1−1g(t)M(σ)+σσ1M(σ)∫0tsσ−1g(s)ds.(5)

Definition 2.6: The function g(t) of order (σ,σ1), for fractal-fractional integral with Mittag-Leffler type kernel is described by:

FFMJ0,tσ,σ1(g(t))=σ1(1−σ)tσ1−1g(t)AB(σ)+σσ1AB(σ)∫0tsσ−1(t−s)g(s)ds.(6)

3  Fractal Fractional Order Model

We suppose the COVID-19 model formulated by Ahmad et al. [39]. In this model, S(t) represents susceptible individuals, H(t) represents resistant or healthy individuals, infected individuals are represented by I(t) and Q(t) represents quarantined individuals. We suppose that the used parameters in the model are non-negative. Hence, N(t)=S(t)+H(t)+I(t)+Q(t). In this model, recruitment rate of susceptible individuals is represented by λ, γ denotes disease transmission rate, Recruitment rate of healthy people is α, Healthy people transmission rate is denoted by β, Cure rate of the infected individuals in the quarantined compartment is θ, δ which represent the rate at which quarantined people get infections, μ represents death rate of suspected or infected individuals due to disease and d denotes natural death rate. We present the COVID-19 classical model [39] using fractal-fractional Atangana–Baleanu derivative. We have the following model:

{FFD0,tσ,σ1S(t)=λ−γS(t)I(t)−(d+μ)S(t)FFD0,tσ,σ1H(t)=α−βH(t)I(t)+θI(t)−(d+μ)H(t)FFD0,tσ,σ1I(t)=γS(t)I(t)+βH(t)I(t)+δQ(t)−(d+μ+η+θ)I(t)FFD0,tσ,σ1Q(t)=ηI(t)−(d+μ+δ)Q(t),(7)

with initial conditions are

S(0) ≥ 0, H(0) ≥ 0, I(0) ≥ 0 and Q(0) ≥ 0.(8)

3.1 Equilibrium Points

In this section, we will discuss the equilibrium points of the given COVID-19 model (7). Equilibrium points have two types namely as disease free equilibrium and endemic equilibrium. We obtained these points by putting the number zero on the right side of the system (7). We suppose that E’ represents disease free equilibrium and endemic equilibrium is represented by E*. If we take both of our equilibriums, we have

E′=(S′,H′,I′,Q′)=(λd+μ,αd+μ,0,0)

S∗=λγI∗+d+μ, H∗=λ+θI∗βI∗+d+μ, I∗=δQ∗(d+μ+η+θ)−γS∗−βH∗, Q∗=ηI∗(d+μ+δ).

We obtain the basic reproductive number R0 by [37], we have

R0=(γλ+βα)(d+μ+δ)(d+μ)[(d+μ+η+θ)(d+μ+δ)−δη].

4  Existence and Stability Theory

4.1 Existence

We consider [40]

{0ABRD0σS(t)=σ1tσ1−1C(t,S,H,I,Q)0ABRD0σH(t)=σ1tσ1−1D(t,S,H,I,Q)0ABRD0σI(t)=σ1tσ1−1E(t,S,H,I,Q)0ABRD0σQ(t)=σ1tσ1−1F(t,S,H,I,Q),(9)

where

C(t,S,H,I,Q)=λ−γS(t)I(t)−(d+μ)S(t),

D(t,S,H,I,Q)=α−βH(t)I(t)+θI(t)−(d+μ)H(t),

E(t,S,H,I,Q)=γS(t)I(t)+βH(t)I(t)+δQ(t)−(d+μ+η+θ)I(t),

F(t,S,H,I,Q)=ηI(t)−(d+μ+δ)Q(t).

We can write system (9) as:

{0ABRDtσΠ(t)=σ1tσ1−1Λ(t,Π(t))Π(0)=Π0,(9’)

By replacing 0ABRD0σ,σ1 by 0ABCD0σ,σ1 and applying fractional integral, we get

Π(t)=Π(0)+σ1tσ1−1(1−σ)AB(σ)Λ(t,Π(t))+σσ1AB(σ)Γ(σ)∫0tτσ1−1(t−τ)σ1−1Λ(τ,Π(τ))dτ,

where

Π(t)={S(t)H(t)I(t)Q(t),Π(0)={S(0)H(0)I(0)Q(0),Λ(t,Π(t))={C(t,S,H,I,Q)D(t,S,H,I,Q)E(t,S,H,I,Q)F(t,S,H,I,Q)

We describe a Banach space B=C×C×C×C, where C=[0,T] under the norm

‖Π‖=maxt∈[0,T]|S(t)+H(t)+I(t)+Q(t)|

Define as operator ℵ:B→B as:

ℵ(Π)(t)=Π(0)+σ1tσ1−1(1−σ)AB(σ)Λ(t,Π(t))+σσ1AB(σ)Γ(σ)∫0tτσ1−1(t−τ)σ1−1Λ(τ,Π(τ))dτ.(10)

We suppose that

•   For each Π∈B,∃ constants CΛ>0 and MΛ such that

|Λ(t,Π(t))|≤CΛ|Π(t)|+MΛ.(11)

•   Considering Π,Π¯∈B,∃ a constant LΛ>0 such that

|Λ(t,Π(t))−Λ(t,Π¯(t))|≤LΛ|Π(t)−Π¯(t)|.(12)

Theorem 4.1: Suppose that the state (11) exists. Let Λ:[0,T]×B→R be a continuous function. The system having at least one solution, the condition given in [41,42].

Proof: First of all, considering the Eq. (10) is completely continuous which is described by operator ℵ. Since Λ and ℵ are continuous operators,

Suppose that H=∖{Π∈B:‖Π‖≤R, R>0∖}. For some Π∈B, we have

|ℵ(Π)|=maxt∈[0,T]|Π(0)+σ1tσ1−1(1−σ)AB(σ)Λ(t,Π(t))+σσ1AB(σ)Γ(σ)∫0tτσ1−1(t−τ)σ1−1Λ(τ,Π(τ))dτ|≤Π(0)+σ1Tσ1−1(1−σ)AB(σ)(CΛ‖Π‖+MΛ)+maxt∈[0,T]σσ1AB(σ)Γ(σ)∫0tτσ1−1(t−τ)σ1−1Λ(τ,Π(τ))dτ≤Π(0)+σ1Tσ1−1(1−σ)AB(σ)(CΛ‖Π‖+MΛ)+σσ1AB(σ)Γ(σ)(CΛ‖Π‖+MΛ)Tσ+σ1−1H(σ,σ1)≤R.

Therefore, we get

|ℵ(Π)(t2)−ℵ(Π)(t1)|=|σ1t2σ1−1(1−σ)AB(σ)Λ(t2,Π(t2))+σσ1AB(σ)Γ(σ)∫0t2τσ1−1(t2−τ)σ1−1Λ(τ,Π(τ))dτ−σ1t1σ1−1(1−σ)AB(σ)Λ(t1,Π(t1))+σσ1AB(σ)Γ(σ)∫0t1τσ1−1(t1−τ)σ1−1Λ(τ,Π(τ))dτ|≤|σ1t2σ1−1(1−σ)AB(σ)(CΛ|Π(t)|+MΛ)+σσ1AB(σ)Γ(σ)(CΛ|Π(t)|+MΛ)t2σ+σ1−1H(σ,σ1)−σ1t1σ1−1(1−σ)AB(σ)(CΛ|Π(t)|+MΛ)+σσ1AB(σ)Γ(σ)(CΛ|Π(t)|+MΛ)t1σ+σ1−1H(σ,σ1)|

When t1→t2 then |ℵ(Π)(t2)−ℵ(Π)(t1)|→0|.

‖ℵ(Π)(t2)−ℵ(Π)(t1)‖→0, as t1→t2.

Thus, ℵ is equicontinuous. Then, by Schauder’s fixed point result, the condition is held.

Theorem 4.2: [38] Suppose that the condition (12) holds. If ρ<1, where

ρ=(σ1Tσ1−1(1−σ)AB(σ)+σσ1AB(σ)Γ(σ)(CΛ‖Π‖+MΛ)Tσ+σ1−1H(σ,σ1))LΛ.

Then the solution of the system is unique.

Proof: For all Π,Π¯∈B, acquire the following:

|ℵ(Π)−ℵ(Π¯)|=maxt∈[0,T]|σ1tσ1−1(1−σ)AB(σ)(Λ(t,Π(t))−Λ(t,Π¯(t)))+σσ1AB(σ)Γ(σ)∫0tτσ1−1(t−τ)σ1−1dτ(Λ(τ,Π(τ))−Λ(τ,Π¯(τ)))

≤[σ1Tσ1−1(1−σ)AB(σ)+σσ1AB(σ)Γ(σ)(CΛ‖Π‖+MΛ)Tσ+σ1−1H(σ,σ1)]‖Π−Π¯‖≤ρ‖Π−Π¯‖

Therefore, ℵ is a contraction. Thus, the solution of the system is unique according to Banach contraction principle [43].

We denote by:

{T1(S(t), H(t), I(t),Q(t))=S(t)−0ABRD0σS(t)+σ1tσ1−1C(t,S,H,I,Q)T2(S(t), H(t), I(t),Q(t))=H(t)−0ABRD0σH(t)+σ1tσ1−1D(t,S,H,I,Q)T3(S(t), H(t), I(t),Q(t))=I(t)−0ABRD0σI(t)+σ1tσ1−1E(t,S,H,I,Q)T4(S(t), H(t), I(t),Q(t))=Q(t)−0ABRD0σQ(t)+σ1tσ1−1F(t,S,H,I,Q).(13)

Let Pcl(R) be the set of all nonempty closed subsets of R.

Theorem 4.3: [44] Let A∈Mm,m(R+). The following are equivalents:

(i) A is a matrix which converges to zero;

(ii) An→0 as n→∞;

(iii) The modulus for every eigen-values of A is lower than 1;

(iv) The matrix I—A is non-singular, with (I−A)−1=I+A+⋯+An+⋯.

We give the following theorem, for the hypothesis that Ti:R→ Pcl(R)  for i∈{1,…,4} are contractions. Pcl(R)  is the set of all nonempty closed subsets of R, where R represents the set of real numbers. Examples of conditions for them to be contractions are, for instance, the cases in which the absolute values of the derivatives are lower than 1. This situation is possible when the variations of functions S(t), H(t), I(t),Q(t) are low.

Theorem 4.4: Let Ti:R→ Pcl(R)  for i∈{1,…,4} be contractions and  0≤akk≤1,k∈{1,…,4}. Let (S(t1), H(t1), I(t1), Q(t1)),(S(t2), H(t2), I(t2), Q(t2))∈R4 where t1 and t2 are moments from a time interval J. If for each yk=Tk(S(t1), H(t1), I(t1), Q(t1)), k∈{1,…,4} there exists zk=Tk(S(t1), H(t1), I(t1), Q(t1)) such that for all k∈{1,…,4} in [44]:

|yk−zk|≤a11|S(t2)−S(t1)|,

|yk−zk|≤a22|H(t2)−H(t1)|,

|yk−zk|≤a33|I(t2)−I(t1)|,

|yk−zk|≤a44|Q(t2)−Q(t1)|,

then, the semi linear inclusion system:

{S(t)∈T1(S(t), H(t), I(t),Q(t))H(t)∈T2(S(t), H(t), I(t),Q(t))I(t)∈T3(S(t), H(t), I(t),Q(t))Q(t)∈T4(S(t), H(t), I(t),Q(t))(14)

has at least one solution in R4.

Proof: The theorem is a particular case of Theorem 3.11 from [44]. It is demonstrated the same as Theorem 3.11 from [44], with (u1, u2, u3,u4)= (S(t1), H(t1), I(t1), Q(t1)) and (v1, v2, v3,v4) = (S(t2), H(t2), I(t2), Q(t2)) and working with ‖u−v‖=(|u1−v1||u2−v2||u3−v3||u4−v4|). In this case, the diagonal matrix A = (aij) converges to 0.

The demonstration of the theorem uses elements from the fixed-point theory and results from the fact that T=(T1, …,T4): R4→ Pcl(R4) is a multivalued operator A-contraction to the left, thus T is an MWP operator. The concept of multivalued weakly Picard operator (briefly MWP operator) was introduced by Rus et al. in [45]. The authors created this concept in connection to the successive approximation technique for the fixed-point set of multivalued operators defined on a complete metric space. As R4 is a Banach space, T has at least one fixed point [44], therefore, the conclusion to this theorem is verified.

4.2 Ulam-Hyers Stability

Definition 4.1: The system is Ulam-Hyers stable if ∃ ℵσ,σ1≥0 such that for any ε>0 and for every Π∈(C[0,T],R) is satisfied the following:

|0FFMDtσ,σ1Π(t)−Λ(t,Π(t))|≤ε, t∈[0,T],

And there exists unique solution Ω∈(C[0,T],R) such that

|Π(t)−Ω(t)|≤ℵσ,σ1ε, t∈[0,T],

•   Φ(t)≤ε for ε>0.

•   0FFMDtσ,σ1Π(t)=Λ(t,Π(t))+Φ(t).

Lemma 4.1: Perturbed solution for the system according the given result in [37].

0FFMDtσ,σ1Π(t)=Λ(t,Π(t))+ϕ(t)

Π(0)=Π0

satisfies the following relation:

|ℵ(t)−(Π(0)+σ1tσ1−1(1−σ)AB(σ)Λ(t,Π(t))+σσ1AB(σ)Γ(σ)∫0tτσ1−1(t−τ)σ1−1Λ(τ,Π(τ))dτ)|≤(σ1Tσ1−1(1−σ)AB(σ)+σσ1AB(σ)Γ(σ)Tσ+σ1−1H(σ,σ1))ε,(15)

We note:

σσ,σ1∗=σ1Tσ1−1(1−σ)AB(σ)+σσ1AB(σ)Γ(σ)Tσ+σ1−1H(σ,σ1).(16)

Lemma 4.2: The solution of the system is Ulam-Hyers stable if ρ<1, in condition (12) along with Lemma (4.1).

Proof: Suppose that Ω∈B and Π∈B is unique and any solution, respectively, we have

|Π(t)−Ω(t)|=|Π(t)−(Ω(0)+σ1tσ1−1(1−σ)AB(σ)Λ(t,Ω(t))+σσ1AB(σ)Γ(σ)∫0tτσ1−1(t−τ)σ1−1Λ(τ,Ω(τ))dτ)|≤|Π(t)−(Π(0)+σ1tσ1−1(1−σ)AB(σ)Λ(t,Π(t))+σσ1AB(σ)Γ(σ)∫0tτσ1−1(t−τ)σ1−1Λ(τ,Π(τ))dτ)|+|Π(0)+σ1tσ1−1(1−σ)AB(σ)Λ(t,Π(t))+σσ1AB(σ)Γ(σ)∫0tτσ1−1(t−τ)σ1−1Λ(τ,Π(τ))dτ|−|Ω(0)+σ1tσ1−1(1−σ)AB(σ)Λ(t,Ω(t))+σσ1AB(σ)Γ(σ)∫0tτσ1−1(t−τ)σ1−1Λ(τ,Ω(τ))dτ|≤σσ,σ1∗ε+(σ1Tσ1−1(1−σ)AB(σ)+σσ1AB(σ)Γ(σ)Tσ+σ1−1H(σ,σ1))LΛ|Π(t)−Ω(t)|≤σσ,σ1∗ε+ρ|Π(t)−Ω(t)|.

Consequently, one can write

||Π(t)−Ω(t)||≤σσ,σ1∗ε+ρ||Π(t)−Ω(t)||.

We can write the above relation is

||Π(t)−Ω(t)||≤ℵσ,σ1ε,

where ℵσ,σ1=σσ,σ1∗1−ρ. Therefore, system is stable.

4.3 Fractal-Fractional Integral with Mittag-Leffler Kernel

Consider:

{S(t)=S(0)+σ1tσ1−1(1−σ)AB(σ)P1(t,S,H,I,Q)+σσ1AB(σ)Γ(σ)∫0tτσ1−1(t−τ)σ−1P1(τ,S,H,I,Q)dτH(t)=H(0)+σ1tσ1−1(1−σ)AB(σ)P2(t,S,H,I,Q)+σσ1AB(σ)Γ(σ)∫0tτσ1−1(t−τ)σ−1P2(τ,S,H,I,Q)dτI(t)=I(0)+σ1tσ1−1(1−σ)AB(σ)P3(t,S,H,I,Q)+σσ1AB(σ)Γ(σ)∫0tτσ1−1(t−τ)σ−1P3(τ,S,H,I,Q)dτQ(t)=Q(0)+σ1tσ1−1(1−σ)AB(σ)P4(t,S,H,I,Q)+σσ1AB(σ)Γ(σ)∫0tτσ1−1(t−τ)σ−1P4(τ,S,H,I,Q)dτ.(17)

We construct the numerical scheme at t=tn+1:

{Sn+1=S0+σ1tnσ1−1(1−σ)AB(σ)P1(tn,Sn,Hn,In,Qn)+σσ1AB(σ)Γ(σ)∫0tτσ1−1(tn+1−τ)σ−1P1(τ,S,H,I,Q)dτHn+1=H0+σ1tnσ1−1(1−σ)AB(σ)P2(tn,Sn,Hn,In,Qn)+σσ1AB(σ)Γ(σ)∫0tτσ1−1(tn+1−τ)σ−1P2(τ,S,H,I,Q)dτIn+1=I0+σ1tnσ1−1(1−σ)AB(σ)P3(tn,Sn,Hn,In,Qn)+σσ1AB(σ)Γ(σ)∫0tτσ1−1(tn+1−τ)σ−1P3(τ,S,H,I,Q)dτQn+1=Q0+σ1tnσ1−1(1−σ)AB(σ)P4(tn,Sn,Hn,In,Qn)+σσ1AB(σ)Γ(σ)∫0tτσ1−1(tn+1−τ)σ−1P4(τ,S,H,I,Q)dτ(18)