The role of noise in the tumor dynamics under chemotherapy treatment

Dynamical systems modeling tumor growth have been investigated to analyze the dynamics between tumor and healthy cells. Recent theoretical studies indicate that these interactions may lead to different dynamical outcomes under the effect of particular cancer therapies. In the present paper, we derive a system of nonlinear differential equations, in order to investigate solid tumors in vivo, taking into account the impact of chemotherapy on both tumor and healthy cells. We start by studying our model only in terms of deterministic dynamics under the variation of a drug concentration parameter. Later, with the introduction of noise, a stochastic model is used to analyze the impact of the unavoidable random fluctuations. As a result, new insights into noise-induced transitions are provided and illustrated in detail using techniques from dynamical systems and from the theory of stochastic processes.


Introduction
The term cancer designates generically a group of diseases involving abnormal cell growth with the potential for spreading to new locations and invade other parts of the body. This process contrasts with benign tumors, which do not spread. Cancer can be detected by certain signs and symptoms or screening tests. Nevertheless, despite the encouraging progresses, we have only a limited understanding of cancer and there is no known effective cure for the disease [1]. As a matter of fact, in 2015, about 90.5 million people had cancer. In very recent years, annually, it caused about 8.8 million deaths (15.7% of deaths) (additional information can be found in [2] and [3]). Metastasis is the spread of cancer to other locations in the body. The dispersed tumors are called metastatic tumors. As a matter of fact, most cancer deaths are due to cancer that has metastasized. In order to control this tragic outcome, different types of therapeutic methodologies have been adopted by oncologists, such as surgery, radiotherapy, a e-mail: jesus.seoane@urjc.es (corresponding author) anti-angiogenic drugs, immunotherapy, chemotherapy, among others. The chance of survival depends on the type of cancer and extent of disease at the start of treatment.
In order to gain revealing insights into important phenomena involved in tumor growth and to predict its future behavior, different cancer models have been used. Complex interactions among different body cells, namely tumor cells, healthy tissue cells and activated immune cells, are considered in the modeling process, which usually includes a treatment strategy. Key factors that influence the treatment choice are: tumor severity, age, sex and immune system state of the patient. It has been shown in the literature that cancer growth models are particularly complex in its dynamics (see, for instance, [1,4,5] and references therein).
Depending on the type of cancer and stage, particularly when cells seem to be resistant to immune elimination ( [6,7]), chemotherapy is often applied against tumor progression. Even with their limitations and side effects, chemotherapy may be useful to reduce symptoms such as pain or to reduce the size of an inoperable tumor in the hope that surgery will become possible in the future. It is a fact that the application of chemotherapy remains an important and sensitive research topic to improve [8].
In order to study the reduction of tumor population in an optimal manner, particular treatment protocols, mainly chemotherapy, immunotherapy or radiotherapy, have been taken into consideration on cancer models. Some recent works have focused their attention on solid tumors in vitro [9].
In the present study, with the purpose of enhancing our understanding of tumor growth in vivo in the presence of traditional therapies, we will be concerned with the effect of chemotherapy on both the healthy tissue and tumor cells. In particular, our main goal is to illustrate in simple terms and to provide insights into the dynamical effect of adding noise to a two-dimensional deterministic model, studying a stochastic system under the usage of chemotherapy. More specifically, these insights correspond to a deeper understanding and a wider perspective into the different dynamical scenarios that can emerge from the suggested numerical simulations in a stochastic context, when noise is added to a deterministic system. Accordingly, as a contribution to the literature, we will discuss the use of bifurcation analysis to design treatment protocols of different stages. More precisely, the values of specific parameters at the bifurcation points are biologically meaningful and will offer us the possibility for recognizing the sufficient drug concentration to be used in order to suppress tumor cells. In fact, throughout our work, the bifurcation values of the drug concentration C, precisely C = C 1 and C = C 2 , allow us to characterize dynamical scenarios corresponding to C < C 1 , C 1 < C < C 2 and C > C 2 that guide the entire discussion on the impact of the considered noise intensity. Biologically speaking, this noise intensity corresponds to the amount of random variability in different quantities (e.g., dynamical variables and parameters) arising in cellular biology.
The subsequent part of this article is organized as follows. In Sect. 2, we describe the basic features and dynamical properties of the deterministic version of the model. Divided into two subsections, Sect. 3 is devoted to the analysis of the main topic of the present work-the stochastic effects on the tumor growth model provided by the noise-induced dynamics. The closing Sect. 4 is devoted to the usual final considerations about the undertaken study.

