Dynamical behaviour of fractional order tumor model with Caputo and conformable fractional derivative

Contents

lists

available

at

ScienceDirect

Chaos,

Solitons

and

Fractals

Nonlinear

Science,

and

Nonequilibrium

and

Complex

Phenomena

journal

homepage:

www.elsevier.com/locate/chaos

Frontiers

Dynamical

behaviour

of

fractional

order

tumor

model

with

Caputo

and

conformable

fractional

derivative

R

Ercan

Balcı

a

,

∗

_,

_˙Ilhan

_Öztürk

a

_,

_Senol

_Kartal

b

a Department of Mathematics, Erciyes University, Kayseri 38039, Turkey

b Department of Science and Mathematics Education, Nevsehir Haci Bektas Veli University, Nevsehir 50300, Turkey

a

r

t

i

c

l

e

i

n

f

o

Article history:

Received 27 January 2019 Revised 20 March 2019 Accepted 28 March 2019 Available online 2 April 2019 Keywords:

Tumor-immune interaction Fractional order

Piecewise-constant arguments Neimark-Sacker bifurcation Conformable fractional derivative

a

b

s

t

r

a

c

t

Inthispaper,tumor-immunesysteminteractionhasbeenconsideredbytwofractionalordermodels.The firstand thesecondmodelconsistofsystemoffractionalorderdifferentialequationswithCaputoand conformablefractional derivativerespectively.Firstofall,the stabilityofthe equilibriumpointsofthe firstmodelisstudied.Then,adiscretizationprocessisappliedtoobtainadiscreteversionofthesecond modelwhereconformablefractionalderivativeistakenintoaccount.Indiscretemodel,weanalyzethe stabilityoftheequilibriumpointsandprovetheexistenceofNeimark-Sackerbifurcationdependingon theparameter

σ

.Moreover,thedynamicalbehavioursofthemodelsarecomparedwitheachotherand weobservethat thediscrete version ofconformable fractional ordermodel exhibits chaotic behavior. Finally,numericalsimulationsarealsopresentedtoillustratetheanalyticalresults.

© 2019ElsevierLtd.Allrightsreserved.

1.

Introduction

Research

on

fractional

calculus

has

gained

much

interest

over

the

past

decades

and

furthermore,

fractional

order

differential

equations

has

been

well

applied

to

many

research

field

such

as

physics

[1,2]

,

chemistry

[3]

,

medicine

[5]

,

finance

[4]

and

engineer-ing

[6]

.

These

equations

are

also

widely

used

to

model

biological

phenomenon

and

there

are

successful

applications

in

this

field.

It

is

also

demonstrated

that

biological

models

constructed

by

fractional

order

differential

equations

exhibit

more

realistic

results

compar-ing

with

integer

order

counterpart

[7–9]

.

This

is

for

the

reason

that

fractional

order

derivatives

involve

memory

and

which

is

quite

fa-vorable

to

work

on

biological

processes.

Trying

to

generalize

notion

of

differentiation

to

arbitrary

or-der

brings

out

several

approaches.

The

most

knowns

are

Riemann-Liouville,

Caputo

and

Grünwald-Letnikov

definitions

[10]

.

In

ad-dition

to

these

definitions,

a

new

definition

called

“conformable

fractional

derivative” has

been

introduced

by

Khalil

et

al.

in

2014

[13]

.

According

to

this

definition,

the

left

fractional

derivative

start-ing

from

a

of

the

function

f

:

[

a

,

∞

)

→

∞

of

order

0

<

α

≤ 1

has

its

R This work was supported by Research Fund of the Erciyes University. Project Number: FDK-2018-8119.

∗ _{Corresponding author.}

E-mail addresses: ercanbalci@erciyes.edu.tr (E. Balcı), ozturki@erciyes.edu.tr ( ˙I. Öztürk), senol.kartal@nevsehir.edu.tr (S. Kartal).

own

limit-based

definition

as

follows:

(

T

a

α

f

)(

t

)

=

lim

_→0

f

