In silico tumor control induced via alternating immunostimulating and immunosuppressive phases

A. I. Reppas, J. C. L. Alfonso and H. Hatzikirou

Center for Advancing Electronics, Technische Universität Dresden, 01062 Dresden, Germany.

Corresponding author: haralambos.hatzikirou@tu-dresden.de


Despite recent advances in the field of Oncoimmunology, the success potential of immunomodulatory therapies against cancer remains to be elucidated. One of the reasons is the lack of understanding on the complex interplay between tumor growth dynamics and the associated immune system responses. Towards this goal, we consider a mathematical model of vascularized tumor growth and the corresponding effector cell recruitment dynamics. Bifurcation analysis allows for the exploration of model’s dynamic behavior and the determination of these parameter regimes that result in immune-mediated tumor control. Here, we focus on a particular tumor evasion regime that involves tumor and effector cell concentration oscillations of slowly increasing and decreasing amplitude, respectively. Considering a temporal multiscale analysis, we derive an analytically tractable mapping of model solutions onto a weakly negatively damped harmonic oscillator. Based on our analysis, we propose a theory-driven intervention strategy involving immunostimulating and immunosuppressive phases to induce long-term tumor control.


The immune system is widely recognized for its capacity to detect and destroy cancer cells, as well as to prevent tumor recurrence maintaining an immunological memory [1]. Indeed, every known innate and adaptive immune effector mechanism has been reported that participates in tumor recognition and rejection [2]. However, experimental observations support that tumorigenic processes themselves can promote immunosuppression or immune tolerant states that facilitates neoplastic growth and progression [3, 4]. Cancer cells employ diverse strategies to inhibit or block immune responses, including tumor-induced impairment of antigen presentation, secretion of immunosuppressive cytokines or expression of surface molecules, as well as production of diverse pro-apopoptic factors [5]. Nevertheless, there are clinical and preclinical evidences supporting that activation of the innate antitumor immunity can result in tumor regression and provide therapeutic benefits [6].

The main goal of oncoimmunology is to strengthen the immune system’s innate ability to combat and kill cancer cells by enhancing the effectiveness of the immune responses. Among the different immunotherapeutic techniques are checkpoint inhibitors, immune response modifiers (cytokines), monoclonal antibodies and vaccines [7]. Passive and active immunotherapy has been successfully applied to the treatment of a wide variety of human cancers [8] and holds promise of a lifelong cure [9]. However, tumor-induced immunosuppression still represents a major obstacle to effective cell-mediated immunity and immunotherapy [5, 3]. Accordingly, more insights into the main mechanisms associated with immune responses based on tumor specific features are required to obtain successful therapeutic outcomes with immunomodulatory strategies. Despite years of research devoted to understand the underlying mechanisms of immune-tumor interactions, there are still many unanswered questions. In particular, those related with the impact of tumor-associated vascularization on immune responses, as well as determination of optimal and effective therapeutic protocols in cancer immunotherapy, are far from being completely elucidated.

Mathematical oncology is a valuable descriptive and predictive analytic framework to address such open questions. Continuum Mechanics concepts have been widely considered to investigate tumor growth and therapy implications [10, 11, 12, 13, 14]. In addition, several mathematical models of tumor growth, where some forms of the immune dynamics are often included, have been also extensively studied in the last years [15, 16, 17, 18, 19, 20, 21, 22, 23, 24]. Clinical data evidence that cancer cells can survive in a undetectable dormant state for extended periods of time [25], which has been also predicted by several models of tumor-immune cell interactions [19, 20, 26, 27, 22]. However, the neoplasm develops diverse strategies to circumvent the anti-tumor action of the immune system [28, 29]. In particular, this equilibrium state can be disrupted by different events affecting the immune system which could result in tumor regrowth [29]. Sustained oscillations by the immune system have been observed both in its healthy state and pathological situations [30, 31, 27]. Therefore, the presence of an immune component in mathematical modeling has been described crucial for reproducing clinically observed phenomena such as tumor dormancy, oscillations in tumor size and spontaneous tumor regression [32, 33, 34, 35, 27, 36]. Among the several reviews on the subject are those covering mathematical models of tumor growth mainly focused on the cancer and immune system interactions [37, 38, 39, 40, 41, 42, 43, 44, 45]. The mathematical modeling of the entire immuno-oncology dynamics is an enormously difficult and complex task. In consequence, models describing interactions between growing tumors and immune dynamics should focus on the crucial factors that are known to allow tumor escape from immuno-surveillance. Tumor-induced angiogenesis is a crucial mechanism for cancer survival and proliferation, allowing a continuous supply of oxygen and nutrients needed for tumor growth and progression [46, 47, 48]. However, the effectiveness of antitumor immune responses is associated with the functional levels of the tumor blood vessels, which allow a wider range of effector cell types to penetrate the tumor bulk and further exterminate cancer cells [49, 50]. These opposing effects demand for a mathematical model of vascularized tumor growth that allows to explore the therapeutic potential of immunomodulatory interventions when innate immune responses are insufficient for long-term tumor control.

The work herein reported intends to yield new insights in the potential of immunomodulatory interventions for cancer therapy. To this end, we propose a tumor-effector cell model based on well-known biological assumptions that combines a model of radially symmetric tumor growth with an immune cell recruitment model [51, 52, 32, 20]. The main feature of our model is the modeling of the interplay between functional tumor-associated vasculature and effector cell dynamics as described in [53]. Model results predict that, depending on the functional tumor vascularization degree and effector cell recruitment rate, long-term tumor control cannot be always reached. We particularly focus in such situations where tumors escape immuno-surveillance to suggest an optimized theory-driven therapeutic strategy against tumor growth. A temporal multiscale approach is then implemented to describe the tumor-immune system interactions, where an analytically tractable approximation of the cancer-effector cell dynamics is derived. We find that an efficient modulation of the immunostimulating and immunosuppressive phases could induce long-term tumor control.

Mathematical model description

The present model describes the interactions between growing tumors and induced immune system dynamics. More precisely, the proposed tumor-effector cell model is a combination of a radial tumor growth model and an effector cell recruitment model originally proposed in [51] and [32] respectively, see also [52, 54, 55, 20]. The main feature is the low dimensional modeling of the complex interplay between tumor-associated functional vasculature and immune recruitment dynamics as described in [53]. This model can be interpreted as the temporal evolution of the average tumor radius, since radially symmetric growth is not a realistic behaviour. The system variables are the average tumor radius and effector cell concentration in the tumor vicinity.

