In silico tumor control induced via alternating immunostimulating and immunosuppressive phases

Abstract

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. Toward 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. In this work, 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.

Keywords:

• tumor-immune system interactions
• immuno-modulatory interventions
• mathematical modeling

Introduction

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.¹ Indeed, every known innate and adaptive immune effector mechanism has been reported that participates in tumor recognition and rejection.² However, experimental observations support that tumorigenic processes by 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-apoptotic factors.⁵ Nevertheless, there are clinical and preclinical evidences supporting that activation of the innate antitumor immunity can result in tumor regression and provide therapeutic benefits.⁶

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.⁷ Passive and active immunotherapy has been successfully applied to the treatment of a wide variety of human cancers and holds promise of a lifelong cure.^8,9 However, tumor-induced immunosuppression still represents a major obstacle to effective cell-mediated immunity and immunotherapy.^3,5 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-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-24 Clinical data evidence that cancer cells can survive in a undetectable dormant state for extended periods of time,²⁵ which has been also predicted by several models of tumor-immune cell interactions.^{19,20,22,26,27} 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 that could result in tumor regrowth.²⁹ Sustained oscillations by the immune system have been observed both in its healthy state and pathological situations.^27,30,31 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.^27,32-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-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-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 exploring 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.^20,32,51,52 The main feature of our model is the modeling of the interplay between functional tumor-associated vasculature and effector cell dynamics.⁵³ 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 refs.^51,32 respectively, see also refs.^20,52,54,55 The main feature is the low dimensional modeling of the complex interplay between tumor-associated functional vasculature and immune recruitment dynamics.⁵³ This model can be interpreted as the temporal evolution of the average tumor radius, since radially symmetric growth is not a realistic behavior. The system variables are the average tumor radius ${\displaystyle R\left(t\right)}$ and effector cell concentration ${\displaystyle E\left(t\right)}$ in the tumor vicinity.

Radial Tumor Growth, R(t)

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 ${\displaystyle {\lambda }_M}$ at the tumor-host interface. The death rate ${\displaystyle {\lambda }_A}$ is uniform within the tumor and constant in time. Moreover, we assume that the death rate ${\displaystyle {\lambda }_A}$ 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 ${\displaystyle B}$, where ${\displaystyle 0\le B\le 1}$, represents the net effect of functional tumor-associated vasculature on the tumor radius evolution. In the limit of avascular tumor growth ${\displaystyle B=0}$, such tumor-effector interactions take place only at the tumor surface. At the other extreme, for ${\displaystyle B=1}$ effector cells can potentially interact with any cancer cell within the tumor bulk. Moreover, we consider an intrinsic length scale ${\displaystyle L_D}$ 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.⁵⁶ We consider this process through a phenomenological scaling function $f\left(R,B\right)=\frac{R^{B-1}}{R^{B-1}+1}\in \left[0,1\right]$, for ${\displaystyle B\in \left[0,1\right]}$, 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 ${\displaystyle c}$. Such scaling functions have been also considered in the classical von Bertalanffy approach, and more recently in allometric models.⁵⁷

Taking these factors into account, we deduce the tumor radius ${\displaystyle R\left(t\right)}$ dynamic under the assumption of radial symmetry according to the following ODE equation:(1) $\begin{array}{l}\frac{dR}{dt}=\frac{1}{3}\left({\lambda }_{{\rm M}}{\rm B}-{\lambda }_{{\rm A}}\right)R\\ +{\lambda }_{{\rm M}}\left(1-{\rm B}\right)L_D\left(\frac{1}{tanh\left(R/L_D\right)}-\frac{L_D}{R}\right)\\ -cER f\left(R,B\right),\end{array}$(1) where the time coordinate ${\displaystyle t}$ has been omitted for notational simplicity, and ${\displaystyle {\lambda }_M, {\lambda }_A, L_D}$ and ${\displaystyle c}$ are non-negative constants.

Effector Cell Concentration, E(t)

We assume that effectors are recruited at a rate ${\displaystyle r}$ depending on tumor cells following the Michaelis-Menten kinetics, where ${\displaystyle K}$ is a positive constant denoting the concentration at which the immune recruitment is half-maximal. Effector cells die at a rate ${\displaystyle d_0}$ and become inactivated at a rate ${\displaystyle d_1}$ due to their antitumor activity. In particular, the inactivation of effectors by tumor cells is modeled through the function ${\displaystyle f\left(R,B\right)}$ 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 a background rate of immune effector recruitment ${\displaystyle \text{σ}}$, even in the absence of tumor cells.

