key: cord-0747238-ybc0it78
authors: Ghanbari, Behzad
title: On the modeling of the interaction between tumor growth and the immune system using some new fractional and fractional-fractal operators
date: 2020-10-19
journal: Adv Differ Equ
DOI: 10.1186/s13662-020-03040-x
sha: af4be4ffb4bc0ec55ee860e2d082e235545ce550
doc_id: 747238
cord_uid: ybc0it78

Humans are always exposed to the threat of infectious diseases. It has been proven that there is a direct link between the strength or weakness of the immune system and the spread of infectious diseases such as tuberculosis, hepatitis, AIDS, and Covid-19 as soon as the immune system has no the power to fight infections and infectious diseases. Moreover, it has been proven that mathematical modeling is a great tool to accurately describe complex biological phenomena. In the recent literature, we can easily find that these effective tools provide important contributions to our understanding and analysis of such problems such as tumor growth. This is indeed one of the main reasons for the need to study computational models of how the immune system interacts with other factors involved. To this end, in this paper, we present some new approximate solutions to a computational formulation that models the interaction between tumor growth and the immune system with several fractional and fractal operators. The operators used in this model are the Liouville–Caputo, Caputo–Fabrizio, and Atangana–Baleanu–Caputo in both fractional and fractal-fractional senses. The existence and uniqueness of the solution in each of these cases is also verified. To complete our analysis, we include numerous numerical simulations to show the behavior of tumors. These diagrams help us explain mathematical results and better describe related biological concepts. In many cases the approximate results obtained have a chaotic structure, which justifies the complexity of unpredictable and uncontrollable behavior of cancerous tumors. As a result, the newly implemented operators certainly open new research windows in further computational models arising in the modeling of different diseases. It is confirmed that similar problems in the field can be also be modeled by the approaches employed in this paper.

The immune system is a real masterpiece that does extraordinary things every day while we do not realize such numerous activities. The job of the immune system is protecting our bodies against invading factors such as bacteria and viruses. Many factors such as low nutrient intake, lack of sleep, and high stress weaken a person's immune system. On the other hand, people with cancer are also at risk for infectious diseases since they usually undergo special chemotherapy treatments using special drugs that weaken their immune system. The main purpose of these special drugs is destroying diseased cells without damaging adjacent tissues. Unfortunately, some healthy cells and tissues in the body may also be affected during this treatment.

Nowadays, with the pandemic outbreak of the dangerous infectious disease COVID-19 all over the world, infectious disease specialists and epidemiologists constantly emphasize the strengthening of the immune system as one of the most important ways to prevent this deadly disease.

In this paper, we provide some novel approximate solutions to a computational model that formulates the interaction between tumor growth and the immune system, including several fractional and fractal operators. The model consists of three state variables, each representing the number of specific cells in the problem. The interaction between these variables is presented by Itik and Banks [26] through following nonlinear differential equation system:

subject to initial conditions (T(0), H(0), E(0)) = (T 0 , H 0 , E 0 ) ≥ 0. In this model, T(τ ) is used to count the number of tumor cells at time τ , H(τ ) represents the number of healthy host cells, and E(τ ) refers to the number of effector immune cells in the single tumor-site compartment. Moreover, τ represents the rate of change in the population of the tumor cells, k 1 describes the growth of the tumor cells, and s 1 is their maximum carrying capacity. The first term of Eq. (1) expresses to the logistic growth of the tumor cells in the absence of any effect from other cell populations. The competition between the host cells H(τ ) and the tumor cells T(τ ), which results in the loss of the tumor cell population, is given by the term β 12 T(τ )H(τ ). Moreover, β 12 refers to the tumor cell killing rate by the effector cells E(τ ). In the second equation of system (1) the healthy tissue cells also grow logistically with the growth rate k 2 and maximum carrying capacity s 2 . We assume that the cancer cells proliferate faster than the healthy cells, and thus k 1 > k 2 . The tumor cells also inactivate the healthy cells at the rate β 21 . Also, the third equation of system (1) illustrates the stimulation of the immune system by the tumor cells with tumor specific antigens. One of the other assumptions considered in the model is that the immune system depends directly on the number of tumor cells a the rate of positive constants k 3 and s 3 . Finally, β 31 is the rate of change corresponding to the inactivated effector cells by the tumor cells, and c 3 is their rate of natural death. All parameters in the model are positive constants.

For simplicity of analysis, we first nondimensionalize system (1) by employing the definitions

In addition, we use the following new parameters:

