# A size and space structured model of tumor growth describes a key role for protumor immune cells in breaking equilibrium states in tumorigenesis

Kevin Atsou, Fabienne Anjuère, Veronique M Braud, Thierry Goudon

► To cite this version:

Kevin Atsou, Fabienne Anjuère, Veronique M Braud, Thierry Goudon. A size and space structured model of tumor growth describes a key role for protumor immune cells in breaking equilibrium states in tumorigenesis. PLoS ONE, 2021, 16 (11), pp.e0259291. 10.1371/journal.pone.0259291. hal-03452282

HAL Id: hal-03452282

https://hal.science/hal-03452282

Submitted on 26 Nov 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
RESEARCH ARTICLE

# A size and space structured model of tumor growth describes a key role for protumor immune cells in breaking equilibrium states in tumorigenesis

Kevin Atsou<sup>1*</sup>, Fabienne Anjuère<sup>2</sup>, Véronique M. Braud<sup>2</sup>, Thierry Goudon<sup>1*</sup>

1 Université Côte d'Azur, Inria, CNRS, LJAD, Parc Valrose, Nice, France, 2 Université Côte d'Azur, CNRS, Institut de Pharmacologie Moléculaire et Cellulaire UMR 7275, Valbonne, France

Current address: Inria-Inserm Team COMPO, Inserm U1068, CNRS UMR7258, Institut Paoli-Calmettes Pharmacy faculty, Aix-Marseille University, Marseille, France

* thierry.goudon@inria.fr

<figure><img src="image_2.png" /></figure>

## Abstract

Switching from the healthy stage to the uncontrolled development of tumors relies on complicated mechanisms and the activation of antagonistic immune responses, that can ultimately favor the tumor growth. We introduce here a mathematical model intended to describe the interactions between the immune system and tumors. The model is based on partial differential equations, describing the displacement of immune cells subjected to both diffusion and chemotactic mechanisms, the strength of which is driven by the development of the tumors. The model takes into account the dual nature of the immune response, with the activation of both antitumor and protumor mechanisms. The competition between these antagonistic effects leads to either equilibrium or escape phases, which reproduces features of tumor development observed in experimental and clinical settings. Next, we consider on numerical grounds the efficacy of treatments: the numerical study brings out interesting hints on immunotherapy strategies, concerning the role of the administered dose, the role of the administration time and the interest in combining treatments acting on different aspects of the immune response. Such mathematical model can shed light on the conditions where the tumor can be maintained in a viable state and also provide useful hints for personalized, efficient, therapeutic strategies, boosting the antitumor immune response, and reducing the protumor actions.

<figure><img src="image_3.png" /></figure>

### OPEN ACCESS

**Citation:** Atsou K, Anjuère F, Braud VM, Goudon T (2021) A size and space structured model of tumor growth describes a key role for protumor immune cells in breaking equilibrium states in tumorigenesis. PLoS ONE 16(11): e0259291.
https://doi.org/10.1371/journal.pone.0259291

**Editor:** Christophe Letellier, Normandie Université, FRANCE

**Received:** June 17, 2021

**Accepted:** October 15, 2021

**Published:** November 22, 2021

**Peer Review History:** PLOS recognizes the benefits of transparency in the peer review process; therefore, we enable the publication of all of the content of peer review and author responses alongside final, published articles. The editorial history of this article is available here: https://doi.org/10.1371/journal.pone.0259291

**Copyright:** © 2021 Atsou et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data Availability Statement:** Relevant data are within the manuscript. Codes and numerical data are available at the URL https://github.com/atsoukevin93/tumorgrowth.

## Introduction

The immune system can both constrain and promote tumor development through several complex processes, encompassed in the concept of cancer immunoediting [1]. The antagonistic effects of the immune response on tumor growth shape the different phases that have been identified to characterize their interaction: elimination, when the immune system is able to detect and eradicate the developing tumors; equilibrium, when the immune system is able to
**Funding:** This work was supported by the French Government (National Research Agency, ANR) through the "Investments for the Future" programs LABEX SIGNALIFE ANR-11-LABX-0028 and IDEX UCAJedi ANR-15-IDEX-01. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

**Competing interests:** The authors have declared that no competing interests exist.

maintain the tumor expansion in a cancer-persistent state; and escape, when the tumor develops in an uncontrolled manner [1, 2].

In this context, the identification of the immune components of the tumor microenvironment (TME) reveals valuable information about the stage of cancer development and helps predict patient outcome. This concept called "the immune contexture" has improved the classification of cancers [3, 4]. The antitumor immune response is characterized by the activation and the recruitment of innate immune cells such as natural killer (NK) cells, tumor-associated neutrophils (TAN-N1), tumor-associated macrophages (TAM-M1) and adaptive immune CD8+ T cells. They migrate to the tumor site where they can eliminate tumor cells. They have been found to be highly active on early-stage tumors and associated with good clinical outcome [5, 6]. While this antitumor immune response can be expected to control tumor growth or maintain their development in a viable equilibrium, later phases are characterized by an uncontrolled tumor growth associated with a shift of the immune response towards protumor functions and the establishment of multiple mechanisms of immunosuppression [7]. Among others, the ratio of effector immune cells/protumor immune cells is considered as a relevant indicator of patient survival, the higher the ratio, the better the patient vital prognostic [8]. The ratio evolves dynamically: tumor cells and other components in the TME can produce inhibitory factors such as anti-inflammatory cytokines, interleukins 10 and 4 (IL-10 and IL-4), Transforming Growth Factor-beta (TGF-β) which favor the polarization of antitumor immune cells into protumor ones. For instance, antitumor neutrophils and macrophages are converted into protumor TAN-N2 and TAM-M2 [9, 10]. They are part of a pool of myeloid-derived suppressor cells (MDSCs) which can also be directly recruited from the bone marrow [11]. They promote tumor growth, tissue remodeling, angiogenesis and suppress adaptive immunity [12]. Moreover, the antigen-presenting cells such as dendritic cells (DC) become tolerogenic which leads to exhausted and tolerant T cells, apoptosis of T cells and to the priming and proliferation of regulatory T cells (Tregs) [12, 13]. Besides, tumors subvert mechanisms of immune self-tolerance and inhibit antitumor immune responses through the use of immune checkpoints. The TME is therefore infiltrated with hyporesponsive exhausted T cells [14]. Blockade of these immune checkpoints as exemplified with anti-CTLA4 and anti-PD-1/PD-L1 restore efficient effector functions and has revolutionized cancer treatment [15].

Mathematical modeling might shed some light on these complex interactions, and, based on numerical simulations, provide useful information to elaborate more efficient therapeutic strategies. Quite intricate ordinary differential systems have been developed so far, see for instance [16–22]. Further references and discussion of the various viewpoints can be found in the review [23]. In [24], we proposed a system of partial differential equations (PDE), describing the earliest stages of the tumor/immune system interactions. The system couples an integro-differential equation for the size-structured population of tumor cells, inspired from [25–28], to a convection-diffusion equation for the space-structured immune cells. The latter accounts for chemotaxis mechanisms that drive the immune cells towards the tumor. This model, which only considers the antitumor actions of the immune system, is able to reproduce the equilibrium phase: the large time behavior of this PDEs system is a state where residual tumor cells and a positive concentration of active immune cells exist in equilibrium. However, the simple model of [24] does not consider the contribution of immune cells with protumor functions and the establishment of numerous mechanisms of immunosuppression. This is the issue addressed in this work. To be more specific, our purpose is two-fold:

* First, we incorporate in the model protumor effects that can both reduce the antitumor immune response and strengthen the factors of tumor growth. We shall see that the
protumor immune response can break the equilibrium and lead to an escape phase characterized by the uncontrolled growth of the tumor.

* Second, we complete the model by discussing the effect of different type of targeted treatments that can act on the immune response, either by restoring the effector functions of anti-tumor cells which became exhausted as a result of chronic activation and protumor factors, or by limiting the recruitment of protumor immune cells. The investigation demonstrates the interest of combining both approaches.

The paper is organized as follows. We introduce the modeling assumptions and we set up the equations in the **Mathematical Model** section. We pay a specific attention to the description of the activation of protumor mechanisms, based on the action of cytokines, in response to the growth of the tumor mass. The modeling assumptions naturally induce some delay mechanisms in the protumor response. We also bring out some mathematical properties of the model. We start by considering a simplified situation which reduces the model to a nonlinear system of ODEs. We identify several stationary solutions, free of tumors, free of protumor cells or with all populations of immune cells, and we discuss their linear stability. This discussion provides some hints on the role of the parameters. Next, we establish the existence of equilibrium states for the full size- and space-dependent problem, extending to the model with protumor activities the results of [24]. The **Results** section is devoted to the numerical investigation of the PDE system. We show that depending on the parameters of the model, the solutions either converge to an equilibrium or describe an escape phase with an unlimited growth of the tumor. These results illustrate the critical role of the protumor immune responses. Next, we address specifically the question of treatments. As detailed below, the immune response might have multi-faceted protumor actions. Among other, effector immune cells can become exhausted, a state where they are hyporesponsive and cannot kill the tumor. We consider the effect of treatments that either restore antitumor activity or reduce the recruitment of protumor immune cells. The investigations demonstrate the interest of combining both approaches and discuss the role of the dose and time of administration.

## Mathematical model

A schematic overview of the geometry and the leading mechanisms that guide the construction of the mathematical model is provided in Figs 1 and 2.

### Modeling assumptions

We take into account three populations of interacting cells:

* the cytotoxic effector cells including CD8+ T cells and NK cells as well as myeloid effector cells, TAN-N1 and TAM-M1 that will be referred to as the “antitumor” immune cells;

* the “protumor” immune cells, including Treg, MDSCs, TAN-N2 and TAM-M2 favoring tumor growth;

* the tumor cells.

For the purpose of the mathematical modeling, we are collapsing into single “averaged” quantities the behavior of several distinct cells, that might have different developmental properties and interactions mechanism, to focus on the outcome of their action on tumor growth. The construction of the model uses the same basis as in [24], to which we incorporate the “protumor” immune cells. The modeling assumptions are as follows.
<figure><img src="image_4.png" /><figcaption>Fig 1. Schematic view of the geometry of the mathematical model.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0259291.g001

A.1. the environmental constraints such as nutrient concentrations, temperature, etc. are assumed to be constant. Nevertheless, in late stages of tumor growth, some phenomena such as hypoxia or difficulties in accessing the nutrients can limit the tumor cell expansion;

A.2. the state of the tumor cells is characterized by their size (with a similar setting, it could be their content of cyclins as well [25, 29]); the dynamics of the tumor cells is driven by two phenomena: each tumor cell grows with a certain rate, possibly depending on its size, and it can divide into daughter cells;

A.3. activated antitumor immune cells are able to destroy the targeted tumor cells;

A.4. activated protumor immune cells suppress the antitumor immune cells by direct contact or by the release of soluble substances (like immunosuppressive cytokines);

A.5. activated protumor immune cells favor the tumor growth by enhancing the growth rate of the tumor cells and by favoring angiogenesis.

Moreover, the tumor cells produce several signals of chemical nature (cytokines and chemokines), which drive the immune response as follows:

A.6. a chemotactic signal, proportional to the tumor mass, induces a potential, the gradient of which drives the anti and protumor cells towards the TME;
<figure><img src="image_5.png" /><figcaption>**Fig 2.** Schematic view of the leading mechanisms that guide the construction of the mathematical model. AIC: antitumor cells, PIC: protumor cells, TC: tumor cells.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0259291.g002

A.7. the tumor antigen-specific CD8⁺ T cells are activated by APCs in draining lymph nodes and recruited to the tumor site. The NK cells as well as the myeloid cells (TAN-N1, TAM-M1) are recruited from the circulation and activated from a bath of non-activated immune cells at the tumor site. The signal that defines the recruitment/activation rate is directly related to the tumor mass.