(

t

+



(

t

− a

)

1−α

₎

_{− f}

₍

_t

₎



provided

the

limit

exist.

The

right

fractional

derivative

of

order

0

_<

α

_{≤ 1}

terminating

at

b

of

f

is

defined

by

(

b

α

T

f

)(

t

)

=

− lim

_→0

f

(

t

+



(

b

− t

)

1−α

₎

_{− f}

₍

_t

₎



.

Note

also

that

if

f

is

differentiable

in

usual

sense,

we

have

the

fol-lowing

equalities:

(

T

_αa

f

)(

t

)

=

(

t

− a

)

1−α

_f



₍

_t

₎

_,

₍

b

α

T

f

)(

t

)

=

(

t

− b

)

1−α

f



(

t

)

.

(1)

Unlike

classical

fractional

derivative

definitions,

conformable

fractional

derivative

share

some

basic

properties

with

integer

order

derivative.

In

[14]

,

Abdeljawad

introduced

conformable

versions

of

exponential

functions,

Gronwall’s

inequality,

integration

by

parts,

Taylor

series

expansion

and

Laplace

transforms.

Moreover,

biolog-ical

and

physical

applications

of

conformable

fractional

derivative

can

be

found

in

[15–19]

.

Together

with

these

definitions,

there

are

numerous

approaches

being

studied

by

applied

mathematicians.

Researchers

trying

to

find

the

most

efficient

approach

while

constructing

or

adjusting

their

models.

Meanwhile,

some

numerical

methods

such

as

Ado-mian

decomposition

[20]

,

Adams-type

predictor-corrector

[21]

and

homotopy

perturbation

method

[22]

are

also

developed

to

solve

fractional

differential

equations

[23]

.

In

addition

to

the

these

nu-merical

methods,

several

papers

have

used

piecewise

constant

ap-proximations

to

discretize

fractional

order

differential

equations

https://doi.org/10.1016/j.chaos.2019.03.032

[27–29]

.

In

[17]

,

Kartal

and

Gurcan

considered

the

conformable

fractional

logistic

equation

with

piecewise

constant

arguments

with

adopting

the

method

presented

by

Gopalsamy

in

[30]

.

The

tumor-immune

interaction

model

studied

in

this

paper

originated

from

the

following

model

presented

by

Kuznetsov

et

al.

[24]

:



dE dt

=

s

+

F

(

E,

T

)

− mET

− dE,

dT dt

=

aT

(

1

− bT

)

− nET

,

(2)

where

F

(

E

,

T

)

=

_g

pE

₊

T

_T

.

In

study

[25]

,

Galach

suppose

that

F

(

E

,

T

)

=

θ

ET

and

thus

model

(2)

takes

the

form



dE

dt

=

s

+

α

1

ET

− dE,

dT

dt

=

aT

(

1

− bT

)

− nET

,

(3)

where

α1

=

θ

− m.

Then,

the

dimensionless

form

of

the

model

(3)

can

be

obtained

as



dx dt

=

σ

+

wxy

−

δ

x

,

dy dt

=

γ

y

(

1

−

β

y

)

− xy,

(4)

where

x

and

y

denote

the

dimensionless

density

of

ECs

and

TCs

re-spectively,

x

=

E

/E

0

,

y

=

T

/T

0

,

,

γ

=

nTa0

β

=

bT

0

,

δ

=

d nT0

,

σ

=

s nE0T0

,

ω

=

α1

/n,

E

0

and

T

0

are

the

initial

conditions.

In

the

study

[26]

,

the

fractional

order

form

of

the

model

(4)

is

considered

with

the

Caputo

sense

as

follows:



_D

_α

_x

₍

_t

₎

₌

_σ

₊

_ω

_xy

₋

_δ

_x,

D

α

y

(

t

)

=

γ

y

(

1

−

β

y

)

− xy

,

(5)

with

initial

conditions

x

(

0

)

=

x

0

≥ 0

and

y

(

0

)

=

y

0

≥ 0

.

Thus,

con-formable

fractional

order

version

of

the

model

(4)

is

given

as

fol-lows:



T

_α

E

(

t

)

=

σ

+

ω

ET

−

δ

E,

T

_α

T

(

t

)

=

γ

T

(

1

−

β

T

)

− ET

(6)

where

we

take

α

∈

(0,

1)

as

conformable

fractional

order.

The

aim

of

this

study

is

to

investigate

the

dynamical

behaviour

of

model

(5)

and

model

(6)

and

to

compare

the

obtained

results.

2.

Dynamical

behaviour

of

fractional

order

tumor

model

Here,

we

adopt

fractional

order

model

(5)

with

Caputo

frac-tional

derivative

which

is

defined

by

D

αa

f

(

t

)

=



t a

f

n

₍

_x

₎

(

t

− x

)

α−n+1

dx.

2.1.

Stability

analysis

Theorem

1

[11

,

12]

.

Consider

the

system

D

α_a

f

(

t

)

=

f

(

t

,

X

(

t

))

,

X

(

t

0

)

=

X

0

.

(7)

Let

J

(

X

∗

)

denote

the

Jacobian

matrix

of

the

system

(7)

evaluated

at

the

equilibrium

point

X

∗

.

(1)

The

equilibrium

point

X

∗

is

locally

asymptotically

stable

if

and

only

if

all

the

eigenvalues

λ

_i

,

i

₌

1

_,

2

_,

_.

_.

_.

_,

n

of

J

(

X

∗

)

satisfy

|

arg

(

λ

i

)

|

>

απ2

,

(2)

The

equilibrium

point

X

∗

is

stable

if

all

the

eigenvalues

λ

i

,

i

₌

1

_,

2

_,

_.

_.

_.

_,

n

of

J

(

X

∗

)

satisfy

|

arg

(

λ

_i

)

|

_≥

απ₂

and

eigenvalues

with

|

arg

(

λ

i

)

|

=

απ2

have

the

same

geometric

and

algebraic

multiplicity,

and

(3)

The

equilibrium

point

X

∗

is

unstable

if

and

only

if

there

exist

eigenvalues

λ

i

for

some

i

=

1

,

2

,

.

.

.

,

n

of

J

(

X

∗

)

satisfy

|

arg

(

λ

i

)

|

<

απ2

.

An

equilibrium

point

of

model

(5)

is

obtained

by

solving

the

following

system:



D

α

x

(

t

)

=

0

,

D

α

y

(

t

)

=

0

.

That

is



σ

+

ω

xy

−

δ

x

=

0

,

γ

y

(

1

−

β

y

)

− xy

=

0

.

Hence,

we

have

two

equilibrium

points:

i.

The

tumor-free

equilibrium

point

E

0

=

(

σ_δ

,

0

)

,

ii.

The

coexistence

equilibrium

point

E

1

=

(

E

,

T

)

=

(

γ (−βδ+ω)+ √  2ω

,

γ (βδ+ω)− √  2γ βω

)

where



=

4

γ βσω

+

γ

2

₍

_βδ

₋

_ω

₎

2

_.

_Under

_the

_condition

_σ

_<

_{γ δ}

_,

the

coexistence

equilibrium

point

is

always

positive.

Theorem

2.

For

the

equilibrium

point

E

0

=

(

σ_δ

,

0

)

of

model

(5)

,

the

following

results

holds

true;

i.

If

σ

>

γ δ

,

then

E

0

is

locally

asymptotically

stable,

ii.

If

σ

<

γ δ

,

then

E

0

is

unstable

and

is

a

saddle

point.

Proof.

The

Jacobian

matrix

of

the

model

(5)

evaluated

at

equilib-rium

point

E

0

is

given

by

J

(

E

0

)

=



−

δ

σ ω_δ

0

γ

−

σ_δ



.

Hence

the

eigenvalues

of

J

(

E

0

)

are

λ1

=

−

δ

and

λ2

=

γ

−

σ_δ

.

Since

λ1

<

0,

we

have

arg

(

λ1

)

₌

π

which

satisfies

|

arg

(

λ1

)

|

_>

απ₂

_.

If

σ

>

γ δ

,

then

λ2

<

0

and

arg

(

λ2

)

=

π

which

results

in

|

arg

(

λ2

)

|

>

απ

2

.

According

to

Theorem

1

,

equilibrium

point

E

0

is

locally

asymp-totically

stable

if

σ

_>

γ δ

.

If

σ

_<

γ δ

,

then

λ2

_>

0.

Hence

arg

(

λ2

)

₌

0

_,

which

always

satisfies

|

arg

(

λ2

)

|

<

απ2

.

By

Theorem

1

,

the

equilib-rium

point

E

0

is

a

saddle

point

so

unstable.



Theorem

3.

Consider

the

coexistence

equilibrium

point

E

1

of

the

model

(5)

.

Under

the

positivity

condition

σ

<

γ δ

,

E

1

is

locally

asymp-totically

stable.

Moreover,

if

σ=βγ2₍δ(β(γ+δ)−ω)

−βγ+ω)2 ,