Deterministic model
The mathematical model that we are using was developed in [10], and it corresponds to the one described in [11], though no constant input of effector immune cells is considered. Each variable represents a cell population, where T (t) corresponds to the tumor cells, H (t) to the healthy host cells, and E(t) to the effector immune cells. The growth of cancer and host cells follows a logistic curve with growth rate r i and carrying capacity k i . As a competitive model, the parameters a i j represent the competitive parameters.
The parameter r 3 represents the immune effector cell production rate related to the response of the presence of tumor cells. On the other hand, k 3 corresponds to the steepness of the response curve; that is, the value of the tumor cells at which the immune response rate is half of the maximum production, threshold for which the response curve saturates. Immune effector cells only compete with cancer cells. In the absence of these cells, they die off with a constant per capita rate d 3 . Thus, the system reads as follows: The nondimensionalization and parameter reduction of this system is thoroughly studied in [5,10], yielding the set of equationṡ where x is the density of population of tumor cells, y is the density of population of healthy cells, z is the density of population of immune cells, and a 12 , a 13 , a 21 , a 31 , r 2 , r 3 and d 3 are the dimensionless parameters of the system. For our purpose, we take the two-dimensional system obtained after making z = 0 from Eq. 2. This means that we do not consider the existence of the immune cells. In this situation, our system can be written as where we fix the parameters a 12 = 0.5, a 21 = 4.8 and r 2 = r = 1.2 following [12].
On the other hand, to take into account the impact of chemotherapy, we follow the fractional cell kill law used in [9] and derive the model where the dynamical variables and parameter values are listed and described in Table 1. Terms associated with the treatment contain six parameters: a, b, ρ 1 , ρ 2 , s 1 and s 2 , which have no units since they are dimensionless quantities. We suppose a >> b, because the tumor must be more sensitive to drugs than the healthy cells. The behavior of Eq. 4 is studied under a variation of the parameter C of the drug concentration. Notice that, in the context of our study, we do not need to pay attention with the pharmacokinetics of the drug itself since here the drug concentration is treated as a parameter. As far as the administration of chemotherapy is concerned, this model allows us to represent the so-called Norton-Simon hypothesis. As Notice that both the variables and the parameters have not units since they are dimensionless quantities pointed out in [9], this hypothesis states that the rate of destruction of a cytotoxic agent is proportional to the rate of growth of the unperturbed tissue (tumor cells/healthy cells). The fractional cell kill by cytotoxic agents can then be described by means of an exponential kill model, which has been tested against experimental results. This model describes the rate of cell kill by using the mathematical functions a(1 − e −ρ 1 C ) and b(1 − e −ρ 2 C ), corresponding to the tumor cells and healthy tissue, respectively. In order to gain direct insights into the dynamics of the system in the absence of noise, we first examine the equilibria and perform a phase plane analysis, having C as the most relevant parameter. Indeed, Eq. 3 possesses several equilibrium points called M n (x, y) where n = 0, 1, 2. The trivial one, M 0 (0, 0), is unstable for the considered choice of parameters. On the other hand, coordinates of equilibria M 1 (x, 0) and M 2 (0, y) can be obtained from the equations: Here, we consider only nonnegative roots; that is, the ones that are biologically meaningful. For s 1 = s 2 = 1, we have: Moreover, the system can possess the equilibrium point M 3 (x, y) in which x, y > 0 and they are governed by the equations: where We evaluate now how the chemotherapy with increasing drug concentration C affects the changes in the density of tumor cells x and in the density of healthy cells y. The key dynamical regimes of Eq. 4 under the variation of C in the interval 0 ≤ C ≤ 1 are well seen in the bifurcation diagram shown in Fig. 1. Here, C 1 = 0.2939 and C 2 = 0.6394 are bifurcation points. The equilibrium M 1 (red color) is stable in the interval 0 ≤ C < C 2 . The equilibrium M 2 (blue color) is stable in the interval C 1 < C ≤ 1. The equilibrium M 1 corresponds, according to [9], to the state of an active tumor (AT), while M 2 stands for the dead tumor (DT). Three dynamical regimes can be described in biological terms: (i) in the interval 0 ≤ C < C 1 , that can be interpreted as weak therapy, Eq. 4 is monostable-independently of the initial state, the system transits to the state of AT; (ii) in the interval C 2 < C ≤ 1 of the strong therapy, the system is also monostable, and all solutions tend to the state of DT; (iii) in the interval C 1 < C < C 2 of intermediate strength of therapy, the system (4) is bistable: situation where both regimes, AT and DT, coexist. Depending on the ratio of the tumor and healthy cells, the solution converges either to M 1 (AT) or to M 2 (DT). These two equilibria, M 1 and M 2 , are biologically meaningful. The active tumor (AT) equilibrium corresponds to a malignant attractor for which the tumor has an unacceptable level of tumor cells and can develop an increased metastatic potential. Actually, the DT equilibrium represents the tumor extinction.
For the three dynamical scenarios above mentioned, phase portraits are shown in Fig. 2 for different values of the drug concentration C , namely for C = 0.2 (AT), C = 0.5 (AT + DT) and C = 0.7 (DT). In the bistability case (Fig. 2b), the saddle equilibrium M 3 (empty point) plays an important role: its stable manifold (green dashed) separates basins of attraction of the stable equilibria M 1 and M 2 .
Based on results of the bifurcation analysis depicted in Fig. 1, the following protocol of two-stage treatment can be proposed. At the first stage, we use drug concentration slightly above C 2 . This allows us to suppress tumor cells independently of the initial values of tumor and healthy cells. In the second stage, to provide this state of DT, it is sufficient to use drugs with concentration C slightly above C 1 . Such a protocol of treatment is valid only from the point of view of the deterministic model, since the presence of random disturbances can significantly change the behavior comparing with the deterministic prediction.
At this moment of our study, and in the sequence of the previous deterministic approach, it is necessary to consider the dynamical behavior of system of Eq. (4) under the effect of

