New approaches to the fractional dynamics of schistosomiasis disease model

Physica A 525 (2019) 373–393

Contents lists available at ScienceDirect

Physica A journal homepage: www.elsevier.com/locate/physa

New approaches to the fractional dynamics of schistosomiasis disease model ∗

Mehmet Yavuz a , , Ebenezer Bonyah b a b

Department of Mathematics-Computer, Faculty of Science, Necmettin Erbakan University, 42090, Konya, Turkey Department of Mathematics Education, University of Education Winneba (Kumasi Campus), Ghana

article

info

Article history: Received 25 November 2018 Received in revised form 10 February 2019 Available online 28 March 2019 Keywords: Schistosomiasis model Mittag-Leffler kernel Exponential law kernel Trematodes Power law distribution

a b s t r a c t In this paper, schistosomiasis fractional order dynamic model is examined via exponential law kernel sense and Mittag-Leffler kernel in Liouville–Caputo sense. Some special solutions for two operators are obtained using the iterative scheme through Laplace transform and Sumudu–Picard integration technique. The uniqueness and existence of solution for both operators are established. The numerical solutions for both operators approve that the desirable results can be obtained when alpha value is less than 1. © 2019 Published by Elsevier B.V.

1. Introduction Schistosomiasis is one of the neglected tropical diseases and second most prevalent among the tropical diseases. This disease is caused by flatworms which belong to a group of trematodes within the phylum, platyhelminthes and genus schistosoma. In some years, the schistosoma infection has been estimated to be 218 million people over the world. The most common geographical regions where the disease is prevalent are Africa, Eastern Asian and South America [1–4]. Mathematical modelling of Schistosomiasis is not new but mostly formulated in ordinary differential equations. Several researchers have developed new schistosomiasis models or modified existing ones using a deterministic approach [3,5–7]. These integer order mathematical models cannot capture the most of the hidden properties of many natural phenomena. In the recent years, the fractional calculus has been intensively used to describe the hereditary properties, dissipative effects, memory, and damage structures in the real world problems. Riemann–Liouville (RL) and Caputo fractional operators are outstanding definitions of fractional calculus defined by the convolution of a given function/its derivative and the power decay function as a kernel [8]. According to this definition, these operators are non-local. On the other hand, the singularity arising from power decay kernel function reveals many significant computational difficulties and hence requires introducing numerical solutions. For several decades mathematicians have applied the concept of fractional differential and integral operators in many areas of interest in science and technology related fields with power law distribution. This power law distribution unfortunately has no interpretations in statistics when the fractional order is less than 1. The class of fractional differential operators hinged on the so-called power law kernel happen to meet some classical criterions such as index law and classical mechanics law and singular kernels [9–14]. A thorough analysis of mechanic indicates that a certain classic of forces such as weak and strong forces interestingly have nothing to do with commutativity or index law rather weak forces have. This implies that those operators having ∗ Corresponding author. E-mail addresses: [email protected] (M. Yavuz), [email protected] (E. Bonyah). https://doi.org/10.1016/j.physa.2019.03.069 0378-4371/© 2019 Published by Elsevier B.V.

374

M. Yavuz and E. Bonyah / Physica A 525 (2019) 373–393

the form of power law distribution are physically weak and do not have strong capabilities to capture more complex phenomenon. Another problem is the singularity which means object cannot be defined in mathematical sense and hence poses a great challenge in describing many natural phenomena with singularity. In recent times, some mathematicians have made the point that singular differential operator is sufficient to present information from the past or to capture real problems with singularities which are highly misleading. Because there are different kinds of singularities in real life that cannot be stood for one operator hinged on power-law singularities [15–17]. Recent developments in the concept of fractional calculus have led to the establishment of two powerful operators which are Caputo–Fabrizio and Atangana–Baleanu. These new derivatives have non-singular kernels and do not operate under the power-distribution. The Atangana–Baleanu operator is based on generalized Mittag-Leffler function which possesses strong forces as their diffusive effects are crossovers that give good predictions. The Caputo–Fabrizio is also based on exponential law which highlighted in many natural processes and also has crossover effect with statistical representation. These operators have shown to be the future for modelling of several scientific processes including chaotic theory, biological processes, heat conduction problem, control theory and financial problem, etc. [18–36]. Especially, Singh et al. [37] proposed a new fractional model of chemical kinetic system and analysed using Atangana–Baleanu derivative. The numerical solution based on new operator was compared with classical derivative which proved to be efficient. In [38] Singh et al. constructed a fractional epidemiological model for computer viruses pertaining to a new fractional derivative and obtained uniqueness solution using fixed point theory. A numerical solution was also carried out to depict impact of the arbitrary order derivative. Kumar et al. [39] investigated a fractional model of the Ambartsumian equation and obtained solution using HATM approach. The solution obtained in the proposed study was represented in the form of series and also graphical drawn to depict the solution. In [40] Singh et al. developed an efficient numerical algorithm for the fractional Drinfeld–Sokolov–Wilson equation and obtained a numerical scheme that produced a result which is very accurate, flexible, effective and easy to use. In [41] Gómez-Aguilar and Atangana proposed a new operator comprised operator power-law-exponential-Mittag-Leffler kernel. Moreover they presented numerical scheme with some examples of physical problem and obtained desirable results. Furthermore, Atangana and Gómez-Aguilar [42] constructed a new differential operator and integral operator by using the normal distribution as kernel. Also they established the relations of new operators and existing integrals transform. Atangana–Baleanu and Caputo–Fabrizio fractional derivatives have the properties of mean squares distribution and their probability distributions are Gaussian to non-Gaussian crossover. The Caputo–Fabrizio derivative seems to be less noisy while Atangana–Baleanu derivative presents an excellent description, as a result of its Mittag-Leffler memory. The stated two operators do not have artificial singularity that is incorporated into a given model. Each of them has different waiting times distribution which is the ideal waiting time distribution as observed in several biological processes including schistosomiasis. The aim of the paper is to explore the role of exponential law distribution and generalized Mittag-Leffler distribution in investigating the transmission dynamics of schistosomiasis. 2. Preliminaries We start with recalling some well-known mathematical tools used in the subsequent sections. Definition 2.1 ([43]). The Caputo–Fabrizio fractional derivative in Liouville–Caputo sense (CFC) is defined as CFC α 0 Dt f

(t ) =

M (α) 1−α

t

∫

[ ] α (t − τ ) f (τ ) exp − dτ , 1−α ′

0

0 < α ≤ 1,

(1)

where M (α) is a normalization function such that M (0) = M (1) = 1. The CF definition can also be applied for the functions which do not belong to H 1 (a, b) and the kernel does not have singularity for t = τ unlike the Caputo fractional derivative. In the present study, we apply the Laplace transform of CF derivative with respect to the time variable t. CFC α 0 Dt f

(t ) is given by: [ ]} {CFC α+n } { (α+n) } 1 αt L 0 Dt f (t ) (p) = L f (t ) L exp − 1−α 1−α pn+1 L {f (t )} − pn f (0) − pn−1 f ′ (0) − · · · − f (n) (0) . = p + α (1 − p)

Definition 2.2 ([43,44]). The Laplace transform of