By involving the introduced changes in equations (2) and (3), we obtain a new dimensionless form for the problem [26] :

Due to the great importance of model (4), many researchers have shown interest in examining various technical and computational aspects of this model. The authors explain the biological relevance of model (7) . In [41] the authors have utilized the localization of compact invariant sets (LCIS) method and Lyapunov stability theory to investigate the global dynamics of model (4) . In [28] the authors have focused on the mathematical points of view, such as the existence of Hopf bifurcations corresponding to the model. In [24] the authors considered the model using the Caputo-Fabrizio-Caputo and new fractional derivative with Mittag-Leffler kernel in the Liouville-Caputo sense. Starko and Coria [40] provided sufficient conditions on model parameters and treatment parameters under which all trajectories in the positive orthant tend to the tumor-free equilibrium point. The authors in [32] have developed an efficient methodology of partial control and applied it to avoid the extinction of the healthy tissue. In [42] , taking impulsive differential equations into account, the authors have presented a mathematical formulation of tumor-immune interaction with periodically pulsed immunotherapy. In [3] , model (4) was modified to include three delay parameters in the problem. Fractional calculus has a relatively long history almost as long as a integer-order differential account. However, in recent decades, the implementations of these concepts was neglected compared to standard calculus. This trend seems to have changed in the past years in general. A defining sign of this change is the increasing use of these tools in the literature. Thanks to the researchers' efforts, many differentiated and integral operators based on different approaches were proposed and then successfully implemented in the past few years [1, 4, 5, 10, 13, 23, 25, 30] . From the perspective of numerical aspects, a wide range of new mathematical methods was successfully applied in various branches of science [6, 13, 16-19, 27, 29, 35, 37, 39] . For example, in [9] the authors have developed an efficient numerical treatment for ordinary fractional and fractal-fractional differential differential equations, which is based on Newton polynomials. In [8] , Atangana and his  collaborator have successfully applied a new numerical algorithm to approximate a modified version of the Chua attractor model with both fractional and fractal-fractional oper- ators. An efficient numerical technique, the Atangana-Seda numerical scheme, based on Newton polynomials, is utilized in [7] to handle a chaotic problem with fractional operators, which include the exponential decay, power law, and Mittag-Leffler kernel. In [23] the author has employed some novel differential and integral operators of fractional order and fractal dimension using he uville-Caputo and Atangana-Baleanu definitions to obtain multiple attractors and periodicity on the Vallis model for El Niño/La Niña-Southern oscillation model. A new fractional-order compartmental SEIRS model with Caputo-type fractional-order derivative was also studied in [25] . Chaotic systems are almost one of the most important and applicable types of nonlinear equations [12, 21, 34, 36, 38] . Therefore in many cases the exact solution is not available for such equations. On the other hand, the use of new derivative operators in structures of chaotic systems has made significant development in this field [2, 22] . In some cases the researchers have obtained desirable attractors, which were not achievable by common integer-order operators. This fact highlights the importance of new derivative operators in other real-world models. Motivated by these achievements, especially following the work [24] , we intend to investigate the model presented in equation (4) using some new efficient fractional and fractionalfractional operators.

To reach this goal, the subsequent parts of the paper are structured as follows. The analysis of model equilibrium points is presented in Sect. 2. This model is examined in the next section via the Liouville-Caputo fractional-order derivative. This section also confirms that under appropriate assumptions, the model always possesses a unique solution. Then a numerical method corresponding to this structure is designed and then utilized. Besides, detailed numerical simulations are presented. Similar processes will be followed in Sects. 4 and 5 of the paper, with the Caputo-Fabrizio-Caputo and Atangana-Baleanu-Caputo fractional derivative operators, respectively. In Sect. 6, we examine the model via several fractal-fractional operators. This section also presents the numerical methods corresponding to each of these operators. To investigate the dynamic behavior of the results, we added several numerical simulations. Finally, in Sect. 7, we present a summary of the results and achievements of the paper.

In this section, we analyze the equilibrium points of the considered model (4). These points are in fact the roots of the nonlinear algebraic system that is contracted on the right-hand side of the model. By solving this system, we determine six possible equilibrium points for the model [28] . 13 ). The first coordinate T * is determined by finding the nonnegative root(s) of the characteristic equation

This equation has the acceptable root

Necessary conditions for the existence of this equilibrium point are

The singular point P 5 = ( k2(β 12 -1) β12β21-k2 , β21-k2 β 12 β 21 -k2 , 0), which implies the coexistence of cancer and host cells. The necessary conditions for the existence of this equilibrium point are