Radial tumor growth,

The temporal evolution of the average tumor radius is considered, where for simplicity invasive and diffusive tumor properties are not taken into account. The tumor is modeled as an incompressible fluid flowing through a porous medium, where tissue elasticity is simplified. The tumor-host interface is assumed to be sharp and cell-to-cell adhesive forces are modeled as a surface tension at that interface. The tumor expands as a mass whose growth is governed by a balance between cell birth (mitosis) and death (apoptosis and necrosis). The mitotic rate within the tumor is assumed to be linearly dependent on the nutrient concentration (oxygen, glucose, etc.) and is characterized by its maximal value at the tumor-host interface. The death rate is uniform within the tumor and constant in time. Moreover, we assume that the death rate reflects the lump effect of apoptotic/necrotic processes and any other cell death inducing factor. The concentration of nutrients (e.g. oxygen or glucose) obeys a reaction-diffusion equation in the tumor volume, where nutrients are supplied from the tumor-associated functional vasculature and consumed by the tumor cells at a uniform consumption rate.

To gain insight on the impact of vascularization in tumor growth and immune responses, we assume that the non-negative and dimensionless parameter , where , represents the net effect of functional tumor-associated vasculature on the tumor radius evolution. In the limit of avascular tumor growth , such tumor-effector interactions take place only at the tumor surface. At the other extreme, for effector cells can potentially interact with any cancer cell within the tumor bulk. Moreover, we consider an intrinsic length scale representing the average length of nutrient gradient, i.e. supply, diffusion and consumption.

The efficacy of immune killing depends on the ability of effector cells to penetrate the tumor bulk via the functional tumor-associated blood vessels. With improved vascularization, the effectors kill tumor cells not only on the surface of the tumor, but also further inside [56]. We consider this process through a phenomenological scaling function , for , that models the penetration of effector cells in the tumor parenchyma through the existing functional vasculature. This function modulates the term related to tumor-effector cell interactions, such as killing of tumor cells due to effectors represented by a rate . Such scaling functions have been also considered in the classical von Bertalanffy approach, and more recently in allometric models [57].

Taking these factors into account, we deduce the tumor radius dynamic under the assumption of radial symmetry according to the following ODE equation:


where the time coordinate has been omitted for notational simplicity, and , , and are non-negative constants.

Effector cell concentration,

We assume that effectors are recruited at a rate depending on tumor cells following the Michaelis-Menten kinetics, where is a positive constant denoting the concentration at which the immune recruitment is half-maximal. Effector cells die at a rate , and become inactivated a rate due to their antitumor activity. In particular, the inactivation of effectors by tumor cells is modeled through the function described above. As functional vascularization increases, effectors can kill cancer cells throughout the tumor bulk and tumor-immune cell interactions increase resulting in more inactivated immune cells. Moreover, innate immunity or base immuno-surveillance is represented as a minimum presence of active effector cells at any time given by an innate immune growth rate , even in the absence of tumor cells.

The resulting ODE equation for the effector cell concentration dynamic is given by:


where the time coordinate has been omitted for sake of simplicity in the notation, and , , , and are non-negative constants.

Model parameterization

Under the small tumor radius assumption, the very early tumor growth is always of exponential nature and does not depend on the vascularization effects, i.e. parameter , [58, 52, 59, 54]. Accordingly, we assume that this initial growth depends exclusively on the net proliferation rate in the absence of adaptive antitumor immune responses at those stages of growth, see the first term of Eq. (1). Thus, considering experimental estimates of the growth rate at early phases of spheroids tumor growth for the mouse colon carcinoma cell line CT26 as day [60, 53] and the physiological plausible value day for CT26 murine cells [60, 61, 53], we have that day. The characteristic nutrient diffusion length has been experimentally estimated that ranges between 0.2 and 0.3 mm [62, 63, 64, 53]. The effector cell characteristic concentration is at the order of magnitude cells. The latter estimate is justified since the characteristic length scale of the system is at the order of 1.0 mm, and given that cells are commonly assumed with a diameter between 10m and 20m [61], then for a volume of the concentration is at the order of magnitude cells. Moreover, we consider cells days as measured from murine CT26 tumor growth experiments in [53], see also [32, 20, 55]. The remaining parameter values day, mm days, mm and cells days are considered from previously reported experimental data [32, 20, 55, 65, 66], and properly rescaled to the magnitudes and units considered in our model. Through a parameter sensitivity analysis, we explore the effects on tumor growth of the effector cell recruitment rate and functional vascularization degree for different choices of the initial tumor radius and concentration of effector cells .

For convenience of the reader, we summarize in Tab. 1 the parameter values used in numerical simulations of the tumor-immune dynamics considering the effect of tumor-associated vasculature.

Description Parameter Value Units Sources
Tumor mitotic rate 1.34 days [60, 61, 53]
Tumor death rate 0.14 days [60, 61, 53]
Characteristic nutrient diffusion length [0.2 - 0.3] mm [62, 63, 64, 53]
Tumor cell kill by effectors 0.03 cells days [32, 20, 55, 53]
Tumor volumen where is half-maximal 2.72 mm [32, 20, 55, 66]
Effectors inactivation rate by tumor cells 0.01 mm days [32, 20, 55, 66]
Effectors death rate 0.37 days [32, 20, 65, 66]
Innate immune growth rate cells days [32, 20, 55, 66]
Table 1: Model parameters. The effects of model parameters and , i.e. functional tumor-associated vasculature and effector cell recruitment rate respectively, are investigated by a model sensitivity analysis.

Dynamical analysis of the model

In this section, we analyze the model’s behaviour with respect to two parameters, namely the effector cell recruitment rate and functional tumor-associated vasculature .

Fixed point analysis