E

1

is asymptotically

sta-ble

under

the

condition

γ<_β2ω_δ2 +δω.

Proof.

The

Jacobian

matrix

of

the

model

(5)

evaluated

at

equilib-rium

point

E

1

is

given

by

J

(

E

1

)

=



−

βγ δ+₂√_βγ−γ ω 1 2

(

−

βγ δ

+

√



+

γ ω

)

√ −γ (βδ+ω) 2βγ ω √ −γ (βδ+ω) 2ω



.

Then,

under

the

positivity

condition

σ

<

γ δ

of

the

coexistence

equilibrium

point,

the

determinant

and

the

trace

of

J

(

E

1

)

are

det

(

J

(

E

1

))

=

√

δ

(

βγ δ

−

√



+

γ ω

)

2

βγ ω

>

0

tr

(

J

(

E

1

))

=

12



−

βγ δ+_βγ√−γ ω

+

√−γ (βδ_ω +ω)

≤ 0

.

Thus,

the

eigenvalues

of

J

(

E

2

)

are

written

as

λ

1

=

12



tr

(

J

(

E

1

))

+

tr

2

(

_J

(

_E

1

))

− 4

det

(

J

(

E

1

))

λ

2

=

1₂



tr

(

J

(

E

1

))

−

tr

2

(

_J

(

_E

₁

))

