Abstract
Analyzing the dynamics of tumor-immune systems can play an important role in the fight against cancer, since it can foster the development of more effective medical treatments. This paper was aimed at making a contribution to the study of tumor-immune dynamics by presenting a new model of cancer growth based on fractional-order differential equations. By investigating the system dynamics, the manuscript highlights the chaotic behaviors of the proposed cancer model for both the commensurate and the incommensurate cases. Bifurcation diagrams, the Lyapunov exponents, and phase plots confirm the effectiveness of the conceived approach. Finally, some considerations regarding the biological meaning of the obtained results are reported through the manuscript.
1. Introduction
In the last fifty years, great research efforts and economic resources have been directed to win the fight against cancer. In order to tackle the problem, one of the key issues is to active and control the immune system in its competition against neoplastic cells [1]. To this purpose, the study of tumor-immune dynamics can play a role of paramount importance, given that the mathematical modeling of cancer growth is considered one of the useful tools for the development of effective medical treatments [2, 3]. Over the years, the study of the tumor-immune dynamics has led to the discovery of remarkable phenomena, including the presence of chaos in the system dynamics. By considering integer-order dynamical systems (i.e., biological systems described by integer-order differential equations), in [4], a simple chaotic model of three competing cell populations (host, immune, and tumor cells) is introduced. Topological analysis and computing observability coefficients are illustrated, with the aim to suggest new trends in understanding the interactions of some tumor cells [4]. The authors of reference [5] have suggested a suitable model for the tumor growth, i.e., a discrete-time system capable of exhibiting periodic and chaotic behaviors. The model, which is validated through experimental data, can explain a number of biologically observed tumor states and dynamics [5]. Another interesting model of tumor growth is proposed in [6], based on the interactions among tumor cells, healthy tissue cells, and activated immune system cells. The study, besides analyzing the stability of the system equilibria, highlights the presence of chaotic behaviors in the system dynamics [6]. Referring to biological systems, it should be noted that the behavior of most of these systems has memory or aftereffects [7]. Moreover, biological systems are usually characterized by hereditary properties and nonlocal distributed behaviors. As a consequence, the modeling of these systems by fractional-order differential equations has more advantages than integer-order modeling, in which such effects are neglected [7]. This explains why fractional calculus has recently emerged as a valuable tool for describing a number of dynamic phenomena in biological systems [8]. Regarding tumor-immune dynamics, in [7], a fractional-order model with two immune effectors interacting with the cancer cells is introduced. The conditions that guarantee the stability of the equilibrium points in the considered fractional cancer model are discussed in details [7]. In [9], a mathematical model of cancer chemotherapy effect involving the Caputo fractional derived is presented. In [10], a fractional model of cancer-immune system with Caputo and Caputo-Fabrizio derivatives is investigated. In particular, after examining the stability of the system with singular kernel, the existence and uniqueness of the numerical solution is discussed [10]. In [11], a novel fractional model for a tumor-immune surveillance mechanism is introduced. The approach, besides analyzing the interactions between various tumor cell populations and immune systems, provides an optimal control strategy for investigating the effects of chemotherapy treatments [11]. While references [9–11] have mainly investigated the properties of the mathematical model under consideration, a number of papers have recently focused on the presence of chaos in fractional-order cancer models [12–16]. For example, reference [12] has been one of the first papers to investigate the presence of chaos in fractional-order cancer models. In particular, in [12], the authors have developed a fractional chaotic dynamical model of cancer growth, which includes the interactions between healthy tissue cells, tumor cells, and activated immune system cells. The existence of chaos for the commensurate and incommensurate fractional cancer systems (with order less than 3) is investigated [12]. In [13], a fractional discrete version of a tumor-immune system interaction is analyzed. This model, derived via a discretization process where conformable fractional derivatives are taken into account, exhibits bifurcations and chaotic behaviors [13]. In [14], a study of tumor and effector cells through fractional tumor-immune dynamical model is conducted. By using the Mittag-Leffler law, the paper highlights the existence of chaos in the considered fractional tumor-immune model for cancer treatment [14]. In [15], three-dimensional cancer models that include the interactions between tumor cells, healthy tissue cells, and activated immune system cells are considered. The systems, which are described via Liouville-Caputo, Caputo-Fabrizio, Atangana-Baleanu, and fractional conformable derivatives, show a number of chaotic attractors with symmetric scrolls, depending on the type of the selected derivative [15]. In [16], a cancer model involving the new fractional derivative with the Mittag-Leffler kernel in Liouville-Caputo sense is investigated. A large variety of chaotic attractors is shown, along with the uniqueness and existence of the solutions in the fractional cancer system [16]. Based on these considerations, this paper was aimed at making a contribution to the study of tumor-immune dynamics by presenting a new model of cancer growth based on fractional-order differential equations. By investigating the system dynamics, the manuscript highlights the chaotic behaviors of the proposed cancer model for both the commensurate and the incommensurate cases. Moreover, some considerations regarding the biological meaning of the obtained results are reported. The paper is organized as follows. In Section 2, a novel fractional-order cancer model based on the Caputo derivative is presented. Moreover, a stability analysis of the system equilibria is conducted. In Section 3, by varying the value of the fractional order as well as the values of the system parameters, the dynamics of the commensurate fractional cancer model are analyzed via bifurcation diagrams, the Lyapunov exponents, and phase plots. When the order of the derivative goes beyond the threshold value, , chaotic behaviors are found, indicating that the number of the tumor cells of the healthy host cells and of the effector cells becomes unpredictable. Finally, in Section 4, the dynamics of the incommensurate fractional cancer model are analyzed in details by varying the value of the fractional order in each system equation. Simulation results reported through the manuscript highlight that the proposed approach can explain many biologically observed tumor states (including stable, periodic, and chaotic behaviors), indicating that under some conditions the interactions between tumor cells, healthy tissue cells, and activated immune system cells could lead to invasive tumor growth.
2. Fractional-Order Cancer Model and Its Equilibrium Points
A three-dimensional integer-order cancer growth model has been studied in [6]. Its dynamic equations are described by [17]: where denotes the number of tumor cells at time , is the number of healthy host cells at time , and refers to the number of effector immune cells at time in the single tumor-site compartment. Here, the parameters , , and are positive real numbers representing the growth rates of populations of , , and (see [6]). Specifically, the parameter represents the growth rate of the tumor cells (measured in sec-1), the parameter is the growth rate of the healthy host cells (measured in sec-1), whereas represents the growth rate of the effector immune cells (measured in sec-1). Generally, the model parameters are chosen such that the system dynamic analogies with clinical evidences reported in literatures [4, 18], where depending on control parameter values and initial conditions, the considered biological cancerous system should also approach different states [19]: stationary equilibrium state where any changes are damped, stable periodic process (a limit cycle), and state of instability with chaotic behavior. As shown in [6], particular values of these growth rates lead to make the behavior of system (1) chaotic. To this purpose, the chaotic attractor of system (1) for parameters and initial conditions is shown in Figure 1.
Herein, the fractional version of system (1) is considered. Namely, the dynamics of the proposed fractional-order cancer model (FOCM) are described by where is -order Caputo differential operator, are the derivative orders of the state variables , and (see appendix A). The fractional-order system (2) is called as commensurate if and incommensurate otherwise.
Using the definitions (Equations (B.3) and (B.4)) in Appendix B, then the numerical solution of the FOCM can be given as in (3) with parameters defined in (4) and (5), where and .
Note that system (2) has five equilibrium points [6], four of them are obtained analytically and can be described as follows: (1)(2) if (3) if (4) if
The last equilibrium corresponding to the case does not possess an analytical expression. It could be obtained by intersecting the three surfaces corresponding to the following equations:
By taking the system parameters and , the fixed points become , , , and the fourth equilibrium point found as , as shown in Figure 2.
Note that the fixed points and have negative coordinates, indicating that the dynamics cannot take place since it is not possible to define negative populations in system (2). The fixed point , which corresponds to a situation where there is no cell at all, is unstable, since its eigenvalues are given by . The fixed point , which is associated with the coexistence of the three different types of cells, represents a saddle-focus equilibrium, since its eigenvalues are given by .
3. Dynamics of the Commensurate Fractional-Order Cancer Model
In this section, the dynamics of the proposed commensurate fractional-order cancer model (2) are studied by varying the fractional-order and the system parameters , , and . The bifurcation diagrams, Lyapunov exponents, time behaviors, and phase plots are illustrated to investigate the system dynamics in detail. Moreover, some considerations regarding the biological meaning of the obtained results are reported.
3.1. Analysis of the System Dynamics by Varying the Fractional-Order
The study of the stability of the equilibria is important to understand the system dynamics in the proposed cancer model. Herein, analytical and numerical analyses are conducted to determine the behavior of the system trajectories when the value of the fractional order is properly varied. To this purpose, a theorem proved in reference [20] is now exploited.
Theorem 1. Given the fractional system (2), a necessary condition to have a chaotic attractor around the equilibrium point is that the eigenvalues of its Jacobian matrix satisfy the condition [20]: By taking the fractional system (2) with parameters , and , the eigenvalues of the Jacobian matrix are evaluated at the equilibrium point are given by . By considering that the application of Theorem 1 to the equilibrium gives it follows that a necessary condition to have a chaotic attractor in the fractional system (2) is to satisfy the condition .
In order to investigate the system dynamics and numerically search for proper values of the fractional-order which is able to generate chaotic behaviors, the bifurcation diagram is plotted in Figure 3 for and initial conditions .
From the bifurcation diagram, it can be seen that system (2) is asymptotically stable when , whereas a number of periodic windows appear for . Moreover the FOCM (2) exhibits chaotic behavior for as confirmed by the positive values of the maximum Lyapunov exponents (see Figure 4). From the biological point of view, this behavior can be explained as follows. When , the system become fractional and, consequently, memory effects and hereditary properties appear in the modeling of the system dynamics. When these effects are not so strong (i.e., ), the system dynamics undertake chaotic behaviors. On the other hand, when these effects become stronger (i.e., ), they overwhelm the system dynamics, which undertake stable behaviors.
By varying the value of the fractional-order , Figure 5 shows the time behaviors of the three state variables: , and (in red, blue, and green color, respectively) along with the corresponding phase portraits in the - plan, for the system parameters and initial conditions . When , it can be observed that the FOCM (2) is asymptotically stable and the system trajectories converge to the equilibrium point (Figure 5(a)). When , the system loses its stability and a scroll begins to appear around the point (Figure 5(b)). By increasing the values of , periodic attractors appear for and (Figures 5(c) and 5(d)). When , the fractional cancer model (2) exhibits a chaotic attractor (Figure 5(e)), which is similar to that one obtained for the integer-order case (Figure 5(f)).
(a)
(b)
(c)
(d)
(e)
(f)
A projection in the 3D space of the chaotic attractor generated by the proposed fractional-order cancer model is plotted in Figure 6 for . The conducted analyses clearly indicate that, in order to get chaos, the theoretical condition expressed by Theorem 1 is numerically fulfilled when . From this results it can be concluded that, when the value of the fractional order decreases, the system becomes stable, indicating that the number of the tumor cells, of the healthy cells, and of the effector cells asymptotically converges to the equilibrium point. On the other hand, when the order of the derivative increases and goes beyond the value of , the dynamics of the proposed FOCM turn to be chaotic, indicating that the number of tumor cells, of the healthy host cells, and of the effector cells becomes unpredictable.
3.2. Analysis of the System Dynamics by Varying the Parameters , , and
Herein, the analysis of the system dynamics is conducted by taking the fractional-order and the initial conditions and by varying the parameters , , and . At first, the parameters and are selected as and , whereas the parameter is varied in the interval . Note that the parameter is related to the growth rate of host cells. Since the best strategy to face the cancer dynamics, from the biological point of view, is to act on the healthy host cells [4], herein the parameter is varied, with the aim to investigate the behavior of the proposed cancer model. The bifurcation diagrams of the three state variables , , and of the FOCM (2) are shown in Figure 7, where is the bifurcation parameter.
(a)
(b)
(c)
By varying the value of the parameter , Figure 8 shows the time behaviors of the three state variables: , and (in red, blue, and green color, respectively) along with the corresponding phase portraits in the - plan for the system parameters and initial conditions . When , it can be observed that the FOCM (2) is asymptotically stable and the system trajectories converge to the equilibrium point (Figure 8(a)). When the parameter increases, a periodic route to chaos appears in the range . In this range of parameter , the system exhibits limit cycles of different periods (see Figures 8(b) and 8(c)). Then, a chaotic attractor appears at (see Figure 8(d)) and the system exhibits a chaotic behavior for .
(a)
(b)
(c)
(d)
Now, the parameter is fixed at the value , whereas the parameter , which represents the growth rate of the tumor cells, is varied in the interval . The corresponding bifurcation diagram for the three state variables , , and is plotted in Figure 9(a). Similarly, by fixing the value , the parameter (i.e, the growth rate of the effector cells) is varied in the interval . The corresponding bifurcation diagram for the three state variables , , and is plotted in Figure 9(b). By analyzing the two bifurcation diagrams, it can be argued that the FOCM (2) loses its stability when the values of the parameters and are increased. Moreover, chaotic behaviors appear in the FOCM (2) when and .
(a)
(b)
Regarding the biological meaning of these results, it should be noted that, for low values of the growth rates, the FOCM (2) has a stable equilibrium point. On the other hand, when the growth rates increase, the system loses its stability. At this stage, the tumor is ready to become invasive and even malignant [5]. With the further increase of the growth rates, the chaotic attractor of the tumor appears, being this higher tumor burden complicated by the presence of several periodic and chaotic dynamics. This is similar to what happens when parameter increases in Figure 8. On the other hand, when the growth rates decrease, the attractor corresponding to high tumor burden disappears. This is similar to what happens when , , and decrease in Figures 7 and 9. These results can help the doctors for controlling the tumor burden, thus giving suggestions regarding the medical treatments.
3.3. Comparison between the Dynamics of Integer-Order and Commensurate Fractional-Order Cancer Models
Now, comparisons between the dynamics of integer-order and commensurate fractional-order cancer models are carried out. The time behaviors of the state variable (representing the tumor population) are plotted in Figure 10 by selecting , , , and and by taking different values of the system parameters , , and ranges. It can be observed that for smaller values of the parameter (Figure 10(a)), of the parameter (Figure 10(c)), and of the parameter (Figure 10(e)), the commensurate fractional derivatives damp the oscillation behavior. Consequently, the three states of tumor, host, and effector cells approach faster the equilibrium point, indicating that the commensurate fractional derivatives enlarge the region of stability. When the values of the parameters , , and increase, the system is stable for small values of the fractional orders (i.e., and ). Namely, by looking at Figures 10(b), 10(d), and 10(e), it can be observed that the system trajectories tend to the equilibrium point for and . On the other hand, chaotic oscillations with different amplitudes appear when for and . Note that the amplitude of the chaotic oscillation reaches the maximum value for the integer-order case ().
(a)
(b)
(c)
(d)
(e)
(f)
4. Dynamics of the Incommensurate Fractional-Order Cancer Model
This section analyzes the dynamics of the incommensurate FOCM (2) by taking the parameters , , and and initial conditions and by selecting different values of the fractional-orders , , and . At first, the bifurcation diagrams of the variable are plotted in Figure 11(a) for three cases: and ; and and . From Figure 11(a), it can be seen that the equilibrium point is asymptotically stable when , , and . When the values of increase, periodic windows appear for , , and , whereas chaotic behaviors are exhibited for , , and . The existence of positive Lyapunov exponents is confirmed by the plot as a function of the fractional-order , as shown in Figure 11(b). Namely, from Figure 11(b), it can be seen that the fractional cancer system (2) is chaotic for , , and . Note that the maximum value of the variable is obtained by varying (see Figure 11(a)).
(a)
(b)
Figure 12 shows the time behaviors of the three state variables , , and (in red, blue, and green color, respectively) along with the corresponding phase portraits in the - plan, for the system parameters and initial conditions . By taking different values of the fractional-order , , and in the incommensurate FOCM (2), some chaotic attractors appear. For example, Figure 12(a) plots the chaotic attractor obtained for and whereas Figures 12(b) and 12(c) illustrate the chaotic attractors obtained for and for , , respectively. By looking at the time behaviors of the state variables, it can be noticed that the maximum amplitudes of the trajectories change from one plot to the other when the incommensurate orders are varied. Specifically, the population of the healthy host cells (i.e., the state variable ) is the largest when (see Figure 12(a)), the population of the effector immune cells (i.e., the state variable ) is the largest when (see Figure 12(b)), whereas the population of the tumor cells (i.e., the state variable ) is the largest when (see Figure 12(c)).
(a)
(b)
(c)
Now, by fixing the system parameters and , the bifurcation diagrams for the variable as a function of the parameter are derived for three cases: , ; , ; and (see Figure 13). It can be noticed that the incommensurate system (2) exhibits chaos in all the three cases when .
(a)
(b)
(c)
Figure 14 presents the chaotic attractors of the incommensurate FOCM (2) in 3D projection by taking the parameter for: (a) , ; (b) , ; and (c) , . By comparing this chaotic range with the range that has been obtained in Section 3 (see Figures 7 and 8), it can be observed that the incommensurate fractional derivatives enlarge the chaotic range of the solution. From the biological point of view, it can be deduced that, when the growth rate decreases, the attractor corresponding to the high tumor burden disappears. This is in accordance with the results in Figure 13, since when decreases the system dynamics go towards stable behaviors.
(a)
(b)
(c)
Now, comparisons between the dynamics of integer-order and incommensurate fractional-order cancer models are carried out. The time behaviors of the state variable (representing the tumor population) are plotted in Figure 15 by selecting different values of the fractional-orders , , and , when the parameter assumes the two values (corresponding to the stable range) and (corresponding to the chaotic range). It can be observed that for (see Figures 15(a), 15(c), and 15(e)), the incommensurate fractional derivatives damp the oscillation behavior. Consequently, the three states of tumor, host, and effector cells approach faster the equilibrium point, indicating that the incommensurate fractional derivatives enlarge the region of stability. When the parameter assumes the value , the system is stable for small values of the fractional orders (i.e., , , , , , ). Namely, by looking at Figures 15(b), 15(d), and 15(e), it can be observed that the system trajectories tend to the equilibrium point for , , and , , . On the other hand, chaotic oscillations with different amplitudes appear when , , and , , . Note that the amplitude of the chaotic oscillation reaches the maximum value for the integer-order case (, , ). The motivation of the manuscript is to provide a complete study of tumor-immune dynamics by presenting a new model of cancer growth. In order to explain the physical meaning of introducing the fractional-order into the model, it is worth noting that biological systems are characterized by memory or aftereffects, hereditary properties, and nonlocal distributed behaviors [7]. Since these features are neglected in integer-order modeling, this has motivated the use of fractional calculus as a tool for accurately describing dynamic phenomena in tumor-immune systems. The main advantage of the results in this paper, compared with others published in the literature, is that our approach represents an exhaustive study of tumor-immune dynamics, since it includes the bifurcation diagrams, Lyapunov exponents, and phase plots for both the commensurate and the incommensurate cases. No paper published in the literature so far (to the best of the authors’ knowledge) includes such a complete analysis of the fractional chaotic dynamics of tumor-immune systems [12–16].
(a)
(b)
(c)
(d)
(e)
(f)
5. Conclusion
This paper has made a contribution to the study of tumor-immune dynamics by presenting a new model of cancer growth based on fractional-order differential equations. By investigating the system dynamics, the manuscript has highlighted the chaotic behaviors of the proposed cancer model for both the commensurate and the incommensurate cases. In particular, by using the bifurcation diagrams, Lyapunov exponents, phase plots, and a necessary condition to get chaos, the paper has shown that, when the order of the derivative goes beyond the threshold value , different chaotic behaviors are found, indicating that the number of the tumor cells, of the healthy host cells, and the effector cells becomes unpredictable. Finally, simulation results reported through the manuscript have highlighted that the proposed approach can explain many biologically observed tumor states, including stable, periodic, and chaotic behaviors. Regarding open research problems, an important issue is related to the development of control techniques for suppressing chaos in fractional-order biological systems. Our future plan is to work on this issue, since we believe that controlling chaos in fractional tumor-immune systems might help biologists in the fight against cancer.
Appendix
A. Basic Concepts
We state certain key preliminaries in regard with the noninteger calculus [21]:
Definition A.1. The integral operator of fractional-order in the sense of Riemann-Liouville of the function is outlined as where , and .
Definition A.2. The differential operator of fractional-order in the sense of Caputo of the function is outlined as where , and .
B. Numerical Method for Solving Fractional Differential Equations
To approximate the fractional-order dynamical system using the ABM and PECE numerical approximation method consider [22] where is the fractional-order, with initial values for . Equation (B.1) can be solved using the Volterra integral equation (VIE) given by
The numerical approximation form of (B.2) can be defined as where
, and with . The error for this method can be estimated as , where and .
Data Availability
The data that support the findings of this study are available upon request.
Disclosure
An earlier version of this manuscript has been presented in 2021 International Conference on Recent Advances in Mathematics and Informatics (ICRAMI) Tebessa, Algeria.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The research was supported by the Taif University Researchers Supporting Project number (TURSP-2020/77), Taif University, Taif, Saudi Arabia.