Abstract
The complexity of cancer has motivated the development of different approaches to understand the dynamics of this large group of diseases. One that may allow us to better comprehend the behavior of cancer cells, in both short- and long-term, is mathematical modelling through ordinary differential equations. Several ODE mathematical models concerning tumor evolution and immune response have been formulated through the years, but only a few may exhibit chaotic attractors and oscillations such as stable limit cycles and periodic orbits; these dynamics are not that common among cancer systems. In this paper, we apply the Localization of Compact Invariant Sets (LCIS) method and Lyapunov stability theory to investigate the global dynamics and the main factors involved in tumor growth and immune response for a chaotic-cancer system presented by Itik and Banks in 2010. The LCIS method allows us to compute what we define as the localizing domain, which is formulated by the intersection of all lower and upper bounds of each cells population in the nonnegative octant, . Bounds of this domain are given by inequalities in terms of the system parameters. Then, we apply Lyapunov stability theory and LaSalle’s invariance principle to establish existence conditions of a global attractor. The latter implies that given any nonnegative initial condition, all trajectories will go to the largest compact invariant set (a stable equilibrium point, limit cycles, periodic orbits, or a chaotic attractor) located either inside or at the boundaries of the localizing domain. In order to complement our analysis, numerical simulations are performed throughout the paper to illustrate all mathematical results and to better explain their biological implications.
1. Introduction
Cancer is a major public health problem; with 8.8 million deaths recorded in 2015, it is considered to be a leading cause of death globally, second only to heart disease, and it should be noticed that several cases are being reported each year. In general, cancer is a complex and large group of diseases with high mortality rates regardless of age, gender, or race. Incidence and mortality statistics from GLOBOCAN and the WHO indicate that breast, colorectal, lung, cervix, and stomach cancer are the most common types among women, while lung, prostate, colorectal, stomach, and liver cancer are the most common in men [1–3]. Cancer cells have the particularity of avoiding the apoptosis process, which regulates the life cycle of each cell; this allows them to grow beyond their normal limits and progressively invade different organs or tissues through the lymphatic or circulatory system; this process is known as metastasis and it is considered to be the most common cause of death in cancer patients [4].
Therefore, the complexity and diversity of cancer demands that specialists of different areas gather together in order to find a feasible way to understand cancer dynamics in the short- and long-term. One approach of particular interest is mathematical modelling by ordinary differential equations which may allow us to describe tumor evolution without the need for clinical trials, which can be quite costly [5] and are always restricted to one case of study at a time [6].
Mathematical models can help us study the dynamics between cancer tumor growth, the immune system, and different treatments such as chemotherapy or immunotherapy in the short- and long-term. Diverse mathematical models have been propounded with this purpose; see, for example, [7] where authors explore the effects of adoptive cellular immunotherapy and describe under what circumstances the tumor can be eliminated. Two systems that describe tumor growth, immune escape, and siRNA treatment of tumors are shown in [8]; these models predict conditions under which siRNA treatment can be successful in returning an aggressive TGF- producing tumor to a passive nonimmune evading state. Paper [9] aims to express the spontaneous regression and progression of a malignant tumor dynamics as a prey-predator-like system that describes interactions between tumor cells, hunting predator cells, and resting predator cells; further, local stability analysis is performed along with numerical simulations to support the analytical findings. In [10] authors design a six-dimensional model to explain tumor evolution and the corresponding immune response under chemotherapy and immunotherapy treatments. In [11] they formulate and analyze a mathematical model of tumor-immune interactions with chemotherapy; their model allows them to test and compare various optimal control strategies related to drug administration. The first mathematical model that describes tumor-immune interactions in the bladder as a result of the Bacillus Calmette–Guérin (BCG) therapy is developed in [12], whose parameters are estimated by using published data taken from in vitro studies of mouse and human patients. In [13] the model for the basic interactions of a vaccine, the immune system, and prostate cancer cells is presented; this mathematical model is validated by the results of clinical trials, testing an allogeneic prostate cancer whole-cell vaccine. In [14] they analyzed a coupled ordinary differential equation model of tumor-immune dynamics, which accounts for biological and clinical factors that regulate the interaction rates of cytotoxic T lymphocytes on the surface of the tumor mass.
Our paper is focused on investigating the global dynamics of a three-dimensional chaotic-cancer mathematical model formulated by Itik and Banks in [15]. It is important to notice that this system was developed based on the de Pillis and Radunskaya tumor growth with an immune response and chemotherapy system shown in [16]. The mathematical model from Itik and Banks differs from the one of de Pillis and Radunskaya as it does not have a constant influx of effector immune cells. A more detailed explanation about this is provided by the authors in [15] and discussed by Letellier et al. in [17]. One characteristic of particular interest about the system under study is the existence of a chaotic attractor, which is assumed to be produced by the unpredictable growth of cancer cells, their competition for space and resources with healthy host cells, and their interaction with immune effector cells. This chaotic and uncontrolled growth of the tumor population may cause them to reach its carrying capacity. The latter will eventually trigger the metastasis process, which is the most common cause of death in cancer patients. Furthermore, the mathematical model was designed on similar principles as those shown in [7], which have been broadly accepted by the scientific community in both mathematical and biology fields. Hence, our interest in investigating the global dynamics of this chaotic-cancer system by means of the Localization of Compact Invariant Sets (LCIS) method [18, 19] and stability theories like Lyapunov’s direct method and LaSalle’s invariance principle [20, 21].
The remainder of this paper proceeds as follows. In Section 2 we provide the mathematical preliminaries about the LCIS method, describe the chaotic-cancer mathematical model, and provide a brief summary of the work already done on this system related to its dynamical analysis. In Section 3 we compute the bounds of a domain containing all compact invariant sets of the chaotic-cancer system, we establish sufficient conditions for the existence of a Bounded Positively Invariant Domain (BPID) in the nonnegative octant, and, we illustrate all our analytical results by means of numerical simulations. In Section 4 we discuss the biological implications of our mathematical results. Finally, conclusions are given in Section 5.
2. Materials and Methods
2.1. Mathematical Preliminaries
The LCIS method is used to determine a domain on the state space where all compact invariant sets are located. As examples of compact invariant sets, we may recall equilibrium points, periodic orbits, homoclinic and heteroclinic orbits, limit cycles, and chaotic attractors. The relevance of this analysis is due to the fact that it is useful to study the long-time dynamics of the system. Let us consider an autonomous nonlinear system of the form:where is a differentiable vector function and is the state vector. Let be a differentiable function; by we denote the restriction of on a set . The function used in this statement is called localizing and we assume that is not the first integral of . By we denote the set , where represents the Lie derivative of and is given by the following . Now, let us defineandthen, the General Theorem concerning the localization of all compact invariant sets of a dynamical system is established as follows.
General Theorem (see [18, 19]). Each compact invariant set of is contained in the localizing domain:
Localizing functions are selected by an heuristic process; this means that one may need to analyze several functions in order to find a proper set that will allow fulfilling the General Theorem. If we consider the location of all compact invariant sets inside the domain , we have the set , with defined in the General Theorem. It is evident that if all compact invariant sets are located in the sets and , with , then they are located in the set as well. Furthermore, a refinement of the localizing domain is realized with the help of the Iterative Theorem stated as follows.
Iterative Theorem (see [18, 19]). Let , , be a sequence of differentiable functions. Setswith contain any compact invariant set of the system and
This method has already been applied combined with other stability theories to analyze mathematical models that describe cancer tumor growth, the response of the immune system, and the effects of chemotherapy and immunotherapy on cancer evolution. For example, in paper [22] one can see the study of the global behavior of a three-dimensional tumor system, where they find the lower and upper bounds for one of the cells populations and determine conditions under which trajectories go to a specific equilibrium point. In [23] the existence of a positively invariant domain in the positive octant for a cancer-immunotherapy mathematical model is demonstrated. The analysis of the dynamics of a superficial bladder cancer model with Bacillus Calmette–Guérin immunotherapy is shown in [24], in which sufficient conditions of global asymptotic stability of the tumor-free equilibrium point are derived by complementing the analysis with the LCIS method. Sufficient conditions to prove that dynamics of a three-dimensional chronic myelogenous leukemia model around the tumor-free equilibrium point is unstable are presented in [25]; further, ultimate upper bounds are computed for all three cell populations and existence conditions of a positively invariant polytope are provided. In [19] authors study the global dynamics of a cancer chemotherapy system and determine sufficient conditions for tumor clearance and global asymptotic stability of the tumor-free equilibrium point by means of the LCIS method and LaSalle’s invariance principle. By combining the LCIS method with Lyapunov’s stability and LaSalle’s invariance principle we can perform a global analysis of these biological systems; for example, we may compute lower and upper bounds for each variable, determine conditions for the existence of a global attractor, and establish sufficient conditions for tumor clearance.
2.2. The Chaotic-Cancer Mathematical Model
Modelling cancer with ordinary differential equations is an approach that allows us to understand the evolution of tumors in diverse scenarios that will be difficult to study in a patient affected by this disease. Hitherto, as we mention in the last section, diverse mathematical models have been propounded with this spirit. In this paper, we analyze the mathematical model presented by Itik and Banks in [15]. This model describes the dynamics of tumor growth and its interaction with healthy cells and the immune system by the following equations:
Equation (8) represents the rate of change in the cancer cells population. The first term refers to the logistic growth of this population, and the competition between tumor cells and healthy and effector cells is given, respectively, by the second and third term. In (9), healthy host cells also grow logistically and they are assumed to be deactivated by their interaction with cancer cells. As presented in (10), the evolution of effector cells is directly related to the presence of tumor cells as they both stimulate and inactivate them at different rates. Effector cells die naturally as given by the third term. The description of all parameters is shown in Table 1 and was obtained from [15].
Further, it should be noticed that dynamics of the nondimensionalized system (8)–(10) is located in the nonnegative octant:
As one can see, all parameters are positive and units of cells were scaled. The maximum carrying capacity for both tumor and healthy cells populations has been normalized to . Nonetheless, one can trace the origins of the values for the chaotic-cancer system following the path of the work of de Pillis and Radunskaya in [26]. In this paper the authors mention that one unit represents the carrying capacity of the normal cells in the region of the tumor which could be considered on the order of cells (although this depends on the type of tumor). Further, authors suggest that their original mathematical model reflects the qualitative behavior of real tumors with respect to treatment response and not quantitative results that reflect real-life quantities. Concerning the competition parameters (, and ) the model assumes a destructive competition, i.e., , with larger since the competition between immune cells and tumor cells is most detrimental to the tumor than its competition to normal cells. About the per unit growth rates, it is assumed that the tumor cell population grows more rapidly than the normal cell population; hence , but depending on the type of cancer and its stage of growth these rates could change [16, 26].
However, we are interested in one particular dynamic of the system; therefore, we used the values shown in Table 1 as it was demonstrated by Itik and Banks that the phase space solutions of (8)–(10) contain a chaotic attractor for these values; see Theorem 4.1 in [15].
The chaotic-cancer mathematical model (8)–(10) provides important information about tumor evolution through different behaviors such as coexisting equilibria, limit cycles, periodic orbits, and chaotic attractors. All these dynamics are directly related to the bifurcation parameter which represents the competition between healthy and cancer cells. Solutions for and are also bounded from above due to the logistic growth in each equation. About the different dynamics of the system, numerical simulations illustrate the following concerning the values shown in Table 1: when , the system exhibits chaotic behavior; if , the trajectories of the system will begin to show oscillations such as stable limit cycles and periodic orbits, whereas a value implies that all solutions will converge to a healthy tumor-free equilibrium point. Although the system exhibits a wide variety of biologically meaningful dynamics, it is important to notice that it does not take into account the evolutionary evasion of tumors from the immune system which is a very important phenomenon that should be considered in the modelling of a disease as complex as cancer [27–30].
Recently, the global dynamics of model (8)–(10) was discussed by [15, 17, 23]. In general, they conclude that for all all trajectories of the system with feasible nonnegative initial conditions are bounded in . In these previous works one can also find the analytical and numerical computation of equilibrium points as well as sufficient conditions to ensure that all trajectories inside the nonnegative domain will lead to the healthy equilibrium point. Nonetheless, a Bounded Positively Invariant Domain (BPID) has not been reported yet; i.e., an upper bound for the effector cells populations has not been calculated; we determine this bound by means of the LCIS method, and by Lyapunov stability theory we solve the BPID problem. The latter implies the existence of a global attractor; this means that all trajectories will go to a compact domain or to its boundaries and remain there for all future time. Further, numerical simulations are performed for better understanding of all mathematical results.
3. Results
3.1. Localization of Compact Invariant Sets
The LCIS method is applied in order to define the localizing domain illustrated in Figure 1, which is located in the nonnegative octant ; this domain contains all compact invariant sets of the chaotic-cancer system (8)–(10). Boundaries are given by inequalities in terms of the system parameters and were derived by means of five localizing functions. Dynamics of the system, i.e., chaotic attractors, limit cycles, periodic orbits, and some equilibrium points, are directly related to the value of the fractional tumor cells killed or inactivated by healthy host cells, which is given by parameter as indicated in both panels in Figure 1.
All simulations illustrated in this section are developed with the following characteristics: the chaotic attractor with and initial conditions , , ; the limit cycle with and , , ; periodic orbit with and , , ; the set of equilibria shown in the right panel of Figure 1 are plotted with . Initial conditions and values for the parameter are chosen to better illustrate the dynamics in the numerical simulations. Values of other parameters are taken from Table 1. Boundaries of the domain are established in the following result.
Theorem 1. All compact invariant sets of the chaotic-cancer system (8)–(10) are located inside the bounded domainwhere
It should be noted that the boundaries of the domain do not depend on the bifurcation parameter . Upper bounds for the effector immune cells population are restricted by a combination of the system parameters. Furthermore, planes , , and are invariant and we consider them to be the lower bound of each equation in the chaotic-cancer system; this means that , , . Additionally, upper bounds of variables and are evident due to the logistic growth in each equation and can be calculated by applying integration by separation of variables in the absence of immune response and healthy cells competition. Nonetheless, we demonstrate how these bounds can be easily estimated by taking into consideration the biological meaning of the system; i.e., all dynamics are located in nonnegative octant ; this assertion allows us to apply linear localizing functions to solve the bounding problem and estimate all upper limits for all solutions of the chaotic-cancer system (8)–(10). Now, by means of the LCIS method we demonstrate Theorem 1.
First, let us take a simple linear localizing function that will allow us to compute the upper bound of the tumor population; such function is given byand this function is trivial as we need to determine a result of the form , so we do not need to explore different combinations of linear terms or even nonlinear ones.
Now, we compute its Lie derivative as follows:and, then, we get the set , as indicated belowfrom which by taking the corresponding operations we getand thereforeand the maximum value of the cancer cells population in absence of immune response is given by the next resultIn Figure 2, we illustrate our result by plotting the solution of (8) with three dynamics given by different values on the parameter . Since the upper bound is given by the maximum tumor carrying capacity, one can see that for any initial condition inside the bound, the corresponding solution will not go beyond this limit. Mathematically, the latter implies that given any initial condition the corresponding solution will stay inside ; i.e., the domain is positively invariant. Biologically, the tumor carrying capacity could be defined as the maximum size that can be sustained by the surrounding tissue without causing organ malfunction; it can also be seen as a threshold on which the tumor can no longer be sustained by its initial environment; therefore, to ensure its survival the tumor begins the angiogenesis process to metastasize through the body [31]. However, a patient might die before the tumor becomes close to its carrying capacity [32] and, as illustrated in Figure 2, the chaotic behavior has some instants in time where the solution reaches values dangerously close to . The latter phenomenon could be controlled with treatments such as chemotherapy or immunotherapy which could be considered as new equations or control parameters.
Further, in order to get the maximum value of the healthy cells population we exploit the next localizing functionand compute the Lie derivative as shown belowso we get the set , as shown belowand by taking the corresponding operations we getand we formulate the next equationfrom which we are able to conclude that the maximum value of the healthy cells population is given by the next resultThe upper bound is a natural limit to this population since it is considered to be the maximum carrying capacity due to the logistic growth in (9), as discussed on the previous result . We illustrate in Figure 3 that whatever the type of compact invariant set is (limit cycle, periodic orbit, or chaotic attractor), the dynamics is located inside the upper bound given by the set . This is the reason why none of the solutions go further than this limit. As contrary to cancer cells, healthy cells have the ability to control their population size by the apoptosis process. This mechanism can be activated by intracellular signals, also known as the mitochondrial pathway, or by extracellular signals which are mainly produced by activated macrophages and in particular cases, as described by (9) by their interaction with cancer cells. A small population of healthy cells implies more resources for a tumor to grow and become malignant invading adjoining tissues and organs; this is the stage where immune effector cells fulfill the main role to avoid the metastasis process.
Now, by applying the Iterative Theorem we can establish an upper bound for the effector cells populations by proposing a localizing function with a combination of linear and nonlinear terms as given belowOne should notice that the upper bound of the effector cells population is not evident, as it was with the healthy and tumor cells populations. Therefore, by analyzing different functions we conclude that was the one that gave us the best results concerning the upper bound to the solution of (10). The Lie derivative is computed and simplified as shown belowand by considering the next conditionthe set can be written as follows:where ; hence, we disregard this term in order to get the next subsetNow, we substitute and formulate the next inequalitywhich implies the following:and, thus, from the latter set we can determine an upper bound for (10) by applying the Iterative Theorem with and determine the following result:where the existence of the boundary also depends on condition (28). According to the numerical values shown in Table 1, condition (28) holds and equals ; hence we illustrate in Figure 4 our result in the temporal plane for three values of . The bound is not intrinsic to the structure of (10) and it can only be computed by applying the Iterative Theorem. That is why is restricted and it is given by an inequality in terms of the system parameters and the carrying capacity of the cancer cells population, which is a typical result of the LCIS method. Although we did not solve (10), we are able to provide a domain where all its compact invariant sets are located. The importance of our result lies on the fact that biologically represents the maximum value at which the effector cells population can grow. In the next section, we provide conditions for the existence of a global attractor inside the domain ; this attractor could be an equilibrium point, a periodic orbit, a limit cycle, or a chaotic attractor, each one with its corresponding biological implications.
The latter three functions allow us to compute all upper bounds. Now, in order to complement our results, we take two nonlinear localizing functions that allow us to improve the localizing domain where all compact invariant sets of the chaotic-cancer system (8)–(10) are located. First, let us compute the Lie derivative of the functionas follows:and therefore we get the setwhich we rewrite by taking into account the fact that as shown belowand, thus, we apply the Iterative Theorem and we get the next expressionand, hence, we formulate the nextand conclude the following result
Another nonlinear localizing function we analyze is given byand by calculating its Lie derivative we getand therefore we get the setwhich we rewrite by considering that as shown belowNow, we apply the Iterative Theorem and get the setfrom the latter we formulate the followingand, thus, if the next condition holdswe conclude the next resultand, hence, conclude the proof to Theorem 1.
Now, in Figure 5 we illustrate results obtained with the nonlinear localizing functions, and , with a chaotic attractor as an example. Results derived by means of these two nonlinear localizing functions provide us a hint of the overall dynamics between cancer cells and healthy and effector cells. In the left panel of Figure 5 we illustrate that the cancer and the healthy cells population cannot grow into their carrying capacity at the same time, which is expected as they compete for the same resources in the tumor site. On the contrary, both cancer and effector cells populations may grow simultaneously to their maximum capacity, since the existence of the effector cells is directly related to the tumor population; this phenomenon is illustrated in the right panel of Figure 5.
The localizing domains are defined by the sets and , which are both open as one can see in Figure 5; the invariant planes , , and are asymptotes. Nonetheless, the intersection of all our results given by the domain is a compact domain as illustrated in Figure 1.
3.2. Existence of a Global Attractor
Lyapunov’s direct method can be applied together with the LCIS method to derive sufficient conditions under which the bounded domain is positively invariant in the nonnegative octant . The latter implies that for any nonnegative initial condition on system (8)–(10) all trajectories will go to the domain and remain either inside it or at its boundaries. We solved the BPID problem by taking the following candidate Lyapunov function:and imposing the following conditions:andand we are able to conclude the next result.
Theorem 2. If conditions (28), (47), (50), and (51) are fulfilled then the chaotic-cancer system (8)–(10) has a global attractor located either inside or at the boundaries of the domain .
Now, in order to illustrate the previous result, we perform numerical simulations with different initial conditions outside the domain as shown in Figure 6.
One can see that the limit set to the trajectory with initial conditions , , and is given by the chaotic attractor and it is located inside the domain . Therefore, the limit set for the trajectory with initial conditions , , and is given by the healthy equilibrium point which is located at the boundary of the domain. Initial conditions are chosen to better illustrate the dynamics in the numerical simulations and values of parameters are taken from Table 1.
Therefore, to demonstrate Theorem 2, let us compute the Lie derivative of the candidate Lyapunov function as follows:and, thus, by the corresponding mathematical operations we get the following: so, we rewrite as given belowand, then, if the next condition holdsand because the next condition on the system parameters should be fulfilledthen in order to fulfill conditions to apply LaSalle’s Theorem, let us define the domain in , and we have in the complement to in denoted by that the next mathematical expression holds Further, we notice that is a compact domain. Hence for large the domain contains the domain Next, on the boundary of we have that the inequality is fulfilled. So it means that is a bounded absorbing domain. Hence for the trajectory passing through any point in its limit set is not empty and it is a compact invariant set [20, 21]. Therefore, we derive the following:and our proof is now complete. Additionally, in Figure 1 we illustrate that given initial conditions inside the domain all trajectories remain inside for all future time. The main biological implication of Theorem 2 goes to the effector cells populations described by (10); our results allow us to ensure that this population will not grow out of control at any given time, which had not been demonstrated before for the chaotic-cancer mathematical model (8)–(10).
4. Discussion
Biomathematics has been a subject of increasing interest through the last years. Mathematical models can be helpful to comprehend complex biological phenomena such as cancer. These models may provide critical insights concerning the short- and long-term evolution of this large group of diseases.
Nonetheless, having a mathematical model that is able to accurately describe tumor growth, the corresponding immune response, the competition between cancer cells with other healthy cells, and even the effects of treatments is only the first step in order to understand cancer evolution on different real-life scenarios. Hence, it is necessary to establish a methodology that allows the scientific community to develop proper treatment schedules that may be able to decrease tumor burden while also diminishing their side effects, thus, improving the overall health and increasing life expectancy of the patient.
Up to date, several approaches have been developed throughout the years to study and control the dynamics of cancer systems described by ordinary differential equations; one may see [33–35] and many others. Nevertheless, with our methodology we can either investigate the overall dynamics as we demonstrate in this work or derive sufficient conditions for global stability and successfully solve the tumor clearance problem as shown in [19, 24].
As future work we intend to apply our approach to develop our own models concerning tumor evolution, the immune response, and the application of treatments as control parameters, such as chemotherapy and immunotherapy either alone [7, 11] or combined [10], as well as other biological processes involved in cancer growth such as angiogenesis [36] and immune evasion from the tumor [8].
In fact, it was proposed to consider a positive influx of immune cells to further investigate if it is possible to destroy the chaotic regime of the system. The latter has been done in many works, e.g., [7, 16, 26, 37–39]. Hence, we consider the influx of effector cells by adding the parameter to (10) and, as demonstrated through this paper, we apply the LCCI method to compute the lower bound with . Now, as shown in [19], one can apply Lyapunov’s direct method to derive sufficient conditions for tumor clearance, which in this case will be given by . This result is illustrated in Figure 7, where the bifurcation parameter and initial conditions , , , the same values we used for the chaotic attractor in the numerical simulations above.
It is evident that all solutions of the system go to the tumor-free equilibrium point; i.e., , , and as illustrated by the dashed line in all panels in Figure 7. The latter was of course to be expected because, in the absence of tumor cells, healthy cells will grow to their maximum carrying capacity and the density of effector cells will be governed only by the value of the influx parameter . Therefore, with a value of the chaotic behavior of the system is now suppressed. The influx of immune effector cells could be considered either as constant input of immunotherapy or as effector cells that are continuously attracted towards the site of the tumor by chemotaxis.
5. Conclusions
We study the global dynamics of a chaotic-cancer system given by a set of three ordinary differential equations by means of the LCIS method, Lyapunov’s direct method, and LaSalle’s invariance principle. Our main contribution consists in the definition of a Bounded Positively Invariant Domain where all compact invariant sets of the system are located; this localizing domain is given by the set . As we illustrate in Figure 1, the dynamics of equations (8)–(10) are directly related to the value of parameter , which represents the strength of the immune response against cancer.
The LCIS method and the Iterative Theorem allow us to determine all upper bounds for the system (8)–(10). Unlike previous works, we were able to establish the upper bound to the effector cells population ; it should be noticed that this bound is positive according to the values shown in Table 1. Bounds defined in sets and are given by nonlinear functions; these results allows us to improve the initial polytope domain by reducing it with additional excisions. Inside this localizing domain we have equilibrium points, periodic orbits, limit cycles, and chaotic attractors, namely, every interesting dynamics of the chaotic-cancer model (8)–(10).
Additionally, by means of Lyapunov stability theory and LaSalle’s invariance principle we were able to derive sufficient conditions for the existence of a global attractor in the nonnegative octant. The latter implies that all trajectories with initial conditions inside the domain remain inside for all future time, as illustrated in Figure 1. Further, all trajectories with nonnegative initial conditions outside the domain will eventually go inside it or to its boundaries, as illustrated in Figure 6. Concerning the restrictions over the boundaries of the domain and its attractivity proof with condition (51), it is biologically feasible and it is given in terms of the system parameters.
About the biological implications of the results we have determined that global stability can only be achieved if the maximum recruitment rate of effector cells is stronger than its inactivation rate by cancer cells. The latter implies that if the immune system is not able to detect and attack cancer cells, a treatment should be considered to control tumor growth as we briefly elaborated in our discussion section.
Our numerical simulations illustrate in Figure 2 that for the three values of all solutions to the cancer cells population, i.e., the chaotic attractor, the periodic orbit, and the stable limit cycle, at some intervals in time, become very small, reaching values very close to zero. In real life, we can consider one cancer cell as a finite critical value below which tumor survival is impossible. Then, if there is a solution such that goes below this threshold, one should conclude the complete elimination of the cancer cells population; see the works of d’Onofrio [40] and Korobeinikov [30].
As future work, we intend to rescale all parameters in order to investigate the time of occurrence of this event (e.g., ) and calculate, at least empirically, the average extinction time of the tumor. In order to do this, we recall that the units of cells were scaled so that one unit represents the carrying capacity for both healthy and cancer cells. According to different sources, the carrying capacity of a tumor population could be considered of cells [7], cells [10], cells [37], [11]. Authors of the original model on which the chaotic-cancer system (8)–(10) was formulated consider the carrying capacity to be in the order of cells [26]. Additionally, the positive influx of immune cells that we consider in the discussion section could be used to explore the case of time periodic immunotherapies that could induce further chaotic dynamics due to nonlinear resonance to limit cycles or strange phenomena as discussed by d’Onofrio in [38].
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This work is fulfilled within the TecNM project number 6575.18-P “Ingeniería aplicada mediante el modelizado matemático para comprender la evolución del cáncer, la respuesta inmunológica y el efecto de algunos tratamientos”, and the TecNM project number 6178.17-P “Modelos matemáticos para cáncer, VIH y enfisema”.