_{− 4}

_det

(

_J

(

_E

₁

))

If

tr

2

₍

_J

₍

_E

1

))

− 4

det

(

J

(

E

1

))

>

0

,

then

the

eigenvalues

becomes

negative

real

numbers;

if

tr

2

₍

_J

₍

_E

1

))

− 4

det

(

J

(

E

1

))

<

0

,

then

we

ob-tain

a

pair

of

complex

conjugate

eigenvalues

λ1

and

λ

2

=

λ

1

.

Since

tr

(

J

(

E

1

))

<

0,

we

have

Re

(

λ1

)

=

Re

(

λ2

)

<

0

and

consequently

we

have

|

arg

(

λ1

_,2

)

|

>

απ2

.

If

σ=βγ2₍δ(β(γ_−βγ++_ω)δ)2−ω),

then

tr

(

J

(

E

1

))

=

0

.

So,

we

obtain

a

pair

of

complex

conjugate

eigenvalues

λ1

and

λ2

=

λ1

.

Since

Re

(

λ1

)

=

Re

(

λ2

)

=

tr

(

J

(

E

1

))

=

0

,

we

have

arg

(

λ1

)

=

π2

and

arg

(

λ2

)

=

−

π2

Fig. 1. Stable dynamical behaviour of the model (5) for the parameter values given in Table 1 with

σ

= 0 . 1181 in (a) and

σ

= 0 . 5 in (b) with initial condition (E, T ) = (3 , 10)

where blue and red curves represent population density of ECs and TCs respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

2.2.

Numerical

simulations

In

this

section,

the

predictor-corrector

method

is

used

for

nu-merical

simulations

of

the

model

(5)

.

This

method

introduced

in

[36,37]

which

is

a

combination

of

some

product

integration

rules,

known

as

fractional

Adams-Bashforthm-Moulton

methods

[38]

.

To

analyze

the

effects

of

the

model

parameters

on

its

dynamics,

it

is

easier

to

make

qualitative

analysis

on

the

dimensionless

form

of

the

model

(3)

.

That

is

why

we

obtained

non-dimensionalized

model

on

which

qualitative

analysis

is

performed.

In

Fig.

1

,

we

get

asymptotically

stable

coexistence

equilibrium

points

for

different

values

of

σ

.

In

Fig.

1

a

for

smaller

values

of

σ

,

we

observed

an

oscillatory

behaviour

for

both

ECs

and

TCs.

Ac-cording

to

Theorem

1

,

fractional

derivatives

enlarges

the

regions

of

stability.

Here,

we

also

observed

that

smaller

fractional

order

derivatives

damps

the

oscillation

behavior

and

for

smaller

frac-tional

derivatives,

both

ECs

and

TCs

approaches

quicker

to

the

equilibrium

point.

In

Fig.

1

b

for

σ

=

0

.

5

,

ECs

are

more

success-ful

but

insufficient

in

eradicating

the

TCs

and

oscillatory

behaviour

disappears

comparing

with

Fig.

1

a.

Both

situations

corresponds

to

the

dormant

tumor

state

[24,25]

.

3.

Dynamical

behavior

of

conformable

fractional

order

tumor

model

3.1.

Discretizations

process

In

this

section,

we

will

discretize

the

model

(6)

by

using

piece-wise

constant

approximation

[17]

.

Consider

the

conformable

frac-tional

order

model

(6)

as



T

_α

E

(

t

)

=

σ

+

ω

E

(

t

)

T

(

[

t h

]

h

)

−

δ

E

(

t

)

,

T

_α

T

(

t

)

=

γ

T

(

t

)(

1

−

β

T

(

t

))

− E

(

[

t h

]

h

)

T

(

t

)

,

(8)

with

E

(

0

)

=

E

0

and

T

(

0

)

=

T

0

,

where

[

t

]

denotes

the

integer

part

of

t

∈

[0,

∞

)

and

h

>

0

is

discretization

parameter.

Appliying

the

property

(1)

of

conformable

fractional

derivative

to

the