It should also be noted that for β 12 = 1, this equilibrium point becomes the equilibrium point P 2 . Point 6: The interior fixed point

where T * is a positive root of β 13 T * 2 + T * (c 3 + s 3 β 31k 3 ) + s 3 c 3 = 0, as mentioned earlier. In this case, all three cell populations are present in the problem. The necessary conditions for the existence of this equilibrium point are

Moreover, the Jacobi matrix corresponding to this system is

In this section, we consider model (4) with the Liouville-Caputo (LC) fractional derivative,

where the LC fractional derivative is defined as [14] LC

Utilizing the Laplace transform on the LC derivative (8), we get

Taking (9) into account and then utilizing the inverse Laplace transform on Eq. (7), we get

From (10) we suggest the following iterative schemes:

where

The desired approximate solutions can be obtained by computing the limits

Let us assume that a Banach space like has a closed convex bounded subset that contains a fixed point of . Moreover, let ω : → be a condensing map. Moreover, let us assume that there exists

Then we have:

1. X * , Y * , and Z * are Lipschitz and bounded.

2. X * * , Y * * , and Z * * are compact and bounded.

With the help of the Riemann-Liouville integral, Eq. (7) is written as

Theorem 1 As soon as

then under two properties 1 and 2, the existence of a solution to the problem is guaranteed.

Hence we obtain that X, Y, Z, and W are condensing, and the fixed points of X, Y, Z, and W are verified.

Similarly, we have

and therefore

II) Let us show that X * , Y * , and Z * are contractions. For m, z ∈ χ , we conclude

These statements confirm that X * , Y * , and Z * are contractions. III) Let us show that X * * , Y * * , and Z * * are compact. For 0 ≤ j 1 ≤ j 2 ≤ T, we have

Similarly, we get

Finally, we have

Now by the Arzelà-Ascoli theorem [15] X * * ( ξ ), Y * * , and Z * * are relatively compact. So X * * , Y * * , and Z * * are compact.

Since X * , Y * , Z * , and W * are contractions and X * * , Y * * and Z * * are compact and therefore continuous, the maps X = X * + X * * , Y = Y * + Y * * , and Z = Z * + Z * * are condensing on ξ . Hence the existence of a fixed point for each point X, Y, and Z is proved.

IV) We will show that the problem has a unique solution. To this end, let us define the map H as follows:

These results indicate that model (7) will always has a unique solution.

In this section, we use the Adams-Bashforth-Moulton (ABM) numerical method. Here we follow the steps to apply this method in solving the following fractional-order problem:

Now taking the Liouville-Caputo fractional integration on (24) yields

The following predictor-corrector form determines an approximate solution to the problem [31] :

where

Utilizing the numerical algorithm presented in (26), we determine an approximate solution to the fractional problem (7) from the formulae 

In this section, we study the following model:

The fractional derivative operator CFC 0 D α t in this model is Caputo-Fabrizio-Caputo (CFC) given by [13] CFC

where

The CFC fractional integral is also defined by [33] CFC

Applying the CFC integral definition (29), we obtain the following relationships:

Now we consider the following kernels:

Theorem 2 The initial value problem and the kernels K 1 (T (t), H(t), E(t)), K 2 (T (t), H(t), E(t)), and K 3 (T (t), H(t), E(t)) satisfy the Lipschitz condition.

Proof See [24] .

Theorem 3 The fractional nonlinear system (29) admits at least one solution.

Proof See [24] .

Proof See [24] .

Now let us focus on determining an approximate solution to the following CFC fractional Cauchy problem:

Using the corresponding fractional integral operator yields

Taking t = t n+1 in (36), we have

and

Inserting Eq. (38) into Eq. (37), we get

where

So we have

As a result, the following recursive relations are determined to approximate the CFC problem (29) as in [22] :

where 1 

Now let us consider the model via the Atangana-Baleanu-Caputo fractional derivative as

where the Atangana-Baleanu fractional integral of order α of a function φ(t) is defined as

where B(α) = 1α + α (α) is a normalization function. The Atangana-Baleanu fractional integral of order α of a function φ(t) is also defined as [10] ABC

Taking the definition of the Atangana-Baleanu fractional integral (45) on both sides of equations in system (43), we get the Volterra integral system

and the following iterative formulas:

As n → ∞, Eq. (47) suggests the exact solution for the model. Proof First, we define

