On the modeling of the interaction between tumor growth and the immune system using some new fractional and fractional-fractal operators

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.


Introduction
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.

Investigation of stability of equilibrium points of models
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

Point 5:
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

The model via the Liouville-Caputo fractional derivative
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

Existence and uniqueness
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.
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.

The numerical method
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

The model via Caputo-Fabrizio-Caputo fractional derivative
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

Existence of the coupled solutions
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.
Theorem 3 The fractional nonlinear system (29) admits at least one solution.

Numerical method
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 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.

Exponential decay-law fractal-fractional derivative
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).

Mittag-Leffler-law fractal-fractional derivative
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.

Conclusion
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.