first

equation

of

the

system

(8)

for

t

∈

[

nh

,

(

n

+

1

)

h

)

gives

(

t

− nh

)

1−α

dE

(

t

)

dt

=

σ

+

ω

E

(

t

)

T

(

nh

)

−

δ

E

(

t

)

.

By

simplifying

this

equation,

we

get

E



(

t

)

+

E

(

t

)



_δ

₋

_ω

_T

₍

_nh

₎

(

t

− nh

)

1−α



=

₍

σ

t

− nh

)

1−α

.

Clearly,

this

is

a

first-order

linear

ordinary

differential

equation.

Solving

this

equation

with

respect

to

t

∈

[

nh,

t

),

we

obtain

E

(

t

)

=

(

δ

−

ω

T

(

nh

))

E

(

nh

)

+

σ

e

(δ−ωT(nh))(t−nh_α)α

_{− 1}



e

(δ−ωT(nh))(t−nh)α

α

(

δ

−

ω

T

(

nh

))

and

by

taking

t

_→

(

n

₊

1

)

h

_,

we

get

the

following

difference

equa-tion

E

((

n

+

1

)

h

)

=

σ

+

[

(

δ

−

ω

T

(

nh

))

E

(

nh

)

−

σ

]

e

(ωT(nh)−δ)

hα

α

δ

−

ω

T

(

nh

)

.

Finally,

adjusting

difference

equation

notation

and

replacing

E

(

nh

)

and

T

(

nh

)

by

E

(

n

)

and

T

(

n

)

yields

E

(

n

+

1

)

=

σ

+

[

(

δ

−

ω

T

(

n

))

E

(

n

)

−

σ

]

e

(ωT(n)−δ)

hα

α

δ

−

ω

T

(

n

)

.

In a

similar fashion,

discretizing the second equation

of the system

(8)

T

_α

T

(

t

)

=

γ

T

(

t

)(

1

−

β

T

(

t

))

− E



t

h



h

T

(

t

)

leads

to

the

following

difference

equation

T

(

n

+

1

)

=

T

(

n

)(

γ

− E

(

n

))

(

γ

− E

(

n

)

−

γ β

T

(

n

))

e

(E(n)−γ )hα

α

+

γ β

T

(

n

)

.

Therefore,

we

get

the

two-dimensional

discrete

system

⎧

⎨

⎩

E

(

n

+

1

)

=

σ+[(δ−ωT(n))E(n)−σ]e(ωT(n)−δ)hα α δ−ωT(n)

T

(

n

+

1

)

=

T(n)(γ−E(n)) (γ−E(n)−γ βT(n))e(E(n)−γ )hα α +γ βT(n)

.

(9)

3.2.

Stability

analysis

Now,

we

analyze

local

asymptotic

stability

of

the

system

(9)

.

We

note

that

system

(9)

and

system

(5)

have

the

same

equilibrium

points

that

is

E

0

,

E

1

.

We

linearize

the

system

(9)

about

the

equilibrium

point

E

0

.

Fig. 2. Stable dynamical behaviours of the model (9) for the parameter values given in Table 1 with

σ

= 0 . 1181 in (a),

σ

= 0 . 5 in (b) and initial condition (E, T ) = (3 , 10)

where blue and red curves represents population density of ECs and TCs respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

λ2

=

e

−hα (σαδ−γ δ)

.

It

is

easy

to

prove

that

E

0

is

locally

asymptotically

stable

if

γ δ

<

σ

,

and

unstable

if

γ δ

>

σ

.

The

equilibrium

point

E

1

has

nonnegative

coordinates

if

the

condition

σ

_<

δγ

is

satisfied.

Now,

the

Jacobian

matrix

obtained

by

linearizing

the

system

(9)

about

the

equilibrium

point

E

1

=

(

E

,

T

)

is

given

by

A

_(E,T)

=



a

11

a

12

a

21

a

22



=

⎛

⎝

e

−2hαβγαc 4β2γ2σ ω(1−e − hαc 2αβγ₎ c2 −1+e− h2αωαd βγ

e

− hαd 2αω

⎞

⎠

where

c

=

βγ δ

+

√



−

γ ω

and

d

=

βγ δ

−

√



+

γ ω

.

The

corresponding

characteristic

polynomial

is

λ

2

₊

_p

1

λ

+

p

0

where

p

0

=

e

− hαc 2αβγ

_e

−h2αωαd

+

4

1−e− h2αβγαc



_1−e− h2αωαd



βγ σ ω c2

,

p

1

=

−

e

−2hαβγαc

₊

_e

−2hααωd



.

Theorem

4.

Assume

that

σ

<

δγ

.

The

coexistence

equilibrium

point

(

E

,

T

)

is

locally

asymptotically

stable

if

and

only

if

σ

>

σ

=

e hα (−√−γ (βδ−ω)) 2αβγ _−ehα (βγ δ− √₊ γ ω) 2αω



c2 4

−1+e hα (−√−γ (βδ−ω)) 2αβγ



₋₁₊ehα (βγ δ− √₊ γ ω) 2αω



βγ ω

.

(10)

Proof.

To

check

asymptotically

stability

of

E

1

,

we

use

the

following

Jury

conditions

[31]

:

i

)

