[BACK]
images Computer Modeling in Engineering & Sciences images

DOI: 10.32604/cmes.2022.019221

ARTICLE

Modeling the Spread of Tuberculosis with Piecewise Differential Operators

Abdon Atangana1,2 and Ilknur Koca3,*

1Institute for Groundwater Studies, Faculty of Natural and Agricultural Sciences, University of the Free State, Bloemfontein, 9301, South Africa
2Department of Medical Research, China Medical University Hospital, China Medical University, Taichung, 008864, Taiwan
3Department of Accounting and Financial Management, Seydikemer High School of Applied Sciences, Mugla Sitki Kocman University, Mugla, 48300, Turkey
*Corresponding Author: Ilknur Koca. Email: ilknurkoca@mu.edu.tr
Received: 09 September 2021; Accepted: 09 November 2021

Abstract: Very recently, a new concept was introduced to capture crossover behaviors that exhibit changes in patterns. The aim was to model real-world problems exhibiting crossover from one process to another, for example, randomness to a power law. The concept was called piecewise calculus, as differential and integral operators are defined piece wisely. These behaviors have been observed in the spread of several infectious diseases, for example, tuberculosis. Therefore, in this paper, we aim at modeling the spread of tuberculosis using the concept of piecewise modeling. Several cases are considered, conditions under which the unique system solution is obtained are presented in detail. Numerical simulations are performed with different values of fractional orders and density of randomness.

Keywords: Spread of tuberculosis; piecewise differentiation; numerical simulation

1  Introduction

Although several studies have been done on behaviors of the tuberculosis virus, its spread, and its effect on the human’s body until today, this virus persists and kills humans around the world each year. It is even believed that the tuberculosis virus has affected about 25 percent of the world population since about one percent of the world population is infected each year according to what is reported in the literatures [14]. Tuberculosis is a seasonal transmissible disease, as the peaks are reached every spring and summer. However, there is no apparent scientific reason recorded that can explain this variation. Nevertheless, it is recorded that the virus spreads more during weather conditions like low temperature, humidity, and low rainfall. Thus tuberculosis incidence rates could be linked to change in the climate. Having peaks that occurred during some period of the year show that the spread had many waves since antiquity. Indeed, each wave has a specific pattern different from others or similar in some cases. It can be concluded that the virus spread follows piecewise patterns. Mathematicians have tried to provide mathematical models to depict the spread behaviors as a function of time. Several studies have been performed in the decades. The reproductive number of this virus has been calculated in many studies. New and modified models have been provided and studied in detail. Several differential and integral operators have been used, for example, fractional differential operators to replicate spread behaviors. Fractional derivative based on power law was introduced to replicate behaviors resembling the power law [511]. Different techniques have been employed, for example, the stochastic process to capture random behaviors. Nevertheless, the problem of different was not really addressed. The concept of piecewise differential and integral operators was recently suggested and employed to model some complex real-world problems, such as chaos and other epidemiological problems [12]. The concept seems to be efficient when modeling problems with crossover behaviors. In this paper, we aim to modify an existing tuberculosis model with the concept of piecewise differentiation.

1.1 Important Definitions of Fractional Modelling

Definition 1: Let α > 0 of a function h:(0,)R and the Riemann-Liouville derivative of fractional order is presented as

Dtαh(t)=1Γ(1α)ddt0t(tx)αh(x)dx,0<α1.(1)

Definition 2: Let h : H1(a, b), b > a, 0 < α < 1 then, the Caputo-Fabrizio derivative of fractional derivative is presented as

aCFDtαh(t)=11αath(x)exp[α(tx)1α]dx.(2)

Definition 3: Let h : H1(a, b), b > a, α ∈ (0, 1) then, the definition of the new fractional derivative (Atangana-Baleanu derivative in Caputo sense) is presented as

aABCDtαh(t)=AB(α)1αath(x)Eα[α(tx)α1α]dx,(3)

where aABCDtα is fractional operator with Mittag-Leffler kernel in the Caputo sense with order α with respect to t and

AB(α)=1α+αΓ(α),(4)

is a normalization function.

Definition 4: Let h be continuous not necessary differentiable in [t1, T]. Thus, the piecewise Riemann-Liouville derivative is presented as