{

(2)

Definition 2.3 ([45]). The Caputo–Fabrizio fractional integral of order α ( 0 < α ≤ 1) of the function f (t ) is defined as CF α 0 It f

(t ) =

2 (1 − α)

(2 − α) M (α)

f (t ) +

2α

(2 − α) M (α)

∫

t

f (s) ds,

t ≥ 0,

(3)

0

where M (α) =

2 2−α

,

0 < α < 1.

(4)

M. Yavuz and E. Bonyah / Physica A 525 (2019) 373–393

375

Definition 2.4 ([12]). The Atangana–Baleanu fractional derivative in Caputo sense (ABC) is defined as: ABC α 0 Dt f

N (α)

(t ) =

t

∫

n−α

[

f (n) (τ ) Eα −

0

α (t − τ )α n−α

]

dτ ,

n − 1 < α ≤ n,

(5)

where N (0) = N (1) = 1 is a normalization function and Eα is the Mittag-Leffler function. Definition 2.5 ([12]). The Laplace transform of L

{ABC 0

ABC α 0 Dt f

(t ) is defined as

α (t − τ )α n−α 1−α 0 N (α) pα L {f (t )} (p) − pα−1 f (0) . = α 1−α pα + 1− p

Dαt f (t ) (p) =

}

N (α)

{∫

t

L

[

f ′ (τ ) Eα −

]

dτ

}

(p) (6)

Definition 2.6 ([12]). The Atangana–Baleanu fractional integral of order α of a function f (t ) is given as AB α 0 It f

(t ) =

1−α N (α)

f (t ) +

α Γ (α) N (α)

t

∫

f (s) (t − s)α−1 ds.

(7)

0

Theorem 2.7 ([12]). The following time fractional ordinary differential equation ABC α 0 Dt f

(t ) = ψ (t )

(8)

has a unique solution considering the inverse LT and the convolution property below: f (t ) =

1−α N (α)

ψ (t ) +

α Γ (α) N (α)

t

∫

ψ (s) (t − s)α−1 ds.

(9)

0

Definition 2.8 ([46]). The Sumudu transform of the AB fractional operator is given by ST

{ABC 0

Dαt f (t ) (s) =

}

N (α) 1−α

(

[ ]) uα α Γ (α + 1) Eα − (ST {f (t )} − f (0)) . 1−α

(10)

3. Mathematical model formulation of the schistosomiasis disease The model sub-divides the entire human population at time t, denoted by N (t ), into the following sub-populations of susceptible individuals who have not infected with schistosomiasis yet S (t ), those that have been infected with schistosomiasis and infectious Ih (t ), those have recovered from the infection R (t ). Thus the total human population is N = S + Ih + R. Furthermore, the total snail vector population, given by Nv , is also partitioned into susceptible snails who have not been infected with schistosomiasis Sv and those who have been infected with schistosomiasis and infectious Iv. This implies that the total vector population is given by Nv = Sv + Iv . The recruitment rate for susceptible human and snail population are denoted by Λ and Λv respectively. Susceptible humans acquire schistosomiasis infection following effective contact with infectious snails at a rate βh . Susceptible snails acquire schistosomiasis infection through effective contact with infected humans at a rate βv . The infected humans progress to the recovery class at a rate δ . The recovered individuals loss immunity at a rate γ . The mortality of humans and snails are denoted by µh and µv respectively. Based on the interrelationship with the compartments the following nonlinear ordinary differential equations are obtained. dSh dt dIh dt dR dt dSv dt dIv dt

= Λ − βh Iv Sh + γ R − µh Sh , = βh Iv Sh − (δ + µh ) Ih , = δ Ih − (γ + µh ) R, = Λv − βv Ih Sv − µv Sv , = βv Ih Sv − µv Iv .

(11)

376

M. Yavuz and E. Bonyah / Physica A 525 (2019) 373–393

3.1. The schistosomiasis model in the exponential law kernel sense The fractional schistosomiasis model with the exponential law kernel is given as CFC α 0 Dt Sh

= Λ − βh Iv Sh + γ R − µh Sh ,

CFC α 0 D t Ih

= βh Iv Sh − (δ + µh ) Ih ,

CFC α 0 Dt R

= δ Ih − (γ + µh ) R,

CFC α 0 Dt Sv

= Λv − βv Ih Sv − µv Sv ,

CFC α 0 D t Iv

= βv Ih Sv − µv Iv .

(12)

CFC α 0 Dt

where represents the fractional operator of type Caputo–Fabrizio–Caputo and α with 0 < α ≤ 1 is the fractional order. This model is subject to initial conditions Ih(0) (t ) = Ih (0) , Iv(0) (t ) = Iv (0) .

Sh(0) (t ) = Sh (0) , Sv(0) (t ) = Sv (0) ,

R0 (t ) = R (0) ,

(13)

By using the fixed-point theorem, we define the existence of the solution of the system (12). Firstly, if we transform the system (12) into an integral equation as follows: Sh (t ) − Sh (0) = Ih (t ) − Ih (0) = R (t ) − R (0) =

CF α 0 It

CF α 0 It

CF α 0 It

Sv (t ) − Sv (0) = Iv (t ) − Iv (0) =

[Λ − βh Iv Sh + γ R − µh Sh ] ,

[βh Iv Sh − (δ + µh ) Ih ] ,

[δ Ih − (γ + µh ) R] ,

CF α 0 It

CF α 0 It

(14)

[Λv − βv Ih Sv − µv Sv ] ,

[βv Ih Sv − µv Iv ] .

In view of the fractional integral of order α given by Eq. (3), we have 2 (1 − α) [Λ − βh Iv (t ) Sh (t ) + γ R (t ) − µh Sh (t )] (2 − α) M (α) ∫ t 2α [Λ − βh Iv (k) Sh (k) + γ R (k) − µh Sh (k)] dk, + (2 − α) M (α) 0 2 (1 − α) [βh Iv (t ) Sh (t ) − (δ + µh ) Ih (t )] Ih (t ) = Ih (0) + (2 − α) M (α) ∫ t 2α [βh Iv (k) Sh (k) − (δ + µh ) Ih (k)] dk, + (2 − α) M (α) 0 2 (1 − α) [δ Ih (t ) − (γ + µh ) R (t )] R (t ) = R (0) + (2 − α) M (α) ∫ t 2α [δ Ih (k) − (γ + µh ) R (k)] dk, + (2 − α) M (α) 0 2 (1 − α) [Λv − βv Ih (t ) Sv (t ) − µv Sv (t )] Sv (t ) = Sv (0) + (2 − α) M (α) ∫ t 2α [Λv − βv Ih (k) Sv (k) − µv Sv (k)] dk, + (2 − α) M (α) 0 2 (1 − α) [βv Ih (t ) Sv (t ) − µv Iv (t )] Iv (t ) = Iv (0) + (2 − α) M (α) ∫ t 2α [βv Ih (k) Sv (k) − µv Iv (k)] dk. + (2 − α) M (α) 0 Sh (t ) = Sh (0) +

(15)

We now consider the following kernels

τ (t , Sh (t )) = Λ − βh Iv Sh + γ R − µh Sh , ϕ (t , Ih (t )) = βh Iv Sh − (δ + µh ) Ih , ξ (t , R (t )) = δ Ih − (γ + µh ) R, η (t , Sv (t )) = Λv − βv Ih Sv − µv Sv , Φ (t , Iv (t )) = βv Ih Sv − µv Iv .

(16)

M. Yavuz and E. Bonyah / Physica A 525 (2019) 373–393

Theorem 3.1.

377

The kernels τ , ϕ, ξ , η and Φ given in Eq. (16) satisfy the Lipschitz condition.

Proof. Let Sh and sh , for the kernel τ , Ih and ih , for the kernel ϕ , R and r, for the kernel ξ , Sv and sv , for the kernel η and Iv and iv , for the kernel Φ be two functions; then we obtain the following:

∥τ (t , Sh (t )) − τ (t , sh (t ))∥ = ∥Λ − βh Iv (Sh (t ) − sh (t )) + γ R − µh (Sh (t ) − sh (t ))∥ , ∥ϕ (t , Ih (t )) − ϕ (t , ih (t ))∥ = ∥βh Iv Sh − (δ + µh ) (Ih (t ) − ih (t ))∥ , ∥ξ (t , R (t )) − ξ (t , r (t ))∥ = ∥δ Ih − (γ + µh ) (R (t ) − r (t ))∥ , ∥η (t , Sv (t )) − η (t , sv (t ))∥ = ∥Λv − βv Ih (Sv (t ) − sv (t )) − µv (Sv (t ) − sv (t ))∥ , ∥Φ (t , Iv (t )) − Φ (t , iv (t ))∥ = ∥βv Ih Sv − µv (Iv (t ) − iv (t ))∥ .

(17)

Using Cauchy’s inequality in Eq. (17), we have

∥τ (t , Sh (t )) − τ (t , sh (t ))∥ ≤ ∥Λ − βh Iv (Sh (t ) − sh (t )) + γ R − µh (Sh (t ) − sh (t ))∥ , ∥ϕ (t , Ih (t )) − ϕ (t , ih (t ))∥ ≤ ∥βh Iv Sh − (δ + µh ) (Ih (t ) − ih (t ))∥ , ∥ξ (t , R (t )) − ξ (t , r (t ))∥ ≤ ∥δ Ih − (γ + µh ) (R (t ) − r (t ))∥ , ∥η (t , Sv (t )) − η (t , sv (t ))∥ ≤ ∥Λv − βv Ih (Sv (t ) − sv (t )) − µv (Sv (t ) − sv (t ))∥ , ∥Φ (t , Iv (t )) − Φ (t , iv (t ))∥ ≤ ∥βv Ih Sv − µv (Iv (t ) − iv (t ))∥ .

(18)

Considering the following recursive formula, we can write Sh (t ) =

2 (1 − α)

( ) τ t , Sh(n−1) +

2α

∫

t

( ) τ k, Sh(n−1) dk,

(2 − α) M (α) (2 − α) M (α) 0 ∫ t ( ( ) ) 2α 2 (1 − α) Ih (t ) = ϕ k, Ih(n−1) dk, ϕ t , Ih(n−1) + (2 − α) M (α) (2 − α) M (α) 0 ∫ t ( ) ( ) 2 (1 − α) 2α ξ k, R(n−1) dk, R (t ) = ξ t , R(n−1) + (2 − α) M (α) (2 − α) M (α) 0 ∫ t ( ) ( ) 2 (1 − α) 2α Sv (t ) = η k, Sv(n−1) dk, η t , Sv(n−1) + (2 − α) M (α) (2 − α) M (α) 0 ∫ t ) ) ( ( 2 (1 − α) 2α Iv (t ) = Φ k, Iv(n−1) dk. Φ t , Iv(n−1) + (2 − α) M (α) (2 − α) M (α) 0

(19)

We now present the difference between the successive terms, applying the norm and the triangular inequality, we obtain

  ) ) ( 2 (1 − α)  ( τ t , Sh(n−1) (t ) − τ t , sh(n−2) (t )  ∥Gn (t )∥ = Sh(n) (t ) − sh(n−1) (t ) ≤  ∫ t (2 − α) M (α)  [ ( ) ( )]  2α  τ k, Sh(n−1) (k) − τ k, sh(n−2) (k) dk + , (2 − α) M (α)  0    ) ) ( ( 2 (1 − α) ϕ t , Ih(n−1) (t ) − ϕ t , ih(n−2) (t )  ∥Kn (t )∥ = Ih(n) (t ) − ih(n−1) (t ) ≤ ∫ t (2 − α) M (α)   [ ( ) ( )]  2α  , + ϕ k , I k − ϕ k , i k dk ( ) ( ) h n − 1 h n − 2 ( ) ( )  (2 − α) M (α)  0   ) ( ) 2 (1 − α)  ( ξ t , R(n−1) (t ) − ξ t , r(n−2) (t )  ∥Pn (t )∥ = R(n) (t ) − r(n−1) (t ) ≤ 2 − α) M ( (α)  ∫ t  [ ( ) ( )]  2α ,  + ξ k , R k − ξ k , r k dk ( ) ( ) (n−1) (n−2)  (2 − α) M (α)  0   ) ( ) 2 (1 − α)  ( η t , Sv (n−1) (t ) − η t , sv (n−2) (t )  ∥Qn (t )∥ = Sv (n) (t ) − sv (n−1) (t ) ≤ 2 − α) M ( (α)  ∫ t  [ ( ) ( )]  2α ,  + η k , S t − η k , s k dk ( ) ( ) v (n−1) v (n−2)  (2 − α) M (α)  0   ) ( ) 2 (1 − α)  ( Φ t , Iv (n−1) (t ) − Φ t , iv (n−2) (t )  ∥Zn (t )∥ = Iv (n) (t ) − iv (n−1) (t ) ≤  ∫ t (2 − α) M (α)  [ ( ) ( )]  2α  + Φ k, Iv (n−1) (k) − Φ k, iv (n−2) (k) dk , (2 − α) M (α)  0

(20)

where Sh(n) (t ) =

∞ ∑

Gm (t ) ,

Ih(n) (t ) =

m=0

Sv(n) (t ) =

∞ ∑ m=0

∞ ∑

Km (t ) ,

m=0

Qm (t ) ,

Iv(n) (t ) =

∞ ∑ m=0

Rn (t ) =

∞ ∑ m=0

Zm (t ) .

Pm (t ) , (21)

378

M. Yavuz and E. Bonyah / Physica A 525 (2019) 373–393

Since the kernels τ , ϕ, ξ , η and Φ satisfy the Lipschitz condition, we can write

    2 (1 − α) ∥Gn (t )∥ = Sh(n) (t ) − sh(n−1) (t ) ≤ ∆1 Sh(n−1) (t ) − sh(n−2) (t ) ∫ t (2 − α) M (α)   2α Sh(n−1) (k) − sh(n−2) (k) dk, ∆2 + (2 − α) M (α) 0     2 (1 − α) ∥Kn (t )∥ = Ih(n) (t ) − ih(n−1) (t ) ≤ ∆3 Ih(n−1) (t ) − ih(n−2) (t ) ∫ t (2 − α) M (α)   2α Ih(n−1) (k) − ih(n−2) (k) dk, ∆4 + (2 − α) M (α) 0     2 (1 − α) ∥Pn (t )∥ = R(n) (t ) − r(n−1) (t ) ≤ ∆5 R(n−1) (t ) − r(n−2) (t ) ∫ t(2 − α) M (α)   2α R(n−1) (k) − r(n−2) (k) dk, ∆6 + (2 − α) M (α) 0     2 (1 − α) ∥Qn (t )∥ = Sv (n) (t ) − sv (n−1) (t ) ≤ ∆7 Sv (n−1) (t ) − sv (n−2) (t ) ∫ t (2 − α) M (α)   2α Sv (n−1) (t ) − sv (n−2) (k) dk, ∆8 + (2 − α) M (α) 0     2 (1 − α) ∥Zn (t )∥ = Iv (n) (t ) − iv (n−1) (t ) ≤ ∆9 Iv (n−1) (t ) − iv (n−2) (t ) ∫ t (2 − α) M (α)   2α Iv (n−1) (k) − iv (n−2) (k) dk. + ∆10 (2 − α) M (α) 0 These results complete the proof.

(22)

□

Theorem 3.2 (Existence of the Solution). The system given by Eq. (12) has a solution. Proof. Taking into consideration the results of Eq. (22) and using the recursive technique, we get the following relation:

∥Gn (t )∥ ≤ ∥Sh (0)∥ + ∥Kn (t )∥ ≤ ∥Ih (0)∥ +

2 (1 − α)

{{

(2 − α) M (α) 2 (1 − α)

{{

(2 − α) M (α)

{{

2 (1 − α)

∆1

∆3

}n

{ +

}n

(2 − α) M (α)

{ +

}n

2α 2α

(2 − α) M (α)

{

2α

∆2 t

∆4 t

}n }

}n }

}n }

,

,

, (2 − α) M (α) (2 − α) M (α) }n { }n } {{ 2α 2 (1 − α) ∥Qn (t )∥ ≤ ∥Sv (0)∥ + , ∆7 + ∆8 t (2 − α) M (α) (2 − α) M (α) {{ }n { }n } 2 (1 − α) 2α ∥Zn (t )∥ ≤ ∥Iv (0)∥ + ∆9 + ∆10 t . (2 − α) M (α) (2 − α) M (α) ∥Pn (t )∥ ≤ ∥R (0)∥ +

∆5

+

∆6 t

(23)

Therefore, Eq. (23) is exist and smooth. Nevertheless, to prove that the functions in Eq. (23) are a system of solutions of Eq. (12), we first assume

Sh (t ) = Sh(n) (t ) − Θ1(n) (t ) , Ih (t ) = Ih(n) (t ) − Θ2(n) (t ) , R (t ) = R(n) (t ) − Θ3(n) (t ) , Sv (t ) = Sv (n) (t ) − Θ4(n) (t ) , Iv (t ) = Iv (n) (t ) − Θ5(n) (t ) ,

(24)

M. Yavuz and E. Bonyah / Physica A 525 (2019) 373–393

379

where Θ1(n) (t ) , Θ2(n) (t ) , Θ3(n) (t ) , Θ4(n) (t ) and Θ5(n) (t ) are reminder terms of series solution. Thus, we have

( ) 2 (1 − α) τ t , Sh − Θ1(n) (t ) (2 − α) M (α) ∫ t ( ) 2α + τ k, Sh − Θ1(n) (k) dk, (2 − α) M (α) 0

Sh (t ) − Sh (n−1) (t ) =

Ih (t ) − Ih(n−1) (t ) =

2 (1 − α)

( ) ϕ t , Ih − Θ2(n) (t )

(2 − α) M (α) ∫ t ( ) 2α + ϕ k, Ih − Θ2(n) (k) dk, (2 − α) M (α) 0

R (t ) − R(n−1) (t ) =

2 (1 − α)

( ) ξ t , R − Θ3(n) (t )

(2 − α) M (α) ∫ t ( ) 2α + ξ k, R − Θ3(n) (k) dk, (2 − α) M (α) 0

(25)

( ) 2 (1 − α) η t , Sv − Θ4(n) (t ) (2 − α) M (α) ∫ t ) ( 2α + η k, Sv − Θ4(n) (k) dk, (2 − α) M (α) 0

Sv (t ) − Sv(n−1) (t ) =

Iv (t ) − Iv(n−1) (t ) =

2 (1 − α)

) ( Φ t , Iv − Θ5(n) (t )

(2 − α) M (α) ∫ t ) ( 2α + Φ k, Iv − Θ5(n) (k) dk. (2 − α) M (α) 0

Applying the norm on both sides and using the Lipschitz condition, we have

  ∫ t   2α Sh (t ) − 2 (1 − α) τ (t , Sh ) − Sh (0) − τ (k, Sh (k)) dk   (2 − α) M (α) (2 − α) M (α) 0 [ { }]   2 (1 − α) 2α ∆1 + ∆2 t , ≤ Θ1(n) (t ) 1 + (2 − α) M (α) (2 − α) M (α)   ∫ t   2α Ih (t ) − 2 (1 − α) ϕ (t , Ih ) − Ih (0) − ϕ (k, Ih (k)) dk   (2 − α) M (α) (2 − α) M (α) 0 [ }] {   2 (1 − α) 2α ≤ Θ2(n) (t ) 1 + ∆3 + ∆4 t , (2 − α) M (α) (2 − α) M (α)   ∫ t   2α R (t ) − 2 (1 − α) ξ (t , R) − R (0) −  ξ k , R k dk ( ( ))   (2 − α) M (α) (2 − α) M (α) 0 [ { }]   2 (1 − α) 2α ≤ Θ3(n) (t ) 1 + ∆5 + ∆6 t , (2 − α) M (α) (2 − α) M (α)   ∫ t   2α Sv (t ) − 2 (1 − α) η (t , Sv ) − Sv (0) −  η k , S k dk ( ( )) v   (2 − α) M (α) (2 − α) M (α) 0 [ { }]   2 (1 − α) 2α ∆7 + ∆8 t , ≤ Θ4(n) (t ) 1 + (2 − α) M (α) (2 − α) M (α)   ∫ t   2α Iv (t ) − 2 (1 − α) Φ (t , Iv ) − Iv (0) − Φ k , I k dk ( v ( ))    (2 − α) M (α) (2 − α) M (α) 0 [ { }]   2 (1 − α) 2α ≤ Θ5(n) (t ) 1 + ∆9 + ∆10 t . (2 − α) M (α) (2 − α) M (α)

(26)

380

M. Yavuz and E. Bonyah / Physica A 525 (2019) 373–393

Taking the limit n → ∞ of Eq. (26), we obtain

2 (1 − α)

Sh (t ) = Sh (0) +

2 (1 − α)

Ih (t ) = Ih (0) +

R (t ) = R (0) +

(2 − α) M (α)

(2 − α) M (α) 2 (1 − α)

(2 − α) M (α)

Sv (t ) = Sv (0) +

Iv (t ) = Iv (0) +

τ (t , Sh (t )) +

ϕ (t , Ih (t )) +

ξ (t , R (t )) +

2 (1 − α)

(2 − α) M (α) 2 (1 − α)

(2 − α) M (α)

2α

(2 − α) M (α) 2α

η (t , Sv (t )) +

Φ (t , Iv (t )) +

t

∫

ϕ (k, Ih (k)) dk, 0

ξ (k, R (k)) dk, 0

2α

(2 − α) M (α)

t

∫

(2 − α) M (α) 2α

(27)

t

∫

(2 − α) M (α)

τ (k, Sh (k)) dk, 0

(2 − α) M (α) 2α

t

∫

η (k, Sv (k)) dk, 0

t

∫

Φ (k, Iv (k)) dk. 0

This says that Eq. (27) is the solution of Eq. (12) which completes the proof.

□

Theorem 3.3 (Uniqueness of the Solution). The solution of system given by Eq. (12) is unique.

Proof. We can get the other solutions for the system (12), say Sh∗ (t ), Ih∗ (t ), R∗ (t ), Sv∗ (t ), and Iv∗ (t ); then

)] ( 2 (1 − α) [ τ (t , Sh (t )) − τ t , Sh∗ (t ) (2 − α) M (α) ∫ t [ ( )] 2α + τ (k, Sh (k)) − τ k, Sh∗ (k) dk, (2 − α) M (α) 0 ( )] 2 (1 − α) [ Ih (t ) − Ih∗ (t ) = ϕ (t , Ih (t )) − ϕ t , Ih∗ (t ) (2 − α) M (α) ∫ t [ ( )] 2α + ϕ (k, Ih (k)) − ϕ k, Ih∗ (k) dk, (2 − α) M (α) 0 ( )] 2 (1 − α) [ R (t ) − R∗ (t ) = ξ (t , R (t )) − ξ t , R∗ (t ) (2 − α) M (α) ∫ t [ ( )] 2α + ξ (k, R (t )) − ξ k, R∗ (k) dk, (2 − α) M (α) 0 ( )] 2 (1 − α) [ Sv (t ) − Sv∗ (t ) = η (t , Sv (t )) − η t , Sv∗ (t ) (2 − α) M (α) ∫ t [ ( )] 2α + η (k, Sv (k)) − η k, Sv∗ (k) dk, (2 − α) M (α) 0 ( )] 2 (1 − α) [ Φ (t , Iv (t )) − Φ t , Iv∗ (t ) Iv (t ) − Iv∗ (t ) = (2 − α) M (α) ∫ t [ ( )] 2α + Φ (k, Iv (k)) − Φ k, Iv∗ (k) dk. (2 − α) M (α) 0

Sh (t ) − Sh∗ (t ) =

(28)

M. Yavuz and E. Bonyah / Physica A 525 (2019) 373–393

381

Applying the norm to both sides of Eq. (28), we obtain 2 (1 − α)

 ) ( τ (t , Sh (t )) − τ t , S ∗ (t )  h (2 − α) M (α) ∫ t )] [ ( 2α τ (k, Sh (k)) − τ k, S ∗ (k)  dk, + h (2 − α) M (α) 0   ) ( 2 (1 − α)  ∗ Ih (t ) − I (t ) ≤ ϕ (t , Ih (t )) − ϕ t , I ∗ (t )  h h (2 − α) M (α) ∫ t )] ( [ 2α ϕ (k, Ih (k)) − ϕ k, I ∗ (k)  dk, + h (2 − α) M (α) 0   ) ( 2 (1 − α)  ∗ R (t ) − R (t ) ≤ ξ (t , R (t )) − ξ t , R∗ (t )  (2 − α) M (α) ∫ t )] [ ( 2α ξ (k, R (t )) − ξ k, R∗ (k)  dk, + (2 − α) M (α) 0   ) ( 2 (1 − α)  ∗ Sv (t ) − S (t ) ≤ η (t , Sv (t )) − η t , S ∗ (t )  v v (2 − α) M (α) ∫ t )] [ ( 2α η (k, Sv (k)) − η k, S ∗ (k)  dk, + v (2 − α) M (α) 0    ) ( Iv (t ) − I ∗ (t ) ≤ 2 (1 − α) Φ (t , Iv (t )) − Φ t , I ∗ (t )  v v (2 − α) M (α) ∫ t )] [ ( 2α Φ (k, Iv (k)) − Φ k, I ∗ (k)  dk. + v (2 − α) M (α) 0   Sh (t ) − S ∗ (t ) ≤ h

(29)

Considering the results of Theorems 3.1 and 3.2, we have 2 (1 − α)

{

2α

}n

∆1 Υ1 + ∆2 Υ2 t , (2 − α) M (α) (2 − α) M (α) { }n   2α Ih (t ) − I ∗ (t ) ≤ 2 (1 − α) ∆3 Υ3 + ∆ Υ t , 4 4 h (2 − α) M (α) (2 − α) M (α) { }n   2α R (t ) − R∗ (t ) ≤ 2 (1 − α) ∆5 Υ5 + ∆6 Υ6 t , (2 − α) M (α) (2 − α) M (α) { }n   2 1 − α) 2α ( Sv (t ) − S ∗ (t ) ≤ ∆7 Υ7 + ∆8 Υ8 t , v (2 − α) M (α) (2 − α) M (α) { }n   2α Iv (t ) − I ∗ (t ) ≤ 2 (1 − α) ∆9 Υ9 + ∆ Υ t . 10 10 v (2 − α) M (α) (2 − α) M (α)

  Sh (t ) − S ∗ (t ) ≤ h

(30)

Eq. (30) holds for any n, hence we can obtain Sh (t ) = Sh∗ (t ) , Ih (t ) = Ih∗ (t ) , R (t ) = R∗ (t ) , Sv (t ) = Sv∗ (t ) , Iv (t ) = Iv∗ (t ) .

(31)

Then this result completes the proof. □ Now, we derive the approximate solution of the system given by Eq. (12) using the Laplace transform operator given by Eq. (2). Applying it to both sides of Eq. (12), we get pL {Sh (t )} − Sh (0) p + α (1 − p) pL {Ih (t )} − Ih (0)

= L {βh Iv Sh − (δ + µh ) Ih } (p) ,

p + α (1 − p) pL {R (t )} − R (0) p + α (1 − p)

= L {δ Ih − (γ + µh ) R} (p) ,

pL {Sv (t )} − Sv (0) p + α (1 − p) pL {Iv (t )} − Iv (0) p + α (1 − p)

= L {Λ − βh Iv Sh + γ R − µh Sh } (p) ,

= L {Λv − βv Ih Sv − µv Sv } (p) , = L {βv Ih Sv − µv Iv } (p) .

(32)

382

M. Yavuz and E. Bonyah / Physica A 525 (2019) 373–393

By applying the inverse Laplace transform of Eq. (32), we get

Sh (t ) = Sh (0) + L−1 Ih (t ) = Ih (0) + L R (t ) = R (0) + L

−1

p

{

p + α (1 − p) p

−1

Sv (t ) = Sv (0) + L

p + α (1 − p)

{

{

p + α (1 − p) p

−1

Iv (t ) = Iv (0) + L−1

{

} L {βh Iv Sh − (δ + µh ) Ih } (p) (t ) ,

} L {δ Ih − (γ + µh ) R} (p) (t ) ,

p + α (1 − p) p

{

p + α (1 − p) p

}

L {Λ − βh Iv Sh + γ R − µh Sh } (p) (t ) ,

(33)

} L {Λv − βv Ih Sv − µv Sv } (p) (t ) , }

L {βv Ih Sv − µv Iv } (p) (t ) .

The following recursive formula is then obtained

Sh(n) (t ) = L

−1

Ih(n) (t ) = L R(n) (t ) = L

p + α (1 − p) p

{

Sv(n) (t ) = L−1 Iv(n) (t ) = L−1

p + α (1 − p) p

{

−1

−1

{

p + α (1 − p) p

{

} } { L βh Iv(n−1) Sh(n−1) − (δ + µh ) Ih(n−1) (p) (t ) , L δ Ih(n−1) − (γ + µh ) R(n−1)

{

p + α (1 − p) p

{

p + α (1 − p) p

} { } L Λ − βh Iv(n−1) Sh(n−1) + γ R(n−1) − µh Sh(n−1) (p) (t ) ,

}

} (p) (t ) ,

(34)

}

L Λv − βv Ih(n−1) Sv(n−1) − µv Sv(n−1) (p) (t ) ,

}

{

}

L βv Ih(n−1) Sv(n−1) − µv Iv(n−1) (p) (t ) ,

{

}

where Sh(0) (t ) = Sh (0) , Ih(0) (t ) = Ih (0) , R(0) (t ) = R (0) , Sv(0) (t ) = Sv (0) , Iv(0) (t ) = Iv (0). The approximate solution is assumed to be obtained as a limit when n tends to infinity

Sh (t ) = lim Sh(n) (t ) , Ih (t ) = lim Ih(n) (t ) , R (t ) = lim R(n) (t ) , n→∞

n→∞

n→∞

(35)

Sv (t ) = lim Sv(n) (t ) , Iv (t ) = lim Iv(n) (t ). n→∞

n→∞

3.2. The numerical scheme and simulation with the Caputo–Fabrizio operator

This section shows the numerical results obtained the following numerical scheme as in [47] for varying fractional derivative order. The parameter values used for the simulation are Λ = 0.5, µh = 0.014, βh = 0.09, βv = 0.6, δ = 0.08, γ = 0.008, µh = 0.008 and the initial conditions employed are Sh (0) = 0.5, Ih (0) = 0.75, R(0) = 0.01, Sv (0) = 0.6, Iv (0) = 0.05. The simulation time used for this model is 800 s and the step size employed in evaluating the approximate solutions is h = 0.002. Fig. 1 shows the numerical simulation results of system model (11) with different fractional order values. It can be observed that as the fractional order value increases from 0.5 to 1 the rate of convergence reduces. The infective compartments of Fig. 1 show that as α = 0.5 increases to α = 1 the number of infected individuals increase. It can also be seen in Fig. 1 that as α = 0.5 increases to α = 1 the number of humans recovered increase. This shows that the fractional order and particular operator has a role in accurate predictions as Caputo–Fabrizio operator is less noisy. This implies that the crossover property of Caputo–Fabrizio is strong to capture complex model and give productive predictions.

M. Yavuz and E. Bonyah / Physica A 525 (2019) 373–393

383

3.2.1. R0 Sensitivity analysis In order to explore the sensitivity of R0 to each of its factors, we find the partial derivative for each parameter in the reproduction number.

∂ R0 βh βv Λv √ = > 0, βh βv Λh Λv ∂Λ 2 µ 2µh (δ + µh ) 2 v µ (δ+µ )µ h

v

h

∂ R0 βh βv Λ √ > 0, = βh βv Λh Λv ∂ Λv µ2 2µh (δ + µh ) µ (δ+µ )µ2 v h

v

h

∂ R0 Λβv Λv √ > 0, = Λβh βv Λv ∂βh 2 µ 2µh (δ + µh ) 2 v µ (δ+µ )µ h

v

h

∂ R0 Λβh Λv √ > 0, =− Λβh βv Λv ∂βv µ2 2µh (δ + µh ) µ (δ+µ )µ2 v

(36)

h h v βh βv Λh Λv

h βv Λh Λv − µβ(δ+µ 2 2 − µ2 (δ+µ )µ2 ∂ R0 h h v h ) µv h √ =− < 0, βh βv Λh Λv ∂µh 2

µh (δ+µh )µ2v

Λβh βv Λv ∂ R0 √ =− < 0, βh βv Λv ∂µv 3 µh (δ + µh ) µΛ(δ+µ µ 2 v )µ h

v

h

∂ R0 βh βv Λh Λv √ =− < 0. βh βv Λh Λv ∂δ 2µh (δ + µh )2 µ2 µ (δ+µ )µ2 v h

h

v

In the sensitivity analysis, it has been established that R0 increases with the following parameters Λ, Λv , βh , βv and decreases with µv , µh , δ . This means that in order to reduce the spread of the disease, those parameters with negative values of the partial derivative must be minimized in the environment. 3.2.2. Stability analysis The system given by Eq. (12) has two equilibrium points E0 = (Λ/µh , 0, 0, Λv /µv , 0) and E1 . E0 is the schistosomiasisfree equilibrium and the Jacobian matrix JE0 of Eq. (12) evaluated at the schistosomiasis-free equilibrium E0 is given by

⎛

−µh

⎜ 0 ⎜ JE0 = ⎜ 0 ⎝ 0 0

γ

0

− (δ + µh ) δ −Λv /µv ∗ βv Λv /µv ∗ βv

0

− (γ + µh )

−Λ/µh ∗ βh Λ/µh ∗ βh −µv 0

−µv

0

0 0

⎞

0 0 0 0

⎟ ⎟ ⎟. ⎠

(37)

The transmission matrix F and transition matrix V are obtained as

( F =

0

βv ∗ Λv /µv

βh ∗ Λ/µh 0

)

( and

V =

δ + µh 0

0

µv

)

.

By the next generation matrix method, the next generation matrix is defined by FV −1 with the basic reproduction number R0 of system (12) given by the spectral radius of FV −1 as

√ R0 =

Λv βh βv Λ . µh (δ + µh ) µ2v

(38)

Theorem 3.4. The disease — free equilibrium E0 is locally asymptotically stable if R0 < 1, and unstable if R0 > 1. Proof. The disease free-steady-state is locally asymptotically stable if all eigenvalues λi , i ∈ {1, 2, 3, 4, 5} of the Jacobian matrix JE0 meet the following condition [48]

⏐ απ ⏐⏐ ⏐ ⏐arg(λi ) > ⏐.

2 The eigenvalues of the Jacobian matrix JE0 are λ1 = −µh , λ2 = −µv , λ3 = −(γ + µh ); the remaining√two roots can be

obtained by using the quadratic equation λ2 + (δ + µh + µv )λ + (δµv + µv µh )(1 − R0 ) = 0 where R0 = the proof is complete. □

Λv βh βv Λ . µh (δ+µh )µ2v

Thus

384

M. Yavuz and E. Bonyah / Physica A 525 (2019) 373–393

Fig. 1. Numerical results for system (12) with respect to the Caputo–Fabrizio derivative.

3.3. The Schistosomiasis model in the Mittag-Leffler kernel sense

The fractional schistosomiasis model with the Mittag-Leffler kernel is given as ABC α 0 Dt Sh ABC α 0 Dt Ih ABC α 0 Dt R ABC α 0 Dt Sv ABC α 0 Dt Iv

= Λ − βh Iv Sh + γ R − µh Sh , = βh Iv Sh − (δ + µh ) Ih , = δ Ih − (γ + µh ) R,

(39)

= Λv − βv Ih Sv − µv Sv , = βv Ih Sv − µv Iv ,

α where ABC 0 Dt represents the fractional operator of type Atangana–Baleanu–Caputo (ABC) and α with 0 < α ≤ 1 is the fractional order. This model is subject to initial conditions

Sh(0) (t ) = Sh (0) , Sv(0) (t ) = Sv (0) ,

Ih(0) (t ) = Ih (0) , Iv(0) (t ) = Iv (0) .

R0 (t ) = R (0) ,

The system in Eq. (39) can be converted to the Volterra-type integral equation by using the ABC fractional integral.

(40)

M. Yavuz and E. Bonyah / Physica A 525 (2019) 373–393

385

By Theorem 2.7, the model can be written as:

(1 − α) {Λ − βh Iv (t ) Sh (t ) + γ R (t ) − µh Sh (t )} N (α) ∫ t α + (t − k)α−1 {Λ − βh Iv (k) Sh (k) + γ R (k) − µh Sh (k)} dk, Γ (α) N (α) 0 (1 − α) {βh Iv (t ) Sh (t ) − (δ + µh ) Ih (t )} Ih (t ) − Ih (0) = N (α) ∫ t α + (t − k)α−1 {βh Iv (k) Sh (k) − (δ + µh ) Ih (k)} dk, Γ (α) N (α) 0 (1 − α) {δ Ih (t ) − (γ + µh ) R (t )} R (t ) − R (0) = N (α) ∫ t α + (t − k)α−1 {δ Ih (k) − (γ + µh ) R (k)} dk, Γ (α) N (α) 0 (1 − α) {Λv − βv Ih (t ) Sv (t ) − µv Sv (t )} Sv (t ) − Sv (0) = N (α) ∫ t α + (t − k)α−1 {Λv − βv Ih (k) Sv (k) − µv Sv (k)} dk, Γ (α) N (α) 0 (1 − α) {βv Ih (t ) Sv (t ) − µv Iv (t )} Iv (t ) − Iv (0) = N (α) ∫ t α + (t − k)α−1 {βv Ih (k) Sv (k) − µv Iv (k)} dk, Γ (α) N (α) 0 Sh (t ) − Sh (0) =

(41)

For simplicity, we can assign the following kernels:

τ2 (t , Sh (t )) = Λ − βh Iv (t ) Sh (t ) + γ R (t ) − µh Sh (t ) , ϕ2 (t , Ih (t )) = βh Iv (t ) Sh (t ) − (δ + µh ) Ih (t ) , ξ2 (t , R (t )) = δ Ih (t ) − (γ + µh ) R (t ) , η2 (t , Sv (t )) = Λv − βv Ih (t ) Sv (t ) − µv Sv (t ) , Φ2 (t , Iv (t )) = βv Ih (t ) Sv (t ) − µv Iv (t ) .

(42)

Theorem 3.5. The kernels τ2 , ϕ2 , ξ2 , η2 and Φ2 given in Eq. (42) satisfy the Lipschitz condition and contraction if the following inequality holds: 0 ≤ ω1 , ω2 , ω3 , ω4 , ω5 < 1.

(43)

Proof. Consider the kernel τ2 (t , Sh (t )) = Λ − βh Iv (t ) Sh (t ) + γ R (t ) − µh Sh (t ). Let Sh and sh be two functions; then we obtain the following:

∥τ2 (t , Sh (t )) − τ2 (t , sh (t ))∥ = ∥− (βh Iv + µh ) (Sh (t ) − sh (t ))∥ ≤ βh ∥Iv ∥ ∥Sh (t ) − sh (t )∥

(44)

≤ βh k1 ∥Sh (t ) − sh (t )∥ = ω1 ∥Sh (t ) − sh (t )∥ , where ω1 = βh k1 , ω2 = δ + µh , ω3 = γ + µh , ω4 = βv k2 , ω5 = µv , k1 = maxt ∈I ∥Iv (t )∥ , k2 = maxt ∈I ∥Ih (t )∥ are bounded functions. Then we get

∥τ (t , Sh (t )) − τ (t , sh (t ))∥ ≤ ω1 ∥Sh (t ) − sh (t )∥ , ∥ϕ2 (t , Ih (t )) − ϕ2 (t , ih (t ))∥ ≤ ω2 ∥Ih (t ) − ih (t )∥ , ∥ξ2 (t , R (t )) − ξ2 (t , r (t ))∥ ≤ ω3 ∥R (t ) − r (t )∥ , ∥η2 (t , Sv (t )) − η2 (t , sv (t ))∥ ≤ ω4 ∥Sv (t ) − sv (t )∥ , ∥Φ2 (t , Iv (t )) − Φ2 (t , iv (t ))∥ ≤ ω5 ∥Iv (t ) − iv (t )∥ .

(45)

386

M. Yavuz and E. Bonyah / Physica A 525 (2019) 373–393

Considering the kernels of the model, Eq. (41) can be rewritten as:

∫ t α (1 − α) Sh (t ) = Sh (0) + τ2 (t , Sh (t )) + τ2 (k, Sh (k)) dk, N (α) Γ (α) N (α) 0 ∫ t α (1 − α) Ih (t ) = Ih (0) + ϕ2 (t , Ih (t )) + ϕ2 (k, Ih (k)) dk, N (α) Γ (α) N (α) 0 ∫ t α (1 − α) ξ2 (t , R (t )) + R (t ) = R (0) + ξ2 (k, R (k)) dk, N (α) Γ (α) N (α) 0 ∫ t α (1 − α) Sv (t ) = Sv (0) + η2 (t , Sv (t )) + η2 (k, Sv (k)) dk, N (α) Γ (α) N (α) 0 ∫ t α (1 − α) Iv (t ) = Iv (0) + Φ2 (t , Iv (t )) + Φ2 (k, Iv (k)) dk. N (α) Γ (α) N (α) 0

(46)

We now present the following recursive formula:

∫ t ) ) ( α (1 − α) ( τ2 k, Sh(n−1) (k) dk, τ2 t , Sh(n−1) (t ) + N (α) Γ (α) N (α) 0 ∫ t ) ( ) α (1 − α) ( Ih(n) (t ) = ϕ2 k, Ih(n−1) (k) dk, ϕ2 t , Ih(n−1) (t ) + N (α) Γ (α) N (α) 0 ∫ t ( ) ) α (1 − α) ( R(n) (t ) = ξ2 k, R(n−1) (k) dk, ξ2 t , R(n−1) (t ) + N (α) Γ (α) N (α) 0 ∫ t ) ) ( α (1 − α) ( Sv(n) (t ) = η2 k, Sv(n−1) (k) dk, η2 t , Sv(n−1) (t ) + N (α) Γ (α) N (α) 0 ∫ t ) ( ) α (1 − α) ( Iv(n) (t ) = Φ2 k, Iv(n−1) (k) dk. Φ2 t , Iv(n−1) (t ) + N (α) Γ (α) N (α) 0 Sh(n) (t ) =

(47)

with the initial conditions Sh(0) (t ) = Sh (0) , Sv(0) (t ) = Sv (0) ,

Ih(0) (t ) = Ih (0) , Iv(0) (t ) = Iv (0) .

R0 (t ) = R (0) ,

(48)

We next get the difference between the iterative terms in the expression

) ( )) (1 − α) ( ( τ2 t , Sh(n−1) (t ) − τ2 t , Sh(n−2) (t ) N (α) ∫ t ( ( ) ( )) α α−1 + τ2 k, Sh(n−1) (k) − τ2 k, Sh(n−2) (k) dk, (t − k) Γ (α) N (α) 0 ) ( )) (1 − α) ( ( Kn∗ (t ) = Ih(n) (t ) − Ih(n−1) (t ) = ϕ2 t , Ih(n−1) (t ) − ϕ2 t , Ih(n−2) (t ) N (α) ∫ t ( ( ) ( )) α α−1 + ϕ2 k, Ih(n−1) (k) − ϕ2 k, Ih(n−2) (k) dk, (t − k) Γ (α) N (α) 0 ) ( )) (1 − α) ( ( Pn∗ (t ) = R(n) (t ) − R(n−1) (t ) = ξ2 t , R(n−1) (t ) − ξ2 t , R(n−2) (t ) N (α) ∫ t ( ( ) ( )) α + (t − k)α−1 ξ2 k, R(n−1) (k) − ξ2 k, R(n−2) (k) dk, Γ (α) N (α) 0 ) ( )) (1 − α) ( ( ∗ Qn (t ) = Sv (n) (t ) − Sv (n−1) (t ) = η2 t , Sv (n−1) (t ) − η2 t , Sv (n−2) (t ) N (α) ∫ t ( ( ) ( )) α + (t − k)α−1 η2 k, Sv (n−1) (k) − η2 k, Sv (n−2) (k) dk, Γ (α) N (α) 0 ) ( )) (1 − α) ( ( ∗ Zn (t ) = Iv (n) (t ) − Iv (n−1) (t ) = Φ2 t , Iv(n−1) (t ) − Φ2 t , Iv(n−2) (t ) N (α) ∫ t ( ( ) ( )) α + (t − k)α−1 Φ2 k, Iv(n−1) (k) − Φ2 k, Iv(n−2) (k) dk, Γ (α) N (α) 0 G∗n (t ) = Sh(n) (t ) − Sh(n−1) (t ) =

(49)

M. Yavuz and E. Bonyah / Physica A 525 (2019) 373–393

387

where Sh(n) (t ) = Sv(n) (t ) =

∞ ∑

G∗m (t ) , Ih(n) (t ) =

∞ ∑

m=0

m=0

∞ ∑

∞ ∑

Qm∗ (t ) , Iv(n) (t ) =

m=0

Km∗ (t ) , Rn (t ) =

∞ ∑

∗ Pm (t ) ,

m=0

∗ Zm (t ) .

m=0

Applying the norm of both sides to Eq. (49) and considering triangular inequality, the equation becomes as the following equation:

 (  ∗    ) ( ) G (t ) = Sh(n) (t ) − Sh(n−1) (t ) ≤ (1 − α) τ2 t , Sh(n−1) (t ) − τ2 t , Sh(n−2) (t )  n N (α) ∫ t   ( ( ) ( ))  α α−1  τ2 k, Sh(n−1) (k) − τ2 k, Sh(n−2) (k) dk + (t − k) , Γ (α) N (α)  0

  ∗    ( ) ( ) K (t ) = Ih(n) (t ) − Ih(n−1) (t ) ≤ (1 − α) ϕ2 t , Ih(n−1) (t ) − ϕ2 t , Ih(n−2) (t )  n N (α)  ∫ t  ))  ( ) ( ( α α−1  + ϕ2 k, Ih(n−1) (k) − ϕ2 k, Ih(n−2) (k) dk (t − k) , Γ (α) N (α)  0

 ∗     ( ) ( ) P (t ) = R(n) (t ) − R(n−1) (t ) ≤ (1 − α) ξ2 t , R(n−1) (t ) − ξ2 t , R(n−2) (t )  n N (α)  ∫ t  ( ( ) ( ))  α α−1  + ξ2 k, R(n−1) (k) − ξ2 k, R(n−2) (k) dk (t − k) , Γ (α) N (α) 

(50)

0

  ∗    ( ) ) ( Q (t ) = Sv (n) (t ) − Sv (n−1) (t ) ≤ (1 − α) η2 t , Sv (n−1) (t ) − η2 t , Sv (n−2) (t )  n N (α)  ∫ t  ))  ( ) ( ( α α−1  + η2 k, Sv (n−1) (k) − η2 k, Sv (n−2) (k) dk (t − k) , Γ (α) N (α)  0

 (  ∗    ) ( ) Z (t ) = Iv (n) (t ) − Iv (n−1) (t ) ≤ (1 − α) Φ2 t , Iv(n−1) (t ) − Φ2 t , Iv(n−2) (t )  n N (α)  ∫ t  ))  ( ) ( ( α α−1  + Φ2 k, Iv(n−1) (k) − Φ2 k, Iv(n−2) (k) dk (t − k) . Γ (α) N (α)  0 Since the kernels satisfy the Lipschitz condition, we get the following

∫ t  ∗  (1 − α)  ∗    α G (t ) ≤   ω1 Gn−1 (t ) + ω1 (t − k)α−1 G∗n−1 (k) dk, n N (α) Γ (α) N (α) 0 ∫ t    ∗  (1 − α)  ∗  α K (t ) ≤ ω2 Kn−1 (t ) + ω2 (t − k)α−1 Kn∗−1 (k) dk, n N (α) Γ (α) N (α) 0 ∫ t    ∗  (1 − α)  ∗  α P (t ) ≤ ω3 Pn−1 (t ) + ω3 (t − k)α−1 Pn∗−1 (k) dk, n N (α) Γ (α) N (α) 0 ∫ t  ∗  (1 − α)  ∗    α Q (t ) ≤ ω4 Qn−1 (t ) + ω4 (t − k)α−1 Qn∗−1 (k) dk, n N (α) Γ (α) N (α) 0 ∫ t    ∗  (1 − α)  ∗  α Z (t ) ≤ ω5 Zn−1 (t ) + ω5 (t − k)α−1 Zn∗−1 (k) dk. n N (α) Γ (α) N (α) 0 Then, the proof is complete.

(51)

□

Theorem 3.6 (Existence of the Solution). The system given by Eq. (39) has a solution under the conditions that we can find tmax satisfying ( ) ωj

N(α)

(1 − α) +

α tmax

Γ (α)

< 1 for j ∈ {1, 2, 3, 4, 5}.

388

M. Yavuz and E. Bonyah / Physica A 525 (2019) 373–393

Proof. Considering the functions Sh (t ) , Ih (t ) , R (t ) , Sv (t ) and Iv (t ) are bounded, we can obtain the following relations from Eq. (51):

( )}n α ω1 tmax 1 − α) + , ( n N (α) Γ (α) ( )}n { α  ∗  tmax K (t ) ≤ ∥Ih (0)∥ ω2 1 − α) + , ( n N (α) Γ (α) ( )}n { α  ∗  tmax P (t ) ≤ ∥R (0)∥ ω3 1 − α) + , ( n N (α) Γ (α) ( )}n { α  ∗  tmax Q (t ) ≤ ∥Sv (0)∥ ω4 1 − α) + , ( n N (α) Γ (α) { ( )}n α  ∗  tmax Z (t ) ≤ ∥Iv (0)∥ ω5 1 − α) + . ( n N (α) Γ (α)

 ∗  G (t ) ≤ ∥Sh (0)∥

{

(52)

Hence, the functions G∗n (t ) , Kn∗ (t ) , Pn∗ (t ) , Qn∗ (t ) and Zn∗ (t ) satisfy Eq. (52) exist and are smooth. Moreover, to show that the functions in Eq. (52) are the solutions of Eq. (39), we assume that

Sh (t ) − Sh (0) = Sh(n) (t ) − Θ1∗(n) (t ) , Ih (t ) − Ih (0) = Ih(n) (t ) − Θ2∗(n) (t ) , (53)

R (t ) − R (0) = R(n) (t ) − Θ3∗(n) (t ) , Sv (t ) − Sv (0) = Sv (n) (t ) − Θ4∗(n) (t ) , Iv (t ) − Iv (0) = Iv (n) (t ) − Θ5∗(n) (t ) ,

∗ where Θ1∗(n) (t ) , Θ2∗(n) (t ) , Θ3∗(n) (t ) , Θ4∗(n) (t ) and Θ we must show  5(n∗) (t ) are  reminder  termsof series  solution. Then,   Θ  → 0, Θ ∗ (t ) → 0, Θ ∗ (t ) → 0, Θ ∗ (t ) → that these terms approach to zero at infinity, that is, t ( ) 1(∞) 2(∞) 3(∞) 4(∞)  ∗  0 and Θ5(∞) (t ) → 0. Thus, we get for the term Θ1∗(n) (t ),

  ∗  ) ( Θ (t ) ≤ (1 − α) τ2 (t , Sh ) − τ2 t , Sh(n−1)  1(n) N (α) ∫ t  ( ) α + (t − k)α−1 τ2 (k, Sh ) − τ2 k, Sh(n−1)  dk Γ (α) N ( (α) 0 )   tα ω1 Sh − Sh(n−1)  . ≤ (1 − α) + N (α) Γ (α)

(54)

Following this step recursively, we have

 ∗  Θ (t ) ≤ W ωn 1(n) 1

{

1 N (α)

( (1 − α) +

tα

Γ (α)

)}n+1

,

(55)

α where W = Sh − Sh(n−1) . In this case, we get the following at tmax :





 ∗  Θ (t ) ≤ W ωn 1(n) 1

{

( )}n+1 tα . (1 − α) + max N (α) Γ (α) 1

(56)

 ∗  When we take  ∗ the limit  of both  ∗sides as n tends to∗ infinity,  we have Θ1∗(∞) (t ) → 0. Following the similar steps, we can obtain Θ2(∞) (t ) → 0, Θ3(∞) (t ) → 0, Θ4(∞) (t ) → 0 and Θ5(∞) (t ) → 0. Hence the proof is complete. □ 



Theorem 3.7 (Uniqueness of the Solution). The model constructed with the Atangana–Baleanu fractional operator and given by Eq. (39) has a unique solution.

M. Yavuz and E. Bonyah / Physica A 525 (2019) 373–393

389

Proof. To prove this, we assume that other solutions for the mentioned system (39) are Sh⊗ (t ) , Ih⊗ (t ) , R⊗ (t ) , Sv⊗ (t ) and Iv⊗ (t ) ; then,

∫ t ( )) ( ( )) α (1 − α) ( τ2 (t , Sh ) − τ2 t , Sh⊗ + (t − k)α−1 τ2 (k, Sh ) − τ2 k, Sh⊗ dk, N (α) Γ (α) N (α) 0 ∫ t ( )) ( )) ( ( 1 − α) α ( ⊗ ⊗ Ih (t ) − Ih (t ) = ϕ2 (t , Ih ) − ϕ2 t , Ih + (t − k)α−1 ϕ2 (k, Ih ) − ϕ2 k, Ih⊗ dk, N (α) Γ (α) N (α) 0 ∫ t )) )) ( ( ( α (1 − α) ( R (t ) − R⊗ (t ) = ξ2 (t , R) − ξ2 t , R⊗ + (t − k)α−1 ξ2 (k, R) − ξ2 k, R⊗ dk, N (α) Γ (α) N (α) 0 ∫ t )) )) ( ( ( ( 1 − α) α ( ⊗ ⊗ Sv (t ) − Sv (t ) = η2 (t , Sv ) − η2 t , Sv + (t − k)α−1 η2 (k, Sv ) − η2 k, Sv⊗ dk, N (α) Γ (α) N (α) 0 ∫ t )) )) ( ( ( α (1 − α) ( Iv (t ) − Iv⊗ (t ) = Φ2 (t , Iv ) − Φ2 t , Iv⊗ + (t − k)α−1 Φ2 (k, Iv ) − Φ2 k, Iv⊗ dk. N (α) Γ (α) N (α) 0

Sh (t ) − Sh⊗ (t ) =

(57)

Applying the norm to both sides of Eq. (57), we obtain

∫ t     ( ) ( ) α Sh (t ) − S ⊗ (t ) ≤ (1 − α) τ2 (t , Sh ) − τ2 t , S ⊗  + (t − k)α−1 τ2 (k, Sh ) − τ2 k, Sh⊗  dk, h h N (α) Γ (α) N (α) 0 ∫ t     ) ( ) ( α Ih (t ) − I ⊗ (t ) = (1 − α) ϕ2 (t , Ih ) − ϕ2 t , I ⊗  + (t − k)α−1 ϕ2 (k, Ih ) − ϕ2 k, Ih⊗  dk, h h N (α) Γ (α) N (α) 0 ∫ t     ( ) ( ) 1 − α) α ( R (t ) − R⊗ (t ) = ξ2 (t , R) − ξ2 t , R⊗  + (58) (t − k)α−1 ξ2 (k, R) − ξ2 k, R⊗  dk, N (α) Γ (α) N (α) 0 ∫ t     ( ) ( ) α Sv (t ) − S ⊗ (t ) = (1 − α) η2 (t , Sv ) − η2 t , S ⊗  + (t − k)α−1 η2 (k, Sv ) − η2 k, Sv⊗  dk, v v N (α) Γ (α) N (α) 0 ∫ t     ( ) ( ) 1 − α) α ( Iv (t ) − I ⊗ (t ) = Φ2 (t , Iv ) − Φ2 t , I ⊗  + (t − k)α−1 Φ2 (k, Iv ) − Φ2 k, Iv⊗  dk. v v N (α) Γ (α) N (α) 0 Taking into consideration Theorems 3.5 and 3.6, we obtain

      tα Sh (t ) − S ⊗ (t ) ≤ (1 − α) ω1 Sh (t ) − S ⊗ (t ) + ω1 Sh (t ) − S ⊗ (t ) , h h h N (α) Γ (α) N (α)      (1 − α)  tα ⊗ ⊗ Ih (t ) − I ⊗ (t ) , Ih (t ) − I (t ) = ω2 Ih (t ) − Ih (t ) + ω2 h h N (α) Γ (α) N (α)       1 − α) tα ( R (t ) − R⊗ (t ) = R (t ) − R⊗ (t ) , ω3 R (t ) − R⊗ (t ) + ω3 N (α) Γ (α) N (α)   (1 − α)     tα ⊗ ⊗ Sv (t ) − S (t ) = Sv (t ) − S ⊗ (t ) , ω4 Sv (t ) − Sv (t ) + ω4 v v N (α) Γ (α) N (α)       tα 1 − α) ( Iv (t ) − I ⊗ (t ) = Iv (t ) − I ⊗ (t ) . ω5 Iv (t ) − Iv⊗ (t ) + ω5 v v N (α) Γ (α) N (α)

(59)

The solution functions in Eq. (59) hold the following inequalities:

( { })   tα Sh (t ) − S ⊗ (t ) 1 − ω1 1 − α) − ≤ 0, ( h N (α) Γ (α) ( { })   tα Ih (t ) − I ⊗ (t ) 1 − ω2 ≤ 0, (1 − α) − h N (α) Γ (α) { }) (   tα R (t ) − R⊗ (t ) 1 − ω3 ≤ 0, (1 − α) − N (α) Γ (α) ( { })   tα Sv (t ) − S ⊗ (t ) 1 − ω4 1 − α) − ≤ 0, ( v N (α) Γ (α) ( { })   tα Iv (t ) − I ⊗ (t ) 1 − ω5 ≤ 0. (1 − α) − v N (α) Γ (α)

(60)

From the last equation, it can be concluded that Sh (t ) = Sh⊗ (t ) , Ih (t ) = Ih⊗ (t ) , R (t ) = R⊗ (t ) , Sv (t ) = Sv⊗ (t ) , Iv (t ) = Iv⊗ (t ) .

(61)

390

M. Yavuz and E. Bonyah / Physica A 525 (2019) 373–393

These results show the uniqueness of the solution of the system (39).

□

We are now going to obtain the approximate solution of the system (39) using the Sumudu transform presented by Eq. (10). Applying it to both sides of Eq. (39), we get

ST {Λ − βh Iv Sh + γ R − µh Sh } (s) =

[ ] α N (α) Γ (α + 1) uα Eα − (ST {Sh (t )} − Sh (0)) , 1−α 1−α

[ ] α N (α) Γ (α + 1) uα ST {βh Iv Sh − (δ + µh ) Ih } (s) = Eα − (ST {Ih (t )} − Ih (0)) , 1−α 1−α [ ] α N (α) Γ (α + 1) uα ST {δ Ih − (γ + µh ) R} (s) = Eα − (ST {R (t )} − R (0)) , 1−α 1−α

(62)

[ ] α N (α) Γ (α + 1) uα ST {Λv − βv Ih Sv − µv Sv } (s) = Eα − (ST {Sv (t )} − Sv (0)) , 1−α 1−α ST {βv Ih Sv − µv Iv } (s) =

[ ] uα α N (α) Γ (α + 1) Eα − (ST {Iv (t )} − Iv (0)) . 1−α 1−α

By considering the last equation in terms of unknown functions, we obtain 1−α [ α ] ST {Λ − βh Iv Sh + γ R − µh Sh } , α N (α) Γ (α + 1) Eα − 1u−α

ST {Sh (t )} = Sh (0) + ST {Ih (t )} = Ih (0) + ST {R (t )} = R (0) +

1−α

[ α ] ST {βh Iv Sh − (δ + µh ) Ih } , α N (α) Γ (α + 1) Eα − 1u−α 1−α

[ α ] ST {δ Ih − (γ + µh ) R} , α N (α) Γ (α + 1) Eα − 1u−α 1−α [ α ] ST {Λv − βv Ih Sv − µv Sv } , α N (α) Γ (α + 1) Eα − 1u−α

ST {Sv (t )} = Sv (0) + ST {Iv (t )} = Iv (0) +

(63)

1−α

[ α ] ST {βv Ih Sv − µv Iv } . α N (α) Γ (α + 1) Eα − 1u−α

Taking the inverse Sumudu transform of Eq. (63), we have

{

1−α

}

[ α ] ST {Λ − βh Iv Sh + γ R − µh Sh } (t ) , α N (α) Γ (α + 1) Eα − 1u−α { } 1−α −1 [ Ih (t ) = Ih (0) + ST α ] ST {βh Iv Sh − (δ + µh ) Ih } (t ) , α N (α) Γ (α + 1) Eα − 1u−α { } 1−α −1 [ R (t ) = R (0) + ST α ] ST {δ Ih − (γ + µh ) R} (t ) , α N (α) Γ (α + 1) Eα − 1u−α { } 1−α −1 [ ] Sv (t ) = Sv (0) + ST ST {Λv − βv Ih Sv − µv Sv } (t ) , α α N (α) Γ (α + 1) Eα − 1u−α { } 1−α −1 [ Iv (t ) = Iv (0) + ST α ] ST {βv Ih Sv − µv Iv } (t ) . α N (α) Γ (α + 1) Eα − 1u−α Sh (t ) = Sh (0) + ST

−1

M. Yavuz and E. Bonyah / Physica A 525 (2019) 373–393

391

The following recursive formula is then proposed Sh(n) (t ) = Sh(n−1) (0)

{ +ST

−1

1−α

} { } ] ST Λ − βh Iv(n−1) Sh(n−1) + γ R(n−1) − µh Sh(n−1) (t ) ,

[ α α N (α) Γ (α + 1) Eα − 1u−α } { } { 1−α −1 [ Ih(n) (t ) = Ih(n−1) (0) + ST (t ) , α ] ST βh Iv(n−1) Sh(n−1) − (δ + µh ) Ih(n−1) α N (α) Γ (α + 1) Eα − 1u−α { } } { 1−α −1 [ R(n) (t ) = R(n−1) (0) + ST (t ) , α ] ST δ Ih(n−1) − (γ + µh ) R(n−1) α N (α) Γ (α + 1) Eα − 1u−α { } { } 1−α −1 [ Sv(n) (t ) = Sv(n−1) (0) + ST (t ) , α ] ST Λv − βv Ih(n−1) Sv(n−1) − µv Sv(n−1) α N (α) Γ (α + 1) Eα − 1u−α { } { } 1−α −1 [ Iv(n) (t ) = Iv(n−1) (0) + ST (t ) , α ] ST βv Ih(n−1) Sv(n−1) − µv Iv(n−1) α N (α) Γ (α + 1) Eα − 1u−α

(64)

where Sh(0) (t ) = Sh (0) , Ih(0) (t ) = Ih (0) , R(0) (t ) = R (0) , Sv(0) (t ) = Sv (0) , Iv(0) (t ) = Iv (0). The approximate solution is assumed to obtain as a limit when n tends to infinity Sh (t ) = lim Sh(n) (t ) , Ih (t ) = lim Ih(n) (t ) , R (t ) = lim R(n) (t ) , n→∞

n→∞

n→∞

(65)

Sv (t ) = lim Sv(n) (t ) , Iv (t ) = lim Iv(n) (t ). n→∞

n→∞

3.4. The numerical scheme and simulation with the Atangana–Baleanu operator As the second part of the simulation, numerical simulation results based recently ABC numerical scheme developed by Atangana and Toufik [11] has been employed with the following parameter values having different fractional order Λ = 0.5, µh = 0.014, βh = 0.09, βv = 0.6, δ = 0.08, γ = 0.008, µh = 0.008. The initial conditions employed for the ABC scheme are Sh (0) = 0.5, Ih (0) = 0.75, R(0) = 0.01, Sv (0) = 0.6, Iv (0) = 0.05. The simulation time used for this work is 800 s and the step size employed in evaluating the approximate solutions is h = 0.002. Fig. 2 depicts the numerical simulation results based on different fractional order values using ABC numerical scheme as in Atangana and Toufik [11]. Similarly, as the fractional order values get closer to 1 then the convergence gets reduced. This indicates that the newly developed operator is very effective and has the crossover waiting time behaviour that guaranteed accurate prediction. A major reason of obtaining accurate prediction hinged on the new scheme is that there is no artificial singularity that is incorporated in Eq. (39). In Fig. 2 (a–d), it is worth noting that as the fractional order derivative increases from α = 0.5 to α = 1, the number of humans infected with schistosomiasis also gets asymptotically zero. The case of snails infected with schistosomiasis is similar to humans. This therefore suggests the accurate prediction depends on the fractional order and the operator used. It can also be seen that there is a new asymptotic behaviour as a result of a strong memory effect of ABC operator and this could not be observed when alpha is 1. 4. Conclusions In this work, the schistosomiasis fractional order dynamical model for the transmission of the disease was analysed. The model was considered in two folds through the concept of Caputo–Fabrizio–Caputo and Atangana–Baleanu–Caputo fractional order derivatives. The model in Caputo–Fabrizio sense was solved using an iterative scheme based on Laplace transform and also the model in Atangana–Baleanu–Caputo sense solution obtained through Sumudu transform. The detailed analysis of iterative techniques and uniqueness of special solutions for the two operators used in the model are presented. The numerical solution for both operators pointed out interesting results when alpha value is 0.9. This is because two operators have crossover property which makes it possible to capture complex phenomena and give accurate predictions. The results obtained in the numerical simulations are also in agreement with the conclusion drawn in the study of [47]. The numerical simulations presented in this paper are in perfect agreement with the conclusion of the study of [17]. This crossover behaviour of these new fractional differential operators is due to their capacity of not obeying the index law imposed in fractional calculus.

392

M. Yavuz and E. Bonyah / Physica A 525 (2019) 373–393

Fig. 2. Numerical results for system (12) with respect to the Atangana–Baleanu derivative.

References [1] B. Oyinloye, F. Adenowo, N. Gxaba, A. Kappo, The promise of antimicrobial peptides for treatment of human schistosomiasis, Current Drug targets 15 (2014) 852–859. [2] A.F. Adenowo, B.E. Oyinloye, B.I. Ogunyinka, A.P. Kappo, Impact of human schistosomiasis in sub-saharan africa, Braz. J. Infec. Dis. 19 (2015) 196–205. [3] E.T. Chiyaka, G. Magombedze, L. Mutimbu, Modelling within host parasite dynamics of schistosomiasis, Comput. Math. Methods Med. 11 (2010) 255–280. [4] H.M. Yang, Comparison between schistosomiasis transmission modelings considering acquired immunity and age-structured contact pattern with infested water, Math. Biosci. 184 (2003) 1–26. [5] T.D. Mangal, S. Paterson, A. Fenton, Predicting the impact of long-term temperature changes on the epidemiology and control of schistosomiasis: a mechanistic model, PLoS one 3 (2008) e1438. [6] K.O. Okosun, R. Smith, Optimal control analysis of malaria–schistosomiasis co-infection dynamics, Math. Biosci. Eng. 14 (2017) 377–405. [7] Z. Chen, L. Zou, D. Shen, W. Zhang, S. Ruan, Mathematical modelling and control of schistosomiasis in hubei province, China, Acta Tropica 115 (2010) 119–125. [8] D. Baleanu, K. Diethelm, E. Scalas, J.J. Trujillo, Fractional Calculus: Models and Numerical Methods, World Scientific, 2012. [9] B.S.T. Alkahtani, Atangana-batogna numerical scheme applied on a linear and non-linear fractional differential equation, Eur. Phys. J. Plus 133 (2018) 111. [10] A. Atangana, K.M. Owolabi, New numerical approach for fractional differential equations, Math. Model. Nat. Phenom. 13 (2018) 3. [11] M. Toufik, A. Atangana, New numerical approximation of fractional derivative with non-local and non-singular kernel: application to chaotic models, Eur. Phys. J. Plus 132 (2017) 444. [12] A. Atangana, D. Baleanu, New fractional derivatives with nonlocal and non-singular kernel: theory and application to heat transfer model, Thermal Sci. 20 (2016) 763–769. [13] I. Koca, Modelling the spread of ebola virus with atangana-baleanu fractional operators, Eur. Phys. J. Plus 133 (2018) 100. [14] A. Yokus, T.A. Sulaiman, H.M. Baskonus, S.P. Atmaca, On the exact and numerical solutions to a nonlinear model arising in mathematical biology, Proc., ITM Web Conf.: EDP Sci. 01061 (2018).

M. Yavuz and E. Bonyah / Physica A 525 (2019) 373–393

393

[15] M.S. Asl, M. Javidi, Novel algorithms to estimate nonlinear fdes: applied to fractional order nutrient-phytoplankton–zooplankton system, J. Comput. Appl. Math. 339 (2018) 193–207. [16] A. Atangana, J. Gómez-Aguilar, Decolonisation of fractional calculus rules: breaking commutativity and associativity to capture more natural phenomena, Eur. Phys. J. Plus 133 (2018) 1–22. [17] A. Atangana, J. Gómez-Aguilar, Fractional derivatives with no-index law property: application to chaos and statistics, Chaos Solitons Fractals 114 (2018) 516–535. [18] D. Baleanu, A. Jajarmi, E. Bonyah, M. Hajipour, New aspects of poor nutrition in the life cycle within the fractional calculus, Adv. Difference Equ. 2018 (2018) 230. [19] A. Atangana, Non validity of index law in fractional calculus: a fractional differential operator with markovian and non-markovian properties, Physica A 505 (2018) 688–706. [20] B.S.T. Alkahtani, A. Atangana, I. Koca, Novel analysis of the fractional zika model using the adams type predictor-corrector rule for non-singular and non-local fractional operators, J. Nonlinear Sci. Appl 10 (2017) 3191–3200. [21] M. Yavuz, N. Ozdemir, H.M. Baskonus, Solutions of partial differential equations using the fractional operator involving mittag-leffler kernel, Eur. Phys. J. Plus 133 (2018) 215. [22] M. Yavuz, N. Özdemir, European vanilla option pricing model of fractional order without singular kernel, Fractal Fract. 2 (2018) 3. [23] F. Evirgen, M. Yavuz, An alternative approach for nonlinear optimization problem with caputo-fabrizio derivative, Proc., ITM Web Conf.: EDP Sci. 01009 (2018). [24] E. Bonyah, K. Badu, S.K. Asiedu-Addo, Optimal control application to an ebola model, Asian Pacific J. Tropical Biomed. 6 (2016) 283–289. [25] J.F. Gómez-Aguilar, M.G. López-López, V.M. Alvarado-Martí nez, D. Baleanu, H. Khan, Chaos in a cancer model via fractional derivatives with exponential decay and mittag-leffler law, Entropy 19 (2017) 681. [26] J. Hristov, Steady-state heat conduction in a medium with spatial non-singular fading memory: derivation of caputo-fabrizio space-fractional derivative with jeffrey’s kernel and analytical solutions, Thermal Sci. 21 (2017) 827–839. [27] A. Coronel-Escamilla, J.F. Gómez-Aguilar, D. Baleanu, T. Córdova-Fraga, R.F. Escobar-Jiménez, V.H. Olivares-Peregrino, M.M.A. Qurashi, Bateman–Feshbach Tikochinsky and Caldirola–Kanai Oscillators With new fractional differentiation, Entropy 19 (2017) 55. [28] A. Atangana, J. Gómez-Aguilar, Hyperchaotic behaviour obtained via a nonlocal operator with exponential decay and mittag-leffler laws, Chaos Solitons Fractals 102 (2017) 285–294. [29] J. Gómez-Aguilar, H. Yepez-Martinez, R. Escobar-Jiménez, C. Astorga-Zaragoza, J. Reyes-Reyes, Analytical and numerical solutions of electrical circuits described by fractional derivatives, Appl. Math. Model. 40 (2016) 9079–9094. [30] F. Jarad, T. Abdeljawad, Z. Hammouch, On a class of ordinary differential equations in the frame of atangana–baleanu fractional derivative, Chaos Solitons Fractals 117 (2018) 16–20. [31] K. Ait Touchent, Z. Hammouch, T. Mekkaoui, F. Belgacem, Implementation and convergence analysis of homotopy perturbation coupled with sumudu transform to construct solutions of local-fractional pdes, Fractal Fractional 2 (2018) 22. [32] S. Uçar, E. Uçar, N. Özdemir, Z. Hammouch, Mathematical analysis and numerical simulation for a smoking model with atangana–baleanu derivative, Chaos Solitons Fractals 118 (2019) 300–306. [33] M. Yavuz, N. Ozdemir, Comparing the new fractional derivative operators involving exponential and mittag-leffler kernel, Discrete Contin. Dyn. Syst. Ser. 13 (2020) in press. [34] D. Avc, A. Yetim, Cauchy and source problems for an advection-diffusion equation with atangana–baleanu derivative on the real line, Chaos Solitons Fractals 118 (2019) 361–365. [35] D. Avc, A. Yetim, Analytical solutions to the advection-diffusion equation with the atangana-baleanu derivative over a finite domain, Balkesir Univ. Fen Bilimleri Enstitüsü Dergisi 20 (2018) 382–395. [36] M. Modanli, Difference scheme to the fractional telegraph model with atangana-baleanu-caputo derivative, Chaos 29 (2019) in press. [37] J. Singh, D. Kumar, D. Baleanu, On the analysis of chemical kinetics system pertaining to a fractional derivative with mittag-leffler type kernel, Chaos 27 (2017) 103113. [38] J. Singh, D. Kumar, Z. Hammouch, A. Atangana, A fractional epidemiological model for computer viruses pertaining to a new fractional derivative, Appl. Math. Comput. 316 (2018) 504–515. [39] D. Kumar, J. Singh, D. Baleanu, S. Rathore, Analysis of a fractional model of the ambartsumian equation, Eur. Phys. J. Plus 133 (2018) 259. [40] J. Singh, D. Kumar, D. Baleanu, S. Rathore, An efficient numerical algorithm for the fractional drinfeld–sokolov–wilson equation, Appl. Math. Comput. 335 (2018) 12–24. [41] J. Gómez-Aguilar, A. Atangana, New insight in fractional differentiation: power, exponential decay and mittag-leffler laws and applications, Eur. Phys. J. Plus 132 (2017) 13. [42] A. Atangana, J. Gómez-Aguilar, A new derivative with normal distribution kernel: theory, methods and applications, Physica A 476 (2017) 1–14. [43] M. Caputo, M. Fabrizio, A new definition of fractional derivative without singular kernel, Progress Fract. Differentiation Appl. 1 (2015) 1–13. [44] A. Atangana, B.S.T. Alkahtani, New model of groundwater flowing within a confine aquifer: application of caputo-fabrizio derivative, Arab. J. Geosci. 9 (2016) 8. [45] J. Losada, J.J. Nieto, Properties of a new fractional derivative without singular kernel, Progr. Fract. Differ. Appl 1 (2015) 87–92. [46] A. Atangana, I. Koca, Chaos in a simple nonlinear system with atangana–baleanu derivatives with fractional order, Chaos Solitons Fractals 89 (2016) 447–454. [47] A. Atangana, S. Jain, The.Role.of.Power. Decay, The role of power decay exponential decay and mittag-leffler function’s waiting time distribution: application of cancer spread, Physica A 512 (2018) 330–351. [48] M. Saeedian, M. Khalighi, N. Azimi-Tafreshi, G. Jafari, M. Ausloos, Memory effects on epidemic evolution: The susceptible-infected-recovered epidemic model, Phys. Rev. E 95 (2017) 022409.