1

+

p

1

+

p

0

>

0

,

ii

)

1

− p

1

+

p

0

>

0

,

iii

)

1

− p

0

>

0

.

Under

assumption

σ

_<

δγ

,

we

always

have

c,

d

_>

0.

These

two

in-equality

assures

that

0

<

e

−2hαβγαc

_,

_e

−h2αωαd

<

1

(11)

and

consequently

we

have

1

+

p

1

+

p

0

=

1

− e

−hαc 2αβγ



₁

_{− e}

−h2αωαd



+

4

1

− e

−2hαβωαc



₁

_{− e}

−h2ααωd



βγ σω

c

2

>

0

,

1

− p

1

+

p

0

=

1

+

e

−2hαβγαc



₁

₊

_e

−h2αωαd



+

4

1

− e

−2hαβωαc



1

− e

−h2ααωd



βγ σω

c

2

>

0

.

In

addition,

under

the

condition

(10)

,

we

have

p

0

=

e

− hαc 2αβγ

_e

−h2ααωd

+

4

1−e− h2αβγαc



_1−e− h2αωαd



βγ σ ω c2

<

1

,

(12)

essentially

1

− p

0

>

0

and

the

proof

is

completed.

Fig.

2

shows

asymptotically

stable

coexistence

equilibrium

points

for

different

values

of

σ

similarly

to

Fig.

1

.

In

Fig.

2

a

for

smaller

values

of

σ

,

we

observed

more

oscillatory

behaviour

for

both

ECs

and

TCs

comparing

with

Fig.

1

a.

This

is

for

the

reason

that,

σ

₌

0

_.

1181

is

very

close

to

the

critical

value

σ

₌

0

_.

075445

in

Theorem

4

which

will

be

the

critical

threshold

for

Neimark-Sacker

bifurcation

as

we

will

see

in

next

section.

In

Fig.

1

b

for

σ

=

0

_.

5

_,

the

system

loses

its

oscillatory

behavior

and

the

tumor

cells

extinct

for

a

while,

then

approaches

to

the

coexistence

equi-librium

point.

This

situation

corresponds

to

the

state

of

a

“return-ing

tumor”[35]

.

Moreover,

in

Fig.

1

b

we

observe

that

for

a

smaller

fractional

derivatives,

both

ECs

and

TCs

approaches

quicker

to

the

equilibrium

point.

3.3.

Neimark-Sacker

bifurcation

analysis

In

this

section,

we

analyze

the

Neimark-Sacker

bifurcation

of

the

system

(9)

at

the

equilibrium

point

E

1

.

First

of

all,

we

convert

the

equilibrium

point

E

1

=

(

E

,

T

)

of

the

system

(9)

into

the

origin

by

change

of

variables

x

1

=

E

− E

and

x

2

=

N

− N

.

Thus,

the

system

(9)

turns

into



x

1

x

2



→

A

(

σ

)



x

1

x

2



+



F

1

(

x

1

,

x

2

,

σ

)

F

2

(

x

1

,

x

2

,

σ

)



(13)

where

A

(

σ

)

=

A

₍

_E,T

₎

and

F

1

(

x

1

,

x

2

,

σ

)

=

h

α

a

11

ω

α

x

1

x

2

+

h

2α

_a

11

ω

2

2

α

2

x

1

x

2 2

−

8

a

11

β

2

γ

2

σ ω

2