where E 1 (t, T (t)), E 2 (t, H(t)), and E 3 (t, E(t)) are contractions respect to T (t), H(t), and E(t), respectively. Moreover, we set

where

Considering the Picard operator, we have : ( a,b 1 , a,b 2 , a,b 3 ) → ( a,b 1 , a,b 2 , a,b 3 ) (52) defined as follows

where (t) = (T (t), H(t), E(t)), 0 (t) = (T (0), H(0), E(0)), and (t, (t)) = E 1 (t, T (t)), E 2 (t, H(t)), E 3 (t, E(t)).

Now we assume that all the solutions are bounded within a period of time:

provided that

The fixed point theorem of a Banach space together with the metric suggests that

Since is a contraction along with η < 1, we must have ωη < 1. Therefore we conclude that is a contraction operator. Finally, it is proved that (43) always possesses a unique solution.

Consider the following fractional initial value problem:

Employing the product integration rule, Ghanbari and Kumar [20] have developed an efficient scheme to obtain the approximate solution of (56), given by

where p n = (n -1) α+1n α (nα -1) (α + 2) ,

, j = 1, 2, . . . , n -1.

(58)

Using this numerical approximation, we get the following iterative scheme: Thus we obtain

-K 1 (t ζ -1 ) (-ζ + 1 + n) α+1 -(-ζ + n) α (-ζ + n + α + 1) ,

where K j , ζ = 1, 2, 3, are given in (63). 

In this part, we replace the derivative in (4) by the fractal-fractional derivative with exponential decay-law:

where the fractal-fractional derivative with exponential decay law of a function φ(t) is defined as [11] FF-E 0

and d dt ψ φ(u) is introduced in (62). Applying the Caputo-Fabrizio integral to Eq. (69), we obtain

Setting t = t n+1 in (71), based on the proposed scheme in [11] , we get

Taking the difference between T n+1 and T n yields

Therefore the approximate solution to the problem can be determined using the following iterative procedures:

where F j , ζ = 1, 2, 3, are given in (63). 

In this subsection, we use the fractal-fractional derivative with Mittag-Leffler law in (4).

So, we achieve

where the fractal-fractional derivative with Mittag-Leffler law of a unction φ(t) is defined as [11] 0 FF-ABC D α,τ t φ(t)

and d dt ψ φ(u) is introduced in (62). Applying the Atangana-Baleanu integral to (75), we obtain

Based on the scheme proposed in [11] , we get

Now using the Lagrange polynomial piecewise interpolation given by Eq. (64), we obtain

(80) Fig. 18 . In Figs. 19 and 20, we also investigate the sensitivity of the state variables when two parameters β 12 and c 3 change, respectively. In these two figures, we can see that changes in each of these parameters cause changes that are somewhat meaningful and constructive in the behavior of variables. As a biological conclusion, we can point to the fact that the stability in the model can only be achieved when the recruitment rate of effector cells is greater than the rate of inactivation by cancer cells. In other words, if the immune system is unable to detect and attack cancer cells, then effective treatment must be used to control the tumor growth.

Recent advances in the presentation of efficient numerical methods in fractional differential calculus are a great help to researchers in modeling real phenomena. It has also been proven that differential calculus of integer order in some cases is incapable of describing the behavior of some phenomena or faces fundamental problems. In this paper, we study the dynamics of a tumor-immune model due to the unpredictable growth of tumor cells described by an attractive form of nonlinear differential equations. The importance of this model has led to many different approaches of studying the problem. Our study in the paper makes it clear how tumor cells interact with the immune system. The main dif-ference of this paper from other papers on this system is that we have used new fractional and fractal-fractional definitions for the derivative in the system structure. Some of these new concepts of differentiation were introduced in 2017 by Atangana and his collaborators. They combine the idea of fractal derivative and fractional differentiation, which takes into account the memory, fractal effect, and nonlocality. The model considers processes like power-law, fading memory, and crossover. It is important to note that the numerical methods and analysis used in this paper are quite different from those presented in [24] . In fact, this work can be considered as a complement to the content presented in this paper. The basis of these mathematical algorithms is applying some fundamental axioms of fractional derivative and the first-order interpolation. The existence and uniqueness of the model solutions were also investigated. The interesting attractors obtained in this paper imply that these new fractional and fractal-fractional operators can describe newer aspects of the behavior of these systems than fractional derivatives. Some of these features cannot be described by conventional integer-order operators. Through numerical simulations, we have confirmed the chaotic dynamics of the model by taking certain parameters in the model and suitable initial conditions. These results indicate that the values considered for the fractional-order derivatives significantly influence the dynamic behavior of the problem. The fractal-fractional operators allow us to describe self-similar problems with power-law, fading memory, and crossover behavior. In addition, the chaotic behavior seen in the results is entirely consistent with the inherent nature of the problem. These systems are recognized as complex real-world problems that cannot be represented with classical or fractional differential operators. Numerical simulations confirm that a change in fractional order exhibits memory properties and very strange complex dynamical behaviors. The results of this paper show that other similar problems in this field can be described and investigated using the numerical methods used in this paper.