A.8. similarly, protumor immune cells can be recruited from a bath of immune cells (that might differ from the source of antitumor immune cells, though) according to a signal directly related to the tumor mass;

A.9. the signal triggers the shift of certain antitumor immune cells into protumor immune cells.

Assumptions **A.1–A.3** appeared in [24] where they are discussed in details. The protumor effects become sensitive in a later stage of the tumor growth, and, as we shall see, play a central role in the transition to the escape phase. Assuming a constant growth rate of the tumor cells becomes questionable in this regime, and we shall model it by means of a Gompertz law, which accounts for size-limitation mechanisms, see (1) below. Assumption **A.4** describes immunosuppression mechanisms mediated by protumor immune cells. In addition to the contact-dependent suppression of antitumor immune cells, the secretion of immunosuppressive cytokines abrogates the effector functions of T cells and NK cells and negatively modifies their proliferation. It triggers the reverse conditioning of DCs and can induce the apoptosis of effector T cells through the depletion of IL-2 from the TME. It is worth bearing in mind that not all the antitumor immune cells are eliminated: they become exhausted and can no longer
kill tumor cells; however, they can be reactivated by the action of some treatment (i. e. anti-PD-1 therapy). Assumption A.5 corresponds to the overexpression of VEGF by protumor immune cells favoring the accumulation of microvessels supplying the tumor in nutrients. Assumption A.6 already appeared in [24]; here the chemotactic mechanisms apply on both type of immune cells. Assumption A.7 corresponds to a rough description of the complex activation and recruitment mechanisms which are related to the tumor mass.

Similarly, assumption A.8 corresponds to the recruitment of MDSC coming from the bone marrow and Tregs from the circulation. Assumption A.9 corresponds to the possible conversion of some immune cells that eliminate the tumor into protumor immune cells (i. e. macrophages and neutrophils becoming TAN M2 and TAM N2, conventional T cells becoming immunosuppressive Treg).

## Construction of the model

Following [24], the model uses two distinct length scales:

* the length scale of the displacement of the immune cells,

* the length scale describing the size of the tumor cells.

The modeling assumes that the former is “infinitely large” compared to the latter. This is motivated by the fact that immune cell displacement (for instance from the lymph nodes toward the tumor site) occurs on several centimeters while the estimated radius of the tumor cells is about a few micrometers. The interactions between the tumor and the immune system are described by the dynamics of the following unknowns:

* **Tumor cell density.** The population of tumor cells is described by $(t, z) \mapsto n(t, z)$, the volumic density of tumor cells. Given $z_2 > z_1 > 0$, the integral $\int_{z_1}^{z_2} n(t, z) dz$ gives the number of tumor cells having a size in the interval $[z_1, z_2]$ at time $t$.

* **Cytotoxic effector cell concentration.** The concentration of antitumor immune cells that are actively fighting against the tumor is $(t, x) \mapsto c(t, x)$.

* **Protumor immune cell concentration.** Similarly, $(t, x) \mapsto c_r(t, x)$ stands for the concentration of the protumor immune cells favoring tumor growth.

* **Chemoattractant potential.** We denote by $(t, x) \mapsto \phi(t, x)$ the concentration of the chemical signal (chemokines) that attracts the immune cells towards the TME.

* **Cytokine concentration.** Finally, let $t \mapsto I(t)$ be the concentration of cytokines in the overall TME.

The dynamics of the population of tumor cells is governed by volume growth and cellular division, see [25–28]. We add to these effects a death rate induced by the activated antitumor immune cells. Let $z \mapsto V(z) \ge 0$ stand for the tumor cell growth rate. We can assume it is a positive constant, but in the present context it is more appropriate to adopt a size dependent model, that incorporates size-limitation effects, like the Gompertz law

$$V(z) = rz \ln(b/z), \qquad (1)$$

with $r > 0$ and $b > 0$, the maximal size. Accordingly, when using (1), the size variable $z$ lies in the interval $[0, b]$. We refer the reader to [30–34] for derivation and use of this law in tumor growth modeling, in particular when taking into account the limited access to nutrients or necrotic mechanisms, see A.1, A.2. The simpler case where $V$ is constant can be useful since it leads to explicit formulae that can be used to check more easily whether or not the conjectured
behaviors hold. The cell division mechanism is described by the operator