0PRLDtαh(t)={h(t),if 0tt1t1RLDtαh(t),if t1tT,(5)

where 0PRLDtα presents classical derivative on 0≤ tt1 and Riemann-Liouville fractional derivative on t1tT.

Definition 5: The piecewise derivative with classical and exponential decay kernel is defined as

0PCFDtαh(t)={h(t),if 0tt1t1CFDtαh(t),if t1tT(6)

and

0PCFDtαh(t)={h(t),if 0tt1t1CFRDtαh(t),if t1tT(7)

where 0PCFDtα presents classical derivative on 0≤ tt1 and Caputo-Fabrizio fractional derivative on t1tT.

Definition 6: The piecewise derivative with classical and Mittag-Leffler kernel is given as

0PABDtαh(t)={h(t),if 0tt1t1ABCDtαh(t),if t1tT(8)

where 0PABDtα presents classical derivative on 0≤ tt1 and Atangana-Baleanu fractional derivative on t1tT.

Definition 7: Let h be continuous and α > 0 then a piecewise integral of h is given as

PPLJtαh(t)={0t1h(τ)dτ,if 0tt11Γ(α)t1t(tτ)α1h(τ)dτ,if t1tT(9)

where PPLJtαh(t) presents classical integral on 0≤ tt1 and the integral with power-law kernel on t1tT.

Definition 8: Let h be continuous and α > 0 then a piecewise integral of h is given as

PCFJtαh(t)={0t1h(τ)dτ,if 0tt11αM(α)h(t)+αM(α)t1th(τ)dτ,if t1tT(10)

where PCFJtαh(t) presents classical integral on 0≤ tt1 and Caputo-Fabrizio integral on t1tT.

2  Tuberculosis Epidemic Model

In this section, we take into account the following piecewise model of tuberculosis:

dS(t)dt=λβ1S(t)I1(t)β2S(t)I2(t)μS(t),dE(t)dt=β1p1S(t)I1(t)β2q1S(t)I2(t)(μ+γ)E(t),dI1(t)dt=pβ1S(t)I1(t)+qβ2S(t)I2(t)+γE(t)(ϕ+μ+δ1)I1(t),dI1(t)dt=ϕ(1r1)I1(t)(μ+δ2)I2(t)φr2I2(t).(11)

The initial conditions are taken as follows:

S(0)=S0,E(0)=E0,I1(0)=A0,I2(0)=I0.(12)

But we noted that the model was considered with its classical version in paper [13] before. Now, we give the meanings of the parameters of model considered in this paper given by Table 1 below:

images

2.1 Second Derivative of Lyapunov Function and Strength Number

Lyapunov function formulation has been used in different analyses in different fields in the last past year. In epidemiology, this function has been used to determine the stability analysis of an epidemiological model. It has been reported that the Lyaponuv can be viewed as energy; therefore, a sign of the first derivative of the function can be useful for the determination of stability. Nevertheless, the sign of the first derivative of a function may not be enough to define whether we have a local maximum or local minimum. On this note, it was suggested that the sign of the second derivative should also be studied. In this section, we shall proceed with such analysis to determine the sign of our model’s second derivative of the associate Lyaponuv function.

In this section, we present the second derivative of Lyapunov function for the model [214].

dS(t)dt=λβ1S(t)I1(t)β2S(t)I2(t)μS(t),dE(t)dt=β1p1S(t)I1(t)β2q1S(t)I2(t)(μ+γ)E(t),dI1(t)dt=pβ1S(t)I1(t)+qβ2S(t)I2(t)+γE(t)(ϕ+μ+δ1)I1(t),dI2(t)dt=ϕ(1r1)I1(t)(μ+δ2)I2(t)φr2I2(t).(13)

Now we find second derivative of Lyapunov function for model with following equality:

L¨=dL˙dt=ddt{(1SS)S˙+(1EE)E˙+(1I1I1)I˙1+(1I2(t)I2(t))I˙2,=(S˙S)2S+(E˙E)2E+(I˙1I1)2I1+(I˙2I2)2I2+(1SS)S¨+(1EE)E¨+(1I1I1)I¨1+(1I2I2)I¨2.(14)

Here second derivatives of classes are given as below:

S¨(t)=β1(S˙(t)I1(t)+I1˙(t)S(t))β2(S˙(t)I2(t)+I2.(t)S(t))μS.(t),E¨(t)=β1p1(S.(t)I1(t)+I˙1(t)S(t))+β2q1(S.(t)I2(t)+I2.(t)S(t))(μ+γ)E.(t),I¨1(t)=β1p(S.(t)I1(t)+I˙1(t)S(t))+β2q(S˙(t)I2(t)+I˙2(t)S(t))+γE˙(t)(ϕ+μ+δ1)I˙1(t),I¨2(t)=ϕ(1r1)I˙1(t)(μ+δ2)I˙2(t)φr2I˙2(t).(15)

Then we have

dL.dt=(S.S)2S+(E.E)2E+(I˙1I1)2I1+(I2.I2)2I2+(1SS){β1(S.(t)I1(t)+I˙1(t)S(t))β2(S.(t)I2(t)+I2.(t)S(t))μS.(t)}+(1EE){β1p1(S.(t)I1(t)+I1.(t)S(t))+β2q1(S.(t)I2(t)+I2.(t)S(t))(μ+γ)E.(t)}+(1I1I1){β1p(S.(t)I1(t)+I˙1(t)S(t))+β2q(S.(t)I2(t)+I2.(t)S(t))+γE.(t)(ϕ+μ+δ1)I1.(t)}+(1I2I2){ϕ(1r1)I˙1(t)(μ+δ2)I2.(t)φr2I2.(t)}.(16)

dL.dt=L.(S,E,I1,I2)+S.(t){β1I1(t)β2I2(t)μEEβ1p1I1(t)EEβ2q1I2(t)I1β1pI1I1β2qI2(t)+β1p1I1(t)+β2q1I2(t)+β1pI1(t)+β2qI2(t)+SSβ1I1(t)+SSβ2I2(t)+SSμ}+E.(t){(μ+γ)I1I1γ+γ}+I˙1(t){β1S(t)(ϕ+μ+δ1)EEβ1p1S(t)I1I1β1pS(t)I2I2ϕ(1r1)+β1p1S(t)+β1pS(t)+ϕ(1r1)+Sβ1+I1I1(ϕ+μ+δ1)}+I2.(t){β2S(t)(μ+δ2)φr2EEβ2q1S(t)I1I1β2qS(t)+β2q1S(t)+β2qS(t)+Sβ2+I2I2(μ+δ2)+I2I2φr2}.(17)

Now replacing S.(t), E.(t), I1.(t), and I2.(t) by their respective formula with their positive and negative parts, we have

d2Ldt2=L.(S,E,I1,I2)+Π+Π,d2Ldt2=Π1L.(S,E,I1,I2)+Π+Π2Π(18)

where

Π+=λ(β1p1I1(t)+β2q1I2(t)+β1pI1(t)+β2qI2(t)+SSβ1I1(t)+SSβ2I2(t)+SSμ)+(β1S(t)I1(t)+β2S(t)I2(t)+μS(t)).(β1I1(t)+β2I2(t)+μ+EEβ1p1I1(t)+EEβ2q1I2(t)+I1β1p+I1I1β2qI2(t))+(β1p1S(t)I1(t)+β2q1S(t)I2(t))γ+((μ+γ)E(t))((μ+γ)+I1I1)+(β1pI1(t)S(t)+β2qI2(t)S(t)+γE(t))(β1p1S(t)+β1pS(t)+ϕ(1r1)+Sβ1+I1I1((ϕ+μ+δ1)))+(ϕ+μ+δ1)I1(t)(β1S(t)+(ϕ+μ+δ1)+EEβ1p1S(t)+I1I1β1pS(t)+I2I2ϕ(1r1))+ϕ(1r1)I1(t)(β2q1S(t)+β2qS(t)+Sβ2+I2I2(μ+δ2)+I2I2φr2+(φr2I2)(β2S(t)+(μ+δ2)+φr2+EEβ2q1S(t)+I1I1β2qS(t))).(19)

Π=λ(β1I1(t)+β2I2(t)+μ+EEβ1p1I1(t)+EEβ2q1I2(t)+β1pI1(t)+I1I1β2qI2(t))+(β1S(t)I1(t)+β2S(t)I2(t)+μS(t)).(β1p1I1(t)+β2q1I2(t)+β1pI1(t)+β2qI2(t)SSβ1I1(t)+SSβ2I2(t)+SSμ)+(β1p1S(t)I1(t)+β2q1S(t)I2(t))((μ+γ)I1I1γ)γ(μ+γ)E(t)+(ϕ+μ+δ1)I1(t)(β1p1S(t)+I1I1(ϕ+μ+δ1)+β1pS(t)+β1S+ϕ(1r1))+(β1pS(t)I1(t)+β2qS(t)I2(t))+γE(t)(β1S(t)+(ϕ+μ+δ1)+EEβ1p1S(t)+I1I1β1pS(t)+I2I2ϕ(1r1))+ϕ(1r1)I1(t)(β2S(t)+(μ+δ2)+φr2+EEβ2q1S(t)+I1I1β2qS(t)+(φr2I2)(β2q1S(t)+qβ2S(t)+Sβ2+I2I2(μ+δ2)+I2I2φr2)).

Now we can easly put following results for obtained results above:

d2Ldt2=Π1Π2.(20)

Then

If Π1>Π2 then d2Ldt2>0,If Π1<Π2 then d2Ldt2<0,If Π1=Π2 then d2Ldt2=0.(21)

So, the interpretation associated the sign of second order.

2.2 Strength Number

Without a doubt, the reproductive number has been utilized as a powerful mathematical tool to the stability of a mathematical model for a given infectious disease. While it has been used with some success, it has also been criticized as an insufficient tool to predict the behavior of the spread. For example, it was pointed out that there are several ways to obtain this value on the other hand. However, it was also argued that this value could not help humans t determine whether the model will determine waves. The concept of strength number has been suggested to further the analysis and will be used in this section.

The component FA is obtained with deriving the nonlinear part of the infected classes. In our model there are two infected classes named by I1 and I2. These infected classes given by

I˙1=pβ1S(t)I1(t)+qβ2S(t)I2(t)+γE(t)(ϕ+μ+δ1)I1(t),I2.=ϕ(1r1)I1(t)(μ+δ2)I2(t)φr2I2(t).(22)

But we only use nonlinear part of infected classes. So we use I2. classes. Nonlinear part of I˙1 classes is given by

=S(t)I1(t)+S(t)I2(t),(23)

I1(S(t)I1(t)+S(t)I2(t))=S(t),2I12(S(t))=0.(24)

In this case, we can have the following:

FA=[00].(25)

Then

det(FAV1λI)=0,(26)

leads to

A0=0.(27)

A0 means there is no strength. Also there are more conlusion when strengh is zero

1) The disease will spread with a constant speed.

2) The disease will not renewal process therefore no new wave will be expected.

3) The magnitude of the spread will be the same at all time until extinction.

3  Applications of Piecewise Derivative

3.1 A Mathematical Model of Tuberculosis Epidemic Model with Piecewise Modeling

In this section, we present some applications of piecewise derivative for tuberculosis epidemic model such as

{dS(t)dt=λβ1S(t)I1(t)β2S(t)I2(t)μS(t),dE(t)dt=β1p1S(t)I1(t)β2q1S(t)I2(t)(μ+γ)E(t),dI1(t)dt=pβ1S(t)I1(t)+qβ2S(t)I2(t)+γE(t)(ϕ+μ+δ1)I1(t),dI2(t)dt=ϕ(1r1)I1(t)(μ+δ2)I2(t)φr2I2(t).if0tW1,(28)

{t1CDtαS(t)=λβ1S(t)I1(t)β2S(t)I2(t)μS(t),t1CDtαE(t)=β1p1S(t)I1(t)β2q1S(t)I2(t)(μ+γ)E(t),t1CDtαI1(t)=pβ1S(t)I1(t)+qβ2S(t)I2(t)+γE(t)(ϕ+μ+δ1)I1(t),t1CDtαI2(t)=ϕ(1r1)I1(t)(μ+δ2)I2(t)φr2I2(t).if W1tW2,(29)

{dS(t)=[λβ1S(t)I1(t)β2S(t)I2(t)μS(t)]dt+σ1SdB1(t),dE(t)=[β1p1S(t)I1(t)β2q1S(t)I2(t)(μ+γ)E(t)]dt+σ2EdB2(t),dI1(t)=[pβ1S(t)I1(t)+qβ2S(t)I2(t)+γE(t)(ϕ+μ+δ1)I1(t)]dt+σ3I1dB3(t),dI2(t)=[ϕ(1r1)I1(t)(μ+δ2)I2(t)φr2I2(t)]dt+σ4I2dB4(t).if W2tW.(30)

Let us give necessary conditions for the existence and uniqueness, we must prove that [0, W1] and [W1, W2] fi(S, E, I1, I2) for i = 1, 2, 3, 4 satisfy

1) Linear growth condition

|fi(xi,t)|2ki(1+|xi|2) for i=1,2,3,4.(31)

and

2) The Lipschitz condition

|fi(xi1,t)fi(xi2,t)|2k¯i|xi1xi2|2  for i=1,2,3,4.(32)

Now we define the norm φ=tDφsup|φ(t)|. Now we put forth the existence and uniqueness of the solution piecewisely for [0,W2]. For [0,W2], there exist 4 positive constant M1, M2,M3 and M4 < such that

S <M1, E <M2, I1 <M3, I2 <M4.(33)

Let us write system as below:

{S.=f1(S,E,I1,I2),E.=f2(S,E,I1,I2),I˙1=f3(S,E,I1,I2),I2.=f4(S,E,I1,I2).if 0tW2.(34)

For proof, we consider the function

|f1(S,E,I1,I2)|2=|λβ1S(t)I1(t)β2S(t)I2(t)μS(t)|2,4λ2+4|β1S(t)I1(t)|2+4|β2S(t)I2(t)|+4|μS(t)|2,4λ2+4|β1I1(t)|2|S(t)|2+4|β2I2(t)||S(t)|2+4μ2|S(t)|2,4λ2+4(β1I1(t))2|S(t)|2+4(β2I2(t))2|S(t)|2+4μ2|S(t)|2,4λ2(1+(β1I1(t))2+(β2I2(t))2+μ2λ2|S(t)|2).(35)

under the condition that

(β1I1(t))2+(β2I2(t))2+μ2λ2<1,(36)

then we have

|f1(S,E,I1,I2)|2k1(1+|S(t)|2).(37)

Using same routine

|f2(S,E,I1,I2)|2=|β1p1S(t)I1(t)β2q1S(t)I2(t)(μ+γ)E(t)|2,3|β1p1S(t)I1(t)|2+3|β2q1S(t)I2(t)|2+3|(μ+γ)E(t)|2,3t[0,T2]sup|β1p1S(t)I1(t)|2+3t[0,T2]sup|β2q1S(t)I2(t)|2+3|(μ+γ)E(t)|2,3(β1p1S(t)I1(t))2+3(β2q1S(t)I2(t))2+3(μ+γ)2|E(t)|2,3(β1p1S(t)I1(t))2+3(β2q1S(t)I2(t))2.(1+3(μ+γ)23(β1p1S(t)I1(t))2+3(β2q1S(t)I2(t))2|E(t)|2),(38)

under the condition

(μ+γ)2(β1p1S(t)I1(t))2+(β2q1S(t)I2(t))2<1, then (39)

|f2(S,E,I1,I2)|2k2(1+|E(t)|2).(40)

For the function f3(S,E,I1,I2)

|f3(S,E,I1,I2)|2=|pβ1S(t)I1(t)+qβ2S(t)I2(t)+γE(t)(ϕ+μ+δ1)I1(t)|2,4|pβ1S(t)I1(t)|2+4|qβ2S(t)I2(t)|2+4|γE(t)|2+4|(ϕ+μ+δ1)I1(t)|2,4|pβ1S(t)|2|I1(t)|2+4|qβ2S(t)I2(t)|2+4|γE(t)|2+4(ϕ+μ+δ1)2|I1(t)|2,4t[0,T2]sup|pβ1S(t)|2|I1(t)|2+4t[0,T2]sup|qβ2S(t)I2(t)|2+4t[0,T2]sup|γE(t)|2+4(ϕ+μ+δ1)2|I1(t)|2,4(pβ1S(t))2|I1(t)|2+4(qβ2S(t)I2(t))2+4(γE(t))2+4(ϕ+μ+δ1)2|I1(t)|2,4(qβ2S(t)I2(t))2+4(γE(t))2(1+(pβ1S(t))2+(ϕ+μ+δ1)2(qβ2S(t)I2(t))2+(γE(t))2|I1(t)|2)(41)

under the condition

(pβ1S(t))2+(ϕ+μ+δ1)2(qβ2S(t)I2(t))2+(γE(t))2<1, then (42)

|f3(S,E,I1,I2)|2k3(1+|I1(t)|2).(43)

Finally for the function f4(S,E,I1,I2)

|f4(S,E,I1,I2)|2=|ϕ(1r1)I1(t)(μ+δ2)I2(t)φr2I2(t)|2,3|ϕ(1r1)I1(t)|2+3|(μ+δ2)I2(t)|2+3|φr2I2(t)|2,3ϕ2(1r1)|I1(t)|2+3(μ+δ2)2|I2(t)|2+(φr2)23|I2(t)|2,3ϕ2(1r1)t[0,T2]sup|I1(t)|2+3(μ+δ2)2|I2(t)|2+(φr2)23|I2(t)|2,3ϕ2(1r1)I12(t)+3(μ+δ2)2|I2(t)|2+(φr2)23|I2(t)|2,3ϕ2(1r1)I12(t)(1+(μ+δ2)2+(φr2)2ϕ2(1r1)I12(t)|I2(t)|2),(44)

under the condition

(μ+δ2)2+(φr2)2ϕ2(1r1)I12(t)<1, then (45)

|f4(S,E,I1,I2)|2k4(1+|I2(t)|2).(46)

Therefore the condition of linear growth is verified if

max{(β1I1(t))2+(β2I2(t))2+μ2λ2,(μ+γ)2(β1p1S(t)I1(t))2+(β2q1S(t)I2(t))2,(pβ1S(t))2+(ϕ+μ+δ1)2(qβ2S(t)I2(t))2+(γE(t))2,(μ+δ2)2+(φr2)2ϕ2(1r1)I12(t)}<1.(47)

Now we have to verify Lipschitz condition for equations.

For the function f1(S,E,I1,I2),

|f1(S,E,I1,I2)f1(S,E,I1,I2)|(β1I1(t)+β2I2(t)+μ)|SS|,k¯1|SS|.(48)

For the function f2(S,E,I1,I2),

|f2(S,E,I1,I2)f2(S,E,I1,I2)|(γ+μ)|EE|,k¯2|EE|.(49)

For the function f3(S,E,I1,I2),

|f3(S,E,I1,I2)f3(S,E,I1,I2)|(pβ1S(t)+(ϕ+μ+δ1))|I1I1|,k¯3|I1I1|.(50)

Finally for the function f4(S,E,I1,I2),

|f4(S,E,I1,I2)f4(S,E,I1,I2)|(μ+δ2+φr2)|I2I2|,k¯4|I2I2|.(51)

We verified the Lipschitz condition which completes the proof.

Let us do proof for last part of piecewise equation. Here we consider for ∀ t[ W2, W]. In the model we take for S(t), E(t), I1(t) I2(t)[W2,τe), where τe shows explosion time. To prove the solution is global, one has to prove that such system solution is global, so we have to prove that τe = .

Now we consider l0R+ is a positive constant such that S(W2), E(W2), I1(W2) I2(W2) lies within [1l0,l0]. We define a stopping time as

τl={t[W2,τe):1lmin{S(t),E(t),I1(t),I2(t)} or max{S(t),E(t),I1(t)I2(t)}l},(52)

for each ll0. While as l, τl is monotonically increasing. llimτl=τ with τeτ. ∀ t≥0, if we show that τ = 0, then we can conclude that τe = and S(t),E(t),I1(t),I2(t)R+4 is solution. So we have to prove that τe = .

If we have contradictory for the conclusion, then there exists 0 < W and ɛ ∈ (0, 1) such that

P{Wτ}>ε.(53)

Now we define a function H(X):R+4R+ in HC2 such that

H¯(X)= dH(X)=j=14(11xj)dxj+j=14σj(xj1)dBj(t),   =j=14(11xj)xj+j=14σj(xj1)dBj(t),(54)

where

x1=S(t),x2=E(t),x3=I1(t),x4=I2(t),σj=(σ1,σ2,σ3,σ4),Bj(t)=(B1(t),B2(t),B3(t),B4(t)).(55)

For our model H¯(X) is obtained by following equality:

H¯(X)=j=14(11xj)xj=(11S)S+(11E)E+(11I1)I1+(11I2)I2+j=14σj22.(56)

H¯(X)=λ+β1I1(t)+β2I2(t)+μ+β1p1S(t)I1(t)+β2q1S(t)I2(t)+(μ+γ)+pβ1S(t)I1(t)+β2qS(t)I2(t)+γE(t)+(ϕ+μ+δ1)+ϕ(1r1)I1(t)+(μ+δ2)I2(t)+φr2I2(t){β1S(t)I1(t)+β2S(t)I2(t)+μS(t)+λS+(μ+γ)E(t)+1Eβ1p1S(t)I1(t)+1Eβ2q1S(t)I2(t)+(ϕ+μ+δ1)I1(t)+pβ1S(t)+1I1β2qS(t)I2(t)+1I1γE(t)+(μ+δ2)I2(t)+φr2I2(t)+1I2%ϕ(1r1)I1(t)+j=14σj22,<λ+3μ+γ+ϕ+δ1=θ¯(57)

and

H¯(X)=θ¯dt+j=14σj(xj1)dBj(t).(58)

By taking integration from 0 to τlW, we have

E[H¯(τlX)]H¯(X(W2))+E[0τlWθ¯],H¯(X(W2))+θ¯W. (59)

Setting Πl={W>τl} for l1l and thus P(Πl)ζ.

Notting that for ∀ΩΠl, there must exist at least one X(τl,w) which is equal to 1l or l. Then l −log l −1 or 1l+logl1 as result

(1l+logl1)E(llogl1)<H¯(X(τl)).(60)

From above, we can write

H¯(X(W2))+θ¯W>E(1ΠlH¯(X(τl))),σ[(llogl1)(1l+logl1)].(61)

Here 1Πl is the indicator function of Π. Thus llim leads

>H¯(X(W2))+θ¯W=0.(62)

It is a contradiction. So under the conditions gived earlier τ = which completes the proof.

4  Numerical Schemes for Model with Four Waves Patterns

In this section, we generate a numerical schemes for spread of infectious (specially for pandemic) disease with four patterns. These schemes will consist of three derivatives with randomness [1, 12].

4.1 Case 1: Classical-Power Law-Exponential Decay Law-Randomness

In this case, we consider a version with four waves which have classical derivative starts from 0 to W1, the power law derivative start from W1 to W2, the exponential decay law derivative start from W2 to W3, and the last from W3 to W. So a piecewise mathematical system that is defined as subsection can be given as

{dyidt=g(t,y),if 0tW1yi(0)=yi,0,i=1,2,,nt1CDtαyi=g(t,y),if W1tW2yi(W1)=yi,1,0<α1,i=1,2,,nt2CFDtαyi=g(t,y),if W2tW3yi(W2)=yi,2,0<α1i=1,2,,ndy(t)=g(t,y)dt+σiyidBi(t),if W3tWyi(W3)=yi,3,i=1,2,,n(63)

where σi are densities of randomness and Bi are the functions of noise.

4.2 Case 2: Classical-Power Law-Mittag-Leffler Law-Randomness

In this case, we consider a version with four waves which have classical derivative starts from 0 to W1, the power law derivative start from W1 to W2, the Mittag-Leffler law derivative start from W2 to W3, and the last from W3 to W. So a piecewise mathematical system that is defined as subsection can be given as

{dyidt=g(t,y),if 0tW1yi(0)=yi,0,i=1,2,,nt1CDtαyi=g(t,y),if W1tW2yi(W1)=yi,1,0<α1,i=1,2,,nt2ABCDtαyi=g(t,y),if W2tW3yi(W2)=yi,2,0<α1i=1,2,,ndy(t)=g(t,y)dt+σiyidBi(t),if W3tWyi(W3)=yi,3i=1,2,,n,(64)

where σi are densities of randomness and Bi are the functions of noise.

4.3 Case 3: Classical-Power Law-Fractal-Fractional Power Law Derivative-Randomness

In this case, we consider a version with four waves which have classical derivative starts from 0 to W1, the power law derivative start from W1 to W2, fractal-fractional power law derivative start from W2 to W3, and the last from W3 to W. So a piecewise mathematical system that is defined as subsection can be given as

{dyidt=g(t,y),if 0tW1yi(0)=yi,0,i=1,2,,nt1CDtαyi=g(t,y),if W1tW2yi(W1)=yi,1,0<α1,i=1,2,,nt2FFPDtα,βyi=g(t,y),if W2tW3yi(W2)=yi,2,0<α1i=1,2,,ndy(t)=g(t,y)dt+σiyidBi(t),if W3tWyi(W3)=yi,3i=1,2,,n,(65)

where σi are densities of randomness and Bi are the functions of noise.

4.4 Case 4: Classical-Exponential Decay Law-Fractal-Fractional Exponential Decay Law Derivative-Randomness

In this case, we consider a version with four waves which have classical derivative starts from 0 to W1, the exponential decay law derivative start from W1 to W2, fractal-fractional exponential decay law derivative start from W2 to W3, and the last from W3 to W. So a piecewise mathematical system that is defined as subsection can be given as

{dyidt=g(t,y),if 0tW1yi(0)=yi,0,i=1,2,,nt1CFDtαyi=g(t,y),if W1tW2yi(W1)=yi,1,0<α1,i=1,2,,nt2FFEDtα,βyi=g(t,y),if W2tW3yi(W2)=yi,2,0<α1i=1,2,,ndy(t)=g(t,y)dt+σiyidBi(t),if W3tWyi(W3)=yi,3i=1,2,,n(66)

where σi are densities of randomness and Bi are the functions of noise.

4.5 Case 5: Classical-Mittag-Leffler Law-Fractal-Fractional Mittag-Leffler Law Derivative-Randomness

In this case, we consider a version with four waves which have classical derivative starts from 0 to W1, the Mittag Leffler law derivative start from W1 to W2, fractal-fractional Mittag-Leffler law derivative start from W2 to W3, and the last from W3 to W. So a piecewise mathematical system that is defined as subsection can be given as

{dyidt=g(t,y),if 0tW1yi(0)=yi,0,i=1,2,,nt1ABCDtαyi=g(t,y),if W1tW2yi(W1)=yi,1,0<α1,i=1,2,,nt2FFMDtα,βyi=g(t,y),if W2tW3yi(W2)=yi,2,0<α1i=1,2,,ndy(t)=g(t,y)dt+σiyidBi(t),if W3tWyi(W3)=yi,3i=1,2,,n(67)

where σi are densities of randomness and Bi are the functions of noise.

5  Numerical Schemes of Piecewise Epidemic Disease Models with Four Waves Patterns

In this section we assumed that those kind of epidemic models satisfy existence and uniqueness. So we can put numerical solutions for them. While putting solution results we use in all cases on the Lagrange polynomial interpolation. First we divide [0, W] in four

0t0t1...tn1=W1tn1+1tn1+2...tn2=W2tn2+1tn2+2...tn3=W3tn3+1tn3+2...tn3=W.(68)

5.1 Numerical Method for Case 1:

Let us consider the first case

{dyidt=g(t,y),if 0tW1yi(0)=yi,0,i=1,2,,nt1CDtαyi=g(t,y),if W1tW2yi(W1)=yi,1,0<α1,i=1,2,,nt2CFDtαyi=g(t,y),if W2tW3yi(W2)=yi,2,0<α1,i=1,2,,ndy(t)=g(t,y)dt+σiyidBi(t),if W3tWyi(W3)=yi,3i=1,2,,n(69)

The numerical solution can then be provided as

{yin1=yi(0)+k1=i+1n1{3Δt2g(tk1,y(tk1))g(tk11,y(tk11))Δt2},  0tW1yin2=yi(W1)+(Δt)αΓ(α+2)k2=0n2g(tk2,y(tk2))[(n2k2+1)α(n2k2+2+α)(n2k2)α(n2k2+2+2α)]  W1tW2 (Δt)αΓ(α+2)k2=0n2g(tk21,y(tk21))×[(n2k2+1)α+1(n2k2)α(n2k2+1+α)],yin3=yi(W2)+1αM(α)k3=0n3[g(tk3,y(tk3))g(tk31,y(tk31))]+αM(α)k3=0n3{3Δt2g(tk3,y(tk3))g(tk31%,y(tk31))Δt2}W2tW3yin4=yi(W3)+%k4=i+1n4{3Δt2g(tk4,y(tk4))g(tk41,y(tk41))Δt2} W3tWk4=i+1n4{σ2(y(tk4+1)+y(tk4))(B(tk4+1)B(tk4))}.(70)

5.2 Numerical Method for Case 2:

We deal with the following problem with second case

{dyidt=g(t,y),if 0tT1yi(0)=yi,0,i=1,2,,nt1CDtαyi=g(t,y),if T1tT2yi(T1)=yi,1,0<α1,i=1,2,,nt2ABCDtαyi=g(t,y),if T2tT3yi(T2)=yi,2,0<α1,i=1,2,,ndy(t)=g(t,y)dt+σiyidBi(t),if T3tTyi(T3)=yi,3i=1,2,,n.(71)

The numerical solution for such problem is given by

{yin1=yi(0)+k1=i+1n1{3Δt2g(tk1,y(tk1))g(tk11,y(tk11))Δt2},  0tW1yin2=yi(W1)+(Δt)αΓ(α+2)k2=0n2g(tk2,y(tk2))[(n2k2+1)α(n2k2+2+α)(n2k2)α(n2k2+2+2α)]  W1tW2 (Δt)αΓ(α+2)k2=0n2g(tk21,y(tk21))×[(n2k2+1)α+1(n2k2)α(n2k2+1+α)],yin3=yi(W2)+1αAB(α)f(tk3,y(tk3))+α(Δt)αAB(α)Γ(α+2)k3=0n3g(tk3,y(tk3))×[(n3k3+1)α(n3k3+2+α)(n3k3)α(n3k3+2+2α)] W2tW3α(Δt)αAB(α)Γ(α+2)k3=0n3g(tk31,y(tk31))×[(n3k3+1)α+1(n3k3)α(n3k3+1+α)],  yin4=yi(W3)+k4=i+1n4{3Δt2g(tk4,y(tk4))g(tk41,y(tk41))Δt2} W3tWk4=i+1n4{σ2(y(tk4+1)+y(tk4))(B(tk4+1)B(tk4))}.(72)

5.3 Numerical Method for Case 3:

Now we deal with the following problem with third case

{dyidt=g(t,y),if 0tW1yi(0)=yi,0,i=1,2,,nt1CDtαyi=g(t,y),if W1tW2yi(W1)=yi,1,0<α1,i=1,2,,nt2FFPDtα,βyi=g(t,y),if W2tW3yi(W2)=yi,2,0<α1,i=1,2,,ndy(t)=g(t,y)dt+σiyidBi(t),if W3tWyi(W3)=yi,3i=1,2,,n.(73)

The numerical solution for such problem is given by

{yin1=yi(0)+k1=i+1n1{3Δt2g(tk1,y(tk1))g(tk11,y(tk11))Δt2},  0tW1yin2=yi(W1)+(Δt)αΓ(α+2)k2=0n2g(tk2,y(tk2))×[(n2k2+1)α(n2k2+2+α)(n2k2)α(n2k2+2+2α)]  W1tW2 (Δt)αΓ(α+2)k2=0n2g(tk21,y(tk21))×[(n2k2+1)α+1(n2k2)α(n2k2+1+α)],yin3=yi(W2)+β(Δt)αΓ(α+2)k3=0n3tk3β1g(tk3,y(tk3))×[(n3k3+1)α(n3k3+2+α)(n3k3)α(n3k3+2+2α)] W2tW3+β(Δt)αΓ(α+2)k3=0n3tk31β1g(tk31,y(tk31))×[(n3k3+1)α+1(n3k3)α(n3k3+1+α)],yin4=yi(W3)+k4=i+1n4{3Δt2g(tk4,y(tk4))gg(tk41,y(tk41))Δt2} W3tWk4=i+1n4{σ2(y(tk4+1)+y(tk4))(B(tk4+1)B(tk4))}. (74)

5.4 Numerical Method for Case 4:

Here, we deal with the following problem with fourth case

{dyidt=g(t,y),if 0tW1yi(0)=yi,0,i=1,2,,nt1CFDtαyi=g(t,y),if W1tW2yi(W1)=yi,1,0<α1,i=1,2,,nt2FFEDtα,βyi=g(t,y),if W2tW3yi(W2)=yi,2,0<α1i=1,2,,ndy(t)=g(t,y)dt+σiyidBi(t),if W3tW,yi(W3)=yi,3i=1,2,,n.(75)

The numerical solution for such problem is given by

{yin1=yi(0)+k1=i+1n1{3Δt2g(tk1,y(tk1))g(tk11,y(tk11))Δt2},  0tW1yin2=yi(W1)+1αM(α)k2=0n2[g(tk2,y(tk2))g(tk21,y(tk21))]  W1tW2 +αM(α)k2=0n2{3Δt2g(tk2,y(tk2))g(tk21,y(tk21))Δt2},yin3=yi(W2)+1αM(α)k3=0n3[βtk3β1g(tk3,y(tk3))βtk31β1g(tk31,y(tk31))] W2tW3+αM(α)k3=0n3{3Δt2βtk3β1g(tk3,y(tk3))βtk31β1g(tk31,y(tk31))Δt2},yin4=yi(W3)+k4=i+1n4{3Δt2g(tk4,y(tk4))g(tk41,y(tk41))Δt2} W3tWk4=i+1n4{σ2(y(tk4+1)+y(tk4))(B(tk4+1)B(tk4))}. (76)

6  Numerical Method for Case 5:

Finally, we give numerical method with the following problem with fifth case:

{dyidt=g(t,y),if 0tW1yi(0)=yi,0,i=1,2,,nt1ABCDtαyi=g(t,y),ifW1tW2yi(W1)=yi,1,0<α1,i=1,2,,nt2FFMDtα,βyi=g(t,y),ifW2tW3yi(W2)=yi,2,0<α1i=1,2,,ndy(t)=g(t,y)dt+σiyidBi(t),ifW3tWyi(W3)=yi,3i=1,2,,n.(77)

The numerical solution for such problem is given by

{yin1=yi(0)+k1=i+1n1{3Δt2g(tk1,y(tk1))g(tk11,y(tk11))Δt2}, 0tW1yin2=yi(W1)+1αAB(α)g(tn2,y(tn2))+α(Δt)αAB(α)Γ(α+2)k2=0n2g(tk2,y(tk2))×[(n2k2+1)α(n2k2+2+α)(n2k2)α(n2k2+2+2α)] W1tW2α(Δt)αAB(α)Γ(α+2)k2=0n2g(tk21,y(tk21))×[(n2k2+1)α+1(n2k2)α(n2k2+1+α)],yin3=yi(W2)+1αAB(α)tn3β1g(tn3,y(tn3)) W2tW3+αβ(Δt)αAB(α)Γ(α+2)k3=0n3tk3β1g(tk3,y(tk3))×[(n3k3+1)α(n3k3+2+α)(n3k3)α(n3k3+2+2α)]+αβ(Δt)αAB(α)Γ(α+2)k3=0n3tk31β1g(tk31,y(tk31))×[(n3k3+1)α+1(n3k3)α(n3k3+1+α)],yin4=yi(W3)+k4=i+1n4{3Δt2g(tk4,y(tk4))g(tk41%,y(tk41))Δt2} W3tW+k4=i+1n4{σ2(y(tk4+1)+y(tk4))(B(tk4+1)B(tk4))},(78)

6.1 Numerical Simulation for Stochastic-Deterministic Model of Tuberculosis

In this section, we give numerical simulation of the Tuberculosis epidemic system of fractional stochastic differential equations. We have made use of the model with the piecewise differential operators and the numerical scheme where the Lagrange polynomial interpolation is used. While modelling with piecewise idea, the first part is classical, the second part is fractional and last part is stochastic. The numerical simulation is performed for different values of fractional orders. So the stochastic-deterministic piecewise tuberculosis model is given as

{dS(t)dt=λβ1S(t)I1(t)β2S(t)I2(t)μS(t),dE(t)dt=β1p1S(t)I1(t)β2q1S(t)I2(t)(μ+γ)E(t),dI1(t)dt=pβ1S(t)I1(t)+qβ2S(t)I2(t)+γE(t)(ϕ+μ+δ1)I1(t),dI2(t)dt=ϕ(1r1)I1(t)(μ+δ2)I2(t)φr2I2(t),S(0)=S0,E(0)=E0,I1(0)=I10,I2(0)=I20,if 0tW1(79)

{t1CDtαS(t)=λβ1S(t)I1(t)β2S(t)I2(t)μS(t),t1CDtαE(t)=β1p1S(t)I1(t)β2q1S(t)I2(t)(μ+γ)E(t),t1CDtαI1(t)=pβ1S(t)I1(t)+qβ2S(t)I2(t)+γE(t)(ϕ+μ+δ1)I1(t),t1CDtαI2(t)=ϕ(1r1)I1(t)(μ+δ2)I2(t)φr2I2(t),S(W1)=S1,E(W1)=E1,I1(W1)=I11,I2(W1)=I21,if W1tW2(80)

{dS(t)=[λβ1S(t)I1(t)β2S(t)I2(t)μS(t)]dt+σ1SdB1(t),dE(t)=[β1p1S(t)I1(t)β2q1S(t)I2(t)(μ+γ)E(t)]dt+σ2EdB2(t),dI1(t)=[pβ1S(t)I1(t)+qβ2S(t)I2(t)+γE(t)(ϕ+μ+δ1)I1(t)]dt+σ3I1dB3(t),dI2(t)=[ϕ(1r1)I1(t)(μ+δ2)I2(t)φr2I2(t)]dt+σ4I2dB4(t),S(W2)=S2,E(W2)=E2,I1(W2)=I12,I2(W2)=I22,if W2tW.(81)

For simplicity we consider right side of system as

{S.=f1(S,E,I1,I2),E.=f2(S,E,I1,I2),I˙1=f3(S,E,I1,I2),I2.=f4(S,E,I1,I2),(82)

Using the numerical scheme presented in this paper with piecewise derivative, the numerical solution of the stochastic-deterministic tuberculosis model is given as follows:

{Sin1=Si(0)+k1=i+1n1{3Δt2f1(tk1,S(tk1))f1(tk11,S(tk11))Δt2},  0tW1Sin2=Si(W1)+(Δt)αΓ(α+2)k2=0n2f1(tk2,S(tk2))×[(n2k2+1)α(n2k2+2+α)(n2k2)α(n2k2+2+2α)]  W1tW2 (Δt)αΓ(α+2)k2=0n2f1(tk21,S(tk21))×[(n2k2+1)α+1(n2k2)α(n2k2+1+α)],Sin3=Si(W2)+k3=i+1n3{3Δt2f1(tk3,S(tk3))f1(tk31,S(tk31))Δt2} W2tWk3=i+1n3{σ2(S(tk3+1)+S(tk3))(B(tk3+1)B(tk3))},  % (83)

{Ein1=Ei(0)+k1=i+1n1{3Δt2f2(tk1,E(tk1))-f2(tk1-1,E(tk1-1))Δt2},  0tW1Ein2=Ei(W1)+(Δt)αΓ(α+2)k2=0n2f2(tk2,E(tk2))×[(n2-k2+1)α(n2-k2+2+α)-(n2-k2)α(n2-k2+2+2α)],  T1tT2 -(Δt)αΓ(α+2)k2=0n2f2(tk2-1,E(tk2-1))×[(n2-k2+1)α+1-(n2-k2)α(n2-k2+1+α)]Ein3=Ei(W2)+k3=i+1n3{3Δt2f2(tk3,E(tk3))-f2(tk3-1,E(tk3-1))Δt2}     W2tW.+k3=i+1n3{σ2(E(tk3+1)+E(tk3))(B(tk3+1)-B(tk3))}                  

{I1in1=I1i(0)+%k1=i+1n1{3Δt2f3(tk1,I1(tk1))f3(tk11,I1(tk11))Δt2},  0tW1,I1in2=I1i(W1)+(Δt)αΓ(α+2)k2=0n2f3(tk2,I1(tk2))×[(n2k2+1)α(n2k2+2+α)(n2k2)α(n2k2+2+2α)],  W1tW2 (Δt)αΓ(α+2)k2=0n2f3(tk21,I1(tk21))×[(n2k2+1)α+1(n2k2)α(n2k2+1+α)],I1in3=I1i(W3)+k3=i+1n3{3Δt2f3(tk3,I1(tk3))f3(tk31,I1(tk31))Δt2} W2tW.k3=i+1n3{σ2(I1(tk3+1)+I1(tk3))(B(tk3+1)B(tk3))},

{I2in1=I2i(0)+k1=i+1n1{3Δt2f4(tk1,I2(tk1))f4(tk11,I2(tk11))Δt2},  0tW1I2in2=I2i(W1)+(Δt)αΓ(α+2)k2=0n2f4(tk2,I2(tk2))×[(n2k2+1)α(n2k2+2+α)(n2k2)α(n2k2+2+2α)]  W1tW2(Δt)αΓ(α+2)k2=0n2f4(tk21,I2(tk21))×[(n2k2+1)α+1(n2k2)α(n2k2+1+α)],I2in3=I2i(W2)+k3=i+1n3{3Δt2f4(tk3,I2(tk3))f4(tk31,I2(tk31))Δt2} W2tWk3=i+1n3{σ2(I2(tk3+1)+I2(tk3))(B(tk3+1)B(tk3))}.

7  Numerical Simulations

In this section, we will deal with numerical simulation of the Tuberculosis epidemic system of fractional stochastic differential equations. in order to demonstrate that the proposed method is effective and accurate. We have made use of the model with the piecewise differential operators and the numerical scheme where the Lagrange polynomial interpolation is used. In the numerical scheme, the first part is classical, the second part is fractional and last part is stochastic. We also present the results obtained from the fractional stochastic model, the numerical simulations are shown in Fig. 1 for alpha=1, Fig. 2 for alpha=0.5, Fig. 3 for alpha=0.6 and finally Fig. 4 for alpha=0.9 with density of randomness given by sigma1=0.01, sigma2=0.015, sigma3=0.012, sigma4=0.010. And with same alpha values but different density of randomness given by sigma1=0.1, sigma2=0.2, sigma3=0.3, sigma4=0.4 we put Figs. 58. Also figures including the initial conditions as

S(1)=180,E(1)=130,I1(1)=160,I2(1)=140.

images

Figure 1: Numerical simulation for alpha=1

images

Figure 2: Numerical simulation for alpha=0.5

images

Figure 3: Numerical simulation for alpha=0.6

images

Figure 4: Numerical simulation for alpha=0.9

images

Figure 5: Numerical simulation for alpha=1

images

Figure 6: Numerical simulation for alpha=0.5

images

Figure 7: Numerical simulation for alpha=0.6

images

Figure 8: Numerical simulation for alpha=0.9

8  Conclusion

The spread of tuberculosis within human settlements and has infected and killed millions of humans in the last 200 years. While researchers from all backgrounds have put their efforts together to combat this virus and try to stop its spread, several studies have been performed; however, the virus is still spreading so far. Mathematical models are used to predict the future development of a given real-world problem. While several techniques and models have been proposed, they have not predicted piecewise behaviors of the spread. In this work, we attempted to present a model with piecewise patterns.

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.

References

 1.  Atangana, A., Araz, S. S. (2021). Modeling and forecasting the spread of COVID-19 with stochastic and deterministic approaches: Africa and Europe. Advances in Difference Equations, 2021(1), 1–107. DOI 10.1186/s13662-021-03213-2. [Google Scholar] [CrossRef]

 2.  Atangana, A. (2021). Mathematical model of survival of fractional calculus, critics and their impact: How singular is our world? https://hal.archives-ouvertes.fr/hal-03182789/document. [Google Scholar]

 3.  Nematollahi, M. H., Vatankhah, R., Sharifi, M. (2020). Nonlinear adaptive control of tuberculosis with consideration of the risk of endogenous reactivation and exogenous reinfection. Journal of Theoretical Biology, 486, 110081. DOI 10.1016/j.jtbi.2019.110081. [Google Scholar] [CrossRef]

 4.  Ogundile, O., Edeki, S., Adewale, S. (2018). A modelling technique for controlling the spread of tuberculosis. International Journal of Biology and Biomedical Engineering, 12, 131–136. [Google Scholar]

 5.  Atangana, A., Baleanu, D. (2016). New fractional derivatives with nonlocal and non-singular kernel: Theory and application to heat transfer model. Thermal Science, 20(2), 763. DOI 10.2298/TSCI160111018A. [Google Scholar] [CrossRef]

 6.  Caputo, M., Fabrizio, M. (2017). On the notion of fractional derivative and applications to the hysteresis phenomena. Meccanica, 52(13), 3043–3052. DOI 10.1007/s11012-017-0652-y. [Google Scholar] [CrossRef]

 7.  Caputo, M. (1967). Linear models of dissipation whose Q is almost frequency independent–II. Geophysical Journal International, 13(5), 529–539. DOI 10.1111/j.1365-246X.1967.tb02303.x. [Google Scholar] [CrossRef]

 8.  Dokuyucu, M. A. (2020). Caputo and atangana-baleanu-caputo fractional derivative applied to garden equation. Turkish Journal of Science, 5(1), 1–7. [Google Scholar]

 9.  Dokuyucu, M. A., Celik, E. (2021). Analyzing a novel coronavirus model (COVID-19) in the sense of caputo-fabrizio fractional operator. Applied and Computational Mathematics, 20(1), 49–69. [Google Scholar]

10. Ilknur, K., AKçetin, E., Yaprakdal, P. (2020). Numerical approximation for the spread of siqr model with caputo fractional order derivative. Turkish Journal of Science, 5(2), 124–139. [Google Scholar]

11. Dokuyucu, M. A. (2020). Analysis of a fractional plant-nectar-pollinator model with the exponential kernel. Eastern Anatolian Journal of Science, 6(1), 11–20. [Google Scholar]

12. Atangana, A., Araz, S. S. (2021). New concept in calculus: Piecewise differential and integral operators. Chaos, Solitons & Fractals, 145, 110638. DOI 10.1016/j.chaos.2020.110638. [Google Scholar] [CrossRef]

13. Qomariyah, N., Herdiana, R., Utomo, R., Permatasari, A. (2021). Stability analysis of a tuberculosis epidemic model with nonlinear incidence rate and treatment effects. Journal of Physics: Conference Series. 1943(1). 012118. [Google Scholar]

14https://en.wikipedia.org/wiki/second_derivative. [Google Scholar]

images This work is licensed under a Creative Commons Attribution 4.0 International License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.