Fractional variational calculus in terms of Riesz fractional derivatives

Chua's circuit model with Atangana-Baleanu derivative with fractional order

Analytical and numerical study of Hopf bifurcation scenario for a three-dimensional chaotic system

Fractal-fractional differentiation and integration: connecting fractal calculus and fractional calculus to predict complex system

Non validity of index law in fractional calculus: a fractional differential operator with Markovian and non-Markovian properties

Can transfer function and Bode diagram be obtained from Sumudu transform

Atangana-Seda numerical scheme for labyrinth attractor with new differ

New numerical approximation for Chua attractor with fractional and fractal-fractional operators

New numerical method for ordinary differential equations: Newton polynomial

New fractional derivatives with nonlocal and non-singular kernel: theory and application to heat transfer model

Modeling attractors of chaotic dynamical systems with fractal-fractional operators

Chaotic Dynamics: An Introduction

A new definition of fractional derivative without singular kernel

A new dissipation model based on memory mechanism

The Analysis of Fractional Differential Equations: An Application-Oriented Exposition Using Differential Operators of Caputo Type

New approach for the model describing the deathly disease in pregnant women using Mittag-Leffler function

A new application of fractional Atangana-Baleanu derivatives: designing ABC-fractional masks in image processing

Some new edge detecting techniques based on fractional derivatives with non-local and non-singular kernels

On fractional predator and prey models with mutualistic predation including non-local and nonsingular kernels

Numerical solution of predator-prey model with Beddington-DeAngelis functional response and fractional derivatives with Mittag-Leffler kernel

Some effective numerical techniques for chaotic systems involving fractal-fractional derivatives with different laws

Chaos and multiple attractors in a fractal-fractional Shinriki's oscillator model

Multiple attractors and periodicity on the Vallis model for El Niño/La Niña-Southern oscillation model

Chaos in a cancer model via fractional derivatives with exponential decay and Mittag-Leffler law

A new fractional-order compartmental disease model

Chaos in a three-dimensional cancer model

A new and efficient numerical method for the fractional modeling and optimal control of diabetes and tuberculosis co-existence

The influence of time delay in a chaotic cancer model

Dynamical study of fractional order mutualism parasitism food web module

Theory and Applications of Fractional Differential Equations

On the fractional Adams method

Avoiding healthy cells extinction in a cancer model

Properties of a new fractional derivative without singular kernel

Chaotic Dynamics. Cambridge Books

Global dynamics of a fractional order model for the transmission of HIV epidemic with optimal control

Practical Numerical Algorithms for Chaotic Systems

Fractal-fractional differentiation for the modeling and mathematical analysis of nonlinear diarrhea transmission dynamics under the use of real data

Chaotic Dynamics of Nonlinear Systems

Sir epidemic model with Mittag-Leffler fractional derivative

Global dynamics of the Kirschner-Panetta model for the tumor immunotherapy

Bounding the dynamics of a chaotic-cancer mathematical model

Periodically pulsed immunotherapy in a mathematical model of tumor-immune interaction

Not applicable.

This research work is not supported by any funding agencies.

Not applicable.

for T (0) = 0.3517, H(0) = 0.1115, E(0) = 0.4951, and β 12 = 0.920, we get periodic orbit trajectories in the solutions as depicted in Fig. 9 .

In this section, we will use several well-known fractal-fractional derivatives in model (4).

In this subsection, we replace the derivative in (4) by the fractal-fractional derivative with power-law:where the fractal-fractional derivative with power-law of function φ(t) is defined as [11] FF-P 0andBy employing the corresponding inverse operator in Eq. (60) and then placing t = t n+1 in the resultant we obtain the following recursive form:

These functions can be interpolated in [t ζ , t ζ +1 ] as

The author declares that he has no competing interests.

All authors read and approved the final manuscript.

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.Received: 7 September 2020 Accepted: 8 October 2020