Computer Modeling in Engineering & Sciences

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

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 [1–4]. 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 [5–11]. 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−α)ddt∫0t(t−x)−α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⁡[−α(t−x)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α[−α(t−x)α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 0≤t≤t1t1RLDtαh(t),if t1≤t≤T,(5)

where 0PRLDtα presents classical derivative on 0≤ t ≤ t1 and Riemann-Liouville fractional derivative on t1≤ t ≤ T.

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

0PCFDtαh(t)={h′(t),if 0≤t≤t1t1CFDtαh(t),if t1≤t≤T(6)

and

0PCFDtαh(t)={h′(t),if 0≤t≤t1t1CFRDtαh(t),if t1≤t≤T(7)

where 0PCFDtα presents classical derivative on 0≤ t ≤ t1 and Caputo-Fabrizio fractional derivative on t1≤ t ≤ T.

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

0PABDtαh(t)={h′(t),if 0≤t≤t1t1ABCDtαh(t),if t1≤t≤T(8)

where 0PABDtα presents classical derivative on 0≤ t ≤ t1 and Atangana-Baleanu fractional derivative on t1≤ t ≤ T.

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

PPLJtαh(t)={∫0t1h(τ)dτ,if 0≤t≤t11Γ(α)∫t1t(t−τ)α−1h(τ)dτ,if t1≤t≤T(9)

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

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

PCFJtαh(t)={∫0t1h(τ)dτ,if 0≤t≤t11−αM(α)h(t)+αM(α)∫t1th(τ)dτ,if t1≤t≤T(10)

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

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=ϕ(1−r1)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:

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 [2–14].

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=ϕ(1−r1)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{(1−S∗S)S˙+(1−E∗E)E˙+(1−I1∗I1)I˙1+(1−I2(t)∗I2(t))I˙2,=(S˙S)2S∗+(E˙E)2E∗+(I˙1I1)2I1∗+(I˙2I2)2I2∗+(1−S∗S)S¨+(1−E∗E)E¨+(1−I1∗I1)I¨1+(1−I2∗I2)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)=ϕ(1−r1)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∗+(1−S∗S){−β1(S.(t)I1(t)+I˙1(t)S(t))−β2(S.(t)I2(t)+I2.(t)S(t))−μS.(t)}+(1−E∗E){β1p1(S.(t)I1(t)+I1.(t)S(t))+β2q1(S.(t)I2(t)+I2.(t)S(t))−(μ+γ)E.(t)}+(1−I1∗I1){β1p(S.(t)I1(t)+I˙1(t)S(t))+β2q(S.(t)I2(t)+I2.(t)S(t))+γE.(t)−(ϕ+μ+δ1)I1.(t)}+(1−I2∗I2){ϕ(1−r1)I˙1(t)−(μ+δ2)I2.(t)−φr2I2.(t)}.(16)

dL.dt=L.(S,E,I1,I2)+S.(t){−β1I1(t)−β2I2(t)−μ−E∗Eβ1p1I1(t)−E∗Eβ2q1I2(t)−I1∗β1p−I1∗I1β2qI2(t)+β1p1I1(t)+β2q1I2(t)+β1pI1(t)+β2qI2(t)+S∗Sβ1I1(t)+S∗Sβ2I2(t)+S∗Sμ}+E.(t){−(μ+γ)−I1∗I1γ+γ}+I˙1(t){−β1S(t)−(ϕ+μ+δ1)−E∗Eβ1p1S(t)−I1∗I1β1pS(t)−I2∗I2ϕ(1−r1)+β1p1S(t)+β1pS(t)+ϕ(1−r1)+S∗β1+I1∗I1(ϕ+μ+δ1)}+I2.(t){−β2S(t)−(μ+δ2)−φr2−E∗Eβ2q1S(t)−I1∗I1β2qS(t)+β2q1S(t)+β2qS(t)+S∗β2+I2∗I2(μ+δ2)+I2∗I2φ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=Π1⁡L.(S,E,I1,I2)+Π+⏟−Π2⁡Π−⏟(18)

where

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

Π−=λ(β1I1(t)+β2I2(t)+μ+E∗Eβ1p1I1(t)+E∗Eβ2q1I2(t)+β1pI1∗(t)+I1∗I1β2qI2(t))+(β1S(t)I1(t)+β2S(t)I2(t)+μS(t)).(β1p1I1(t)+β2q1I2(t)+β1pI1(t)+β2qI2(t)S∗Sβ1I1(t)+S∗Sβ2I2(t)+S∗Sμ)+(β1p1S(t)I1(t)+β2q1S(t)I2(t))(−(μ+γ)−I1∗I1γ)−γ(μ+γ)E(t)+(ϕ+μ+δ1)I1(t)(β1p1S(t)+I1∗I1(ϕ+μ+δ1)+β1pS(t)+β1S∗+ϕ(1−r1))+(β1pS(t)I1(t)+β2qS(t)I2(t))+γE(t)(β1S(t)+(ϕ+μ+δ1)+E∗Eβ1p1S(t)+I1∗I1β1pS(t)+I2∗I2ϕ(1−r1))+ϕ(1−r1)I1(t)(β2S(t)+(μ+δ2)+φr2+E∗Eβ2q1S(t)+I1∗I1β2qS(t)+(φr2I2)(β2q1S(t)+qβ2S(t)+S∗β2+I2∗I2(μ+δ2)+I2∗I2φ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.=ϕ(1−r1)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),∂2∂I12(S(t))=0.(24)

In this case, we can have the following:

FA=[00].(25)

Then

det(FAV−1−λ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=ϕ(1−r1)I1(t)−(μ+δ2)I2(t)−φr2I2(t).if0≤t≤W1,(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)=ϕ(1−r1)I1(t)−(μ+δ2)I2(t)−φr2I2(t).if W1≤t≤W2,(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)=[ϕ(1−r1)I1(t)−(μ+δ2)I2(t)−φr2I2(t)]dt+σ4I2dB4(t).if W2≤t≤W.(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)|2≤ki(1+|xi|2) for i=1,2,3,4.(31)

and

2) The Lipschitz condition

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

Now we define the norm ‖φ‖∞=t∈Dφ⁡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 0≤t≤W2.(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)|2≤k1(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)