$$Q(n)(t, z) = -a(z)n(t, z) + \int_{z}^{b} a(z')k(z|z')n(t, z') dz', \quad (2)$$

where $a(z')$ stands for the frequency of the division of cells having size $z'$, and $k(z|z')$ gives the size-distribution produced from the division of a tumor cell with size $z'$. Again a simplified framework assumes that these coefficients are constant, but it is likely relevant to make use of more intricate laws, for instance prohibiting any division below a certain threshold. What is crucial for modeling purposes is the requirement

$$\int_{0}^{z} z'k(z'|z) dz' = z,$$

which is related to the principle of mass conservation. Indeed, it implies that cell-division does not change the total mass

$$\int_{0}^{b} zQ(n) dz = 0.$$

However, the total number of cells increases since $\int_{0}^{b} n(t, z) dz \ge 0$. A relevant example is provided by the binary division operator

$$Q(n)(t, z) = 4a(2z)n(t, 2z) - a(z)n(t, z), \quad (3)$$

which describes the situation where cells with size $2z$ split into two daughter cells, both with size $z$. Further relevant examples of division kernels can be found in [35]. The equation will be completed by the boundary condition $n(t, 0) = 0$, which means that there is no production of cells with size 0. For further purposes, let us introduce the following quantities

$$\text{total number of tumor cells : } \mu_0(t) = \int_0^b n(t, z) dz, \quad (4)$$

$$\text{total mass of tumor cells : } \mu_1(t) = \int_0^b zn(t, z) dz. \quad (5)$$

The displacement of the anti- and protumor immune cells is driven by convection and diffusion, over a domain $\Omega \subset \mathbb{R}^N$. For the sake of simplicity we assume they have the same diffusion coefficient $D$ and, here, we shall work with a constant coefficient. The convection is defined by the chemotactic potential $\phi$, which depends on the total mass of the tumor. It obeys the diffusion equation

$$-\nabla_x \cdot (K \nabla_x \phi) = f(\mu_1)\sigma, \quad K \nabla_x \phi \cdot v|_{\partial\Omega} = 0 \quad (6)$$

endowed with the homogeneous Neumann boundary condition. In (6), $x \mapsto \sigma(x)$ is a given form function with zero-mean, $K > 0$ is a positive coefficient. The coefficient $K$ could be matrix-valued as well, taking into account further details of the vasculature or the properties of the tissues neighboring the tumor, that govern the supply of immune cells. The strength of the potential depends on the total mass of the tumor through the function $\mu_1 \mapsto f(\mu_1) \ge 0$. It is natural to assume that $f(0) = 0$ and $f$ is non decreasing. A typical example of $f$ is the following
Michealis-Menten functional response:

$$f(\mu_1) = \frac{\mu_1}{\eta + \mu_1}, \quad \eta > 0. \qquad (7)$$

We suppose that c and $c_r$ have the same chemotactic sensitivity $\chi > 0$, and they both satisfy homogeneous Dirichlet boundary condition on $\partial\Omega$: the immune cells far from the tumor are non-activated.

Let us now describe the zeroth order terms of the equations, that differ depending on the considered type of cells. Both type of immune cells is subjected to a death rate $\gamma$, $\gamma_r > 0$. The antitumor immune cells are recruited from a source of naive immune cells $(t, x) \mapsto S(t, x)$. The activation process is described through a rate $\mu_1 \mapsto g(\mu_1) \ge 0$, which depends on the total mass of the tumor. Again, it is natural to assume $g(0) = 0$ and $g$ is non decreasing. There are two other mechanisms that lead to a loss of antitumor immune cells. First, according to assumption **A.4**, the protumor immune cells suppress effector cells; this is traduced by a loss term

$$-k_c c c_r$$

where $k_c > 0$ is the rate of this reaction. Second, according to assumption **A.9**, certain activated effector cells can be converted into protumor immune cells under the action of cytokines in the TME. This leads to the loss term

$$-k_r I \theta c$$

where $k_r > 0$ is the rate of this conversion, and $x \mapsto \theta(x)$ is a given form function, say a peaked Gaussian, indicating that such processes occur in the vicinity of the tumor. This loss of antitumor immune cells contributes to the gain term for the population of protumor immune cells. Cytokines also activate protumor immune cells from a distant source denoted by $S_r$ according to assumption **A.8**.

The effector cells release cytotoxic substances in the TME. This effect is described by the death term

$$m(c, n)(t, z) = n(t, z) \times \int_{\Omega} \delta(y) c(t, y) dy \qquad (8)$$

in the tumor growth equation. It involves a non negative space-dependent weight $x \mapsto \delta(x)$, which incorporates both the strength of the immune response and a radius of interaction. According to assumption **A.5**, recruited protumor immune cells favor the tumor growth. Therefore the growth rate of the tumor cells is enhanced by the presence of protumor immune cells and it becomes

$$V(z) \left( 1 + \int_{\Omega} b_1(y) c_r(t, y) dy \right) \qquad (9)$$

with a certain non negative, radially symmetric and compactly supported kernel $b_1$.

Finally, we turn to the dynamics of the tumor-secreted cytokines, which promote the protumor reactions. The production of such cytokines occurs beyond a certain critical mass, denoted by $m$. Moreover, the cytokine concentration is naturally damped with a constant rate $\tau > 0$. This leads to the ODE

$$\partial_t I = \psi(\mu_1) - \tau I \qquad (10)$$

where $\psi$ is a threshold function, non negative and non decreasing. For instance, given a
constant $\bar{\psi} > 0$, it can be defined by:

$$ \psi(\mu_1) = \bar{\psi} \begin{cases} (\mu_1 - m), & \mu_1 > m \\ 0, & \mu_1 \le m \end{cases} \qquad (11) $$

We shall need the technical assumptions $\psi(0) = \psi'(0) = 0$, which clearly holds when $m > 0$. The threshold $m$ can be used to describe the degree of inflammation of the tumor environment: the smaller $m$, the more inflamed the environment, and a reduced $m$ can correspond to altered soil, due to systemic effects caused by a primary tumor. Eventually, we arrive at the following system describing the interactions between the tumor cells, the effector and protumor immune cells:

$$ \partial_t n + \partial_z \left( V(z) \left( 1 + \int_{\Omega} b_1(y) c_r(t, y) dy \right) n \right) = Q(n) - m(n, c), \qquad (12a) $$

$$ \partial_t c + \nabla_x \cdot (c\chi\nabla_x\phi - D\nabla_x c) = g(\mu_1)S - \gamma c - k_r I \theta c - k_c c c_r, \qquad (12b) $$

$$ \partial_t c_r + \nabla_x \cdot (c_r \chi \nabla_x \phi - D \nabla_x c_r) = I(S_r + k_r \theta c) - \gamma_r c_r, \qquad (12c) $$

$$ \partial_t I = \psi(\mu_1) - \tau I, \qquad (12d) $$

$$ -\nabla_x \cdot (K \nabla_x \phi) = f(\mu_1)\sigma, \qquad (12e) $$

$$ n(t, 0) = 0, \quad c|_{\partial\Omega} = 0, \quad c_r|_{\partial\Omega} = 0, \quad K\nabla_x\phi \cdot v(\cdot)|_{\partial\Omega} = 0, \qquad (12f) $$

$$ n(t = 0, z) = n_0(z), \quad c(t = 0, x) = c_0(x), \quad c(t = 0, x) = c_r^0(x), \quad I(t = 0) = I_0. \qquad (12g) $$

We remind the reader that the cell division operator $Q(n)$ and the immune cell-tumor interaction term $m(c, n)$ are defined in (2) and (8) respectively. Table 1 recapitulates the biological meaning of the parameters of the model. We refer the reader to [24] for details on the units and typical values of these quantities; further parameter identification from experimental data and sensitivity analysis is detailed in [36].

## A simplified model to assess the role of the parameters on damping and escape

In order to shed some light on the possible behavior of the solutions and on the role of the parameters, we temporarily adopt a set of simplifying assumptions. Let us consider the very specific case where

* $V$ and $a$ are constant,

* the source $S$ of immune cell is constant,

* the source $S_r = 0$ vanishes and the other parameters $D, K$ are constant and positive,

* the coupling function is linear: $g(\mu_1) = \mu_1$,

* the space variation of the concentrations of immune cells is neglected,

* we consider the binary division model (3) with a constant frequency $a$.
<table><caption>Table 1. Recap of the main definitions and notations.</caption><thead><tr><td>variable</td><td>description</td></tr></thead><tbody><tr><td>z</td><td>volume of tumor cells</td></tr><tr><td>t</td><td>time variable</td></tr><tr><td>x</td><td>space variable</td></tr><tr><td>n</td><td>size-density of tumor cells with a volume z</td></tr><tr><td>V</td><td>tumor cells growth rate</td></tr><tr><td>a</td><td>tumor cells division rate</td></tr><tr><td>μ<sub>0</sub></td><td>total number of tumor cells</td></tr><tr><td>μ<sub>1</sub></td><td>total volume of tumor cells</td></tr><tr><td>c</td><td>concentration of antitumor cells</td></tr><tr><td>c<sub>r</sub></td><td>concentration of protumor cells</td></tr><tr><td>χ</td><td>chemotactic coefficient</td></tr><tr><td>φ</td><td>chemotactic potential</td></tr><tr><td>D</td><td>diffusion coefficient of the immune cells</td></tr><tr><td>S</td><td>source of antitumor immune cells</td></tr><tr><td>S<sub>r</sub></td><td>source of protumor immune cells</td></tr><tr><td>K</td><td>diffusion coefficient of the chemotactic signal</td></tr><tr><td>σ</td><td>chemotactic signal</td></tr><tr><td>I</td><td>cytokine concentration</td></tr><tr><td>γ</td><td>death rate of the antitumor immune cells</td></tr><tr><td>γ<sub>r</sub></td><td>death rate of the protumor immune cells</td></tr><tr><td>τ</td><td>damping rate of the cytokine concentration</td></tr><tr><td>k<sub>c</sub></td><td>suppression rate of antitumor cells by the protumor cells</td></tr><tr><td>k<sub>r</sub></td><td>conversion rate of antitumor cells into protumor cells</td></tr><tr><td>θ</td><td>form function of the conversion of antitumor cells into protumor cells</td></tr><tr><td>δ</td><td>form function of the antitumor action in the TME</td></tr><tr><td>b<sub>1</sub></td><td>tumor growth rate induced by the protumor cells</td></tr><tr><td>m</td><td>threshold on the tumor volume driving the cytokine activation</td></tr></tbody></table>

https://doi.org/10.1371/journal.pone.0259291.t001

These assumptions, that are used only in this Section, completely disregard the space dependence of the unknowns and certainly lack of biological relevance; the ambition of the discussion of this simplified framework is to provide an intuition on the role of the parameters. In this simple situation, the dynamics is described by the following system of ordinary differential equations for $\mu_0, \mu_1$, given by (4) and (5), and the time-dependent concentrations of immune cells and of cytokines:

$$ 
\left\{
\begin{aligned}
\frac{d}{dt}\mu_0 &= \mu_0(a - \delta c), \\
\frac{d}{dt}\mu_1 &= V(1 + b_1 c_r)\mu_0 - \delta\mu_1 c, \\
\frac{d}{dt}c &= \mu_1 S - \gamma c - k_r cI - k_c c c_r, \\
\frac{d}{dt}c_r &= k_r cI - \gamma_r c_r, \\
\frac{d}{dt}I &= \psi(\mu_1) - \tau I.
\end{aligned}
\right.
\qquad (13) $$
The state

$$ \begin{pmatrix} \mu_0^H \\ \mu_1^H \\ c^H \\ c_r^H \\ I^H \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{pmatrix} \qquad (H) $$

is a trivial equilibrium solution to (13) which corresponds to a healthy state. However, we can also find equilibrium states with residual tumor cells, effector immune cells and even protumor immune cells. Indeed, let

$$ \begin{pmatrix} \mu_0^{NP} \\ \mu_1^{NP} \\ c^{NP} \\ c_r^{NP} \\ I^{NP} \end{pmatrix} = \begin{pmatrix} \frac{\gamma a^2}{\delta VS} \\ \frac{\gamma a}{\delta S} \\ \frac{a}{\delta} \\ 0 \\ 0 \end{pmatrix}. \qquad (NP) $$

If $\mu_1^{NP} \le m$, the threshold for the activation of cytokines, this defines an equilibrium solution with a residual tumor and free of protumor immune cells. Next, let us introduce

$$ Q = \frac{k_r a}{\delta \tau}, $$

and let

$$ \mu_1^P = \frac{\frac{\gamma a}{\delta} - Qm - \frac{k_c a}{\gamma_r \delta} Qm}{S - Q - \frac{k_c a}{\gamma_r \delta} Q}. \qquad (14) $$

be the tumor mass at equilibrium with the presence of protumor immune cells. Indeed, if $\mu_1^P > m$ another equilibrium solution is given by

$$ \begin{pmatrix} \mu_0^P \\ \mu_1^P \\ c^P \\ c_r^P \\ I^P \end{pmatrix} = \begin{pmatrix} \frac{a}{V(1 + b_1 c_r^P)} \mu_1^P \\ \mu_1^P \\ \frac{a}{\delta} \\ Q \frac{\psi(\mu_1^P)}{\gamma_r} \\ \psi(\mu_1^P) / \tau \end{pmatrix} \qquad (P) $$
That the definition of this unhealthy state makes sense requires that the right hand side in (14) is positive. It means that

$$ \begin{array}{l} \text{either } S > Q\left(1 + \frac{k_c a}{\gamma_r \delta}\right) > 0 \text{ and } \frac{\gamma a}{\delta} > Qm\left(1 + \frac{k_c a}{\gamma_r \delta}\right) \\ \text{or } S < Q\left(1 + \frac{k_c a}{\gamma_r \delta}\right) \text{ and } \frac{\gamma a}{\delta} < Qm\left(1 + \frac{k_c a}{\gamma_r \delta}\right), \end{array} \qquad (15) $$

imposing constraints on the parameters. Let us discuss the possibility of obtaining the different equilibrium states $\mu_1^{NP}$, $\mu_1^p$, depending on the ratio $\frac{a}{\delta}$ between the tumor division rate $a$ and the strength of the immune response $\delta$. It measures the competitiveness between the tumor and the antitumor immune cells. We thus study respectively the sign of $\mu_1^{NP} - m$ and $\mu_1^p - m$. On the one hand,

$$ \mu_1^{NP} - m < 0 \text{ if and only if } \frac{a}{\delta} < \frac{mS}{\gamma}, \qquad (16) $$

and on the other hand

$$ \mu_1^p - m = \frac{\gamma \frac{a}{\delta} - mS}{S - Q_1 \frac{a}{\delta} - Q_2 \left(\frac{a}{\delta}\right)^2}, \qquad (17) $$

where

$$ Q_1 = \frac{k_r}{\tau} \quad \text{and} \quad Q_2 = \frac{k_c k_r}{\gamma_r \tau}. $$

Let us denote by

$$ X_2 = \frac{-\frac{Q_1}{Q_2} + \sqrt{\left(\frac{Q_1}{Q_2}\right)^2 + \frac{4S}{Q_2}}}{2}, \qquad (18) $$

the non-negative root of the denominator in (17). By analyzing the sign of (17), we distinguish two cases

* if $m < \frac{\gamma X_2}{S}$ (relatively small critical mass),

$$ \mu_1^p - m > 0 \text{ if and only if } \frac{a}{\delta} \in \left(\frac{mS}{\gamma}, X_2\right) \qquad (19) $$

and

$$ \mu_1^p - m < 0 \text{ if and only if } \frac{a}{\delta} \in \left(0, \frac{mS}{\gamma}\right) \cup (X_2, +\infty). \qquad (20) $$
This is summarized in the following table

<table><thead><tr><th>a/δ</th><th>0</th><th>mS/γ</th><th>X<sub>2</sub></th><th>+∞</th></tr></thead><tbody><tr><td>μ<sub>1</sub><sup>P</sup> - m</td><td>&lt; 0</td><td>0</td><td>&gt; 0</td><td>&lt; 0</td></tr><tr><td></td><td></td><td>μ<sub>1</sub><sup>P</sup> admissible</td><td></td><td></td></tr><tr><td>μ<sub>1</sub><sup>NP</sup> - m</td><td>&lt; 0</td><td>0</td><td>&gt; 0</td><td>&gt; 0</td></tr><tr><td></td><td>μ<sub>1</sub><sup>NP</sup> admissible</td><td></td><td></td><td></td></tr></tbody></table>

There is no admissible equilibrium when $a/δ \in [X_2, +\infty)$: the aggressiveness of the tumor is strong and the tumor mass certainly blows up.

* if $m > \frac{\gamma X_2}{S}$, we have

<table><thead><tr><th>a/δ</th><th>0</th><th>X<sub>2</sub></th><th>mS/γ</th><th>+∞</th></tr></thead><tbody><tr><td>μ<sub>1</sub><sup>P</sup> - m</td><td>&lt; 0</td><td>&gt; 0</td><td>0</td><td>&lt; 0</td></tr><tr><td></td><td></td><td>μ<sub>1</sub><sup>P</sup> admissible</td><td></td><td></td></tr><tr><td>μ<sub>1</sub><sup>NP</sup> - m</td><td>&lt; 0</td><td>&lt; 0</td><td>0</td><td>&gt; 0</td></tr><tr><td></td><td>μ<sub>1</sub><sup>NP</sup> admissible</td><td></td><td></td><td></td></tr></tbody></table>

Again, there is no admissible equilibrium when $a/δ$ becomes large.

**Discussion on the stability of the equilibrium points.** The Jacobian matrix evaluated at the healthy state (H) reads

$$J^H = \begin{pmatrix} a & 0 & 0 & 0 & 0 \\ V & 0 & 0 & 0 & 0 \\ 0 & S & -\gamma & 0 & 0 \\ 0 & 0 & 0 & -\gamma_r & 0 \\ 0 & 0 & 0 & 0 & -\tau \end{pmatrix}$$

Since $a > 0$, $J^H$ has a positive eigenvalue and the healthy state is linearly unstable. The equilibrium state (NP), which is free of protumor immune cells, corresponds to the unhealthy state in [24]. The Jacobian matrix at this state reads

$$J^{NP} = \begin{pmatrix} 0 & 0 & -\frac{\gamma a^2}{VS} & 0 & 0 \\ V & -a & -\frac{\gamma a}{S} & 0 & 0 \\ 0 & S & -\gamma & -k_c \frac{a}{\delta} & -k_r \frac{a}{\delta} \\ 0 & 0 & 0 & -\gamma_r & k_r \frac{a}{\delta} \\ 0 & 0 & 0 & 0 & -\tau \end{pmatrix}$$
and its characteristic polynomial is

$$p(\lambda) = -(\gamma_r + \lambda)(\tau + \lambda)(\lambda^3 + (\gamma + a)\lambda^2 + 2a\gamma\lambda + \gamma a^2).$$

As in [24], we distinguish two cases, which depends on the ratio $\frac{\gamma}{a}$. The ratio compares the death rate of the antitumor immune cells to the tumor cells division rate. We get

* if $\gamma > 4a$, the eigenvalues of $J^{NP}$ are given by

$$\lambda_1 = -\gamma_r, \quad \lambda_2 = -\tau, \quad \lambda_3 = -a, \quad \lambda_4 = \frac{1}{2p} \left( -\sqrt{\gamma(\gamma - 4a)} - \gamma \right),$$

$$\lambda_5 = \frac{1}{2} \left( \sqrt{\gamma(\gamma - 4a)} - \gamma \right).$$

They are all real and negative.

* if $\gamma < 4a$, the eigenvalues are given by

$$\lambda_1 = -\gamma_r, \quad \lambda_2 = -\tau, \quad \lambda_3 = -a, \quad \lambda_4 = \frac{1}{2} \left( -i \sqrt{\gamma(\gamma - 4a)} - \gamma \right),$$

$$\lambda_5 = \frac{1}{2} \left( i \sqrt{\gamma(\gamma - 4a)} - \gamma \right).$$

They all have a negative real part.

Therefore, when admissible ($\mu_1^{NP} < m$), the unhealthy state with no protumor immune cells is always linearly stable. In addition, the ratio $\gamma/a$ discriminates between a damped behavior, and an oscillatory behavior. As observed in [24], the greater the cell division, the faster the oscillations of the tumor mass $\mu_1$. For the equilibrium with protumor cells, the Jacobian matrix becomes full and we are not able to find such explicit formula for the eigenvalues.

We observe on numerical grounds that in the case where the equilibrium states (NP) and (P) coexist ($X_2 < \frac{a}{\delta} < \frac{mS}{\gamma}$), there is either the formation of an equilibrium free of protumor cells, namely as in (NP), or the blow up of the tumor mass due to the activation of protumor immune cells; when $\frac{a}{\delta} > \frac{mS}{\gamma}$, the tumor mass always blows up. Typical results are depicted in Fig 3: when the control occurs, the concentration of antitumor immune cells tends to the equilibrium value $a/\delta$ (Fig 3a)); otherwise, it reaches a constant state below the equilibrium value (Fig 3b–3d). The simulations also show that the control of the tumor is very sensitive to the strength of the source of naive immune cells and to the division rate of the tumor cells. It is worth remarking that the damping of the tumor mass towards an equilibrium can be restored by strengthening the activation law of effector immune cells for large tumor masses, for instance by using $g(\mu_1) = \mu_1^2$, see Fig 4.

## Existence of equilibrium phases for the PDE model

We now go back to the full PDE model (12a)–(12g) and we explain why an equilibrium can be expected in certain circumstances. The analysis and simulations carried out in [24] for the model without protumor immune cells reveal the existence and stability of a cancer-persistent equilibrium. It turns out that the equilibrium phase corresponds to the situation where the
<figure><img src="image_6.png" /><figcaption>**Fig 3. Typical behavior of the solutions of the simplified ODE system (13).** Data: $V = 0.616$, $\delta = 1$, $S = 1.5$, $k_r = 1.25$, $k_c = 0.1$, $m = 2$. x-axis: time, y-axis: $\mu_1$, mass of the tumor (plain, left axis), and $c$, strength of the active immune cells (dashed, right axis); the expected equilibrium value is $a/\delta$ (dotted line in fig. (a)). When $a$ is small enough, the tumor is controlled and the concentration of immune cells tends to the equilibrium value (fig. (a)). For larger $a$, the concentration of immune cells tends to a constant smaller than the expected equilibrium and the tumor mass blows up (fig. (b)-(d)). Mind that the time scale differs: in (b)-(d) the tumor mass blows rapidly, the larger $a$, the faster the blow up.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0259291.g003

death rate induced by the effector immune cells precisely counterbalances the natural exponential growth of the tumor cell population.

Indeed, it is known that the growth-fragmentation operator admits an eigenpair $(\lambda, N)$ satisfying

$$ \left\{ \begin{array}{ll} \partial_z(VN) - Q(N) + \lambda N = 0 & \text{for } z \ge 0, \\ N(0) = 0, \quad N(z) > 0 & \text{for } z > 0, \quad \int_0^{+\infty} N(z) \, dz = 1, \quad \lambda > 0. \end{array} \right. \qquad (21) $$

We refer the reader to [35] for a detailed analysis of this eigenproblem. When the action of the immune system is neglected, namely $m(n, c) = 0$ and $b_1 = 0$ in (12a), the population of tumor cells grows exponentially fast and its size-distribution is governed by the eigenfunction: $n(t, z) \sim t \to \infty e^{\lambda t} N(z)$, see [26–28]. Equilibrium occurs when the death rate due to the effector cells reaches the eigenvalue. Namely, the concentration $C$ of cytotoxic effector cells at
<figure><img src="image_7.png" /><figcaption>Fig 4. Typical behavior of the solutions of the simplified ODE system (13) with g(μ₁) = μ₁². Data: V = 0.616, δ = 1, S = 1.5, kᵣ = 1.25, kₑ = 0.1, m = 2. Fig. (a)-(d): x-axis: time, y-axis: μ₁, mass of the tumor (plain, left axis), and c, strength of the active immune cells (dashed, right axis), the expected equilibrium value is a/δ (dotted line in fig. (a) and (b)). Fig. (e)-(f): x-axis: time, y-axis: c strength of the active immune cells (dashed, right axis), cᵣ, concentration of protumor cells (square, left axis). The equilibrium is restored with a strengthened recruitment of antitumor immune cells for large tumor masses. The behavior of the concentration of protumor cells follows the behavior of the tumor mass (compare (c)-(e) or (d)-(f)).</figcaption></figure>

https://doi.org/10.1371/journal.pone.0259291.g004

equilibrium should satisfy

$$ \int_{\Omega} \delta(x) C(x) \, dx = \lambda. \qquad (22) $$
In turn, the definition of the concentration of activated immune cells by means of stationary convection-diffusion-reaction equations defines implicitly the total tumor mass at equilibrium. This intuition is made precise by the following statement which extends the results of [24] to the situation with protumor immune activities and inflammatory signals.

**Theorem 1** Let $\Phi$ be the solution of

$$\nabla_x \cdot (\kappa \nabla_x \Phi) = \sigma,$$

endowed with the homogeneous Neumann boundary condition. There exists $\ell_* > 0$ such that for any $0 < \ell < \ell_*$, there exists a unique $\bar{\mu}_1(\ell) > 0$ and $(C_{\bar{\mu}_1(\ell)}, C_{r,\bar{\mu}_1(\ell)}, I_{\bar{\mu}_1(\ell)})$, solution of the stationary equations

$$\gamma C + k_r I \theta C + k_c C_r - f(\bar{\mu}_1) \nabla_x \cdot (C \chi \nabla_x \phi) - \nabla_x \cdot (D \nabla_x C) = g(\bar{\mu}_1) S, \quad (23a)$$

$$\gamma_r C_r - f(\bar{\mu}_1) \nabla_x \cdot (C_r \chi \nabla_x \phi) - \nabla_x \cdot (D \nabla_x C_r) = I(S_r + k_r \theta C), \quad (23b)$$

$$I = \frac{\psi(\bar{\mu}_1)}{\tau}, \quad (23c)$$

$$C|_{\partial\Omega} = 0, \quad C_r|_{\partial\Omega} = 0, \quad (23d)$$

satisfying $\int_{\Omega} \delta C dx = \ell$.

**Proof.** We adapt the arguments from [24], taking into account protumor immune cells and cytokines. We start by introducing the mapping

$$\mathcal{F} : (\ell, \mu_1) \in [0, \infty) \times [0, \infty) \mapsto \int_{\Omega} \delta C_{\mu_1} dx - \ell$$

where $C_{\mu_1}$ is the solution of (23a) associated to $\mu_1$. We are searching for the zeroes of $\mathcal{F}$. Clearly, when $\mu_1 = 0, C_0 = 0, I_0 = 0, C_{r,0} = 0$ satisfies (23a)-(23d), together with the constraint $\int \delta C_0 dx = 0$, so that $\mathcal{F}(0, 0) = 0$. Next, we have $\partial_{\mu_1} \mathcal{F}(\ell, \mu_1) = \int_{\Omega} \delta C'_{\mu_1} dx$, where $C'_{\mu_1}$ is defined by the system

$$\begin{align*}
& (\gamma + k_r I \theta + k_c C_r) C' + (k_r I' \theta + k_c C'_r) C - \nabla_x \cdot (D \nabla_x C') - f(\mu_1) \nabla_x \cdot (C' \nabla_x \Phi) \\
& = g'(\mu_1) S + f'(\mu_1) \nabla_x \cdot (C_{\mu_1} \nabla_x \Phi) \\
& \gamma_r C'_r - \nabla_x \cdot (D \nabla_x C'_r) - f(\mu_1) \nabla_x \cdot (C'_r \nabla_x \Phi), \\
& = I'(S_r + k_r \theta C) + I k_r \theta C' + f'(\mu_1) \nabla_x \cdot (C_{\mu_1} \nabla_x \Phi), \\
& I' = \frac{\psi'(\mu_1)}{\tau}.
\end{align*}$$

With $\mu_1 = 0$, the right hand side of the equation for $C'_r$ vanishes and we get $C'_r = 0$. In the right hand side of the equation for $C'$, $g'(0)S \neq 0$ is non negative and the maximum principle for elliptic equations tells us that $C'_0 > 0$. It follows that $\partial_{\mu_1} \mathcal{F}(0, 0) = \int_{\Omega} \delta C'_0 dx > 0$. We can thus apply the implicit function theorem: there exists $\ell_* > 0$ and a mapping $\bar{\mu}_1 : \ell \in [0, \ell_*) \mapsto \bar{\mu}_1(\ell)$ such that $\mathcal{F}(\ell, \bar{\mu}_1(\ell)) = 0$ holds for any $\ell \in [0, \ell_*)$. We have

$$\partial_{\ell} \mathcal{F}(\ell, \bar{\mu}_1(\ell)) + \bar{\mu}_1'(\ell) \partial_{\mu_1} \mathcal{F}(\ell, \bar{\mu}_1(\ell)) = -1 + \bar{\mu}_1'(\ell) \partial_{\mu_1} \mathcal{F}(\ell, \bar{\mu}_1(\ell)) = 0$$

with $\partial_{\mu_1} \mathcal{F}(0, 0) > 0$. Hence, $\ell \mapsto \bar{\mu}_1(\ell)$ is increasing on the neighborhood of $\ell = 0$, and it thus takes positive values.
Theorem 1 applies directly when the action of protumor immune cells is neglected on the cell-division equation (namely assuming $b_1 = 0$): the eigenstate $(\lambda, N)$ is defined a priori and the statement directly defines the equilibrium phase with the constraint (22). The statement involves a smallness assumption and it justifies the existence of equilibria with small tumor masses. In the case of the binary division model with a constant division rate $a$ and a constant growth rate $V$, the smallness assumption is equivalent to a smallness assumption on the division rate $a$. The numerical simulations indicate a quite robust property [24, 36] and the smallness assumption could be only technical.

The analysis of the full model accounting for a protumor action on the tumor growth rate is much more involved: $V(z)$ is multiplied by the factor $(1 + \beta(t))$ with

$$ \beta(t) = \int b_1(y)c_r(t, y) dy. $$

As $t \to \infty$, we expect that $c_r(t, x)$ and $c(t, x)$ admits limits so that $\int \delta c(t, x) dx$ tends to some $\lambda > 0$ and $\beta(t)$ tends to an asymptotic value $\beta_\infty$ while the size-distribution of the tumor cells is described by an eigenpair $(\lambda, N)$. However, in general the eigenvalue $\lambda$ depends on the value of $\beta_\infty$, which induces a stronger coupling between the unknowns. The case where both $a$ and $V$ are constant is specific and allows us to strengthen the intuition. In this case, the leading eigenvalue $\lambda = a$ does not change when $V$ is replaced by $V(1 + \beta_\infty)$; only the profile is rescaled into

$$ N_{\beta_{\infty}}(z) = N_1\left(\frac{z}{1 + \beta_{\infty}}\right), $$

where the profile $N_1$ is known. As $\beta_\infty$ increases, the asymptotic size-distribution of tumor cells contains larger cells, see [24, Fig. 3]. For general coefficients, the leading eigenstate of the growth-fragmentation operator is not explicitly known but it can be computed numerically using the power method designed in [36]: Fig 5 gives the value of the eigenvalue $\lambda$ and shows the profiles of the eigenfunctions when the growth rate $z \mapsto V(z)$ follows the Gompertz law (1) and the division rate is given by $z \mapsto a(z) = a1_{z_0 \le z < \infty}$ for some $z_0 > 0$. Changing $V$ now modifies both the profile and the eigenvalue.

For the same reasons, when the protumor cells modify the growth rate of the tumor cells, we cannot use as such the method designed in [36] to predict the equilibrium state. Nevertheless, the numerical simulations of the initial-boundary value problem highlight the following features:

* when the critical mass $m$ is positive, either the steady state is free of protumor immune cells ($c_r \equiv 0$) and the tumor is controlled by the antitumor immune cells or both the tumor mass $\mu_1$ and the concentration of protumor immune cells $c_r$ blow up. The former occurs for small division rates, the latter is observed with more aggressive tumors. Furthermore, when the tumor growth is controlled, we can check that the asymptotic concentration of antitumor immune cells ($x \mapsto C(x)$) satisfies (22) (see the numerical experiments detailed below)

* when the critical mass is equal to zero ($m = 0$), either we observe an equilibrium state containing residual protumor immune cells or the tumor mass $\mu_1$ and the concentration of protumor immune cells $c_r$ blow up.
<figure><img src="image_8.png" /><figcaption>**Fig 5.** Shapes of the leading eigen-function, solution of the growth-division Eq (21). x-axis: z, size of the tumor cells, y-axis: number of tumor cells at the final time. The profiles are obtained for the binary cell division operator (3) with a constant division rate *a* > 0. and the Gompertz law (1) for the growth rate. The shape is driven by the ratio *a/z* where r = (1 + β∞) is the intrinsic growth rate of the tumor cells.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0259291.g005

Results

This section is devoted to numerical investigations of the PDE system (12a)–(12g). In the
absence of protumor immune effects, the establishment of an equilibrium seems to be a quite
robust feature of the model [24, 36]. This raises the issue to determine whether or not the pro-
tumor effects can break the equilibrium and lead to the escape phase and what are the key
parameters to such a transition. We start by discussing this issue, dealing with quite complex
growth-division coefficients, and equally with simpler constant coefficients, for which the
details of the equilibrium are explicitely known, in order to fully validate the intuition on
numerical grounds. The next challenge is to incorporate in the model the effects of treatments
in order to understand how they can help in avoiding the escape phase and restore an equilib-
rium. To this end, we shall need to introduce additional equations describing exhausted
immune cells the antitumor activity of which can be restored by the treatments. Another possi-
bility is to reduce the protumor effects of the cytokine signals. We are going to discuss the
effects of both strategy, bringing out the role of the administration dose and the starting time
of the treatment. Finally, we show that combining the two approaches is highly beneficial, pro-
viding a better control of the tumor mass, with reduced dose and later administration.
<table><caption>Table 2. Data for the simulations.</caption><thead><tr><th>A</th><th>ξ<sup>2</sup></th><th>A<sub>σ</sub></th><th>ξ<sub>σ</sub><sup>2</sup></th><th>a</th><th>V</th><th>r</th><th>b</th><th>χ</th><th>S</th><th>γ</th></tr></thead><tbody><tr><td>1</td><td>0.02</td><td>0.002</td><td>0.05</td><td>0.8</td><td>0.616</td><td>0.616</td><td>10</td><td>0.864</td><td>5</td><td>0.18</td></tr></tbody></table>

<table><thead><tr><th>A<sub>b<sub>1</sub></sub></th><th>ξ<sub>b<sub>1</sub></sub></th><th>A<sub>θ</sub></th><th>ξ<sub>θ</sub></th><th>k<sub>r</sub></th><th>k<sub>c</sub></th><th>τ</th><th>η</th><th>m</th></tr></thead><tbody><tr><td>10<sup>-6</sup></td><td>0.03</td><td>0.2</td><td>0.3</td><td>1</td><td>1</td><td>1</td><td>1</td><td>2</td></tr></tbody></table>

https://doi.org/10.1371/journal.pone.0259291.t002

The mathematical model (12a)–(12g) accounts for the action of both anti- and protumor immune cells shaping tumor growth kinetics. Numerical experiments were used to challenge the model. We perform the numerical simulations considering the binary division operator (3) with a constant division rate $a > 0$. For details on the numerical methods, which are based on finite volume discretizations, we refer the reader to [24, 36]. According to the framework in [24], we assume that the tumor is located at the origin of the computational domain $\Omega$, which here is the unit ball of $\mathbb{R}^2$, and we use the following definitions

$$ \delta(x) = \frac{A}{\xi \sqrt{2\pi}} \exp\left(-\frac{|x|^2}{2\xi^2}\right), \quad b_1(x) = \frac{A_{b_1}}{\xi_{b_1} \sqrt{2\pi}} \exp\left(-\frac{|x|^2}{2\xi_{b_1}^2}\right). \quad (24) $$

For defining the source term of the chemoattractant potential and the form function $\theta$ we also use the following Gaussian profiles:

$$ \sigma(x) = \frac{A_{\sigma}}{\xi_{\sigma} \sqrt{2\pi}} \exp\left(-\frac{|x|^2}{2\xi_{\sigma}^2}\right), \quad \theta(x) = \frac{A_{\theta}}{\xi_{\theta} \sqrt{2\pi}} \exp\left(-\frac{|x|^2}{2\xi_{\theta}^2}\right) \quad (25) $$

In what follows, we use the Michaelis-Menten law (7) and we simply set $g(\mu_1) = \mu_1$. For the simulations, we shall use the following data, otherwise explicitly stated: the initial data are $(c_0(x), c_{r,0}(x)) = (0, 0)$ and $n_0(z) = 1_{0.125 \le z \le 5}$. The parameters are given in Table 2. We use dimensionless equations, without addressing precisely the issue of parameter calibration (see [24, 36] for further details on this issue). For the numerical experiments, we consider the source of antitumor cells, that contains T cells recruited in specific sites like the lymph nodes, as well as NK, N1, M1 taken from the circulation, as space-homogeneous. In contrast, we work with an heterogeneous source $S_r$ of protumor cells, distant from the tumor site (see Fig 6; in practice, we use 1/5 of the source depicted in this figure). This assumption corresponds to a privileged recruitment of protumor cells in specific sites, like the bone marrow.

## Emergent qualitative features: Promotumor cells are determinant for the shift to the escape phase

In order to assess the qualitative behavior of the system (12a)–(12g), we perform a set of simulations with assumptions that incorporate several biological effects, namely

* $z \mapsto V(z)$ obeys the Gompertz law (1);

* the division rate $z \mapsto a(z) = a1_{z_0 \le z < \infty}$ for some $z_0 > 0$ (for the numerical tests we set $z_0 = 1$) vanishes for the smallest cells;

* the presence of the protumor immune cells promotes tumor growth with $b_1 \neq 0$.
<figure><img src="image_9.png" /><figcaption>**Fig 6.** Source *S*, of protumor cells. x-y axis correspond to the space coordinates.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0259291.g006

We generically observe two behaviors: either an equilibrium state establishes, with a residual tumor and free of protumor immune cells, or the immune system fails in controlling the tumor, with a significant concentration of protumor immune cells at the center of the domain (like in Fig 7 which also clearly illustrates how the antitumor resources is stemmed in the vicinity of the tumor) and the tumor mass blows up. This rough conclusion should be nuanced: the threshold *m* certainly plays a critical role. With *m* = 0, the worst situation since protumor immune cells are immediately activated, we can find equilibria with the three types of cells (tumor cells, antitumor and protumor immune cells). Such equilibria occur with quite small values of the cell division rate *a*; increasing *a* leads to an escape phase. It is likely that similar phenomena occur with positive threshold *m* and very small *a*'s. When a steady regime establishes, we check on numerical ground, by evaluating the eigenpair of the growth-division equation, that the asymptotic antitumor cell concentration is consistent with the constraint (22).

We make the parameters vary in order to discuss the influence of their value on the behavior of the system. We only modify one quantity at a time, the others being kept as in Table 2.

* *Tumor aggressiveness.* As indicated in [24], by increasing the division rate *a*, we make the tumor more aggressive. The results with the PDE system are consistent with this intuition and the observations made above on the simplified ODE system: for small division rates, the mass of the tumor is rapidly damped, while the tumor escapes the control of the immune
<figure><img src="image_10.png" /><figcaption>Fig 7. Space distribution of immune cells. Antitumor cells c (left) and protumor cells c_r (right). x-y axis: space coordinates. The concentration of protumor cells is higher in the vicinity of the tumor, where the concentration of antitumor cells has been depleted. The picture has been obtained from the simulation of the system (12a)–(12g), at final time t = 2.23, with a = 4.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0259291.g007

system as a increases, see Fig 8. When control occurs, protumor immune cells can be activated in the transient states, but insufficiently to counterbalance the effector immune response. Therefore, the concentration of protumor immune cells decreases to zero, while

$$ \bar{\mu}_c(t) = \int_{\Omega} \delta(x) c(t, x) dx \qquad (26) $$

tends to the leading eigenvalue λ of the growth-division equation, see Figs 8 and 14a. When the tumor is more aggressive, it recruits more protumor immune cells: in turn, the action of these cells restrains the concentration of antitumor immune cells which remains below the expected equilibrium value, eventually favoring the tumor escape. The specific value of a for the bifurcation from a controlled tumor growth to the escape state depends on the critical mass m: the smaller the critical mass m, the smaller the critical division rate a. Figs 9 and 10 illustrate the space organization of the immune response in a situation where the control of the tumor is lost (the data are the same as in Fig 8(c)). Antitumor cells are activated and recruited from a homogeneous bath while protumor cells come from specific spots, as it appears clearly in the earliest pictures. Both are directed towards the tumor, at the center of the domain, by the chemotactic potential. In the situation depicted in these figures, the protumor action becomes strong in the vicinity of the tumor where it inhibits the antitumor response.

* *Efficiency of the immune response.* The immune response is enhanced by increasing A, that measures the strength of the immune cells against the tumor cells, see (24), or the source S of effector antitumor cells. For small values of these parameters, the tumor escapes the control of the immune response, see Figs 11 and 12. While the tumor mass µ_1 keeps growing, the immune strength µ_c remains limited and cannot balance the growth rate of the tumor. For large S or A, the equilibrium establishes, with µ_c(t) tending asymptotically to counterbalance the tumor growth rate λ.

**Further numerical validation.** The case where the growth rate V and the division rate a are constant is less relevant biologically, but it can be used to check the properties of the model
<figure><img src="image_11.png" /><figcaption>**Fig 8.** Simulation of the PDE model (12a)–(12g) for several values of the division rate *a*. Evolution of the tumor mass $\mu_1$ (plain, left axis), and of the immune strength $\bar{\mu}_c$ (defined in (26), dashed, right axis); expected equilibrium value $\lambda$ (dotted line in fig. (a) and (b)). When *a* is small enough an equilibrium is reached with $\bar{\mu}_c$, tending to the eigenvalue $\lambda$ of the cell-division equation, and a residual tumor mass (fig. (a) and (b)). For larger *a*'s the tumor escapes the control of the immune system and its mass bows up (fig. (c) and (d)): mind the different time scales, since the blows up is faster as *a* increases.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0259291.g008

<figure><img src="image_12.png" /><figcaption>**Fig 9.** Simulation of the PDE model (12a)–(12g): Evolution of the size repartition of the antitumor immune cells. The test corresponds to Fig 8(c), with *a* = 4. x-y axis: space coordinates. The antitumor cells are recruited from a homogeneous bath of non-activated cells, and then are directed towards the tumor by chemotactism. However, the protumor activities lower their concentration in the vicinity of the tumor (mind that the scale of the color map changes with time).</figcaption></figure>

https://doi.org/10.1371/journal.pone.0259291.g009
<figure><img src="image_13.png" /><figcaption>**Fig 10.** Simulation of the PDE model (12a)–(12g): Evolution of the size repartition of the protumor immune cells.
The test corresponds to Fig 8(c), with *a* = 4. The protumor cells are activated in 3 specific sites, and then are directed
towards the tumor by chemotactism. Eventually there is a high concentration of protumor cells in the vicinity of the
tumor (mind that the scale of the color map changes with time).</figcaption></figure>

https://doi.org/10.1371/journal.pone.0259291.g010

<figure><img src="image_14.png" /><figcaption>(a) $S = 0.002$</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(c) $S = 20$</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(d) $S = 200$</figcaption></figure>

<caption>**Fig 11.** Simulation of the PDE model (12a)–(12g) for several values of the source of effector immune cells S.
Evolution of the tumor mass $\mu_1$ (plain, left axis), and of the immune strength $\bar{\mu}_c$ (defined in (26), dashed, right axis).
When S is large enough an equilibrium is reached with $\bar{\mu}_c$ tending to the eigenvalue $\lambda$ of the cell-division equation
(dotted line in fig. (c) and (d)), and a residual tumor mass (fig. (c) and (d)). For smaller S's the tumor escapes the
control of the immune system and its mass blows up (fig. (a) and (b)).</caption>

https://doi.org/10.1371/journal.pone.0259291.g011
<figure><img src="image_17.png" /><figcaption>**Fig 12.** Simulation of the PDE model (12a)–(12g) for several values of the immune strength A. Evolution of the tumor mass $\mu_1$ (plain, left axis), and of the immune strength $\bar{\mu}_c$ (defined in (26), dashed, right axis); expected equilibrium value $\lambda$ (dotted line in fig. (c) and (d)). When $A$ is large enough an equilibrium is reached with $\bar{\mu}_c$ tending to the eigenvalue $\lambda$ of the cell-division equation, and a residual tumor mass (fig. (c) and (d)). For smaller $A$'s the tumor escapes the control of the immune system and its mass bows up (fig. (a) and (b)).</figcaption></figure>

https://doi.org/10.1371/journal.pone.0259291.g012

and of the numerical procedures since the eigenpair $(\lambda, N)$ is explicitly known in this case, see [37, 38] and [28, Lemma 4.1]: in fact we have $\lambda = a$. Still with the purpose of assessing the model and the numerical method on simple basis, it is relevant to perform simulations by assuming also that the protumor immune cells do not enhance the growth rate of the tumor cells ($b_1 = 0$). It means that the protumor effect is limited to the suppression of antitumor capacities. This case is biologically questionable, but, as explained above, we have an intuition a priori on the details of the possible equilibrium state, and we can directly check whether or not the large time behavior is described by this expected equilibrium.

Beyond the validation, it is remarkable that the qualitative conclusions do not substantially change when comparing the results to the more complete situation dealt before: compare Figs 8 and 11–15, respectively. The evolution of the concentration of protumor immune cells in Fig 16 is equally quite generic: we clearly see the difference between an equilibrium case where this concentration vanishes as time grows, and an escape case where it blows up. These observations show the robustness of the model in describing the equilibrium vs escape phenomena and this is very reassuring for further investigations with clinical data, as in [36]. In this direction, identifying the parameters of the equations is a critical issue. It can be interesting, based
<figure><img src="image_18.png" /><figcaption>**Fig 13.** Validation of the PDE model ($b_1 = 0$, $V$ and $a$ constant). Evolution of the tumor mass $\mu_1$ (plain, left axis), and of the immune strength $\bar{\mu}_c$ (defined in (26), dashed, right axis) for several values of the division rate $a$; expected equilibrium value $\lambda$ (dotted line). When $a$ is small enough an equilibrium is reached with $\bar{\mu}_c$ tending to the eigenvalue $\lambda$ of the cell-division equation, and a residual tumor mass (fig. (a) and (b)). For larger $a$'s the tumor escapes the control of the immune system and its mass bows up (fig. (c) and (d)): mind the different time scales, since the blow up is faster as $a$ increases.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0259291.g013

on the present observations, to neglect some phenomena which can only marginally affect the dynamics, while potentially introducing a set of unknown parameters.

## Effect of immunotherapy strategies

Now that we have validated a robust mathematical model of tumor growth which takes into account the contribution of anti and protumor immune cells, we use it to compare the effects of two immunotherapy treatments targeting these immune cells with opposed functions. To this end, we bear in mind that a proportion of the effector cells are just inhibited by immunosuppressive mechanisms; in other words they are not destroyed: they become hyporesponsive. However, they can be re-activated by specific treatments. The restoration of the antitumoral activity of effector T cells can be obtained by using Immune Checkpoint Inhibitors, like anti-PD-1 or anti-CTLA4 antibodies [14, 15]. The infusion of CAR-T and CAR-NK cells can also mimic such rescue [39]. A second strategy is to reduce the recruitment of protumor immune cells by blocking infiltration of MDSCs (anti-CXCR2, cMet) [40, 41] and Tregs (anti-CD25) [13]. We discuss the effect of these approaches, illustrated in Fig 17, independently and we also consider the combination of the two treatments compared to the mono-therapies.
<figure><img src="image_19.png" /><figcaption>**Fig 14.** Validation of the PDE model ($b_1 = 0$, $V$ and $a$ constant). Evolution of the tumor mass $\mu_1$ (plain, left axis), and of the immune strength $\bar{\mu}_c$ (dashed, right axis) for several values of the source of immune cells $S$; expected equilibrium value $\lambda$ (dotted line in fig. (c) and (d)). When $S$ is large enough an equilibrium is reached with $\bar{\mu}_c$ tending to the eigenvalue $\lambda$ of the cell-division equation, and a residual tumor mass (fig. (c) and (d)). For smaller $S$'s the tumor escapes the control of the immune system and its mass bows up (fig. (a) and (b)).</figcaption></figure>

https://doi.org/10.1371/journal.pone.0259291.g014

<figure><img src="image_20.png" /><figcaption>**Fig 15.** Validation of the PDE model ($b_1 = 0$, $V$ and $a$ constant). Evolution of the tumor mass $\mu_1$ (plain, left axis), and of the immune strength $\bar{\mu}_c$ (dashed, right axis) for several values of the source of immune strength $A$; expected equilibrium value $\lambda$ (dotted line in fig. (c) and (d)). When $A$ is large enough an equilibrium is reached with $\bar{\mu}_c$ tending to the eigenvalue $\lambda$ of the cell-division equation, and a residual tumor mass (fig. (c) and (d)). For smaller $A$'s the tumor escapes the control of the immune system and its mass bows up (fig. (a) and (b)).</figcaption></figure>

https://doi.org/10.1371/journal.pone.0259291.g015
<figure><img src="image_21.png" /><figcaption>**Fig 16. Validation of the PDE model (b₁ = 0, V and a constant).** Evolution of the averaged protumor cells concentration $\mu_{c_r} = \int_{\Omega} c_r dx$ (square, left axis), and of the immune strength $\bar{\mu}_c$ (dashed, right axis) for several values of the division rate a. For small a (fig. (a)), the concentration of antitumor immune cells reaches the equilibrium, and the concentration of protumor cells is damped; for large a (fig. (b)), the equilibrium does not establish and the concentration of protumor immune cells keeps growing.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0259291.g016

**Therapy based on the reactivation of exhausted antitumor immune cells.** To consider the action of treatments boosting the immune response against the tumor, we introduce the concentration $c_a$ of exhausted cells. In order to describe the restoration mechanism, we add the following equation to the model:

$$ \partial_t c_a + \nabla_x \cdot (c_a \chi \nabla_x \phi - D \nabla_x c_a) = \alpha k_c c_c - \gamma_a c_a, \quad (27) $$

where the parameter $0 < \alpha < 1$ describes the proportion of effector T cells that become hyporesponsive under the action of the protumor cells (see (28a)). Note that the death rate $\gamma_a$ could be significantly larger than the original death rate $\gamma$: it is believed that exhausted T cells have a shorter life time. Next, the effect of treatments able to restore the antitumor activity of the exhausted immune cells is described by a time-dependent function $t \mapsto T_1(t)$. It is assumed to be proportional to the drug concentration in the TME. Consequently, the dynamic is governed
<figure><img src="image_22.png" /><figcaption>Fig 17. Schematic overview of the action of treatments on the immune response. AIC: antitumor cells, PIC: protumor cells, EIC: exhausted immune cells.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0259291.g017

by the following system, which extends (12a)-(12g)

$$ \partial_t c + \nabla_x \cdot (c\chi\nabla_x\phi - D\nabla_x c) = g(\mu_1)S + T_1c_a - \gamma c - k_r I\theta c - k_c cc_r, \quad (28a) $$

$$ \partial_t c_r + \nabla_x \cdot (c_r \chi \nabla_x \phi - D \nabla_x c_r) = I(S_r + k_r \theta c) - \gamma_r c_r \quad (28b) $$

$$ \partial_t c_a + \nabla_x \cdot (c_a \chi \nabla_x \phi - D \nabla_x c_a) = \alpha k_c cc_r - \gamma_a c_a - T_1 c_a, \quad (28c) $$

$$ \partial_t I = \psi(\mu_1) - \tau I \quad (28d) $$

$$ -\nabla_x \cdot (K \nabla_x \phi) = f(\mu_1)\sigma, \quad (28e) $$

$$ c|_{\partial\Omega} = 0, \quad c_r|_{\partial\Omega} = 0, \quad \nabla_x \phi \cdot \nu(\cdot)|_{\partial\Omega} = 0, \quad (28f) $$

$$ c(t=0, x) = c_0(x), \quad c(t=0, x) = c_r^0(x), \quad I(t=0) = I_0. \quad (28g) $$

The kinetic of the drug effect is described by the following equation

$$ \partial_t T_1 = \kappa(t) - d_{T_1} T_1, \quad (29) $$

where $t \mapsto \kappa(t)$ describes the drug administration protocol and $d_{T_1}$ is the degradation rate of this drug. For the numerical tests, we set

$$ \kappa(t) = \begin{cases} 0, & \forall 0 \le t \le t_0 \\ \sum_{k \ge 0} \kappa_k (t - kT_2), & \forall t \ge t_0 \end{cases} \quad (30) $$
where,

$$ \kappa_k(t) = \begin{cases} q, & t_0 + kT_2 < t \le t_0 + kT_2 + T_1 \\ 0, & t_0 + kT_2 + T_1 \le t < t_0 + (k+1)T_2. \end{cases} \qquad (31) $$

The model depends on

* the time $t_0$ when the treatment starts,

* the duration $T_2$ between two drug administrations,

* the duration $T_1$ of the drug administration,

* the administered drug concentration $q$.

For the numerical tests, we place ourselves in the same configuration as in Fig 8(c) where the tumor escapes the immune control due to the effects of protumor immune cells. We fix $\alpha$, the proportion of effector T cells that become exhausted to 0.5 and we keep the other parameters as in Table 2. We set

$$ T_1 = 1, \quad T_2 = 7, \quad d_{T_1} = 0.05 $$

and we make the starting time $t_0$ and the dose $q$ vary. We indeed observe that these parameters have a critical role on the treatment efficacy.

When the treatment is given early (for instance, when $0 \le t_0 \le 5$), the control of the tumor can be obtained with relatively low drug doses (see Fig 18), in comparison to the cases where the treatment is administered later (Fig 19). At these early administration of the treatment, the tumor growth is controlled with a residue of dormant tumor cells and activated effector immune cells. Reducing the treatment dose reduces the drug efficacy with smaller tumor masses reached over longer period of time. For very small doses, the escape can occur.

<figure><img src="image_23.png" /><figcaption>**Fig 18. Reactivation of exhausted antitumor cells: Early administration of the treatment.** Simulation performed on the system (28a)–(28g), with (29)–(31). Evolution of the tumor mass $\mu_1$ (plain, left axis), and of the immune strength $\bar{\mu}_c$ (dashed, right axis) for several values of the treatment dose $q$. The dash-dotted line represents the time at which the treatment starts.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0259291.g018
<figure><img src="image_24.png" /><figcaption>**Fig 19. Reactivation of exhausted antitumor cells: Late administration of the treatment.** Simulation performed on the system (28a)–(28g), with (29)–(31). Evolution of the tumor mass $\mu_1$ (plain, left axis), and of $\bar{\mu}_c$ (dashed, right axis) for several values of the treatment dose $q$. The dash-dotted line represents the time at which the treatment starts.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0259291.g019

When the treatment is given later (for instance, when $10 \le t_0 \le 15$, see Fig 19): the tumor growth is slowed down by the treatment, but the tumor continues to grow exponentially fast. Increasing the drug dose increases the treatment efficacy. However, this observation raises the issue of the toxicity of the administered drug.

These observations are in line with experimental data using immune checkpoint inhibitors or CART/NK cells. Indeed, syngeneic CMS5 fibrosarcomas allowed to grow for 3 days in vivo were easily eradicated by adoptive transferred tumor-specific T cells while a 100-fold larger number of transferred tumor specific T-cell was mandatory to eradicate tumors that have been grown for an additional 48 hours. The same tumors that have been grown for 7 days before transferring adoptive tumor-specific T cells were not eradicated [42].

**Therapy based on reducing cytokine signals recruiting protumor immune cells.** Treatments based on blocking cytokine signals can help reducing the recruitment of protumor immune cells. A possible strategy uses cytokine traps [43, 44], by means of molecules that inhibit signal transduction from T cell cytokine receptors. Therefore, the treatment acts by down-regulating the effect of the tumor induced cytokines. We denote by $\mathcal{T}_2$, the effect of treatments which are able to block those cytokines. It obeys a kinetic similar to (29)

$$ \partial_t \mathcal{T}_2 = \kappa^{(2)}(t) - d_{\mathcal{T}_2} \mathcal{T}_2, \qquad (32) $$

where

$$ \kappa^{(2)}(t) = \begin{cases} 0, & \forall 0 \le t \le t_0 \\ \sum_{k \ge 0} \kappa_k^{(2)}(t - kT_2), & \forall t \ge t_0 \end{cases} \qquad (33) $$

and

$$ \kappa_k^{(2)}(t) = \begin{cases} q_2, & t_0 + kT_2 < t \le t_0 + kT_2 + T_1 \\ 0, & t_0 + kT_2 + T_1 \le t < t_0 + (k+1)T_2. \end{cases} \qquad (34) $$
The effect on the cytokines is described by modifying in (12b), (12c) the terms related to the cytokine-dependent recruitment of protumor immune cells. Therefore, the equations on the immune response become

$$ \partial_t c + \nabla_x \cdot (c\chi\nabla_x\phi - D\nabla_x c) = g(\mu_1)S - \gamma c - k_r I[1 - T_2]_+ \theta c - k_c cc_r, \quad (35a) $$

$$ \partial_t c_r + \nabla_x \cdot (c_r \chi \nabla_x \phi - D \nabla_x c_r) = I[1 - T_2]_+ (S_r + k_r \theta c) - \gamma_r c_r, \quad (35b) $$

$$ \partial_t c_a + \nabla_x \cdot (c_a \chi \nabla_x \phi - D \nabla_x c_a) = \alpha k_c cc_r - \gamma_a c_a, \quad (35c) $$

$$ \partial_t I = \psi(\mu_1) - \tau I \quad (35d) $$

$$ -\nabla_x \cdot (K \nabla_x \phi) = f(\mu_1)\sigma, \quad (35e) $$

$$ c|_{\partial\Omega} = 0, \quad c_r|_{\partial\Omega} = 0, \quad \nabla_x \phi \cdot v(\cdot)|_{\partial\Omega} = 0, \quad (35f) $$

$$ c(t=0, x) = c_0(x), \quad c(t=0, x) = c_r^0(x), \quad I(t=0) = I_0. \quad (35g) $$

For the numerical tests, we set

$$ T_1 = 1, \quad T_2 = 7, \quad d_{T_2} = 0.0105. $$

For the cytokine-blockade based treatment we observe a similar behavior as with the treatment based on the reactivation of the exhausted immune cells. The efficacy of the treatment is particularly sensitive to the starting time $t_0$, see Fig 20.

**Combination of the two immunotherapy strategies.** When we combine the two treatments described above, acting on both the reactivation of antitumor immune cells and the blockade of the recruitment of protumor immune cells, we observe that this combination is

<figure><img src="image_25.png" /><figcaption>**Fig 20.** Reduction of the protumor recruitment: Early administration of the treatment. Simulation performed on the system (35a)–(35g), with (32)–(34). Evolution of the tumor mass $\mu_1$ (plain, left axis), and of the immune strength $\mu_c$ (dashed, right axis) for several values of the treatment dose $q$. The dash-dotted line represents the time at which the treatment starts.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0259291.g020
<figure><img src="image_26.png" /><figcaption>**Fig 21.** Reduction of the protumor recruitment: Late administration of the treatment. Simulation performed on the system (35a)–(35g), with (32)–(34). Evolution of the tumor mass $\mu_1$ (red curves, left axis), and of the immune strength $\bar{\mu}_c$ (blue curve, right axis) for several values of the treatment dose $q$. The dash-dotted line represents the time at which the treatment starts.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0259291.g021

more efficient than the mono-therapies. Indeed a suitable combination of the treatment doses is able to control the tumor growth. For instance, the treatment based on the reactivation of exhausted immune cells fails in controlling the tumor when given at $t_0 = 10$ with a dose $q = 2$ see Fig 19(b), and the treatment based on cytokine/chemokine blockade fails with a dose $q_2 = 0.12$ at $t_0 = 10$, see Fig 21. However, the combination of the two treatments controls the tumor. Again, we observe that giving the treatments later requires to readjust the doses in order to control the tumor growth, see Fig 22. We notice that the controlled state contains residual tumor cells and activated immune cells, see Fig 23(b) and 23(c).

Conclusion

This work introduces a mathematical model describing the interactions between tumor cells and the immune system that regulate tumor growth, taking into account the antagonistic effects of antitumor and protumor immune cells. While the antitumor action aims at eliminating tumor cells, the protumor effects favor its growth. The later can take different forms: elimination of antitumor cells, conversion of antitumor cells into protumor cells, or enhancement

<figure><img src="image_27.png" /><figcaption>**Fig 22.** Administration of the combined treatments at $t_0 = 15$. Evolution of the tumor mass $\mu_1$ (plain, left axis), and of the immune strength $\bar{\mu}_c$ (dashed, right axis) for several values of the treatment dose $q$. The dash-dotted line represents the time at which the treatment starts.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0259291.g022
<figure><img src="image_28.png" /><figcaption>**Fig 23.** Administration of the combined treatments at $t_0 = 10$. Evolution of the tumor mass $\mu_1$ (plain, left axis), and of the immune strength $\bar{\mu}_c$ (dashed, right axis) for several values of the treatment dose $q$. The dash-dotted line represents the time at which the treatment starts.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0259291.g023

of tumor growth. Mechanisms that dictate the balance between these two conflicting functions with the TME are still not clear, mainly because of their complexity and heterogeneity. Mathematical modeling can help capture such complexity. Our model based on partial differential equations is remarkably able to reproduce equilibrium and escape phases, depending on the value of the biological parameters.

Compared to the situation free of protumor activities we previously modeled [24], where the equilibrium seems to always occur (possibly on very long scale of time, though), the addition of protumor immune cells leads to uncontrolled tumor growth when tumor aggressiveness overcomes the efficiency of the antitumor immune responses. Our model thus pinpoints the critical role of the protumor immune response in the establishment of the escape phase. This is consistent with experimental and clinical data. Indeed, as reported in [41] the blockade of protumor neutrophil recruitment by cMet inhibitor decreases tumor growth and potentiates anti-PD1 immunotherapy. Similarly, anti-CD25 depleting antibody optimized to deplete Treg within tumors enhances anti-tumor immune responses and synergizes with anti-PD1 treatment [45]. Based on these findings, we used our model to study the effects on tumor progression of two common cancer therapeutic strategies, either the reactivation of hyporesponsive antitumor cells or the reduction of the recruitment of protumor cells. Such therapies boost the immune response and restore the equilibrium that maintains
the tumor in a viable state. Importantly, the numerical investigation brings out the influence
of the starting time of the treatment and of the administrated dose. We also show on numeri-
cal grounds that combining the two approaches clearly improves the efficacy of the treat-
ment. Such information is highly valuable and together with clinical observations, can
comfort clinical decisions.

This preliminary study opens challenging perspectives. First, mathematical analysis can provide useful information, starting with further simplified equations, to understand the driving mechanisms of the equilibrium/escape phenomena and the effect of treatments. In particular this raises the practical issue of defining criteria that characterize the efficiency of the immune response, in line with RECIST recommendations [46]. Indeed, not only the residual mass of the tumor can be used as a relevant criterion, as in [36], but one has also to consider the time necessary to reach an equilibrium, as well as the features of the transient states. For treatments, the analysis, that should additionally consider toxicity effects, can help in understanding the optimal balance between dose and time of administration. Second, as simplified as it is, the model contains many parameters. Most of them are not known or even not easily accessible to experiments. Hence, based on optimisation techniques, an important work of parameter calibration should be performed from clinical data, with the two-fold difficulty that available data are rarely structured in time and space, and that data fitting techniques are far less developed for PDE than for ODE. Nevertheless, new mass cytometry imagery techniques open perspectives to address this issue [47]. Such investigation will permit us to determine relevant ranges for the parameters, which, in turn, will allow us to perform a detailed sensitivity analysis, beyond the attempt in [36]. This will be a decisive step to address in details the effects and the optimization of targeted treatments.

**Acknowledgments**

The authors acknowledge the support of UCAncer, an incentive Université Côte d'Azur network, which has permitted and encouraged this collaboration.

**Author Contributions**

**Conceptualization:** Kevin Atsou, Fabienne Anjuère, Véronique M. Braud, Thierry Goudon.

**Formal analysis:** Kevin Atsou, Véronique M. Braud, Thierry Goudon.

**Funding acquisition:** Fabienne Anjuère, Thierry Goudon.

**Investigation:** Kevin Atsou, Véronique M. Braud, Thierry Goudon.

**Methodology:** Kevin Atsou, Thierry Goudon.

**Software:** Kevin Atsou, Thierry Goudon.

**Supervision:** Véronique M. Braud, Thierry Goudon.

**Validation:** Kevin Atsou, Thierry Goudon.

**Visualization:** Kevin Atsou, Véronique M. Braud, Thierry Goudon.

**Writing – original draft:** Kevin Atsou, Fabienne Anjuère, Véronique M. Braud, Thierry Goudon.

**Writing – review & editing:** Kevin Atsou, Fabienne Anjuère, Véronique M. Braud, Thierry Goudon.
References

1. Dunn GP, Bruce AT, Ikeda H, Old LJ, Schreiber RD. Cancer immunoediting: from immunosurveillance to tumor escape. Nature Immunology. 2002; 3:991-998. https://doi.org/10.1038/ni1102-991 PMID: 12407406

2. Dunn GP, Old LJ, Schreiber RD. The Immunobiology Review of Cancer Immunosurveillance and Immunoediting. Immunity. 2004; 21:137-148.

3. Bruni D, Angel HK, Galon J. The immune contexture and immunoscore in cancer prognosis and therapeutic efficacy. Nature Reviews Cancer. 2020; 20:662-680. https://doi.org/10.1038/s41568-020-0285-7 PMID: 32753728

4. Galon J, Costes A, Sanchez-Cabo F, Kirilovsky A, Mlecnik B, Lagorce-Pagès C, et al. Type, density, and location of immune cells within human colorectal tumors predict clinical outcome. Science. 2006; 313:1960-1964. https://doi.org/10.1126/science.1129139 PMID: 17008531

5. Faget J, Groeneveld S, Boivin G, Sankar M, Zangger N, Garcia M, et al. Neutrophils and Snail Orchestrate the Establishment of a Pro-tumor Microenvironment in Lung Cancer. Cell Report. 2017; 21:3190-3204. https://doi.org/10.1016/j celrep.2017.11.052

6. Gentles AJ, Newman AM, Liu CL, Bratman SV, Feng W, Kim D, et al. The prognostic landscape of genes and infiltrating immune cells across human cancers. Nat Med. 2015; p. 938-945. https://doi.org/10.1038/nm.3909 PMID: 26193342

7. Chen DA, Mellman I. Elements of cancer immunity and the cancer-immune set point. Nature. 2017; 541:321-330. https://doi.org/10.1038/nature21349 PMID: 28102259

8. Sato E, Olson SH, Ahn J, Bundy B, Nishikawa H, Qian F, et al. Intraepithelial CD8+ tumor-infiltrating lymphocytes and a high CD8+/regulatory T cell ratio are associated with favorable prognosis in ovarian cancer. Proc Nat Acad Sc. 2005; 102(51):18538-18543. https://doi.org/10.1073/pnas.0509182102

9. Allavena P, Sica A, Garlanda C, Mantovani A. The Yin-Yang of tumor-associated macrophages in neoplastic progression and immune surveillance. Immunological Rev. 2008; 222(1):155-161. https://doi.org/10.1111/j.1600-065X.2008.00607.x PMID: 18364000

10. Fridlender ZG, Sun J, Kim S, Kapoor V, Cheng G, Ling L, et al. Polarization of Tumor-Associated Neutrophil Phenotype by TGF-beta: "N1" versus "N2" TAN. Cancer Cell. 2009; 16(3):183-194. https://doi.org/10.1016/j.ccr.2009.06.017 PMID: 19732719

11. Gabrilovich DI, Nagaraj S. Myeloid-derived suppressor cells as regulators of the immune system. Nat Rev Immunol. 2009; 9:162-174. https://doi.org/10.1038/nri2506 PMID: 19197294

12. Munn DH, Bronte V. Immune suppressive mechanisms in the tumor microenvironment. Current Opinion in Immunology. 2016; 39:1-6. https://doi.org/10.1016/j.coi.2015.10.009 PMID: 26609943

13. Quezada SA P KS, Simpson TR, Allison JP. Shifting the equilibrium in cancer immunoediting: from tumor tolerance to eradication. Immunol Rev. 2011; 241(1):104-118. https://doi.org/10.1111/j.1600-065X.2011.01007.x PMID: 21488893

14. Pardoll DM. The blockade of immune checkpoints in cancer immunotherapy. Nat Rev Cancer. 2012; 12(4):252-264. https://doi.org/10.1038/nrc3239 PMID: 22437870

15. Sharma P, Allison JP. Dissecting the mechanisms of immune checkpoint therapy. Nat Rev Immunol. 2020; 20(2):75-76. https://doi.org/10.1038/s41577-020-0275-8 PMID: 31925406

16. de Pillis LG, Radunskaya AE, Wiseman CL. A Validated Mathematical Model of Cell-Mediated Immune Response to Tumor Growth. Cancer Research. 2005; 65(17):7950-7958. https://doi.org/10.1158/0008-5472.CAN-05-0564 PMID: 16140967

17. d'Onofrio A. A general framework for modeling tumor-immune system competition and immunotherapy: mathematical analysis and biomedical inferences. Physica D. 2005; 208(3-4):220-235. https://doi.org/10.1016/j.physd.2005.06.032

18. Eftimie R, Bramson JL, Earn DJD. Interactions Between the Immune System and Cancer: A Brief Review of Non-spatial Mathematical Models. Bulletin of Mathematical Biology. 2011; 73(1):2-32. https://doi.org/10.1007/s11538-010-9526-3 PMID: 20225137

19. Eftimie R, Hamam H. Modelling and investigation of the CD4+ T cells. Macrophages paradox in melanoma immunotherapies. J Theor Biol. 2017; 420(7):82-104. https://doi.org/10.1016/j.jtbi.2017.02.022 PMID: 28219660

20. Kirschner D, Panetta JC. Modeling immunotherapy of the tumor-immune interaction. J Math Biol. 1998; 37:235-252. https://doi.org/10.1007/s002850050127 PMID: 9785481

21. Kuznetsov VA. Basic models of tumor-immune system interactions: identification, analysis and predictions. In: A Survey of Models for Tumor-Immune System Dynamics. Modeling and Simulation in Science, Engineering & Technology. Birkhäuser; 1997. p. 237-294.
22. Kuznetsov VA, Knott GD. Modeling tumor regrowth and immunotherapy. Math Comput Mod. 2001; 33:1275-1287. https://doi.org/10.1016/S0895-7177(00)00314-9

23. Eftimie R, Gillard JJ, Cantrell DA. Mathematical Models for Immunology: Current State of the Art and Future Research Directions. Bull Math Biol. 2017; 78:2091-2134. https://doi.org/10.1007/s11538-016-0214-9

24. Atsou K, Anjuère F, Braud VM, Goudon T. A size and space structured model describing interactions of tumor cells with immune cells reveals cancer persistent equilibrium states in tumorigenesis. J Theor Biol. 2020; 490:110163. https://doi.org/10.1016/j.jtbi.2020.110163 PMID: 31981572

25. Bekkal Brikci F, Clairambault J, Ribba B, Perthame B. An age-and-cyclin-structured cell population model for healthy and tumoral tissues. J Math Biol. 2008; 57(1):91-110. https://doi.org/10.1007/s00285-007-0147-x PMID: 18064465

26. Michel P. Existence of a solution to the cell division eigenproblem. Models Math Meth App Sci. 2006; 16(Suppl. issue 1):1125-1153. https://doi.org/10.1142/S0218202506001480

27. Michel P, Mischler S, Perthame B. General relative entropy inequality: an illustration on growth models. J Math Pures et Appl. 2005; 84(9):1235-1260. https://doi.org/10.1016/j.matpur.2005.04.001

28. Perthame B. Transport equations in biology. Frontiers in Math. Birkhauser; 2007.

29. Basse B, Baguley BC, Marshall ES, Joseph WR, van Brunt B, Wake G, et al. A mathematical model for analysis of the cell cycle in cell lines derived from human tumors. J Math Biol. 2003; 47(4):295-312. https://doi.org/10.1007/s00285-003-0203-0 PMID: 14523574

30. Alzahrani EO, Kuang Y. Nutrient limitations as an explanation of Gompertzian tumor growth. Disc Cont Dyn Syst B. 2016; 21(2):357-372. https://doi.org/10.3934/dcdsb.2016.21.357

31. Gyllenberg M, Webb GF. Quiescence as an explanation of Gompertzian tumor growth. Growth, Development, and Aging. 1989; 53(1-2):25-33. PMID: 2807642

32. Gyllenberg M, Webb GF. A nonlinear structured population model of tumor growth with quiescence. J Math Biol. 1990; 28:671-694. https://doi.org/10.1007/BF00160231

33. Iwata K, Kawasaki K, Shigesada N. A dynamical model for the growth and size distribution of multiple metastatic tumors. J Theor Biol. 2000; 203(2):177-186. https://doi.org/10.1006/jtbi.2000.1075 PMID: 10704301

34. Kozusko F, Bajzer Z. Combining Gompertzian growth and cell population dynamics. Math Biosci. 2003; 185(2):153-167. https://doi.org/10.1016/S0025-5564(03)00094-4 PMID: 12941534

35. Doumic-Jauffret M, Gabriel P. Eigenelements of a general aggregation-fragmentation model. Math Models Methods Appl Sci. 2010; 20(5):757-783. https://doi.org/10.1142/S021820251000443X

36. Atsou K, Anjuère F, Braud VM, Goudon T. Analysis of the equilibrium phase in immune-controlled tumor predicts best strategies for anti-tumor therapies. Univ. Côte d'Azur, Inria, CNRS, LJAD; 2020.

37. Baccelli F, McDonald DR, Reynier J. A mean field model for multiple TCP connections through a buffer implementing RED. Performance Evaluation. 2002; 49(1-4):77-97. https://doi.org/10.1016/S0166-5316(02)00136-0

38. Perthame B, Ryzhik L. Exponential decay for the fragmentation or cell-division equation. J Differential Equations. 2005; 210:155-177. https://doi.org/10.1016/j.jde.2004.10.018

39. Bailey SR, Maus MV. Gene editing for immune cell therapies. Nat Biotechnol. 2019; 37(12):1425-1434. https://doi.org/10.1038/s41587-019-0137-8 PMID: 31160723

40. Steele CW, Karim SA, Leach JDG, Bailey P, Upstill-Goddard R, Rishi Lea. CXCR2 inhibition profoundly suppresses metastases and augments immunotherapy in pancreatic ductal adenocarcinoma. Cancer Cell. 2016; 29(6):832-845. https://doi.org/10.1016/j.ccell.2016.04.014 PMID: 27265504

41. Glodde N, Bald T, van den Boorn-Konijnenberg D, Nakamura K, O'Donnell JS, Szczepanski S, et al. Reactive Neutrophil Responses Dependent on the Receptor Tyrosine Kinase c-MET Limit Cancer Immunotherapy. Immunity. 2017; 47(4):789-802.e9. https://doi.org/10.1016/j.immuni.2017.09.012 PMID: 29045907

42. Hanson HL, Donermeyer DL, Ikeda H, White JM, Shankaran V, Old LJ, et al. Eradication of Established Tumors by CD8+ T Cell Adoptive Immunotherapy. Immunity. 2000; 13(2):265-276. https://doi.org/10.1016/S1074-7613(00)00026-1 PMID: 10981969

43. Berraondo P, Sanmamed MF, Ochoa MC, Etxeberria I, Aznar MA, Pérez-Gracia JL, et al. Cytokines in clinical cancer immunotherapy. British Journal of Cancer. 2019; 120(1):6-15. https://doi.org/10.1038/s41416-018-0328-y PMID: 30413827

44. Economides AN, Carpenter LR, Rudge JS, Wong V, Koehler-Stec EM, Hartnett C, et al. Cytokine traps: Multi-component, high-affinity blockers of cytokine action. Nature Medicine. 2003; 9(1):47-52. https://doi.org/10.1038/nm811 PMID: 12483208
45. Vargas FA, Furness AJS, et al. Fc-optimized anti-CD25 depletes tumor-infiltrating regulatory T Cells and synergizes with PD-1 blockade to eradicate established tumors. Immunity. 2017; 46:577–586.
https://doi.org/10.1016/j.immuni.2017.03.013

46. Schwartz LH, Litière S, de Vries E, Ford R, Gwyther S, Mandrekar S, et al. RECIST 1.1—Update and Clarification: From the RECIST Committee. Eur J Cancer. 2016; 62:132–137. https://doi.org/10.1016/j.ejca.2016.03.081 PMID: 27189322

47. Elaldi R, Hemon P, Petti L, Cosson E, Desrues B, Sudaka A, et al. High dimensional imaging mass cytometry panel to visualize the tumor immune microenvironment contexture. Front Immunol. 2021; 12: Article #666233. https://doi.org/10.3389/fimmu.2021.666233 PMID: 33936105