((

2

− 2

/a

11

)

αβγ

+

h

α

c

2

α

c

3

x

2 2

−

12

a

11

β

2

γ

2

σω

3

6

α

2

_c

4

×

((

8

− 8

/a

11

)

α

2

β

2

γ

2

+

4

h

α

αβγ

c

+

h

2α

c

2

)

x

32

+

O

(

|

X

|

4

)

,

F

2

(

x

1

,

x

2

,

σ

)

=

2

a

22

((

−2

+

2

a

22

)

αω

+

h

α

d

)

2

αβγ

d

x

2 1

+

4

(

−a

22

+

a

222

)

βγ ω

2

d

x

2 2

+

a

22

((

−4

+

4

a

22

)

αω

+

h

α

d

)

α

d

x

1

x

2

+

3

a

22

6

α

2

βγ

_d

2

(

−8

a

2 22

(

−1

+

1

/a

22

)

α

2

ω

2

+

4

(

−1

+

2

a

22

)

h

α

αω

d

+

h

2α

d

2

)

x

31

+

a

22

2

α

2

_d

2

(

8

α

2

_ω

2

+

24

a

2 22

α

2

ω

2

− 8

h

α

αω

d

+

h

2α

d

2

− 16

a

22

αω

(

2

αω

− h

α

d

))

x

21

x

2

+

4

a

3 22

βγ ω

2

α

d

2

(

6

αω

+

1

/a

2 22

(

4

αω

− h

α

d

)

−

(

2

/a

22

)(

5

αω

− h

α

d

))

x

1

x

22

+

24

a

322

(

−1

+

1

/a

22

)

2

β

2

γ

2

ω

2

6

d

2

x

3 2

+

O

(

|

X

|

4

)

.

For

σ

=

σ

,

the

eigenvalues

of

A

(

σ

)

are

λ

1,2

(

σ

)

=

−p

1

(

σ

)

± i

p

1

(

σ

)

2

− 4

p

0

(

σ

)

2

where

p

1

(

σ

)

=

−e

− hα (βγ δ+√4γ βσ ω+γ2_(βδ−ω)2_{−γ ω)} 2αβγ

_{− e}

hα ( √ 4γ βσ ω+γ2_(βδ−ω)2_{−γ (βδ+ω))} 2αω

and

p

0

(

σ

)

=

1

.

Moreover,

we

obtain

|

λ1

_,2

(

σ

)

|

=

1

.

So,

the

eigenvalues

of

A

(

σ

)

are

complex

conjugates

with

modulus

1

as

required.

In

addi-tion,

p

1

(

σ

)

=

0

,

1

and

we

can

conclude

that

λ

k1,2

(

σ

)

=

1

for

k

=

1

,

2

,

3

,

4

.

Hence,

non-strong

resonance

condition

is

also

satisfied.

The

transversality

condition

yields

to

the

inequality

d|λ1,2(σ )| dσ





σ=σ

=

2βγ ω

e hα (βγδ+√) 2αβγ _−eh2α ωαβ



ehα √ 2αω −ehα γ (βδ+2αωω)



e hα (β2_γ2_{δ+γ ω(γ+δ)+}√_ω) 2αβγ ω _{(βγ δ+}√_{−γ ω)}2

=

0

which

is

satisfied

for

positive

parameters

and

σ

<

βγ

.

Now,

we

calculate

multilinear

functions:

B

1

(

x,

y

)

=

2



j,k=1

∂

2

_F

1

(

ψ

,

σ

)

∂

ψ

j

∂

ψ

k





ψ=0

x

j

y

k

=

h

α

a

11

ω

α

(

x

1

y

2

+

x

2

y

1

)

−

(

8

a

11

β

2

γ

2

σω

2

(

2

− 2

/a

11

)

αβγ

+

h

α

c

)

α

c

3

x

2

y

2

,

B

2

(

x,

y

)

=

2



j,k=1

∂

2

_F

2

(

ψ

,

σ

)

∂

ψ

j

∂

ψ

k





ψ=0

x

j

y

k

=

2

a

22

((

−2

+

2

a

22

)

αω

+

h

α

d

)

αβγ

d

x

1

y

1

+

4

(

−a

22

+

a

222

)

βγ ω

d

x

2

y

2

+

a

22

((

−4

+

4

a

22

)

αω

+

h

α

d

)

α

d

(

x

1

y

2

+

x

2

y

1

)

,

C

1

(

x

,

y

,

u

)

=

2



j,k,l=1

∂

3

_F

1

(

ψ

,

σ

)

∂

ψ

j

∂

ψ

k

∂

ψ

l





ψ=0

x

j

y

k

u

l

=

h

2α

a

11

ω

2

α

2

(

x

2

y

2

u

1

+

x

2

y

1

u

2

+

x

1

y

2

u

2

)

−

12

a

11

β

2

γ

2

σω

3

α

2

_c

4

x

2

y

2

u

2

((

8

− 8

a

22

)

α

2

_β

2

_γ

2

+4

h

α

αβγ

c

+

h

2α

_c

2

₎

_,

C

2

(

x,

y,

u

)

=

2



j,k,l=1

∂

3

_F

2

(

ψ

,

σ

)

∂

ψ

j

∂

ψ

k

∂

ψ

l





ψ=0

x

j

y

k

u

l

=

3

a

22

α

2

βγ

_d

2

(

8

a

22

(

−1

+

a

22

)

α

2

_ω

2

+

4

(

−1

+

2

a

22

)

h

α

αω

d

+

h

2α

d

2

)

x

1

y

1

u

1

+

a

22

α

2

_d

2

(

8

α

2

_ω

2

₊

₂₄

_a

2 22

α

2

ω

2

− 8

h

α

αω

d

+

h

2α

d

2

− 16

a

22

αω

(

2

αω

− h

α

d

))(

x

2

y

1

u

1

+

x

1

y

2

u

1

+

x

1

y

1

u

2

)

+

4

a

22

βγ ω

α

d

2

(

2

(

2

− 5

a

22

+

3

a

2 22

)

αω

+

(

−1

+

2

a

22

)

h

α

d

)(

x

2

y

2

u

1

+

x

2

y

1

u

2

+

x

1

y

2

u

2

)

+

24

a

22

(

−1

+

a

22

)

2

β

2

γ

2

ω

2

d

2

x

2

y

2

u

2

.

Now,

let

q

∼

(

ζ1

+

i

ζ2

,

1

)

∈

C

2

_be

_an

_eigenvector

_of

_A

₍

_σ

₎

corre-sponding

to

λ1

(

σ

)

and

let

¯p

∼

(

ξ1

+

i

ξ2

,

1

)

∈

C

2

_be

_an

_eigenvector

of

A

T

₍

_σ

₎

_{corresponding}

_to

_λ2

₍

_σ

₎

_.

_Afterwards,

_we

_normalize

_¯p

_with

respect

to

q

.

Hence,

q

and

the

normalized

vector

p

are

computed:

p

=



ξ1+iξ2 1+(ξ1+iξ2)(ζ1−iζ2)

,

1 1+(ξ1+iξ2)(ζ1−iζ2)

,

q

=

(

ζ

1

+

i

ζ

2

,

1

)

where

ξ

1

=

−1+e2hααωd



−1+e hα (βγ (βγ δ−√)+(βγ (γ−δ)−√)ω+γ ω2 2αβγ ω



2βγ

ehαωαd−ehα(βγ (βγ δ− √ )+(βγ (γ−δ)−√)ω+γ ω2 2αβγ ω



ξ

2

=

−1+eh2ααωd



eh2ααωd 2βγ

ehαωαd−e hα (βγ (βγ δ−√)+(βγ (γ−δ)−√)ω+γ ω2 2αβγ ω



×



4

−

(

e

−2hαβγαc

)

2

_{− 2}

_e

−2hαβγαc

_e

−2hααωd

−

(

e

−h2αωαd

)

2

_,

ζ

1

=

βγe− h αc 2αβγ

e hαc 2αβγ_−ehαd 2αω



2

−1+e2hααωd



,

ζ

2

=

βγ  4−(e− h αc 2αβγ)2_−2e− h αc 2αβγ_e− h2ααωd−(e− h2ααωd)2 2−2e− h2αωαd

.

Moreover,

for

σ

sufficently

close

to

σ

_,

we

can

express

any

vec-tor

V

∈

R

2

_as

_V

₌

_zq

₊

_¯z

_¯q

_,

_where

_z

_is

_a

_complex

_number.

Accord-ingly,

for

sufficently

small

|

σ

|

(near

σ

),

the

system

(13)

can

be

ex-pressed

in

the

following

form:

z

→

λ

1

z

+

g

(

z,

z

,

σ

)

,

where

λ1

(

σ

)

=

(

1

+

ψ

(

σ

))

e

iθ(σ )

_with

_ψ

₍

_σ

₎

_is

_a

_smooth

func-tion

satisfying

ψ

(

q

)

=

0

,

g

is

a

complex-valued

smooth

function

of

z

,

¯z

_,

σ

_,

whose

Taylor

expression

with

respect

to

(

z

,

¯z

)

contains

quadratic

and

higher

order

terms

g

(

z,

¯z

,

σ

)

=



k+l≥2

1

k

!

l!

g

k,l

(

σ

)

z

k

_¯z

l

_,

_g

kl

∈

C

,

k,

l

=

0

,

1

,

2

,

.

.

.

By

using

multilinear

functions,

we

can

determine

the

Taylor

co-efficients

g

_kl

through

the

following

formulas:

g

20

(

σ

)

=

<

p,

B

(

q,

q

)

>,

g

11

(

σ

)

=

<

p,

B

(

q,

q

)

>,