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 .

In this model, is used to count the number of tumor cells at time τ, represents the number of healthy host cells, and 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, describes the growth of the tumor cells, and 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 and the tumor cells , which results in the loss of the tumor cell population, is given by the term . Moreover, refers to the tumor cell killing rate by the effector cells . In the second equation of system (1) the healthy tissue cells also grow logistically with the growth rate and maximum carrying capacity . We assume that the cancer cells proliferate faster than the healthy cells, and thus . The tumor cells also inactivate the healthy cells at the rate . 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 and . Finally, is the rate of change corresponding to the inactivated effector cells by the tumor cells, and 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 operators. 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 fractional-fractional 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].

Point 1: the no “living cell” singular point .

Point 2: The tumor-free fixed point .

Point 3: The fixed point , which implies the existence of tumor cells in the model.

Point 4: The fixed point . The first coordinate 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 , 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 , this equilibrium point becomes the equilibrium point .

Point 6: The interior fixed point , where is a positive root of , 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] 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 such that . The exist functions such that , , and . Then we have: With the help of the Riemann–Liouville integral, Eq. (7) is written as

, , and are Lipschitz and bounded.

, , and are compact and bounded.

.

Theorem 1

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

Proof

Let us take χ such that , and let be a closed set in the Banach space equipped with the norm of .

Then, by the definition , , we have Hence we obtain that , , , and are condensing, and the fixed points of , , , and are verified.

I) Let us show that . For , we have Similarly, we have and and therefore .

II) Let us show that , , and are contractions. For , we conclude where These statements confirm that , , and are contractions.

III) Let us show that , , and are compact. For , we have Similarly, we get Finally, we have for .

Now by the Arzelà–Ascoli theorem [15] , , and are relatively compact. So , , and are compact.

Since , , , and are contractions and , and are compact and therefore continuous, the maps , , and are condensing on . Hence the existence of a fixed point for each point , , and is proved.

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

For , , , , , , we obtain Besides, we have and 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

Numerical simulations

Figures 1–3 demonstrate the variation of state variables in model (7) when the scheme (26) is utilized for different values of . In this simulations, we have considered the following values in the model: , , , , , , , and . In our performed numerical simulations, we take and . In Fig. 1, we take . In this case the model exhibits chaotic attractor behavior. Also, by taking into consideration , and the model shows the limit cycle behavior in Fig. 2, whereas for , , , and , we get periodic orbit trajectories in the solutions as depicted in Fig. 3.

Figure 1

Simulations for solving (7) using (26) along with

Figure 3

Simulations for solving (7) using (26) along with and

Figure 2

Simulations for solving (7) using (26) along with and

The model via Caputo–Fabrizio–Caputo fractional derivative

In this section, we study the following model: The fractional derivative operator in this model is Caputo–Fabrizio–Caputo (CFC) given by [13] where The CFC fractional integral is also defined by [33]

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 , , and satisfy the Lipschitz condition.

See [24]. □

Theorem 3

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

See [24]. □

Theorem 4

The fractional nonlinear system (29) always admits a unique solution.

See [24]. □

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 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 , , and .

Figures 4–6 are plotted to demonstrate the variation of state variables in model (29) when the scheme (42) is employed for different values of . In these simulations, we considered the following values in the model: , , , , , , , and . In our performed numerical simulations, we considered and . In Fig. 4, we use . In this case the model exhibits chaotic attractor behavior. Also, starting from and , the model shows the limit cycle behavior in Fig. 5. Further, by taking , , , and , we get periodic orbit trajectories in the solutions as depicted in Fig. 6.

Figure 4

Simulations for solving (29) using (42) along with

Figure 6

Simulations for solving (29) using (42) along with and

Figure 5

Simulations for solving (29) using (42) along with and

The model via Atangana–Baleanu–Caputo fractional derivative

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 is defined as [10] where is a normalization function.

The Atangana–Baleanu fractional integral of order of a function is also defined as [10] 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: where As , Eq. (47) suggests the exact solution for the model.

Theorem 5

The initial value problem given by Eq. (43) possesses at least one solution in the interval .

Proof

First, we define where , , and are contractions respect to , , and , respectively.

Moreover, we set where Considering the Picard operator, we have defined as follows where , , and .

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 , we must have . 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 Using this numerical approximation, we get the following iterative scheme:

Figures 7–9 are plotted to demonstrate the variation of state variables in model (43) when the scheme (59) is used for different values of . In this simulations, we have considered the following values in the model: , , , , , , , and . In our performed numerical simulations, and . In Fig. 7, we take . In this case the model exhibits chaotic attractor behavior. Also, taking the initial conditions and the model shows the limit cycle behavior in Fig. 8. Moreover, for , , , and , we get periodic orbit trajectories in the solutions as depicted in Fig. 9.

Figure 7

Simulations for solving (43) using (59) along with

Figure 9

Simulations for solving (43) using (59) along with and

Figure 8

Simulations for solving (43) using (59) along with and

The model via fractal-fractional derivative involving different laws

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

Power-law fractal-fractional derivative

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 is defined as [11] and By employing the corresponding inverse operator in Eq. (60) and then placing in the resultant we obtain the following recursive form: Now let us define the functions These functions can be interpolated in as Thus we obtain where , , are given in (63).

Figures 10–12 demonstrate the variation of state variables in model (60) while applying scheme (65) for different values of , and . In this simulations, we considered the following values in the model: , , , , , , , and . In our numerical simulations, and . In Fig. 10, we put . In this case the model exhibits chaotic attractor behavior. Also, by imposing the initial conditions and the model shows the limit cycle behavior in Fig. 11. Further, for , , , and , Fig. 12 confirms the periodic orbit trajectories in the solutions.

Figure 10

Simulations for solving (60) using (65) along with

Figure 12

Simulations for solving (60) using (65) along with and

Figure 11

Simulations for solving (60) using (65) along with and

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 is defined as [11] and is introduced in (62).

Applying the Caputo–Fabrizio integral to Eq. (69), we obtain Setting in (71), based on the proposed scheme in [11], we get Taking the difference between and yields Therefore the approximate solution to the problem can be determined using the following iterative procedures: where , , are given in (63).

Figures 13–15 demonstrate the variation of state variables in model (69) when scheme (74) is applied for different values of and . In this simulations, we considered the following values in the model: , , , , , , , and . In our numerical simulations, and . In Fig. 13, we take . In this case the model exhibits chaotic attractor behavior. Further, for and , the model shows the limit cycle behavior in Fig. 14. Also, letting , , , and , Fig. 15 confirms the periodic orbit trajectories in the solutions.

Figure 13

Simulations for solving (69) using (74) along with

Figure 15

Simulations for solving (69) using (74) along with and

Figure 14

Simulations for solving (69) using (74) along with and

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 is defined as [11] and 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 where

Figures 16–18 demonstrate the variation of state variables in model (75) when scheme (79) is utilized for different values of and . In this simulations, we considered the following values in the model: , , , , , , , and . In our numerical simulations, and . In Fig. 16, we take . In this case the model exhibits chaotic attractor behavior. Also, by taking into consideration and the model shows the limit cycle behavior in Fig. 17, whereast for , , , and , we get periodic orbit trajectories in the solutions as depicted in Fig. 18.

Figure 16

Simulations for solving (75) using (79) along with

Figure 18

Simulations for solving (75) using (79) along with and

Figure 17

Simulations for solving (75) using (79) along with and

In Figs. 19 and 20, we also investigate the sensitivity of the state variables when two parameters and 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.

Figure 19

Simulations for solving (43) using (59) along with and for different

Figure 20

Simulations for solving (43) using (59) along with and for different

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