Stochastic model
To study the impact of random fluctuations, we will use the following stochastic model Here, ξ 1 (t) and ξ 2 (t) are uncorrelated Gaussian white noises with values of noise intensity ε 1 and ε 2 . In what follows, we consider ε 1 = ε 2 = ε and use a Euler-Maruyama scheme [13] with the time step 10 −3 for numerical simulation of the random solutions. With the purpose of preserving both the physical and the biological sense, we have to keep the variables x and y nonnegative and use a corresponding truncation procedure. Different dynamical scenarios can emerge in a stochastic context, as a result of adding noise to a deterministic system. In this sense, some recent works have been focused on the influence of random fluctuations on the recruitment of effector cells toward a tumor (please see [5,14] and references therein). For the sake of clarity, the following subsections are devoted to noise-induced dynamical scenarios of different nature. Near the right point C 2 = 0.6394 (see Fig. 5 for C = 0.63 ), before the onset of stochastic mixing, the system demonstrates noise-induced transition to the state of M 2 (DT).

Noise-induced excitement in the monostability zone
Now, let us consider the influence of noise in the parameter zone C > C 2 , where the system described in Eq. 4 possesses a single stable equilibrium M 2 corresponding to the state of DT.
For that purpose, we plot Figs. 4, 5 and 6. In Fig. 6, random states (blue) of the stochastic system starting from M 2 are plotted versus the noise intensity along with the mean values of the dynamical variables x and y (light blue).
Even in the monostable system with a singular equilibrium, large-amplitude stochastic oscillations can be induced by noise with increasing intensity. An example of a such excitement is shown in Fig. 7, where time series of stochastic solutions starting from M 2 are plotted for C = 0.7 and ε = 0.01 (blue), ε = 0.1 (green). As we can observe, for ε = 0.01 the system shows small-amplitude noisy oscillations near M 2 . For a larger intensity noise, ε = 0.1, the system exhibits large-amplitude spiking oscillations of complex form.
The transition from small-to large-amplitude stochastic oscillations with increasing noise are accompanied by the significant changes of the form of the probability density. In Fig. 8, one can see these changes for C = 0.64 and C = 0.7 under increasing values of ε.
Here, it should be noted that the form of the probability density function of the variable y, p(y), changes its modality: one peak → two peaks → one peak. Such a qualitative deformation of the form of p(y) is commonly interpreted as an stochastic P-bifurcation. Finally, we consider how these amplitude changes are connected with the frequency characteristics. When studying spikes (which represent quantitative changes in the biological dynamical variables), the mean values of random interspike intervals (ISI), τ , constitute basic statistics that are worth taking in consideration. In practical terms, a short interspike interval means that it has occurred a rapid change and, similarly, a wide interspike interval means that a significant change only occurred after a long period of time. Having established that, in Fig. 9 we plot mean values of ISI, τ , for three values of C versus the noise intensity ε. For weak noises, i.e., low values of ε, the mean values of random interspike intervals, τ , are higher, which indicates that spikes (biological changes) are rare. Therefore, in biological terms, sudden significant changes of the dynamical variables are only expected to occur for increasing values of the noise intensity. The ε-interval, where plots of τ sharply decrease, marks the onset of the noise-induced generation of spikes. The larger the C, the stronger the noise necessary to generate the spiking regime. Biologically, this means that increasing values of the drug concentration require higher values of the noise intensity for significant changes (spikes) of the dynamical variables take place (please observe the three curves of

Concluding remarks
In order to shed light into the study of in vivo tumors, we have designed a two-dimensional mathematical model of nonlinear differential equations, taking into account the impact of chemotherapy on both the tumor and the healthy cells. Studying the model adopting the deterministic and the stochastic approaches allowed us to provide valuable insights into noteworthy and eye-catching features of the dynamics. Interestingly, primary results of the deterministic bifurcation analysis lead us to immediately propose a protocol of a two-stage treatment, involving two distinct levels of drug administration. Taking into account that the tumor growth does not usually take place in a deterministic manner, the addition of noise brings more realism to our study. As a consequence, the discussion of the noise-induced dynamics of the system became particularly relevant. As a matter of fact, it is well known that stochasticity can be found at early stages of tumorigenesis.
Striking features of the dynamical behavior have been pointed out and illustrated within the framework of both the nonlinear dynamics and the theory of stochastic processes using mainly bifurcation diagrams, phase portraits, time series and probability density plots, among others. We found that the dynamics is particularly sensitive to the variation of the drug concentration (chemotherapy), as well as to the noise intensities that have been used in the problem. As a consequence, it has a very important effect on the overall dynamics of the system. In this sense, our theoretical results suggest that a careful control of these parameters from a medical point of view could have a significant role in the fate of tumor cell populations.
This study provides another illustration of how an integrated approach, involving numerical evidences and theoretical reasoning, within the theory of deterministic and stochastic systems, can contribute to our understanding of important biological models. In addition, our work aims to offer a trustworthy explanation of complex phenomena (series of chemical reactions, among other possible events, that result in a transformation) witnessed in biological systems.