The resulting ODE equation for the effector cell concentration ${\displaystyle E\left(t\right)}$ dynamic is given by:(2) ${\displaystyle \frac{dE}{dt}=r\frac{R^3}{K+R^3}E-d_{1}E{R}^3 f\left(R,B\right)-d_{0}E+\sigma ,}$(2) where the time coordinate ${\displaystyle t}$ has been omitted for sake of simplicity in the notation, and ${\displaystyle r, K, d_0, d_1}$ and ${\displaystyle \text{σ}}$ 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 ${\displaystyle B}$.^52,54,58,59 Accordingly, we assume that this initial growth depends exclusively on the net proliferation rate ${\displaystyle {\lambda }_p=\left({\lambda }_M-{\lambda }_A\right)}$ in the absence of adaptive antitumor immune responses at those stages of growth, see the first term of Eq. (1(1) $\begin{array}{l}\frac{dR}{dt}=\frac{1}{3}\left({\lambda }_{{\rm M}}{\rm B}-{\lambda }_{{\rm A}}\right)R\\ +{\lambda }_{{\rm M}}\left(1-{\rm B}\right)L_D\left(\frac{1}{tanh\left(R/L_D\right)}-\frac{L_D}{R}\right)\\ -cER f\left(R,B\right),\end{array}$(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 ${\displaystyle {\lambda }_p\approx 1.20 {\text{days}}^{-1}}$,^53,60 and the physiological plausible value ${\displaystyle {\lambda }_M=1/18 h=1.34 {\text{days}}^{-1}}$ for CT26 murine cells,^53,60,61 we have that ${\displaystyle {\lambda }_A\approx 0.14 {\text{days}}^{-1}}$. The characteristic nutrient diffusion length has been experimentally estimated that ranges between 0.2 and 0.3 mm.^53,62-64 The effector cell characteristic concentration is at the order of magnitude ${\displaystyle {10}^5}$ cells. The latter estimate is justified since the characteristic length scale of the system is at the order of ${\displaystyle 1.0\text{ mm}}$, and given that cells are commonly assumed with a diameter between ${\displaystyle 10\text{ μm}}$ and ${\displaystyle 20\text{ μm}}$,⁶¹ then for a volume of ${\displaystyle 1 {\text{mm}}^3}$ the concentration is at the order of magnitude ${\displaystyle {10}^5}$ cells. Moreover, we consider ${\displaystyle c=0.03 {\text{cells}}^{-1} {\text{days}}^{-1}}$ as measured from murine CT26 tumor growth experiments,⁵³ see also refs.^20,32,55 The remaining parameter values ${\displaystyle d_0=0.37 {\text{days}}^{-1}}$, ${\displaystyle d_1=0.01 {\text{mm}}^{-3} {\text{days}}^{-1}}$, ${\displaystyle K=2.72 {\text{mm}}^3 }$ and ${\displaystyle \sigma =0.13 \times {10}^5 \text{cells} {\text{days}}^{-1}}$ are considered from previously reported experimental data,^{20,32,55,65,66} and properly rescaled to the magnitudes and units considered in our mathematical model. Through a parameter sensitivity analysis, we explore the effects on tumor growth of the effector cell recruitment rate ${\displaystyle r}$ and functional vascularization degree ${\displaystyle B}$ for different choices of the initial tumor radius ${\displaystyle R_0}$ and concentration of effector cells ${\displaystyle E_0}$.

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

Table 1. Model parameters. The effects of model parameters B and r, i.e. functional tumor-associated vasculature and effector cell recruitment rate respectively, are investigated by a model sensitivity analysis.

Description	Parameter	Value	Units	Sources
Mitotic rate of tumor cells	${\displaystyle {\lambda }_M}$	1.34	${\displaystyle {\text{days}}^{-1}}$	53,60,61
Death rate of tumor cells	${\displaystyle {\lambda }_A}$	0.14	${\displaystyle {\text{days}}^{-1}}$	53,60,61
Characteristic nutrient diffusion length	${\displaystyle L_D}$	[0.2 – 0.3]	${\displaystyle \text{mm}}$	53,62-64
Killing rate of tumor cells by effectors	${\displaystyle c}$	0.03	${\displaystyle {\text{cells}}^{-1} {\text{days}}^{-1}}$	20,32,53,55
Tumor volume where effector recruitment rate is half-maximal	${\displaystyle K}$	2.72	${\displaystyle {\text{mm}}^3}$	20,32,55,66
Inactivation rate of effectors by tumor cells	${\displaystyle d_1}$	0.01	${\displaystyle {\text{mm}}^{-3} {\text{days}}^{-1}}$	20,32,55,66
Death rate of effectors	${\displaystyle d_0}$	0.37	${\displaystyle {\text{days}}^{-1}}$	20,32,65,66
Background rate of effector recruitment	${\displaystyle \sigma }$	0.13 × 10^5	${\displaystyle \text{cells} {\text{days}}^{-1}}$	20,32,55,66

Dynamical Analysis of the Model

In this section, we analyze the model's behavior with respect to 2 parameters, namely the effector cell recruitment rate ${\displaystyle r}$ and functional tumor-associated vasculature ${\displaystyle B}$.

Fixed Point Analysis

The first step toward analyzing the model dynamics is the identification of the fixed points along with their stability classification. Figure 1A–D depict the phase portraits of the system of equations (1(1) $\begin{array}{l}\frac{dR}{dt}=\frac{1}{3}\left({\lambda }_{{\rm M}}{\rm B}-{\lambda }_{{\rm A}}\right)R\\ +{\lambda }_{{\rm M}}\left(1-{\rm B}\right)L_D\left(\frac{1}{tanh\left(R/L_D\right)}-\frac{L_D}{R}\right)\\ -cER f\left(R,B\right),\end{array}$(1) )–(2(2) ${\displaystyle \frac{dE}{dt}=r\frac{R^3}{K+R^3}E-d_{1}E{R}^3 f\left(R,B\right)-d_{0}E+\sigma ,}$(2) ) for different values of the model parameter ${\displaystyle B}$, while keeping ${\displaystyle r}$ constant and equal to ${\displaystyle {0.57 \text{days}}^{-1}}$. The black curves represent the nullclines, i.e. curves along which ${\displaystyle dR/dt = 0}$ and ${\displaystyle dE/dt = 0}$, 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 2 fixed points corresponding to a low tumor radius ${\displaystyle \left(L=\left(R_L,E_L\right)\right)}$ and a high tumor radius ${\displaystyle \left(H=\left(R_H,E_H\right)\right)}$.

Figure 1. Phase portraits and long-term characteristic dynamics. (A–D) Phase portraits of the tumor radius R versus the concentration of effector cells E for different functional levels of the tumor-associated vasculature B. Trajectories starting with the initial conditions IC₁ = (R₀, E₀) = (2.0 mm, 20×10⁵ cells) and IC₂ = (R₀, E₀) = (2.0 mm, 10×10⁵ cells) are resented. The nullclines (zero-growth isoclines) of the dynamical system are also plotted. (E–H) Temporal evolution of R and E, that corresponds to the trajectories in panels (A–D) starting at the initial condition IC₁. (I–L) Temporal evolution of R and E, that correspond to the trajectories in panels (A–D), starting at the initial condition IC₂. The immune recruitment rate is r = 0.57 days^-1 and the remaining parameters are as in Table 1.

Figure 1E–L show the evolution of the tumor radius and effector cells of the corresponding trajectories on the phase plain presented in Figure 1A–D. A quick glance at Figure 1 reveals that the system trajectories initiated near the ${\displaystyle L}$ fixed point follow an oscillatory behavior with a slowly varying amplitude that can be either increasing Figure 1F, G or decreasing Figure 1E, H. This behavior indicates that, depending on the model parameter values, the ${\displaystyle L}$ 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. 1I–L).

To gain insight in the system behavior near the fixed points we construct the bifurcation diagrams with respect to the immune recruitment rate ${\displaystyle r}$ for different functional levels of vascularization ${\displaystyle B}$ (see Fig. 2A–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

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

For sufficiently small values of model parameter ${\displaystyle r}$, we obtain that there are no fixed points, which consequently implies that the tumor grows indefinitely. However, for a critical tumor radius ${\displaystyle r_{cr}}$, a homoclinic bifurcation occurs giving rise to 2 states: a lower branch that corresponds to the ${\displaystyle L}$ fixed point, and an upper branch that corresponds to the ${\displaystyle H}$ fixed point.^69,70 The local stability analysis shows that the ${\displaystyle H}$ fixed point is a saddle point, and therefore unstable. On the other hand, the ${\displaystyle L}$ fixed point is a spiral point that, depending on the value of the parameter ${\displaystyle B}$, can be unstable (spiral source) or stable (spiral sink).^71,72 Figure 2D, E show the local stability analysis of ${\displaystyle L}$ and ${\displaystyle H}$ for ${\displaystyle B=0.60}$, which corresponds to the bifurcation diagram in Figure 2B.

In the case of ${\displaystyle L}$ being a spiral sink, a system trajectory initiated inside the homoclinic orbit, i.e., the closed orbit that starts and returns to the saddle point in Figure 3A, 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 ${\displaystyle L}$ is a spiral source, any initialization of the system close to ${\displaystyle L}$ 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 ${\displaystyle H}$, which drives the system toward an exponentially uncontrollable tumor evolution (see Fig. 3B).

Figure 3. Classification of the system trajectories. (A) Homoclinic orbit of the saddle point H (star) when the spiral point L (dot) is stable. (B) The unstable manifold of the saddle point H (star) when the spiral point L (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 ${\displaystyle L}$. Figure 4A, B illustrate the stability of ${\displaystyle L}$ when the parameter ${\displaystyle r}$ is kept fixed and vasculature ${\displaystyle B}$ is allowed to vary. The eigenvalues of the Jacobian estimated at ${\displaystyle L}$ are ${\displaystyle \lambda =\mu \pm i\beta }$. When ${\displaystyle \mu <0}$, ${\displaystyle L}$ is a spiral sink, while for ${\displaystyle \mu >0}$, ${\displaystyle L}$ is a spiral source. Figure 4C, D show how ${\displaystyle \mu }$ changes with respect to vasculature ${\displaystyle B}$. For all values of the parameters considered, we find that when ${\displaystyle L}$ is a spiral source the real part of the eigenvalues is ${\displaystyle \mu <<1}$ with a maximum value of ${\displaystyle {\mu }_{max}\cong 5.5 \times {10}^{-3}}$. Therefore, by setting ${\displaystyle \mu =\epsilon <<1}$, we can identify 2 different time-scales describing the system's behavior: a fast time scale ${\displaystyle t}$ where the regular oscillations occur and a slow time scale ${\displaystyle T=\epsilon t}$ which describes the slowly varying amplitude. Figure 5A depicts the variation of ${\displaystyle \beta }$, which determines the oscillations frequency, with respect to the model parameters ${\displaystyle r}$ and ${\displaystyle B}$. Notice that ${\displaystyle \beta \in \left(0,1\right)}$, i.e. ${\displaystyle \beta \propto O\left(1\right)}$.

Figure 4. Stability analysis of the spiral point ${\displaystyle L}$ with respect to vascularization B. (A, B) Stability analysis of L with respect to the functional tumor-associated vasculature B for the immune recruitment rate ${\displaystyle r}$ equal to 0.54 and 0.06 days⁻¹, respectively. The labels S₁ and S₂ 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 L with respect to B for r = 0.54 days⁻¹. (D) Real part of the eigenvalues of the Jacobian estimated at the spiral point L with respect to B for r = 0.60 days⁻¹.

Figure 5. Frequency dependence on the immune recruitment rate ${\displaystyle r}$ and vascularization B. Imaginary part of the eigenvalues of the Jacobian estimated at the spiral point L with respect to the immune recruitment rate r and the functional tumor-associated vasculature B. The labels S and U represent the regimes where the spiral point L 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 ${\displaystyle L}$ is a spiral sink, the homoclinic orbit creates a trapping region where the tumor remains controlled. However, when ${\displaystyle L}$ 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 ${\displaystyle L}$ 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 ${\displaystyle L}$ is a spiral source.

Mapping to a Negatively Damped Harmonic Oscillator via Multiscale Analysis

The original model given by Eqs. (1(1) $\begin{array}{l}\frac{dR}{dt}=\frac{1}{3}\left({\lambda }_{{\rm M}}{\rm B}-{\lambda }_{{\rm A}}\right)R\\ +{\lambda }_{{\rm M}}\left(1-{\rm B}\right)L_D\left(\frac{1}{tanh\left(R/L_D\right)}-\frac{L_D}{R}\right)\\ -cER f\left(R,B\right),\end{array}$(1) )–(2(2) ${\displaystyle \frac{dE}{dt}=r\frac{R^3}{K+R^3}E-d_{1}E{R}^3 f\left(R,B\right)-d_{0}E+\sigma ,}$(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 ${\displaystyle L}$. 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 or stochastic nonlinear systems.^7-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.⁷⁷

The system's behavior near the spiral source ${\displaystyle L}$ can be described by linearizing the Eqs. (1(1) $\begin{array}{l}\frac{dR}{dt}=\frac{1}{3}\left({\lambda }_{{\rm M}}{\rm B}-{\lambda }_{{\rm A}}\right)R\\ +{\lambda }_{{\rm M}}\left(1-{\rm B}\right)L_D\left(\frac{1}{tanh\left(R/L_D\right)}-\frac{L_D}{R}\right)\\ -cER f\left(R,B\right),\end{array}$(1) )–(2(2) ${\displaystyle \frac{dE}{dt}=r\frac{R^3}{K+R^3}E-d_{1}E{R}^3 f\left(R,B\right)-d_{0}E+\sigma ,}$(2) ) around ${\displaystyle R=R_L}$ and ${\displaystyle E=E_L}$. Then, we obtain a new system of equations approximating the evolution of the perturbations ${\displaystyle u}$ and ${\displaystyle v}$ given by:(3) ${\displaystyle u=R-R_L, v=E-E_L,}$(3) for ${\displaystyle \left|u\right|<<1}$ and ${\displaystyle \left|v\right|<<1}$. Then, the linearized system takes the following form:(4) ${\displaystyle \frac{d}{dt}\left(\begin{array}{l}u\\ v\end{array}\right)\approx J{|}_L \left(\begin{array}{l}u\\ v\end{array}\right) = \left[\begin{array}{ll}\frac{\partial F}{\partial R}{|}_L& \frac{\partial F}{\partial E}{|}_L\\ \frac{\partial G}{\partial R}{|}_L& \frac{\partial G}{\partial E}{|}_L\end{array}\right] \left(\begin{array}{l}u\\ v\end{array}\right)= \left[\begin{array}{ll}{a}_1& a_2\\ a_3& a_4\end{array}\right] \left(\begin{array}{l}u\\ v\end{array}\right),}$(4) where ${\displaystyle J{|}_L}$ represents the Jacobian at the fixed point ${\displaystyle L}$ and ${\displaystyle F,G}$ correspond to the right-hand side of Eqs. (1(1) $\begin{array}{l}\frac{dR}{dt}=\frac{1}{3}\left({\lambda }_{{\rm M}}{\rm B}-{\lambda }_{{\rm A}}\right)R\\ +{\lambda }_{{\rm M}}\left(1-{\rm B}\right)L_D\left(\frac{1}{tanh\left(R/L_D\right)}-\frac{L_D}{R}\right)\\ -cER f\left(R,B\right),\end{array}$(1) )–(2(2) ${\displaystyle \frac{dE}{dt}=r\frac{R^3}{K+R^3}E-d_{1}E{R}^3 f\left(R,B\right)-d_{0}E+\sigma ,}$(2) ) respectively. The eigenvalues of the Jacobian are $\text{λ}=\frac{Tr\left(J{|}_L\right)}{2}\pm \frac{\sqrt{Tr\left(J{|}_L\right)-4D}}{2}=\mu \pm i\beta$. Moreover, ${\displaystyle Tr\left(J{|}_L\right)=a_1+a_4}$ and ${\displaystyle D=a_1{a}_4-a_2{a}_3}$ are the trace and determinant of the Jacobian matrix, respectively. Since ${\displaystyle \epsilon =\mu }$ and ${\displaystyle \beta \propto O\left(1\right)}$, it follows that ${\displaystyle a_1\propto O\left(\epsilon \right)}$, ${\displaystyle a_4\propto O\left(\epsilon \right)}$, whereas ${\displaystyle a_2\propto O\left(1\right)}$ and ${\displaystyle a_3\propto O\left(1\right)}$. Moreover, we observe that ${\displaystyle a_1{a}_4\propto O\left({\epsilon }^2\right)}$ which, consequently, results in $\beta =\sqrt{{\left(a_1+a_4\right)}^2-4\left(a_1{a}_4-a_2{a}_3\right)}/2\approx \sqrt{-a_2{a}_3}$. Thus, the eigenvalues can be approximated by ${\displaystyle \lambda =\hat{\lambda }=\epsilon \pm i\omega }$, where ${\displaystyle \omega =\sqrt{-a_2{a}_3}}$. The condition ${\displaystyle \epsilon <<1}$ holds for each case where ${\displaystyle L}$ is a spiral source and ${\displaystyle a_2}$ is always negative since it represents the death rate of tumor cells. Therefore, we can approximate the linearized system (4) by the following one that incorporates the different order-terms:(5) ${\displaystyle \begin{array}{l}\frac{d\hat{u}}{dt}=\epsilon \hat{u}+{\alpha }_2\hat{v},\\ \frac{d\hat{v}}{dt}={\alpha }_3\hat{u}+\epsilon \hat{v},\end{array}}$(5) where ${\displaystyle \hat{u}\approx u}$ and ${\displaystyle \hat{v}\approx v}$. 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 of equations (5(5) ${\displaystyle \begin{array}{l}\frac{d\hat{u}}{dt}=\epsilon \hat{u}+{\alpha }_2\hat{v},\\ \frac{d\hat{v}}{dt}={\alpha }_3\hat{u}+\epsilon \hat{v},\end{array}}$(5) ) represents a perturbed harmonic oscillator. The unperturbed problem (i.e., when ${\displaystyle \epsilon =0}$) is a linear harmonic oscillator with frequency ${\displaystyle \text{ω}}$. 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 of equations (5(5) ${\displaystyle \begin{array}{l}\frac{d\hat{u}}{dt}=\epsilon \hat{u}+{\alpha }_2\hat{v},\\ \frac{d\hat{v}}{dt}={\alpha }_3\hat{u}+\epsilon \hat{v},\end{array}}$(5) ) can be expressed as:(6) $\left(\begin{array}{l}\hat{u}\\ \hat{v}\end{array}\right)=A(T)\left(\begin{array}{l}\frac{a_2}{\omega }\sin (\omega t)\\ \cos (\omega t)\end{array}\right)+B(T)\left(\begin{array}{l}\frac{a_2}{\omega }\cos (\omega t)\\ -\sin (\omega t)\end{array}\right),$(6) where ${\displaystyle A\left(T\right)}$ and ${\displaystyle B\left(T\right)}$ contain the information regarding the slowly varying amplitude depending on ${\displaystyle T}$, where ${\displaystyle T=\epsilon t}$ for ${\displaystyle \epsilon <<1}$. The multiscale assumption is that the functions ${\displaystyle A\left(T\right)}$ and ${\displaystyle B\left(T\right)}$ evolve at the slow time scale ${\displaystyle T}$ and are considered to be constant with respect to the oscillations with frequency${\displaystyle \text{ ω}}$, evolving on the fast time scale ${\displaystyle t}$. Notice that the variables ${\displaystyle T=\epsilon t\propto O\left(\epsilon \right)}$ and ${\displaystyle t\propto O\left(1\right)}$ are considered to be independent. The constant term ${\displaystyle a_2/\omega }$ in the expression of ${\displaystyle \hat{u}}$ has been used to simplify the upcoming calculations.

In order to estimate the functions ${\displaystyle A\left(T\right)}$ and ${\displaystyle B\left(T\right)}$, we assume that they follow evolution equations of the form:(7) ${\displaystyle dA=f_{1}dT, dB=f_{2}dT.}$(7)

Then, by taking the differentials of the system of equations (6(6) $\left(\begin{array}{l}\hat{u}\\ \hat{v}\end{array}\right)=A(T)\left(\begin{array}{l}\frac{a_2}{\omega }\sin (\omega t)\\ \cos (\omega t)\end{array}\right)+B(T)\left(\begin{array}{l}\frac{a_2}{\omega }\cos (\omega t)\\ -\sin (\omega t)\end{array}\right),$(6) ), we have that:(8) ${\displaystyle \begin{array}{lll}d\hat{u}\hfill & =\hfill & \frac{\partial \hat{u}}{\partial t}dt+\frac{\partial \hat{u}}{\partial A}dA+\frac{\partial \hat{u}}{\partial B}dB\hfill \\ \hfill & =\hfill & a_2\hat{v}dt+\left(f_1\frac{a_2}{\omega }\sin \left(\omega t\right)+f_2\frac{a_2}{\omega }cos\left(\omega t\right)\right)dT,\hfill \\ d\hat{v}\hfill & =\hfill & \frac{\partial \hat{v}}{\partial t}dt+\frac{\partial \hat{v}}{\partial A}dA+\frac{\partial \hat{v}}{\partial B}dB\hfill \\ \hfill & =\hfill & a_3\hat{u}dt+\left(f_{1}cos\left(\omega t\right)-f_2\sin \left(\omega t\right)\right)dT.\hfill \end{array}}$(8)

Notice that the system of equations (5(5) ${\displaystyle \begin{array}{l}\frac{d\hat{u}}{dt}=\epsilon \hat{u}+{\alpha }_2\hat{v},\\ \frac{d\hat{v}}{dt}={\alpha }_3\hat{u}+\epsilon \hat{v},\end{array}}$(5) ) and (8(8) ${\displaystyle \begin{array}{lll}d\hat{u}\hfill & =\hfill & \frac{\partial \hat{u}}{\partial t}dt+\frac{\partial \hat{u}}{\partial A}dA+\frac{\partial \hat{u}}{\partial B}dB\hfill \\ \hfill & =\hfill & a_2\hat{v}dt+\left(f_1\frac{a_2}{\omega }\sin \left(\omega t\right)+f_2\frac{a_2}{\omega }cos\left(\omega t\right)\right)dT,\hfill \\ d\hat{v}\hfill & =\hfill & \frac{\partial \hat{v}}{\partial t}dt+\frac{\partial \hat{v}}{\partial A}dA+\frac{\partial \hat{v}}{\partial B}dB\hfill \\ \hfill & =\hfill & a_3\hat{u}dt+\left(f_{1}cos\left(\omega t\right)-f_2\sin \left(\omega t\right)\right)dT.\hfill \end{array}}$(8) ) are equivalent. Therefore, by equating the terms of the same order, we obtain that:(9) $\begin{array}{l}{f}_1\frac{a_2}{\omega }sin\left(\omega t\right)+f_2\frac{a_2}{\omega }cos\left(\omega t\right)\\ =A\left(T\right)\frac{a_2}{\omega }sin\left(\omega t\right)+B\left(T\right)\frac{a_2}{\omega }cos\left(\omega t\right),\end{array}$(9) (10) ${\displaystyle f_{1}cos\left(\omega t\right)-f_{2}sin\left(\omega t\right)=A\left(T\right)cos\left(\omega t\right)-B\left(T\right)sin\left(\omega t\right).}$(10)

To estimate the functions ${\displaystyle f_1}$ and ${\displaystyle f_2}$, we project the Eqs. (9(9) $\begin{array}{l}{f}_1\frac{a_2}{\omega }sin\left(\omega t\right)+f_2\frac{a_2}{\omega }cos\left(\omega t\right)\\ =A\left(T\right)\frac{a_2}{\omega }sin\left(\omega t\right)+B\left(T\right)\frac{a_2}{\omega }cos\left(\omega t\right),\end{array}$(9) )–(10(10) ${\displaystyle f_{1}cos\left(\omega t\right)-f_{2}sin\left(\omega t\right)=A\left(T\right)cos\left(\omega t\right)-B\left(T\right)sin\left(\omega t\right).}$(10) ) onto the fast dynamics. This is an important step of the calculations in order to isolate the amplitude of functions ${\displaystyle f_1}$ and ${\displaystyle f_2}$. ^75,76 For instance, multiplying Eq. (9(9) $\begin{array}{l}{f}_1\frac{a_2}{\omega }sin\left(\omega t\right)+f_2\frac{a_2}{\omega }cos\left(\omega t\right)\\ =A\left(T\right)\frac{a_2}{\omega }sin\left(\omega t\right)+B\left(T\right)\frac{a_2}{\omega }cos\left(\omega t\right),\end{array}$(9) ) by ${\displaystyle \text{sin}\left(\omega t\right)}$ and integrating over a period of time ${\displaystyle \left[0, 2\pi /\omega \right]}$, we have that:${\displaystyle \begin{array}{ll}\underset{0}{\overset{2\pi /\omega }{{\displaystyle \int }}}\left(f_1\frac{a_2}{\omega }sin\left(\omega t\right)+f_2\frac{a_2}{\omega }cos\left(\omega t\right)\right)sin\left(\omega t\right)dt& =\underset{0}{\overset{2\pi /\omega }{{\displaystyle \int }}}\left(A\left(T\right)\frac{a_2}{\omega }sin\left(\omega t\right)+B\left(T\right)\frac{a_2}{\omega }cos\left(\omega t\right)\right)sin\left(\omega t\right)dt,\end{array}}$ which results in ${\displaystyle f_1=A}$. In a similar way, multiplying Eq. (9(9) $\begin{array}{l}{f}_1\frac{a_2}{\omega }sin\left(\omega t\right)+f_2\frac{a_2}{\omega }cos\left(\omega t\right)\\ =A\left(T\right)\frac{a_2}{\omega }sin\left(\omega t\right)+B\left(T\right)\frac{a_2}{\omega }cos\left(\omega t\right),\end{array}$(9) ) by ${\displaystyle \text{cos}\left(\omega t\right)}$ and following the same procedure, we find that ${\displaystyle f_2=B}$. Consequently, the solution of the approximate linearized system of equations (5(5) ${\displaystyle \begin{array}{l}\frac{d\hat{u}}{dt}=\epsilon \hat{u}+{\alpha }_2\hat{v},\\ \frac{d\hat{v}}{dt}={\alpha }_3\hat{u}+\epsilon \hat{v},\end{array}}$(5) ) is given by:(11) $\left(\begin{array}{l}\hat{u}\\ \hat{v}\end{array}\right)=A_0{e}^T\left(\begin{array}{l}\frac{a_2}{\omega }\sin (\omega t)\\ \cos (\omega t)\end{array}\right)+B_0{e}^T\left(\begin{array}{l}\frac{a_2}{\omega }\cos (\omega t)\\ -\sin (\omega t)\end{array}\right),$(11) where the positive parameters ${\displaystyle A_0}$ and ${\displaystyle B_0}$ are defined by the initial conditions of the original model (1)–(2). In addition, Eq. (11(11) $\left(\begin{array}{l}\hat{u}\\ \hat{v}\end{array}\right)=A_0{e}^T\left(\begin{array}{l}\frac{a_2}{\omega }\sin (\omega t)\\ \cos (\omega t)\end{array}\right)+B_0{e}^T\left(\begin{array}{l}\frac{a_2}{\omega }\cos (\omega t)\\ -\sin (\omega t)\end{array}\right),$(11) ) can be further simplified to take the following form:(12) $\left(\begin{array}{l}\hat{u}\\ \hat{v}\end{array}\right)=\left(\begin{array}{l}M{e}^{T}sin(\omega t+\varphi )\\ N{e}^{T}sin(\omega t+\delta )\end{array}\right),$(12) where ${\displaystyle M cos\left(\phi \right)=A_0{a}_2/\omega }$, ${\displaystyle M sin\left(\phi \right)=B_0{a}_2/\omega }$, ${\displaystyle N sin\left(\delta \right)=A_0}$ and ${\displaystyle N cos\left(\delta \right)=-B_0}$. We show that the functions ${\displaystyle \hat{u}}$ and ${\displaystyle \hat{v}}$ are orthogonal since:${\displaystyle \begin{array}{l}\cos \left(\delta -\phi \right)=\cos \left(\delta \right)cos\left(\phi \right)+sin\left(\delta \right)sin\left(\phi \right)\\ =\frac{-B_0}{N}\frac{A_0}{M}\frac{a_2}{\omega }+\frac{A_0}{N}\frac{B_0}{M}\frac{a_2}{\omega }=0\\ \iff \delta =\phi +\frac{\pi }{2}.\end{array}}$ Then, the system of the approximate solutions becomes:(13) $\left(\begin{array}{l}\hat{u}\hfill \\ \hat{v}\hfill \end{array}\right)=\left(\begin{array}{l}M{e}^{T}sin(\omega t+\varphi )\hfill \\ N{e}^{T}cos(\omega t+\varphi )\hfill \end{array}\right),$(13)

The advantage of using the proposed approach is that the solutions of the approximate model (13) 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 ${\displaystyle \hat{u}}$ and ${\displaystyle \hat{v}}$. To that end, we construct a Lyapounov functional as:(14) ${\displaystyle V\left(\hat{u},\hat{v}\right)=\frac{1}{2}\left(a_3{\hat{u}}^2-a_2{\hat{v}}^2\right),}$(14) where ${\displaystyle V}$ is always positive definite since ${\displaystyle a_2<0}$, and the time-derivative of ${\displaystyle V}$ is given by:(15) ${\displaystyle \begin{array}{lll}\frac{dV}{dt}\hfill & =\hfill & a_3\hat{u}\frac{d\hat{u}}{dt}-a_2\hat{v}\frac{d\hat{v}}{dt}\hfill \\ \hfill & =\hfill & \epsilon \left(a_3{\hat{u}}^2-a_2{\hat{v}}^2\right)+\left(a_3{a}_2\hat{u}\hat{v}-a_3{a}_2\hat{u}\hat{v}\right)\hfill \\ \hfill & =\hfill & \epsilon \left(a_3{\hat{u}}^2-a_2{\hat{v}}^2\right)>0.\hfill \end{array}}$(15) Since ${\displaystyle dV/dt>0}$ the system gains energy. The average energy ${\displaystyle ⟨V⟩}$ over a period is estimated by using the form of the approximate solutions (12) and considering the multiscale assumption:(16) $\begin{array}{l}⟨V⟩=\frac{\omega }{2\pi }\underset{0}{\overset{2\pi /\omega }{{{\displaystyle \int }}^{\text{}}}}Vdt=\frac{\omega }{2\pi }\underset{0}{\overset{2\pi /\omega }{{{\displaystyle \int }}^{\text{}}}}\frac{1}{2}\left(a_3{\hat{u}}^2-a_2{\hat{v}}^2\right)\\ \times dt\approx e^{2T}\frac{1}{4}\left(a_3{M}^2-a_2{N}^2\right).\end{array}$(16) The term $\frac{1}{4}\left(a_3{M}^2-a_2{N}^2\right)$ is constant and depends on time-invariant parameter values and the initial model conditions. Consequently, the average energy gain rate over a period is equal to:(17) ${\displaystyle \frac{d⟨V⟩}{dt}=2\epsilon ⟨V⟩.}$(17)

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 ${\displaystyle 2\epsilon }$ according to relation (16). Thus, if we introduce an external immune-modulatory term ${\displaystyle h\left(\hat{v}\right)}$ to intervene in the tumor-effector cell interactions, the system of equations (5(5) ${\displaystyle \begin{array}{l}\frac{d\hat{u}}{dt}=\epsilon \hat{u}+{\alpha }_2\hat{v},\\ \frac{d\hat{v}}{dt}={\alpha }_3\hat{u}+\epsilon \hat{v},\end{array}}$(5) ) becomes:(18) ${\displaystyle \begin{array}{l}\frac{d\hat{u}}{dt}=\epsilon \hat{u}+{\alpha }_2\hat{v},\\ \frac{d\hat{v}}{ dt}=a_3\hat{u}+\epsilon \hat{v}+h\left(\hat{v}\right).\end{array}}$(18) According to Eq. (15(15) ${\displaystyle \begin{array}{lll}\frac{dV}{dt}\hfill & =\hfill & a_3\hat{u}\frac{d\hat{u}}{dt}-a_2\hat{v}\frac{d\hat{v}}{dt}\hfill \\ \hfill & =\hfill & \epsilon \left(a_3{\hat{u}}^2-a_2{\hat{v}}^2\right)+\left(a_3{a}_2\hat{u}\hat{v}-a_3{a}_2\hat{u}\hat{v}\right)\hfill \\ \hfill & =\hfill & \epsilon \left(a_3{\hat{u}}^2-a_2{\hat{v}}^2\right)>0.\hfill \end{array}}$(15) ), the rate ${\displaystyle dV/dt}$ is at most of the order of ${\displaystyle O\left(\epsilon \right)}$. Therefore, the function ${\displaystyle h\left(\hat{v}\right)}$ should be of the order of ${\displaystyle O\left(\epsilon \right)}$, for instance:(19) ${\displaystyle h\left(\hat{v}\right)=\epsilon \hat{z}\left(t\right).}$(19) In order to calculate the function ${\displaystyle \hat{z}\left(t\right)}$, we require that the energy ${\displaystyle V_h}$ of the system (18) should have a negative rate, that is:(20) ${\displaystyle \begin{array}{lll}\frac{d{V}_h}{dt}\hfill & =\hfill & a_3\hat{u}\frac{d\hat{u}}{dt}-a_2\hat{v}\frac{d\hat{v}}{dt}\hfill \\ \hfill & =\hfill & \epsilon \left(a_3{\hat{u}}^2-a_2{\hat{v}}^2-a_2\hat{v}\hat{z}\left(t\right)\right)+\left(a_3{a}_2\hat{u}\hat{v}-a_3{a}_2\hat{u}\hat{v}\right)\hfill \\ \hfill & =\hfill & \epsilon \left(a_3{\hat{u}}^2-a_2{\hat{v}}^2-a_2\hat{v}\hat{z}\left(t\right)\right)\le 0.\hfill \end{array}}$(20)

The negative energy rate means that a trajectory progresses toward a stable fixed point. In the limit ${\displaystyle d{V}_h/dt=0}$, we find an expression for the zero-rate energy function ${\displaystyle \hat{z}\left(t\right)}$ in terms of the solutions ${\displaystyle \hat{u}}$ and ${\displaystyle \hat{v}}$, given by:(21) ${\displaystyle \hat{z}\left(t\right)=\frac{a_3{\hat{u}}^2-a_2{\hat{v}}^2}{a_2\hat{v}}}$(21) Substituting the solutions of ${\displaystyle \hat{u}}$ and ${\displaystyle \hat{v}}$ from Eq. (13(13) $\left(\begin{array}{l}\hat{u}\hfill \\ \hat{v}\hfill \end{array}\right)=\left(\begin{array}{l}M{e}^{T}sin(\omega t+\varphi )\hfill \\ N{e}^{T}cos(\omega t+\varphi )\hfill \end{array}\right),$(13) ), we have an analytical expression of the zero-rate energy function ${\displaystyle \hat{z}\left(t\right)}$, that is:(22) ${\displaystyle \hat{z}\left(t\right)=\frac{a_3}{a_2}\frac{M^2}{N}{e}^T\tan \left(\omega t+\phi \right)\sin \left(\omega t+\phi \right)-M{e}^{{\rm T}}\cos \left(\omega t+\phi \right).}$(22) The function ${\displaystyle \hat{z}(t)}$ has singularities at the points where the function $\text{cos}(\omega t+\phi )$ is equal to zero, i.e. $t_{sing}=\frac{(\kappa \pi +\frac{\pi }{2}-\phi )}{\omega }, \kappa =0, 1, 2,\dots$ . However, the effect of ${\displaystyle \hat{z}\left(t\right)}$, i.e., the integral of the function within a small time interval ${\displaystyle \left[t_1,t_2\right]}$, is finite.

Figure 6 illustrates the effect of ${\displaystyle \hat{z}\left(t\right)}$ by integrating between ${\displaystyle t=0}$ and ${\displaystyle t=\pi /\omega }$, for values of the model parameters ${\displaystyle r}$ and ${\displaystyle B}$ where ${\displaystyle L}$ is unstable. In each case, the system was initialized near the spiral point ${\displaystyle L}$.

Figure 6. Dependence of the zero-rate energy function ${\displaystyle \hat{z}\left(t\right)}$ effect on the functional tumor-associated vasculature B and immune recruitment rate r. Dependence of the effect of ${\displaystyle \hat{z}\left(t\right)}$ on the model parameters r and B in the unstable regime of ${\displaystyle L}$ denoted by U. The label (S)stands for the regime where the spiral point is stable and tumor control is reached.

Using the function ${\displaystyle \hat{z}\left(t\right)}$ is not feasible for practical purposes. However, we can design meaningful therapy functions ${\displaystyle {\hat{z}}_{ext}\left(t\right)}$ that induce the same or greater effects in a certain time interval ${\displaystyle \left[t_1,t_2\right]}$ as the function ${\displaystyle \hat{z}\left(t\right)}$, that is:(23) ${\displaystyle \underset{t_1}{\overset{t_2}{{\displaystyle \int }}}\hat{z}\left(t\right) dt\le \underset{t_1}{\overset{t_2}{{\displaystyle \int }}}{\hat{z}}_{ext}\left(t\right) dt.}$(23)

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 ${\displaystyle \hat{z}\left(t\right)}$ 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 ${\displaystyle \hat{z}\left(t\right)}$, we suggest an efficient external immune-modulatory term ${\displaystyle {\hat{z}}_{ext}\left(t\right)}$ to fulfill the relation (23) that 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. (14(14) ${\displaystyle V\left(\hat{u},\hat{v}\right)=\frac{1}{2}\left(a_3{\hat{u}}^2-a_2{\hat{v}}^2\right),}$(14) ) 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. (22(22) ${\displaystyle \hat{z}\left(t\right)=\frac{a_3}{a_2}\frac{M^2}{N}{e}^T\tan \left(\omega t+\phi \right)\sin \left(\omega t+\phi \right)-M{e}^{{\rm T}}\cos \left(\omega t+\phi \right).}$(22) ).

To illustrate the simulation results and without loss of generality, we select the following parameter values: ${\displaystyle r=0.6 day{s}^{-1}}$ and ${\displaystyle B=0.8,}$ which results in ${\displaystyle \epsilon =5.85·{10}^{-4}}$ and $\text{ω}=\sqrt{-a_2{a}_3}=0.336 rad/days$. Notice that similar results can be obtained for each parameter set where ${\displaystyle L}$ is a spiral source, since ${\displaystyle \epsilon <<1}$ always holds. For numerical integration, we consider a 4^th order Runge-Kutta method. The time step was set equal to ${\displaystyle dt=0.001}$ days, which is sufficiently small.

Validity of Approximate Solutions

At this point, we need to validate the performance of the approximate solutions. Figure 7A depicts the evolution of the approximate perturbation ${\displaystyle \hat{u}}$ compared to the perturbation ${\displaystyle u=R-R_L}$ of the initial model given by Eqs. (1(1) $\begin{array}{l}\frac{dR}{dt}=\frac{1}{3}\left({\lambda }_{{\rm M}}{\rm B}-{\lambda }_{{\rm A}}\right)R\\ +{\lambda }_{{\rm M}}\left(1-{\rm B}\right)L_D\left(\frac{1}{tanh\left(R/L_D\right)}-\frac{L_D}{R}\right)\\ -cER f\left(R,B\right),\end{array}$(1) )–(2(2) ${\displaystyle \frac{dE}{dt}=r\frac{R^3}{K+R^3}E-d_{1}E{R}^3 f\left(R,B\right)-d_{0}E+\sigma ,}$(2) ), referred in what follows as original model perturbation. The approximate solution ${\displaystyle \hat{u}}$ was estimated by performing a numerical integration of the system (5), while ${\displaystyle u=R-R_L}$ by numerically integrating the original model (1)–(2). Both systems of equations were initialized close to the spiral point ${\displaystyle L}$ and the simulation time was set ${\displaystyle t_{sim}=1000 days}$ (about 3 years). A quick glance at Figure 7A 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, Figure 7A 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 ${\displaystyle {10}^{-3}}$.

Figure 7. Validation of approximate solutions. (A) Evolution of the approximate perturbation ${\displaystyle \hat{u}}$ compared to the original model perturbation ${\displaystyle u=R-R_L}$. (B) Evolution of the energy function of the approximate linearized model defined in Eq. (14(14) ${\displaystyle V\left(\hat{u},\hat{v}\right)=\frac{1}{2}\left(a_3{\hat{u}}^2-a_2{\hat{v}}^2\right),}$(14) ) against the energy estimated by using the original model (1)–(2). (C) The zero-rate energy function in Eq. (21(21) ${\displaystyle \hat{z}\left(t\right)=\frac{a_3{\hat{u}}^2-a_2{\hat{v}}^2}{a_2\hat{v}}}$(21) ) estimated by using the approximate solution of the system (13). (D) The immuno-modulatory function ${\displaystyle {\hat{z}}_{ext}\left(t\right)}$ for ${\displaystyle \alpha =-0.6}$ and ${\displaystyle \gamma =100}$ compared with the zero-rate energy function ${\displaystyle \hat{z}\left(t\right)}$.

Figure 7B compares the evolution of the approximate energy function defined in Eq. (14(14) ${\displaystyle V\left(\hat{u},\hat{v}\right)=\frac{1}{2}\left(a_3{\hat{u}}^2-a_2{\hat{v}}^2\right),}$(14) ) 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.

Proposal of a Therapeutic Term

Figure 7C depicts the zero-rate energy function in Eq. (21(21) ${\displaystyle \hat{z}\left(t\right)=\frac{a_3{\hat{u}}^2-a_2{\hat{v}}^2}{a_2\hat{v}}}$(21) ) estimated by using the approximate solutions ${\displaystyle \hat{u}}$ and ${\displaystyle \hat{v}}$ from the system (18). 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 ${\displaystyle u=R-R_L}$, i.e. the function follows the monotonicity of tumor evolution. Hence, an external immuno-modulatory therapy ${\displaystyle {\hat{z}}_{ext}\left(t\right)}$ that satisfies the condition (23) 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:(24) ${\displaystyle {\hat{z}}_{ext}\left(t\right)=\alpha tanh\left(\gamma \hat{v}\right),}$(24) where ${\displaystyle \gamma >>1}$ expresses the sigmoidal steepness, and the selection of ${\displaystyle \alpha }$ should fulfill the relation (23). Figure 7D depicts the immuno-modulatory function ${\displaystyle {\hat{z}}_{ext}\left(t\right)}$ for ${\displaystyle \alpha =-0.6}$ and ${\displaystyle \gamma =100}$, compared with the estimated zero-rate energy function ${\displaystyle \hat{z}\left(t\right)}$.

Performance of the Proposed Therapeutic Strategy

The application of the external therapy in Eq. (24(24) ${\displaystyle {\hat{z}}_{ext}\left(t\right)=\alpha tanh\left(\gamma \hat{v}\right),}$(24) ) implies the addition of a new term in Eqs. (1(1) $\begin{array}{l}\frac{dR}{dt}=\frac{1}{3}\left({\lambda }_{{\rm M}}{\rm B}-{\lambda }_{{\rm A}}\right)R\\ +{\lambda }_{{\rm M}}\left(1-{\rm B}\right)L_D\left(\frac{1}{tanh\left(R/L_D\right)}-\frac{L_D}{R}\right)\\ -cER f\left(R,B\right),\end{array}$(1) )–(2(2) ${\displaystyle \frac{dE}{dt}=r\frac{R^3}{K+R^3}E-d_{1}E{R}^3 f\left(R,B\right)-d_{0}E+\sigma ,}$(2) ), that is:(25) $\begin{array}{l}\frac{dR}{dt}=\frac{1}{3}\left({\lambda }_{{\rm M}}{\rm B}-{\lambda }_{{\rm A}}\right)R+{\lambda }_{{\rm M}}\left(1-{\rm B}\right)L_D\left(\frac{1}{tanh\left(R/L_D\right)}-\frac{L_D}{R}\right)\\ -cER f\left(R,B\right),\end{array}$(25) (26) ${\displaystyle \frac{dE}{dt}=r\frac{R^3}{K+R^3}E-d_{1}E{R}^3 f\left(R,B\right)-d_{0}E+\sigma +\epsilon \alpha \tanh \left(\gamma v\right),}$(26) where ${\displaystyle v=E-E_L}$.

Figure 8A,B respectively show the evolution of the tumor radius and effector cells by numerically integrating Eqs. (25(25) $\begin{array}{l}\frac{dR}{dt}=\frac{1}{3}\left({\lambda }_{{\rm M}}{\rm B}-{\lambda }_{{\rm A}}\right)R+{\lambda }_{{\rm M}}\left(1-{\rm B}\right)L_D\left(\frac{1}{tanh\left(R/L_D\right)}-\frac{L_D}{R}\right)\\ -cER f\left(R,B\right),\end{array}$(25) )–(26(26) ${\displaystyle \frac{dE}{dt}=r\frac{R^3}{K+R^3}E-d_{1}E{R}^3 f\left(R,B\right)-d_{0}E+\sigma +\epsilon \alpha \tanh \left(\gamma v\right),}$(26) ) that include the proposed external immuno-modulatory function (24). The system is initialized near the spiral point ${\displaystyle L}$ and parameter values of the external function were selected equal to ${\displaystyle \alpha =-0.6}$ and ${\displaystyle \gamma =100}$. In doing so, we obtain that the system reaches a stable steady state, i.e., the tumor remains controlled. Figure 8C illustrates how the energy of the system of Eqs. (25(25) $\begin{array}{l}\frac{dR}{dt}=\frac{1}{3}\left({\lambda }_{{\rm M}}{\rm B}-{\lambda }_{{\rm A}}\right)R+{\lambda }_{{\rm M}}\left(1-{\rm B}\right)L_D\left(\frac{1}{tanh\left(R/L_D\right)}-\frac{L_D}{R}\right)\\ -cER f\left(R,B\right),\end{array}$(25) )–(26(26) ${\displaystyle \frac{dE}{dt}=r\frac{R^3}{K+R^3}E-d_{1}E{R}^3 f\left(R,B\right)-d_{0}E+\sigma +\epsilon \alpha \tanh \left(\gamma v\right),}$(26) ) evolves. The energy decreases with time and becomes zero as the system reaches the steady-state. It is worth pointing out that ${\displaystyle \alpha }$ was selected not only to satisfy the relation (23), 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 (24) with parameters ${\displaystyle \alpha =-0.6}$ and ${\displaystyle \gamma =100}$ (C) Energy of the system of equations (25(25) $\begin{array}{l}\frac{dR}{dt}=\frac{1}{3}\left({\lambda }_{{\rm M}}{\rm B}-{\lambda }_{{\rm A}}\right)R+{\lambda }_{{\rm M}}\left(1-{\rm B}\right)L_D\left(\frac{1}{tanh\left(R/L_D\right)}-\frac{L_D}{R}\right)\\ -cER f\left(R,B\right),\end{array}$(25) )–(26(26) ${\displaystyle \frac{dE}{dt}=r\frac{R^3}{K+R^3}E-d_{1}E{R}^3 f\left(R,B\right)-d_{0}E+\sigma +\epsilon \alpha \tanh \left(\gamma v\right),}$(26) ) with ${\displaystyle \alpha =-0.6}$ and ${\displaystyle \gamma =100}$. Evolution of the tumor radius (D) and concentration of effector cells (E) by initializing the system away from the spiral point L and considering the external immuno-modulatory function with ${\displaystyle \alpha =-25}$ and ${\displaystyle \gamma =100}$. (F) Energy of the system of equations (25(25) $\begin{array}{l}\frac{dR}{dt}=\frac{1}{3}\left({\lambda }_{{\rm M}}{\rm B}-{\lambda }_{{\rm A}}\right)R+{\lambda }_{{\rm M}}\left(1-{\rm B}\right)L_D\left(\frac{1}{tanh\left(R/L_D\right)}-\frac{L_D}{R}\right)\\ -cER f\left(R,B\right),\end{array}$(25) )–(26(26) ${\displaystyle \frac{dE}{dt}=r\frac{R^3}{K+R^3}E-d_{1}E{R}^3 f\left(R,B\right)-d_{0}E+\sigma +\epsilon \alpha \tanh \left(\gamma v\right),}$(26) ), with ${\displaystyle \alpha =-25}$ and ${\displaystyle \gamma =100}$, initialized away from the spiral point L.

Figure 8D, E respectively show the evolution of the tumor radius and concentration of effector cells by numerically integrating Eqs. (25(25) $\begin{array}{l}\frac{dR}{dt}=\frac{1}{3}\left({\lambda }_{{\rm M}}{\rm B}-{\lambda }_{{\rm A}}\right)R+{\lambda }_{{\rm M}}\left(1-{\rm B}\right)L_D\left(\frac{1}{tanh\left(R/L_D\right)}-\frac{L_D}{R}\right)\\ -cER f\left(R,B\right),\end{array}$(25) )–(26(26) ${\displaystyle \frac{dE}{dt}=r\frac{R^3}{K+R^3}E-d_{1}E{R}^3 f\left(R,B\right)-d_{0}E+\sigma +\epsilon \alpha \tanh \left(\gamma v\right),}$(26) ) when the system is initialized away from the spiral point ${\displaystyle L}$. In this case, the approximate solution is expected to deviate from the solutions of the original model. Figure 8F presents the temporal evolution of the corresponding energy. The parameter ${\displaystyle \alpha }$ was selected to be equal to −25 to fulfill the relation (23), 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 ${\displaystyle L}$.

Discussion

In this article, we investigate the therapeutic potential of immunomodulatory interventions against tumor growth. To that end, we consider a model 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.⁷⁸ Therefore, a potential realization of our proposed strategy could be mediated by a combination of state-of-the-art immunotherapies and chemotherapeutic modules.^1,4,8,49 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. Here, we assumed that tumor-associated vascularization is considered time invariant. Therefore, dynamic modeling of vascularization dynamics is required to capture, for instance, potential hypoxic effects due to vasculature perturbations. Moreover, the immune system is much more complicated than the current model description, since it involves more cell types and their corresponding interactions. In particular, immune system is usually regarded as acting against tumor growth. However, recently became clear that it can be both stimulatory and inhibitory, as in the case of tumor-associated macrophages.^79,80 Furthermore, our model cannot be applied to invasive/diffusive tumors, since it requires further modifications. Nevertheless, simple mathematical models allow for a profound understanding of the crucial dynamic phenomena involved in immune-tumor interactions that are essential for generating novel therapeutic ideas.

Disclosure of Potential Conflicts of Interest

No potential conflicts of interest were disclosed.

Funding

AI Reppas, JCL 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. JCL Alfonso and H Hatzikirou also acknowledge the German Federal Ministry of Education and Research (BMBF) for the eMED project SYSIMIT (01ZX1308D).