The first step towards analyzing the model dynamics is the identification of the fixed points along with their stability classification. Fig. 1(A-D) depicts the phase portraits of the system of equations (1)-(2) for different values of the model parameter , while keeping constant and equal to  days. The black curves represent the nullclines, i.e. curves along which and , and the colored curves the system trajectories for different initial conditions. The fixed points of the system are located at the intersection of the nullclines. In each case, we identify the existence of two fixed points corresponding to a low tumor radius ( and a high tumor radius (.

Figure 1: Phase portraits and long-term characteristic dynamics. (A-D) Phase portraits of the tumor radius versus the concentration of effector cells for different functional levels of the tumor-associated vasculature . The colored lines depict trajectories starting with the initial conditions (green) and (purple), respectively. The nullclines (zero-growth isoclines) of the dynamical system (black lines) are also plotted. (E-H) Temporal evolution of (red lines) and (blue lines) which correspond to the trajectories in panels (A-D) starting at the initial condition . (I-L) Temporal evolution of (red lines) and (blue lines) which correspond to the trajectories in panels (A-D) starting at the initial condition . The immune recruitment rate is  days and the remaining parameters are as in Tab. 1.

Fig. 1(E-L) shows the evolution of the tumor radius and effector cells of the corresponding trajectories on the phase plain, presented in Fig. 1(A-D). A quick glance at Fig. 1 reveals that the system trajectories initiated near the fixed point follow an oscillatory behavior with a slowly varying amplitude that can be either increasing Fig. 1(F,G) or decreasing Fig. 1(E,H). This behavior indicates that, depending on the model parameter values, the fixed point can be an attractor or a repellor. On the other hand, the trajectories initiated away from this point follow an exponential behavior which results in a boundless tumor growth while the amount of the effector cells fades out (see Fig. 1(I-L)).

To gain insight in the system behavior near the fixed points we construct the bifurcation diagrams with respect to the model parameter for different values of (see Fig. 2(A-C)). The bifurcation diagrams were calculated by performing an arc-length continuation method. This method is a special case of numerical fixed-point continuation methods that ensures the continuation of solution branches at turning points [67, 68].

For sufficiently small values of , we obtain that there are no fixed points, which consequently implies that the tumor grows indefinitely. However, for a critical tumor radius , a homoclinic bifurcation occurs [69, 70] giving rise to two states: a lower branch which corresponds to the fixed point, and an upper branch which corresponds to the fixed point. The local stability analysis shows that the fixed point is a saddle point, and therefore unstable. On the other hand, the fixed point is a spiral point which, depending on the value of the parameter , can be unstable (spiral source) or stable (spiral sink) [71, 72]. Fig. 2(D,E) shows the local stability analysis of and for , which corresponds to the bifurcation diagram in Fig. 2(B).

Figure 2: Bifurcation diagrams with respect to and local stability analysis. (A-C) One-parameter bifurcation diagrams with respect to the effector cell recruitment rate for different values of the functional tumor-associated vasculature (0.40, 0.60 and 0.95, respectively). The upper branches correspond to the saddle point , whereas the lower branches to the spiral point . Solid lines depict stable fixed points and dotted lines the unstable fixed points. (D) Eigenvalues of the Jacobian estimated at the saddle point with respect to the parameter for . (E) Real part of the eigenvalues of the Jacobian estimated at the spiral point with respect to the parameter for .

In the case of being a spiral sink, a system trajectory initiated inside the homoclinic orbit, i.e. the closed orbit which starts and returns to the saddle point in Fig. 3(A), will follow regular oscillations with a slowly decreasing amplitude until it reaches the fixed point. This implies that the tumor radius will stay in a control bounded state. The homoclinic orbit defines the basin of attraction for the spiral sink [69, 70]. Thus, any trajectory initiated outside the homoclinic orbit will result in an uncontrollable tumor growth while the concentration of effector cells fades out. On the other hand, if is a spiral source, any initialization of the system close to will produce regular oscillations with a slowly increasing amplitude. This behavior persists until the trajectory of the system reaches the unstable manifold of the saddle point which drives the system towards an exponentially uncontrollable tumor evolution (see Fig. 3(B)).

Figure 3: Classification of the system trajectories. (A) Homoclinic orbit of the saddle point (star) when the spiral point (dot) is stable. (B) The unstable manifold of the saddle point (star) when the spiral point (dot) is unstable

Interestingly, the time required for reaching the unstable manifold is significantly large. This system’s behavior can be explained by performing a local stability analysis at the spiral point . Fig. 4(A,B) illustrates the stability of when is kept fixed and is allowed to vary. The eigenvalues of the Jacobian estimated at are . When , is a spiral sink, while for , is a spiral source. Fig. 4(C,D) shows how changes with respect to parameter . For all values of the parameters considered, we find that when is a spiral source the real part of the eigenvalues is with a maximum value of . Therefore, by setting , we can identify two different time-scales describing the system’s behavior: a fast time scale where the regular oscillations occur and a slow time scale which describes the slowly varying amplitude. Fig. 5(A) depicts the variation of , which determines the oscillations frequency, with respect to the parameters and . Notice that , i.e. .

Figure 4: Stability analysis of the spiral point with respect to parameter . (A,B) Stability analysis of with respect to the functional tumor-associated vasculature for the immune recruitment rate equal to 0.54 and 0.60 days, respectively. The labels and represent the bifurcation points, while the solid and dashed lines depict the stable and unstable solution branches, respectively. (C) Real part of the eigenvalues of the Jacobian estimated at the spiral point with respect to for days. (D) Real part of the eigenvalues of the Jacobian estimated at the spiral point with respect to for days.

Figure 5: Frequency dependence on model parameters and . Imaginary part of the eigenvalues of the Jacobian estimated at the spiral point with respect to the immune recruitment rate and the functional tumor-associated vasculature . The labels and represent the regimes where the spiral point is stable and unstable, respectively.

A plausible question relates to the possibility of controlling tumor growth and with the conditions under which this control can be achieved. The local stability analysis reveals that, when is a spiral sink, the homoclinic orbit creates a trapping region where the tumor remains controlled. However, when is a spiral source, the tumor will finally escape innate immune control, albeit the fact that it will stay to a certain region for a long period of time. Therefore, a potential strategy aiming at controlling the tumor growth would be to keep the tumor near the fixed point. In the next sections, we show how an external modulation of the immune system dynamics could be designed to limit the uncontrolled tumor growth in the case that is a spiral source.

Mapping to a negatively damped harmonic oscillator via multiscale analysis

The original model given by Eqs. (1)-(2) cannot be solved analytically and no control strategy can be easily designed. For this reason, we employ a multiscale approach to analyze the system’s behavior near the spiral point . This approach, that efficiently describes the dynamics near a Hopf bifurcation point, reveals the inherent multiscale structure of the problem by capturing the regular oscillations and constructing an envelope of the slowly varying amplitude in deterministic [73] or stochastic nonlinear systems [74, 75, 76]. Thus, adopting such an approach we are able to map the initial system to a simplified negatively damped harmonic oscillator which allows to describe the self-sustained oscillations and the system’s energy gain [77].

The system’s behavior near the spiral source can be described by linearizing the Eqs. (1)-(2) around and . Then, we obtain a new system of equations approximating the evolution of the perturbations and given by:


for and . Then, the linearized system takes the following form:


where represents the Jacobian at the fixed point and correspond to the right-hand side of Eqs. (1)-(2), respectively. The eigenvalues of the Jacobian are . Moreover, and are the trace and determinant of the Jacobian matrix, respectively. Since and , it follows that , , whereas and . Moreover, we observe that which, consequently, results in . Thus, the eigenvalues can be approximated by , where . The condition holds for each case where is a spiral source and is always negative since it represents the death rate of tumor cells. Therefore, we can approximate the linearized system (4) by the following one which incorporates the different order-terms:


where, and . The approximate linearized system (5) explicitly captures the different-order terms and allows to draw analogies with the field of Mechanics. More precisely, the system (5) represents a perturbed harmonic oscillator. The unperturbed problem (i.e. when ) is a linear harmonic oscillator with frequency . The higher-order terms (or correction terms) insert small perturbations which result in a weakly negatively damped harmonic oscillator.

The solution of the approximate system (5) can be expressed as:


where and contain the information regarding the slowly varying amplitude depending on , where for . The multiscale assumption is that the functions and evolve at the slow time scale and are considered to be constant with respect to the oscillations with frequency , evolving on the fast time scale . Notice that the variables and are considered to be independent. The constant term in the expression of has been used to simplify the upcoming calculations.

In order to estimate the functions and , we assume that they follow evolution equations of the form:


Then, by taking the differentials of the system of equations in (6), we have that:


Notice that the system of equations in (5) and (8) are equivalent. Therefore, by equating the terms which are of the same order, we obtain that:


To estimate the functions and , we project the Eqs. (9)-(10) onto the fast dynamics. This is an important step of the calculations in order to isolate the amplitude of functions and [76, 75]. For instance, multiplying Eq. (9) by and integrating over a period of time , we have that:

which results in . In a similar way, multiplying Eq. (9) by and following the same procedure, we find . Consequently, the solution of the approximate linearized system (5) is given by:


where the positive parameters and are defined by the initial conditions of the original model (1)-(2). In addition, Eq. (11) can be further simplified to take the following form:


where , , and . We show that the functions and are orthogonal since:


Then, the system of the approximate solutions becomes:


The advantage of using the proposed approach is that the solutions of the approximate model (16) depend directly on experimentally accessible parameters and the initial conditions. We focus now on the stability properties of the system (5) with respect to the time evolution of the variables and . To that end, we construct a Lyapounov functional as:


where is always positive definite since , and the time-derivative of is given by:


Since the system gains energy. The average energy over a period is estimated by using the form of the approximate solutions (12) and considering the multiscale assumption:


The term is constant and depends on time-invariant parameter values and the initial conditions. Consequently, the average energy gain rate over a period is equal to:


The previous results suggest that a therapeutic strategy, which influences the effector cell dynamics, should be designed in a way to “pump out energy” with an average rate which is greater than the system’s gain energy. In particular, the per period gain rate is equal to according to relation (19). Thus, if we introduce an external immune-modulatory term to intervene in the tumor-effector cell interactions, the system of equations (5) becomes:


According to Eq. (18), the rate is at most of the order of . Therefore, the function should be of the order of , for instance:


In order to calculate the function , we require that the energy of the system (21) should have a negative rate, that is:


The negative energy rate means that a trajectory will flow “downhill” towards a stable fixed point. In the limit , we find an expression for the zero-rate energy function in terms of the solutions and , given by:


Substituting the solutions of and from Eq. (16), we have an analytical expression of the zero-rate energy function , that is:


The function has singularities at the points where the function is equal to zero, i.e. , . However, the effect of , i.e. the integral of the function within a small time interval , is finite.

Fig. 6 illustrates the effect of by integrating between and , for values of the model parameters and where is unstable. In each case, the system was initialized near the spiral point .

Figure 6: Dependence of the zero-rate energy function effect on model parameters and . Dependence of the effect of on the immune recruitment rate and functional tumor-associated vasculature in the unstable regime of denoted by . The label stands for the regime where the spiral point is stable and tumor control is reached.

Using the function is not feasible for practical purposes. However, we can design meaningful therapy functions that induce the same or greater effects in a certain time interval [, ] as the function , that is:


This will result in “pumping out energy” at a rate greater than the system’s energy gain. In the next section, we estimate the function which represents the zero-rate energy scenario, as well as we present numerical simulation results for specific values of the model parameters. Furthermore, based on the behavior of the function , we suggest an efficient external immune-modulatory term to fulfill the relation (26) which results in long-term tumor control.

Results: theory-driven therapeutic design

In this section, we design an immuno-therapeutic proposal derived from our model analysis. We first compare the approximate solutions with those obtained from the original model (1)-(2). Then, we show how the defined energy function in Eq. (17) for the approximate linearized system can be used to describe the system’s energy gain. Finally, we design an external immuno-modulatory strategy following the behavior of the zero-rate energy function given in Eq. (25).

To illustrate the simulation results and without loss of generality, we select the following parameter values: days and , which result in and rad/days. Notice that similar results can be obtained for each parameter set where is a spiral source, since always holds. For numerical integration, we consider a order Runge-Kutta method. The time step was set equal to days which is sufficiently small.

Validity of approximate solution

At this point, we need to validate the performance of the approximate solutions. Fig. 7(A) depicts the evolution of the approximate perturbation compared to the perturbation of the initial model given by Eqs. (1)-(2), referred in what follows as original model perturbation. The approximate solutions were estimated by performing a numerical integration of the system (5), while by numerically integrating the original model (1)-(2). Both systems of equations were initialized close to the spiral point and the simulation time was set days. A quick glance at Fig. 7(A) reveals that the solution derived by the approximate model (5) is very close to the original one, which demonstrates the validity of the approach. Moreover, Fig. 7(A) shows a comparison of the original and approximate model solutions by zooming in the time interval between days 400 and 425. In this interval, the maximum error between such solutions was found to be of the order of .

Figure 7: Validation of approximate solutions. (A) Evolution of the approximate perturbation compared to the original model perturbation . (B) Evolution of the energy function of the approximate linearized model defined in Eq. (17) against the energy estimated by using the original model (1)-(2). (C) The zero-rate energy function in Eq. (24) estimated by using the approximate solution of the system (13). (D) The immuno-modulatory function for and compared with the zero-rate energy function .

Fig. 7(B) compares the evolution of the approximate energy function defined in Eq. (17) against the energy of the original model (1,2). As in the previous case, the energy function of the approximation is in a good agreement with that obtained from the original model. Therefore, we expect that an external immune-modulatory strategy, which results in a negative energy rate, can be used with the same efficacy to the original model.

Construction of the therapy term

Fig. 7(C) depicts the zero-rate energy function in Eq. (24) estimated by using the approximate solutions and from the system (21). As we observe, this function diverges at the singularity points, which means that it cannot be used directly to induce tumor control.

It should be noted that, the zero-rate energy function changes sign according to the evolution of the variable , i.e. the function follows the monotonicity of tumor evolution. Hence, an external immuno-modulatory therapy that satisfies the condition (26) should follow the evolution of the tumor by increasing the recruitment rate of effector cells when the tumor radius increases and decreasing the immune recruitment rate when the tumor radius decreases. A simple approximation of the therapy function would be a step function having a constant value and the same sign to the zero-rate energy function at the same time interval:


where expresses the sigmoidal steepness, and the selection of should fulfill the relation (26). Fig. 7(D) depicts the immuno-modulatory function for and , compared with the zero-rate energy function estimated.

Performance of the suggested therapy

The application of the external therapy in Eq. (27) implies the addition of a new term in Eqs. (1)-(2), that is:


where .

Fig. 8(A,B) respectively shows the evolution of the tumor radius and effector cells by numerically integrating Eqs. (28)-(29) which include the proposed external immuno-modulatory function (27). The system is initialized near the spiral point and parameter values of the external function were selected equal to and . In doing so, we obtain that the system reaches a stable steady state, i.e. the tumor remains controlled. Fig. 8(C) illustrates how the energy of the system of Eqs. (28)-(29) evolves. The energy decreases with time and becomes zero as the system reaches the steady state. It is worth pointing out that was selected not only to satisfy the relation (26), but also to result in the system’s fast convergence to a steady state.

Figure 8: Performance of the proposed therapy term. Evolution of the tumor radius (A) and concentration of effector cells (B) by considering the external immuno-modulatory function (27) with parameters and . (C) Energy of the system of Eqs. (28)-(29) with and . Evolution of the tumor radius (D) and concentration of effector cells (E) by initializing the system away from the spiral point and considering the external immuno-modulatory function with and . (F) Energy of the system of Eqs. (28)-(29), with and , initialized away from the spiral point .

Fig. 8(D,E) respectively shows the evolution of the tumor radius and concentration of effector cells by numerically integrating Eqs. (28)-(29) when the system is initialized away from the spiral point . In this case, the approximate solution is expected to deviate from the solutions of the original model. Fig. 8(F) represents the temporal evolution of the corresponding energy. The parameter was selected to be equal to to fulfil the relation (26), as well as to provide a fast convergence to a steady state. Interestingly enough, in this case the system also reaches a steady state, even though the approximate system is not accurate enough. Consequently, the proposed external immuno-modulatory function is shown to be adequate in controlling tumor growth, even when the system is initialized far to the spiral point .


In this article, we investigate the therapeutic potential of immunomodulatory interventions against tumor growth. To that end, we consider a model introduced in [53] that describes the dynamic interplay between vascularized tumor growth and effector cell responses. Our goal is to identify an external modification of effector cell dynamics that allows for controlling tumor growth. With the help of bifurcation analysis, we identified a unstable oscillatory regime that induces tumor evasion from immuno-surveillance. The characteristic feature of this regime is the oscillations occur at a faster time scale than the amplitude dynamics. Exploiting this time scale separation and via temporal multiscale analysis, we map our model onto a weakly negatively damped harmonic oscillator. This approximation allowed us to identify an analytical expression for the additive control term to the effector cell dynamics. This term acts as an external immunomodulatory therapy where the numerical simulations evidence its efficacy in controlling tumor growth.

The crucial question concerns the translational potential of our theory-driven therapeutic proposal into clinical practice. Our suggested immunomodulatory strategy is relevant to small enough, non-invasive tumors that are initially controlled by the immune system but eventually evade immuno-surveillance. The latter occurs when tumors exceed a critical size where immune responses are unable to confer any control. As stated above, our model suggests that such tumor evasion may take place in the form of oscillations of slowly increasing amplitude. The proposed therapeutic strategy is based on the synchronization of immuno-stimulating and -suppressive phases with tumor growth dynamics. Although immuno-stimulating therapies seem to be expected and plausible [6, 8], immuno-suppression sounds counter-intuitive and dangerous. However, the latter occurs in clinical practice during chemotherapeutic interventions [78]. Therefore, a potential realization of our proposed strategy could be mediated by a combination of state-of-the-art immunotherapies [4, 1, 49, 8] and chemotherapeutic modules. The latter not only will play the role of immuno-suppressor, but also will slow down tumor growth dynamics. Needless to state that such a therapeutic suggestion requires experimental validation and further investigation.

Finally, we conclude by pointing out the limitations of the present work. Paramount among them, introducing vascularization dynamics should be not only a natural, but also an insightful extension of the proposed model. At this state tumor vascularization is considered as a constant parameter in time. Thus, making vascularization dynamic would make sense to model further hypoxic effects, such as necrosis or hypoxic-induced invasion. Moreover, the immune system is much more complicated than its description in the current model, involving much more cell types and interactions. In particular, immune system is often regarded as acting to suppress tumor growth, however it is now clear that it can be both stimulatory and inhibitory, as in the case of tumor-associated macrophages [79, 80]. Moreover, our model needs be extended for invasive tumors. Nevertheless, mathematical modeling can help to improve our understanding of the interplay between the several competing factors that have complex implications for cancer therapy. We hope that the model presented here could provide a step towards this direction.


A. I. Reppas, J. C. L. Alfonso and H. Hatzikirou gratefully acknowledge the support of the German Excellence Initiative via the Cluster of Excellence EXC 1056 Center for Advancing Electronics Dresden (cfAED) at the Technische Universität Dresden. J. C. L. Alfonso and H. Hatzikirou also acknowledge the German Federal Ministry of Education and Research (BMBF) for the eMED project SYSIMIT (01ZX1308D).


  • [1] O. Finn, “Immuno-oncology: understanding the function and dysfunction of the immune system in cancer,” Ann. Oncol., vol. 23, no. 8, pp. 6–9, 2012.
  • [2] G. Dranoff, “Cytokines in cancer pathogenesis and cancer therapy,” Nat. Rev. Cancer, vol. 4, no. 1, pp. 11–22, 2004.
  • [3] T. Stewart and S. Abrams, “How tumours escape mass destruction,” Oncogene, vol. 27, no. 45, pp. 5894–903, 2008.
  • [4] W. Zou, “Immunosuppressive networks in the tumour environment and their therapeutic relevance,” Nat. Rev. Cancer, vol. 5, no. 4, pp. 263–74, 2005.
  • [5] G. Rabinovich, D. Gabrilovich, and E. Sotomayor, “Immunosuppressive strategies that are mediated by tumor cells,” Annu. Rev. Immunol., vol. 25, pp. 267–96, 2007.
  • [6] C. Ghirelli and T. Hagemann, “Targeting immunosuppression for cancer therapy,” J. Clin. Invest., vol. 123, no. 6, pp. 2355–7, 2013.
  • [7] M. Cheever, J. Schlom, L. Weiner, H. Lyerly, M. Disis, A. Greenwood, O. Grad, and W. Nelson, “Translational research working group developmental pathway for immune response modifiers,” Clin. Cancer Res., vol. 14, no. 18, pp. 5692–9, 2008.
  • [8] S. Rosenberg, “Decade in review-cancer immunotherapy: entering the mainstream of cancer treatment,” Nat. Rev. Clin. Oncol., vol. 11, no. 11, pp. 630–2, 2014.
  • [9] O. Finn, “Cancer immunology,” N. Engl. J. Med., vol. 358, no. 25, pp. 2704–15, 2008.
  • [10] D. Ambrosi and F. Mollica, “On the mechanics of a growing tumor,” International journal of engineering science, vol. 40, no. 12, pp. 1297–1316, 2002.
  • [11] H. Byrne and L. Preziosi, “Modelling solid tumour growth using the theory of mixtures,” Mathematical Medicine and Biology, vol. 20, no. 4, pp. 341–366, 2003.
  • [12] R. K. Jain, J. D. Martin, and T. Stylianopoulos, “The role of mechanical forces in tumor growth and therapy,” Annual review of biomedical engineering, vol. 16, p. 321, 2014.
  • [13] A. Ramírez-Torres, R. Rodríguez-Ramos, J. Merodio, J. Bravo-Castillero, R. Guinovart-Díaz, and J. Alfonso, “Action of body forces in tumor growth,” International Journal of Engineering Science, vol. 89, pp. 18–34, 2015.
  • [14] A. Ramírez-Torres, R. Rodríguez-Ramos, J. Merodio, J. Bravo-Castillero, R. Guinovart-Díaz, and J. Alfonso, “Mathematical modeling of anisotropic avascular tumor growth,” Mechanics Research Communications, vol. 69, pp. 8–14, 2015.
  • [15] J. Arciero, T. Jackson, and D. Kirschner, “A mathematical model of tumor-immune evasion and sirna treatment,” Discrete Continuous Dyn Syst Ser B, vol. 4, pp. 39–58, Feb. 2004.
  • [16] N. Bellomo, B. Firmani, and L. Guerri, “Bifurcation analysis for a nonlinear system of integro-differential equations modelling tumor-immune cells competition,” Appl Math Lett, vol. 12, pp. 39–44, Mar. 1999.
  • [17] N. Bellomo and N. Preziosi, “Modelling and mathematical problems related to tumor evolution and its interaction with the immune system,” Math. Comput. Modelling, vol. 32, pp. 413–452, Aug. 2000.
  • [18] M. Galach, “Dynamics of the tumor-immune system competition - the effect of time delay,” Int J Appl Math Comput Sci, vol. 13, no. 3, pp. 395–406, 2003.
  • [19] A. Matzavinos, M. Chaplain, and V. Kuznetsov, “Mathematical modelling of the spatio-temporal response of cytotoxic t-lymphocytes to a solid tumour,” Math Med Biol, vol. 21, pp. 1–34, Mar. 2004.
  • [20] A. d’Onofrio, “A general framework for modeling tumor-immune system competition and immunotherapy: Mathematical analysis and biomedical inferences,” Physica D: Nonlinear Phenomena, vol. 208, pp. 220–235, Sept. 2005.
  • [21] D. Mallet and L. De Pillis, “A cellular automata model of tumor-immune system interactions,” J Theor Biol, vol. 239, pp. 334–350, Apr. 2006.
  • [22] M. Al-Tameemi, M. Chaplain, and A. d’Onofrio, “Evasion of tumours from the control of the immune system: consequences of brief encounters,” Biol Direct, vol. 7, p. 31, Sept. 2012.
  • [23] M. Robertson-Tessi, A. El-Kareh, and A. Goriely, “A mathematical model of tumor-immune interactions,” J Theor Biol, vol. 294, pp. 56–73, Feb. 2012.
  • [24] F. Rihan, D. Abdel Rahman, S. Lakshmanan, and A. Alkhajeh, “A time delay model of tumour–immune system interactions: Global dynamics, parameter estimation, sensitivity analysis,” Appl Math Comput, vol. 232, pp. 606–623, Apr. 2014.
  • [25] C. Koebel, W. Vermi, J. Swann, N. Zerafa, S. Rodig, L. Old, M. Smyth, and R. Schreiber, “Adaptive immunity maintains occult cancer in an equilibrium state,” Nature, vol. 450, pp. 903–907, 2007.
  • [26] A. d’Onofrio, “Tumor evasion from immune control: Strategies of a miss to become a mass,” Chaos Soliton. Fract., vol. 31, no. 2, pp. 261–268, 2007.
  • [27] G. Caravagna, A. d’Onofrio, P. Milazzo, and R. Barbuti, “Tumour suppression by immune system through stochastic oscillations,” J Theor Biol, vol. 265, pp. 336–345, Aug. 2010.
  • [28] D. Pardoll, “Does the immune system see tumors as foreign or self?,” Annu. Rev. Immunol., vol. 21, pp. 807–39, 2003.
  • [29] G. Dunn, L. Old, and R. Schreiber, “The three es of cancer immunoediting,” Annu. Rev. Immunol., vol. 22, pp. 329–60, 2004.
  • [30] J. Stark, C. Chan, and A. George, “Oscillations in the immune system,” Immunol Rev, vol. 216, pp. 213–231, Apr. 2007.
  • [31] B. Coventry, M. Ashdown, M. Quinn, S. Markovic, S. Yatomi-Clarke, and A. Robinson, “Crp identifies homeostatic immune oscillations in cancer patients: a potential treatment targeting tool?,” J. Transl. Med., vol. 7, p. 102, 2009.
  • [32] V. Kuznetsov, I. Makalkin, M. Taylor, and A. Perelson, “Nonlinear dynamics of immunogenic tumors: Parameter estimation and global bifurcation analysis,” Bull Math Biol, vol. 56, pp. 295–321, Mar. 1994.
  • [33] D. Kirschner and J. Panetta, “Modeling immunotherapy of the tumor-immune interaction,” J Math Biol, vol. 37, pp. 235–52, Sept. 1998.
  • [34] A. d’Onofrio, “Tumour-immune system interaction: Modeling the tumour-stimulated proliferation of effectors and immunotherapy,” Math Models Methods Appl Sci, vol. 16, pp. 1375–1401, Aug. 2006.
  • [35] L. de Pillis, W. Gu, and A. Radunskaya, “Mixed immunotherapy and chemotherapy of tumors: modeling, applications and biological interpretations,” J Theor Biol, vol. 238, pp. 841–862, Feb. 2006.
  • [36] A. d’Onofrio, F. Gatti, P. Cerrai, and L. Freschi, “Delay-induced oscillatory dynamics of tumour–immune system interaction,” Math Comput Model, vol. 51, pp. 572–591, Mar. 2010.
  • [37] R. Araujo and D. McElwain, “A history of the study of solid tumour growth: the contribution of mathematical modelling,” Bull Math Biol, vol. 66, pp. 1039–1091, Sept. 2004.
  • [38] J. Nagy, “The ecology and evolutionary biology of cancer: a review of mathematical models of necrosis and tumor cell diversity,” Math Biosci Eng, vol. 2, pp. 381–418, Apr. 2005.
  • [39] H. Byrne, T. Alarcon, M. Owen, S. Webb, and P. Maini, “Modelling aspects of cancer dynamics: a review,” Philos Trans A Math Phys Eng Sci, vol. 364, pp. 1563–1578, June 2006.
  • [40] T. Roose, S. Chapman, and P. Maini, “Mathematical models of avascular tumor growth,” SIAM Rev, vol. 49, pp. 179–208, May 2007.
  • [41] M. Martins, S. Ferreira Jr., and M. Vilela, “Multiscale models for the growth of avascular tumors,” Phys Life Rev, vol. 4, pp. 128–156, June 2007.
  • [42] N. Bellomo and M. Delitala, “From the mathematical kinetic, and stochastic game theory to modelling mutations, onset, progression and immune competition of cancer cells,” Phys Life Rev, vol. 5, pp. 183–206, Dec. 2008.
  • [43] M. Chaplain, “Modelling aspects of cancer growth: Insight from mathematical and numerical analysis and computational simulation,” in Multiscale Problems in the Life Sciences (V. Capasso and M. Lachowicz, eds.), vol. 1940 of Lecture Notes in Mathematics, pp. 147–200, Springer Berlin Heidelberg, 2008.
  • [44] R. Eftimie, J. Bramson, and D. Earn, “Interactions between the immune system and cancer: a brief review of non-spatial mathematical models,” Bull Math Biol, vol. 73, pp. 2–32, Jan. 2011.
  • [45] K. Wilkie, “A review of mathematical models of cancer-immune interactions in the context of tumor dormancy,” Adv Exp Med Biol, vol. 734, pp. 201–234, Aug. 2013.
  • [46] E. Ruoslahti, “Specialization of tumour vasculature,” Nat Rev Cancer, vol. 2, pp. 83–90, Feb. 2002.
  • [47] N. Ferrara and R. Kerbel, “Angiogenesis as a therapeutic target,” Nature, vol. 438, pp. 967–974, Dec. 2005.
  • [48] R. Farnsworth, M. Lackmann, M. Achen, and S. Stacker, “Vascular remodeling in cancer,” Oncogene, vol. 33, pp. 3496–3505, Aug. 2014.
  • [49] W. Fridman, F. Pagès, C. Sautès-Fridman, and J. Galon, “The immune contexture in human tumours: impact on clinical outcome,” Nat Rev Cancer, vol. 12, pp. 298–306, Mar. 2012.
  • [50] M. Junttila and F. de Sauvage, “Influence of tumour micro-environment heterogeneity on therapeutic response,” Nature, vol. 501, pp. 346–354, Sept. 2013.
  • [51] H. Greenspan, “On the growth and stability of cell cultures and solid tumors,” J Theor Biol, vol. 56, pp. 229–242, Jan. 1976.
  • [52] V. Cristini, J. Lowengrub, and Q. Nie, “Nonlinear simulation of tumor growth,” J Math Biol, vol. 46, pp. 191–224, Mar. 2003.
  • [53] H. Hatzikirou, J.C.L. Alfonso, S. Mühle, C. Stern, S. Weiss, and M. Meyer-Hermann, “Cancer therapeutic potential of combinatorial immuno- and vaso-modulatory interventions,” arXiv:1505.05670v1, 2015.
  • [54] H. Hatzikirou, A. Chauvière, J. Lowengrub, J. De Groot, and V. Cristini, “Effect of vascularization on glioma tumor growth,” in Modeling Tumor Vasculature (T. L. Jackson, ed.), pp. 237–259, Springer New York, 2012.
  • [55] L. de Pillis, A. Radunskaya, and C. Wiseman, “A validated mathematical model of cell-mediated immune response to tumor growth,” Cancer Res, vol. 65, pp. 7950–8, Sept. 2005.
  • [56] Y. Huang, S. Goel, D. Duda, D. Fukumura, and R. Jain, “Vascular normalization as an emerging strategy to enhance cancer immunotherapy,” Cancer Res, vol. 73, pp. 2943–2948, May 2013.
  • [57] A. Herman, V. Savage, and G. West, “A quantitative theory of solid tumor growth, metabolic rate and vascularization,” PLoS ONE, vol. 6, no. 9, p. e22973, 2011.
  • [58] A. Brú, S. Albertos, J. Subiza, J. García-Asenjo, and I. Brú, “The universal dynamics of tumor growth,” Biophys J, vol. 85, pp. 2948–61, Nov. 2003.
  • [59] A. Brú, S. Albertos, J. García-Asenjo, and I. Brú, “Pinning of tumoral growth by enhancement of the immune response,” Phys Rev Lett, vol. 92, pp. 1–4, June 2004.
  • [60] K. Alessandri, B. Sarangi, V. Gurchenkov, B. Sinha, T. Kieling, L. Fetler, F. Rico, S. Scheuring, C. Lamaze, A. Simon, S. Geraldo, D. Vignjevic, H. Doméjean, L. Rolland, A. Funfak, J. Bibette, N. Bremond, and P. Nassoy, “Cellular capsules as a tool for multicellular spheroid production and for investigating the mechanics of tumor progression in vitro,” Proc Natl Acad Sci U S A, vol. 110, pp. 14843–14848, Aug. 2013.
  • [61] M. Delarue, F. Montel, O. Caen, J. Elgeti, J. Siaugue, D. Vignjevic, J. Prost, J. Joanny, and G. Cappello, “Mechanical control of cell flow in multicellular spheroids,” Phys Rev Lett, vol. 110, p. 138103, Mar. 2013.
  • [62] H. Acker, Spheroids in cancer research: methods and perspectives. Springer-Verlag Berlin; New York, 1984.
  • [63] H. Frieboes, X. Zheng, C. Sun, B. Tromberg, R. Gatenby, and V. Cristini, “An integrated computational/experimental model of tumor invasion,” Cancer Res, vol. 66, pp. 1597–1604, Feb. 2006.
  • [64] V. Cristini, H. Frieboes, X. Li, J. Lowengrub, P. Macklin, S. Sanga, S. Wise, and X. Zheng, “Nonlinear modeling and simulation of tumor growth,” in Selected topics in cancer modeling: Genesis, evolution, immune competition, and therapy. Modelling and Simulation in Science, Engineering, and Technology (N. Bellomo, M. Chaplain, and E. de Angelis, eds.), ch. 6, pp. 113–82, Boston, MA USA: Birkhäuser, 2008.
  • [65] B. Su, W. Zhou, K. Dorman, and D. Jones, “Mathematical modelling of immune response in tissues,” Comput Math Methods Med, vol. 10, no. 1, pp. 9–38, 2009.
  • [66] A. d’Onofrio, U. Ledzewicz, and H. Schättler, “On the dynamics of tumor-immune system interactions and combined chemo- and immunotherapy,” in New Challenges for Cancer Systems Biomedicine (A. d’Onofrio, P. Cerrai, and A. Gandolfi, eds.), SIMAI Springer Series, pp. 249–266, Springer Milan, 2012.
  • [67] C. Kelley, Iterative methods for optimization, vol. 18. SIAM, 1999.
  • [68] H. B. Keller, “Numerical solution of bifurcation and nonlinear eigenvalue problems,” Applications of bifurcation theory, pp. 359–384, 1977.
  • [69] A. Nayfeh and B. Balachandran, Applied nonlinear dynamics: analytical, computational and experimental methods. John Wiley and Sons, 2008.
  • [70] S. Strogatz, Nonlinear dynamics and chaos: with applications to physics, biology, chemistry and engineering. Perseus Books Group, 2001.
  • [71] W. E. Boyce, R. C. DiPrima, and C. W. Haines, Elementary differential equations and boundary value problems, vol. 9. Wiley New York, 1992.
  • [72] G. Iooss and D. D. Joseph, Elementary stability and bifurcation theory. Springer Science & Business Media, 2012.
  • [73] J. Cole and J. Kevorkian, “Multiple scale and singular perturbation methods,” Appl. Math. Sci, vol. 114, 1996.
  • [74] M. Klosek and R. Kuske, “Multiscale analysis of stochastic delay differential equations,” Multiscale Model Simul., vol. 3, no. 3, pp. 706–729, 2005.
  • [75] R. Kuske, “Multi-scale analysis of noise-sensitivity near a bifurcation,” in IUTAM Symposium on Nonlinear Stochastic Dynamics, pp. 147–156, Springer Netherlands, 2003.
  • [76] R. Kuske, L. Gordillo, and P. Greenwood, “Sustained oscillations via coherence resonance in sir,” J Theor Biol, vol. 245, no. 3, pp. 459–469, 2007.
  • [77] A. Jenkins, “Self-oscillation,” Physics Reports, vol. 525, no. 2, pp. 167–222, 2013.
  • [78] L. Zitvogel, L. Apetoh, F. Ghiringhelli, and G. Kroemer, “Immunological aspects of cancer chemotherapy,” Nat Rev Immunol, vol. 8, no. 1, pp. 59–73, 2008.
  • [79] K. de Visser, A. Eichten, and L. Coussens, “Paradoxical roles of the immune system during cancer development,” Nat Rev Cancer, vol. 6, pp. 24–37, Jan. 2006.
  • [80] A. Mantovani, P. Romero, A. Palucka, and F. Marincola, “Tumour immunity: effector response to tumour and role of the microenvironment,” Lancet, vol. 371, pp. 771–783, Mar. 2008.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description