A mathematical model to the melanoma dynamics
involving CAR T-cells

Guilherme Rodrigues, Jairo G Silva, Mostafa Adimy, Paulo F A Mancera

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

To cite this version:

Guilherme Rodrigues, Jairo G Silva, Mostafa Adimy, Paulo F A Mancera. A mathematical model to the melanoma dynamics involving CAR T-cells. Journal of Computational and Applied Mathematics, 2025, 44 (101), pp.30. 10.1007/s40314-024-03060-3. hal-04864443

HAL Id: hal-04864443

https://inria.hal.science/hal-04864443v1

Submitted on 5 Jan 2025

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
A mathematical model to the melanoma
dynamics involving CAR T-cells

Guilherme Rodrigues<fnref n="1" />*, Jairo G. Silva<fnref n="2" />†, Mostafa Adimy<fnref n="3" />†
and Paulo F. A. Mancera<fnref n="4" />†

Abstract

Melanoma is one of the most aggressive types of cancer. Although it
has a low percentage of incidence in the population, a high degree
of lethality is observed due to its rapid metastasis. As melanoma
is a highly immunogenic cancer, it has been used as an experimen-
tal model in several studies aimed at developing therapies, such as
immunotherapy with Chimeric Antigen Receptor (CAR) T-cells. We pro-
pose a mathematical model of three ordinary differential equations to
describe the dynamics of melanoma in the presence of Tumor-Associated
Macrophages (TAMs) and CAR T-cell therapy, to assess the role of
TAMs cells in the failure of this melanoma therapy. We examine the

<footnote n="1">¹Research Program in Biometry, Institute of Biosciences, São Paulo State University (UNESP), Botucatu, 18618-689, São Paulo, Brazil.</footnote>

<footnote n="2">²Federal Institute of Mato Grosso, Barra do Garças, 78607-899, Mato Grosso, Brazil.</footnote>

<footnote n="3">³Inria, ICJ UMR 5208, CNRS, Ecole Centrale de Lyon, INSA Lyon, Universite Claude Bernard Lyon 1, Universite Jean Monnet, Villeurbanne, F-69603, France.</footnote>

<footnote n="4">⁴Department of Biodiversity and Biostatistics, Institute of Biosciences, São Paulo State University (UNESP), Botucatu, 18618-689, São Paulo, Brazil.</footnote>

<footnote>*Corresponding author(s). E-mail(s): g.rodrigues2@unesp.br;
Contributing authors: jairo.gomes@ifmt.edu.br;
mostafa.adimy@inria.fr; paulo.mancera@unesp.br;
†These authors contributed equally to this work.</footnote>
existence and asymptotic stability of equilibrium points of this sys-
tem, giving a biological interpretation to each of them. Based on our
theoretical and numerical results, we conclude that immunosuppres-
sion has a negative impact on CAR T-cell immunotherapy and that
increasing the immunotherapy dose can improve tumor control. Fur-
thermore, an increase in the action of the TAMs population on tumor
proliferation can induce oscillations that eventually become periodic.

**Keywords:** Immunotherapy - Tumor-Associated Macrophages (TAMs) - Cancer modeling - Stability analysis - Sensitivity analysis

1 Introduction

Melanoma skin cancer is considered one of the most aggressive types of cancer
due to its capacity to invade other tissues and organs in the process of metasta-
sis. Projections of new cases and number of deaths from cutaneous melanoma
indicate that, based on global population changes, more than 500,000 new cases
and nearly 100,000 deaths due to melanoma will occur by 2040 [1]. In addi-
tion, 100,640 new cases and 8,290 deaths from melanoma were estimated for
the year 2024 [2]. Melanoma occurs when melanocytes, the cells that make up
the skin and produce the melanin responsible for skin pigmentation, undergo
genetic mutations that cause them to multiply in a disordered manner, pass-
ing on their characteristics to subsequent generations of cells, thus forming
a malignant tumor. Melanoma is known to be highly resistant to chemother-
apy, even when administered in a targeted manner, which compromises the
long-term success of treatments. When identified at an early stage, approxi-
mately 80% of cases are treated by surgical resection, which greatly increases
the chances of cure. However, if the disease is detected late, survival rates drop
too dramatically [3].

Adoptive Cell Therapy (ACT) has emerged as a promising approach to
cancer treatment. It is a therapeutic technique in which cells from a patient’s
immune system are harvested, modified, or stimulated in the laboratory, and
then transplanted into the patient to fight cancer. The ACT using Chimeric
Antigen Receptor (CAR) T-cells has been employed as a new strategy to
enhance cancer treatment, particularly against hematological malignancies,
exhibiting great efficacy [4, 5]. Due to its immunogenicity, melanoma is
used as an experimental model in several experiments aimed at developing
immunotherapeutic drugs against solid cancers, such as CAR T-cell therapy
[6]. There are many challenges to the applicability of this immunotherapy
against melanoma, such as the difficulty related to the adequate selection
of target antigens, and the resistance encountered by CAR T-cells to break
the barriers that prevent their infiltration in the immunosuppressive Tumor
Micro-Environment (TME) in the solid tumors [6].
One of the mechanisms that induce the rapid progression of melanoma is the recruitment of immune system cells [7]. Among all immune system cells recruited to tumor sites, macrophages are in the majority and act in tumor progression through various pathways. Macrophages that infiltrate TME of solid tumors are defined as TAMs or M2-type macrophages. This class of pro-tumor macrophages favors tumor progression by secreting growth factors, pro-angiogenic molecules, and immunosuppressive factors such as, for example, interleukin-10 (IL-10). TAMs are anti-inflammatory, increase cell motility and, by producing proteolytic enzymes that digest the extracellular matrix, pave the way for the dissemination of tumor cells from the primary site, contributing to metastasis [8–12]. Some studies indicate that TAMs density in the tumor microenvironment correlates with a poor prognosis of melanoma and acts to prevent the cytotoxic action of CAR T-cells against malignant cells due to the immunosuppressive factors secreted in the tumor microenvironment [9, 13–15].

Immunotherapy is one of the most promising approaches to cancer treatment. Some mathematical models have been developed to study the interactions between tumors and immune cells, to elucidate the various underlying mechanisms that are not yet fully understood. For example, Kuznetsov et al. [16], developed a mathematical model to study the response of cytotoxic T-lymphocytes to the growth of an immunogenic tumor, specifically B1-cell lymphoma (BCL1) in mice. Comparison of the results with experimental data allowed to estimate the kinetic parameters intrinsic to the dynamics. The results underlined the ability of their mathematical model to adequately describe cell population kinetics over a wide range of initial tumor cell concentrations.

In another paper, Eftimie and Hamam [17] investigated potential mechanisms behind the elimination of B16 melanoma cells in mice by the immune system, including Type 1 and Type 2 helper T-cells, as well as the role played by Type 1 and Type 2 macrophages in tumor growth and elimination. Their results indicate that the tumor can be eliminated in the presence of both Type 1 and Type 2 immune responses, while tumor growth is systematically associated with the presence of a Type 2 immune response.

Many other studies have addressed the dynamics of tumor and immune cells, including investigations conducted by de Pillis et al. [18–20], Eftimie [21–23], and Shu et al. [24]. Mathematical models also have been developed to investigate the efficacy of CAR T-cell immunotherapy against tumors. One such model, proposed by Sahoo et al. [25], employed a prey-predator dynamics-based approach to analyze the kinetics of CAR T-cell death during glioma treatment. Given the limited effectiveness of CAR-T immunotherapy in solid tumors, the authors conducted analyses of treatment success and failure, aligning the results with experimental data. The findings suggested that the CAR T-cell mortality rate is directly correlated with proliferation and exhaustion rates, implying that lower immunotherapeutic doses may be more effective, despite an increase in CAR T-cell exhaustion compared to higher doses.
Other models developed by Barros *et al.* [26, 27] explored the interaction between tumor cells, effector CAR T-cells, and memory CAR T-cells. Memory CAR T-cells are introduced into the model to study their ability to combat tumor cells after infusion of the immunotherapeutic dose. The results indicate that the proposed model has the potential to contribute to the reduction and optimization of *in vivo* experiments through computational simulations.

Additionally, León-Triana *et al.* [28] conducted *in silico* studies on the utilization of CAR T-cells targeting CD19 antigens expressed in B-cells and an antigen associated with glioblastoma (GBM). The results underscore that the expansion of the CAR T-cell population in the blood and lymphopoietic systems may lead to the massive production of CAR T-cells, directed towards solid tumors, potentially surpassing their immunosuppressive capacities. These mathematical contributions offer valuable perspectives for understanding and enhancing CAR T-cell immunotherapy against various tumor types [28, 29].

Other models to study the dynamics of CAR T-cells have also been developed, such as those by Santurio *et al.* [30], Mostolizadeh *et al.* [31], and Perez-Garcia *et al.* [32].

The main goal of our work is to develop a simple and innovative mathematical model, based on ordinary differential equations (ODEs), to elucidate the subtleties of the dynamics of advanced-stage cutaneous melanoma during immunotherapeutic intervention employing CAR T-cells. This model introduces a paradigm shift by explicitly incorporating the TAMs population. This new improvement enables us to perform comprehensive *in silico* analyses, highlighting the essential role played by the TAMs population in the efficacy and potential failures of the immunotherapeutic approach.

The paper is structured as follows: In Section 2, we present the proposed ODE-based mathematical model and describe the variables and parameters used. In Section 3, we introduce the equilibrium points, analyze their asymptotic stability, and show the existence of oscillations. In Section 4, we present the sensitivity analysis of the parameters. In Section 5, we complete our study with numerical simulations and discussions of several scenarios. Finally, in Section 6, we conclude this work with some remarks.

## 2 Mathematical model

The study of the effectiveness of CAR T-cell immunotherapy against cutaneous melanoma is based on three time-dependent state variables: *T*, the melanoma cell population; *M*, the TAMs population; and *C*, the CAR T-cell population (immunotherapeutic cell population). The proposed model is given by the
following system of ordinary differential equations

$$ \left\{ \begin{aligned} \frac{dT}{dt} &= \alpha_1 T \left(1 - \frac{T}{K_1}\right) (1 + \beta_1 M) - \mu_1 CT, \\ \frac{dM}{dt} &= \alpha_2 MT \left(1 - \frac{M}{K_2}\right) - \kappa_2 M, \\ \frac{dC}{dt} &= \frac{\alpha_3 CT}{K_3 + T} - \mu_3 CT - \kappa_3 C. \end{aligned} \right. \quad (1) $$

The equation associated with the variable *T* is based on the work of Sahoo *et al.* [25] and Eftimie and Hamam [17]. The equation related to the variable *M* is derived from the research by Eftimie and Hamam [17]. The equation associated with the variable *C* is based on the studies of Kuznetsov *et al.* [16] and León-Triana *et al.* [28]. In addition, we make the following assumptions:

(i) Tumor cells naturally grow according to the logistic law, at a rate $\alpha_1$, with a carrying capacity, $K_1$ [25]. The presence of TAMs in the tumor microenvironment significantly contributes to the proliferation of tumor cells due to the production of pro-tumoral cytokines, leading to a more rapid expansion of the tumor population [33]. To model this effect, we multiply the logistic term in the equation of *T* by the factor $1 + \beta_1 M$. CAR T-cells destroy tumor cells at a rate of $\mu_1$, through selective recognition of specific tumor-associated antigens. In the case of melanoma, these antigens include, for example, GD2 (dialoganglioside), VEGF-2 (Vascular Endothelial Growth Factor), MCSP (Melanoma-associated Chondroitin Sulfate Proteoglycan), among others [25, 34].

(ii) Based on the work of Eftimie and Hamam, in the second equation of the model, we assume that TAMs are recruited to the tumor site at a rate of $\alpha_2 T$ and infiltrate the tumor microenvironment in response to cytokines and chemokines produced by tumor cells. In contrast to Eftimie and Hamam, we do not consider the presence of other immune system cells such as Type 1 and 2 helper T-cells (Th1 and Th2), nor Type 1 macrophages (M1). Furthermore, we assume that TAMs population dynamics is characterized by the logistic law with a carrying capacity $K_2$ and that cells undergo apoptosis at a rate $\kappa_2$ [17, 33].

(iii) Considering the equation of *C*, CAR T-cell population expansion depends on their contact with tumor antigen, at a rate represented by $\alpha_3$. We formulated this growth process based on the works of Kuznetsov *et al.* [16] and León-Triana *et al.* [28], employing a Michaelis-Menten term in which $K_3$ determines the point at which the growth rate of CAR T-cells reaches half of its maximum value ($T = K_3$). León-Triana *et al.* postulate that $K_3$ denotes the inflection point from which the stimulation rate of CAR T-cells increases, correlating with the level of antigen expressed by tumor cells. The parameter $\kappa_3$ is employed to model the natural mortality of CAR T-cells.

(iv) Immunosuppression is a reduction in the effectiveness of the immune system. It represents a critical obstacle to CAR T-cell therapy [15]. We incorporate
the immunosuppressive effects on CAR T-cells into our model through the term $-\mu_3 CT$. This term represents the immunosuppression resulting from the recruitment of pro-tumor cells to the TME. The degree of immunosuppression depends on the value of the parameter $\mu_3$. Higher values of $\mu_3$ indicate stronger immunosuppressive activity [35, 36].

A compartmental diagram of the proposed model is exhibited in Figure 1.

<figure><img src="image_3.png" /><figcaption>**Fig. 1** Compartmental model of melanoma cell dynamics, TAMs and CAR T-cells. The figure was created by the authors and partially generated using Servier Medical Art, provided by Servier, under Creative Commons Attribution 3.0 unported license.</figcaption></figure>

Table 1 shows the model parameters and their respective meanings and units.

<table><caption>**Table 1** Descriptions, values, and units of the parameters.</caption><thead><tr><th>Parameter</th><th>Description</th><th>Value</th><th>Unit</th><th>Reference</th></tr></thead><tbody><tr><td>α<sub>1</sub></td><td>Rate of natural proliferation of tumor cells</td><td>[0.69,0.97]</td><td>day<sup>-1</sup></td><td>[17]</td></tr><tr><td>α<sub>2</sub></td><td>Rate of macrophage proliferation promoted by tumor cells</td><td>10<sup>-7</sup></td><td>(cells-day)<sup>-1</sup></td><td>[17]</td></tr><tr><td>α<sub>3</sub></td><td>Rate of proliferation of CAR T-cells due to contact with tumor antigen</td><td>0.9</td><td>day<sup>-1</sup></td><td>[28]</td></tr><tr><td>β<sub>1</sub></td><td>Rate of tumor proliferation due to pro-tumor action of TAMs</td><td>[0, 3, 10<sup>-6</sup>]</td><td>cells<sup>-1</sup></td><td>Assumed</td></tr><tr><td>K<sub>1</sub></td><td>Carrying capacity of tumor cells</td><td>10<sup>9</sup></td><td>cells</td><td>[21]</td></tr><tr><td>K<sub>2</sub></td><td>Carrying capacity TAMs population</td><td>10<sup>9</sup></td><td>cells</td><td>[21, 37]</td></tr><tr><td>K<sub>3</sub></td><td>Saturation coefficient of CAR T-cell population</td><td>2, 10<sup>7</sup></td><td>cells</td><td>[16]</td></tr><tr><td>μ<sub>1</sub></td><td>Rate of death of tumor cells due to cytotoxic action of CAR T-cells</td><td>10<sup>-7</sup></td><td>(cells-day)<sup>-1</sup></td><td>[27]</td></tr><tr><td>μ<sub>3</sub></td><td>Rate of death of CAR T-cells due to tumor immunosuppressive effects in the TME</td><td>[0, 5, 10<sup>-9</sup>]</td><td>(cells-day)<sup>-1</sup></td><td>Assumed</td></tr><tr><td>κ<sub>2</sub></td><td>Rate of apoptosis of TAMs population</td><td>0.34</td><td>day<sup>-1</sup></td><td>[22, 38]</td></tr><tr><td>κ<sub>3</sub></td><td>Rate of apoptosis of CAR T-cells</td><td>[0,0,0.6]</td><td>day<sup>-1</sup></td><td>Assumed</td></tr><tr><td>T(0)</td><td>Initial condition of tumor cells</td><td>2, 10<sup>6</sup></td><td>cells</td><td>Assumed</td></tr><tr><td>M(0)</td><td>Initial condition of TAMs population</td><td>10<sup>6</sup></td><td>cells</td><td>[9]</td></tr><tr><td>C(0)</td><td>Initial condition of CAR T-cells</td><td>[5, 10<sup>3</sup>, 9, 10<sup>7</sup>]</td><td>cells</td><td>Assumed</td></tr></tbody></table>

## 3 Analysis of the model

Firstly, it is necessary to examine the third equation of the model to better understand the dynamics of CAR T-cells in the presence of immunosuppressive
factors. Let us consider the following function

$$f_3(T) := \frac{\alpha_3 T}{K_3 + T}. \qquad (2)$$

Then, the last equation of the model (1) becomes

$$\frac{dC}{dt} = (f_3(T) - (\mu_3 T + \kappa_3)) C.$$

To analyze the interaction between CAR T-cell growth and natural death, we also introduce the following function

$$g_3(T) = \frac{1}{T}(f_3(T) - \kappa_3), \quad T > 0. \qquad (3)$$

The function $g_3$ is specifically designed to measure the net growth/decay of CAR T-cells per tumor cell, independent of the tumor-dependent removal term $\mu_3 T$. By excluding the term $\mu_3 T$, we isolate the effects of natural CAR T-cell dynamics, capturing the balance between the tumor-driven proliferation, represented by $f_3(T)$, and the natural death rate $\kappa_3$. This allows us to assess how the CAR T-cell population evolves as a function of tumor load when the influence of direct interaction with tumor cells is controlled. We have the limits

$$\lim_{T \to 0} g_3(T) = -\infty \quad \text{and} \quad \lim_{T \to +\infty} g_3(T) = 0.$$

The derivative is given by

$$g_3'(T) = \frac{1}{T^2}(f_3'(T)T - (f_3(T) - \kappa_3)).$$

In other words

$$g_3'(T) = \frac{1}{T^2(K_3 + T)^2} ((\kappa_3 - \alpha_3)T^2 + 2\kappa_3 K_3 T + \kappa_3 K_3^2).$$

If $\kappa_3 \ge \alpha_3$, then $g_3'(T) > 0$ for all $T > 0$. Suppose now that $\kappa_3 < \alpha_3$. The discriminant of the second-degree polynomial is

$$\Delta = 4[\kappa_3^2 K_3^2 - (\kappa_3 - \alpha_3)\kappa_3 K_3^2] = 4\kappa_3\alpha_3 K_3^2 > 0.$$

Then, the unique positive solution is given by

$$\bar{T}_0 = \frac{\kappa_3 K_3 + K_3 \sqrt{\kappa_3 \alpha_3}}{\alpha_3 - \kappa_3} = \frac{\sqrt{\kappa_3}}{\sqrt{\alpha_3} - \sqrt{\kappa_3}} K_3.$$
We conclude that $g_3$ is increasing on $(0, \bar{T}_0)$ and decreasing on $(\bar{T}_0, +\infty)$. We define a threshold for the parameter $\mu_3$ as follows

$$ \mu_{3\max} := \max\{g_3(T) : T > 0\} = g(\bar{T}_0), \quad \text{with} \quad \bar{T}_0 = \frac{\sqrt{\kappa_3}}{\sqrt{\alpha_3} - \sqrt{\kappa_3}} K_3. $$

We proved the following proposition.

**Proposition 1** Let $\kappa_3$, $\alpha_3$ and $K_3$ be positive.

* Suppose that $\kappa_3 \ge \alpha_3$. Then, for all $T > 0$, $g_3(T) < 0$.

* Suppose that $\kappa_3 < \alpha_3$.

- If $\mu_3 > \mu_{3\max}$, then $f_3(T) = \mu_3 T + \kappa_3$ has no solution.

- If $\mu_3 < \mu_{3\max}$, then $f_3(T) = \mu_3 T + \kappa_3$ has two solutions $\bar{T}_1 < \bar{T}_2$.

Figure 2 illustrates the Proposition 1.

<figure><img src="image_4.png" /><figcaption>**Fig. 2** Comparison between $f_3(T)$ and $\mu_3 T + \kappa_3$ for different parameter values. Using $\alpha_3 = 0.9$, $\kappa_3 = 0.35$ and $K_3 = 2 \cdot 10^7$, we obtained $\mu_{3\max} \approx 6.37 \cdot 10^{-9}$. Three scenarios are represented: the red dashed line where the growth and decay terms balance out at a single point; the green dashed line where the growth term dominates, leading to two solutions $\bar{T}_1 < \bar{T}_2$; and the orange dashed line where the decay term exceeds the growth term, leading to no solution.</figcaption></figure>

Figure 2 (see also Proposition 1) illustrates the interaction between the growth term $f_3(T)$, which represents the proliferation of CAR T-cells in response to tumor cells, and the linear decay term $\mu_3 T + \kappa_3$, which models the removal of CAR T-cells due to both their interaction with tumor cells (proportional to $\mu_3 T$) and natural cell death ($\kappa_3$).
## 3.1 Steady-states and their local asymptotic stability

The equilibrium points $\bar{E} := (\bar{T}, \bar{M}, \bar{C})$ of the system (1) are the non-negative solutions of

$$ egin{cases} 
\left( \alpha_1 \left( 1 - \frac{\bar{T}}{K_1} \right) (1 + \beta_1 \bar{M}) - \mu_1 \bar{C} \right) \bar{T} = 0, \\ 
\left( \alpha_2 \bar{T} \left( 1 - \frac{\bar{M}}{K_2} \right) - \kappa_2 \right) \bar{M} = 0, \\ 
(f_3(\bar{T}) - (\mu_3 \bar{T} + \kappa_3)) \bar{C} = 0. 
\end{cases} 
\quad (4) $$

We denote by $F_1$, $F_2$, and $F_3$ the three components of the nonlinear function on the right-hand side of System (1). Then, the partial derivatives at an equilibrium point $\bar{E} := (\bar{T}, \bar{M}, \bar{C})$ are given by

$$ \frac{\partial F_1}{\partial \bar{T}} = \alpha_1 \left( 1 - \frac{\bar{T}}{K_1} \right) (1 + \beta_1 \bar{M}) - \frac{\alpha_1 \bar{T} (1 + \beta_1 \bar{M})}{K_1} - \mu_1 \bar{C}, $$

$$ \frac{\partial F_1}{\partial \bar{M}} = \alpha_1 \bar{T} \left( 1 - \frac{\bar{T}}{K_1} \right) \beta_1, $$

$$ \frac{\partial F_1}{\partial \bar{C}} = -\mu_1 \bar{T}, $$

$$ \frac{\partial F_2}{\partial \bar{T}} = \alpha_2 \bar{M} \left( 1 - \frac{\bar{M}}{K_2} \right), $$

$$ \frac{\partial F_2}{\partial \bar{M}} = \alpha_2 \bar{T} \left( 1 - \frac{\bar{M}}{K_2} \right) - \frac{\alpha_2 \bar{M} \bar{T}}{K_2} - \kappa_2, $$

$$ \frac{\partial F_2}{\partial \bar{C}} = 0, $$

$$ \frac{\partial F_3}{\partial \bar{T}} = (f'_3(\bar{T}) - \mu_3) \bar{C}, $$

$$ \frac{\partial F_3}{\partial \bar{M}} = 0, $$

$$ \frac{\partial F_3}{\partial \bar{C}} = f_3(\bar{T}) - (\mu_3 \bar{T} + \kappa_3). $$

We begin by analyzing the local asymptotic stability of trivial steady-states. We introduce notations that will be the main thresholds of the two parameters $\alpha_2$ and $\alpha_3$, for the existence and stability of equilibrium points:

$$ \alpha_{2 \min} := \frac{\kappa_2}{K_1}, \quad \alpha_{3 \min} := \kappa_3 \left( 1 + \frac{K_3}{K_1} \right) \quad \text{and} \quad \alpha_{3 \max}(\alpha_2) := \kappa_3 \left( 1 + \frac{\alpha_2}{\kappa_2} K_3 \right). \quad (5) $$

We can notice that we have $\alpha_{3 \max}(\alpha_{2 \min}) = \alpha_{3 \min}$.
**Proposition 2** The trivial equilibrium $E_0 = (0,0,0)$ always exists and is unstable (saddle point). The equilibrium $E_1 = (K_1, 0, 0)$ always exists and is locally asymptotically stable if and only if $\alpha_2 < \alpha_{2\min}$ and $\mu_3 > g_3(K_1)$, where $\alpha_{2\min}$ is given by (5) and the function $g_3$ is given by (3).

*Proof* It is not difficult to check directly that $E_0$ and $E_1$ are trivial solutions of the system (4). The Jacobian matrix associated with the equilibrium point $E_0$ is

$$J_0 = \begin{pmatrix} \alpha_1 & 0 & 0 \\ 0 & -\kappa_2 & 0 \\ 0 & 0 & -\kappa_3 \end{pmatrix}.$$

The eigenvalues are $\alpha_1, -\kappa_2$ and $-\kappa_3$. Then, $E_0$ is unstable (saddle point).

The Jacobian matrix associated with the equilibrium point $E_1$ is

$$J_1 = \begin{pmatrix} -\alpha_1 & 0 & -\mu_1 K_1 \\ 0 & K_1\alpha_2 - \kappa_2 & 0 \\ 0 & 0 & f_3(K_1) - (\mu_3 K_1 + \kappa_3) \end{pmatrix}.$$

Then, the eigenvalues are $\lambda_1 = -\alpha_1 < 0$, $\lambda_2 = K_1\alpha_2 - \kappa_2$ and $\lambda_3 = f_3(K_1) - (\mu_3 K_1 + \kappa_3)$. We have

* $\lambda_2 < 0$ if $\alpha_2 < \alpha_{2\min}$.

* $\lambda_3 < 0$ if $f_3(K_1) < \mu_3 K_1 + \kappa_3$. This last inequality is equivalent to

$$\mu_3 > \frac{1}{K_1}(f_3(K_1) - \kappa_3) = g_3(K_1).$$

□

To satisfy the case $\mu_3 > g_3(K_1)$, we have two scenarios:

* Suppose that $\alpha_3 \le \kappa_3$. Then, we always have $\mu_3 > g_3(K_1)$.

* Suppose that $\alpha_3 > \kappa_3$. If $\mu_3 \ge \mu_{3\max}$, then we always have $\mu_3 > g_3(K_1)$. If $\mu_3 < \mu_{3\max}$, then $\mu_3 > g_3(K_1)$ if and only if $K_1 \in (0, \bar{T}_1) \cup (\bar{T}_2, +\infty)$, where $\bar{T}_1 < \bar{T}_2$ are the zeros of the function $g_3$ (see Figure 2).

Another equilibrium with the component $\bar{C} = 0$ may exist. More precisely, we have the following result.

**Proposition 3** $E_2 = (K_1, K_2 \left(1 - \frac{\kappa_2}{\alpha_2 K_1}\right), 0)$ is an equilibrium point if and only if $\alpha_2 > \alpha_{2\min}$. It is locally asymptotically stable if and only if $\alpha_2 > \alpha_{2\min}$ and $\mu_3 > g_3(K_1)$.
*Proof* It is not difficult to verify using the system (4) that $E_2$ is an equilibrium if and only if $\alpha_2 > \alpha_{2\min}$. Furthermore, the Jacobian matrix associated with this point is

$$J_2 = \begin{pmatrix} -\alpha_1 \left( K_2 \left( 1 - \frac{\kappa_2}{\alpha_2 K_1} \right) \beta_1 + 1 \right) & 0 & -\mu_1 K_1 \\ \frac{K_2(K_1\alpha_2 - \kappa_2)\kappa_2}{K_1^2\alpha_2} & -K_1\alpha_2 + \kappa_2 & 0 \\ 0 & 0 & f_3(K_1) - (\mu_3 K_1 + \kappa_3) \end{pmatrix}$$

Then, the eigenvalues are

$$\begin{aligned} \lambda_1 &= -\alpha_1 \left( K_2 \left( 1 - \frac{\kappa_2}{\alpha_2 K_1} \right) \beta_1 + 1 \right), \\ \lambda_2 &= -K_1 \alpha_2 + \kappa_2, \\ \lambda_3 &= f_3(K_1) - (\mu_3 K_1 + \kappa_3). \end{aligned}$$

We conclude that $E_2$ is locally asymptotically stable if and only if

$$\alpha_2 > \alpha_{2\min} \text{ and } \mu_3 > g_3(K_1).$$

□

We also have two equilibrium points with the component $\bar{M} = 0$. More precisely, we have the following result.

**Proposition 4** $E_3 = (\bar{T}, 0, \bar{C})$ is an equilibrium point if and only if $\mu_3 \le \mu_{3\max}$ and $0 < \bar{T} < K_1$. This point is locally asymptotically stable if and only if $\mu_3 \le \mu_{3\max}$, $0 < \bar{T} < \min\{K_1, \frac{\kappa_2}{\alpha_2}\}$ and $f'_3(\bar{T}) > \kappa_3$.

*Proof* For $\bar{M} = 0$ the system (4) becomes

$$\begin{cases} g_3(\bar{T}) = \mu_3, \\ \bar{C} = \frac{\alpha_1}{\mu_1} \left( 1 - \frac{\bar{T}}{K_1} \right). \end{cases}$$

The necessary and sufficient condition for the existence of at least one positive solution is

$$\mu_3 \le \mu_{3\max} \text{ and } 0 < \bar{T} < K_1.$$

There are two such equilibrium points, denoted by $E_{3i} := (\bar{T}_i, 0, \bar{C}_i)$, $i = 1, 2$, with

$$g_3(\bar{T}_i) = \mu_3 \text{ and } \bar{C}_i = \frac{\alpha_1}{\mu_1} \left( 1 - \frac{\bar{T}_i}{K_1} \right).$$

The Jacobian matrix is given by

$$J_3 = \begin{pmatrix} a_{11} & a_{12} & a_{13} \\ 0 & a_{22} & 0 \\ a_{31} & 0 & 0 \end{pmatrix},$$

with

$$a_{11} = -\alpha_1 \frac{\bar{T}}{K_1},$$
$$a_{12} = \alpha_1 \bar{T} \left( 1 - \frac{\bar{T}}{K_1} \right) \beta_1,$$

$$a_{13} = -\mu_1 \bar{T},$$

$$a_{22} = \alpha_2 \bar{T} - \kappa_2,$$

$$a_{31} = \frac{\alpha_1}{\mu_1} \left( 1 - \frac{\bar{T}}{K_1} \right) (f_3'(\bar{T}) - \kappa_3).$$

The characteristic equation is then given by

$$(\lambda - a_{22})(\lambda^2 - a_{11}\lambda - a_{13}a_{31}) = 0. \qquad (6)$$

According to the first factor of (6), we have the eigenvalue $\lambda_1 = a_{22} = \alpha_2 \bar{T} - \kappa_2$. It is negative if and only if

$$0 < \bar{T} < \min \left\{ K_1, \frac{\kappa_2}{\alpha_2} \right\}.$$

For the second term of (6), we use the Routh-Hurwitz criterion for the equation $\lambda^2 + p_1\lambda + p_2 = 0$, with

* $p_1 := -a_{11} = \alpha_1 \frac{\bar{T}}{K_1} > 0,$

* $p_2 := -a_{13}a_{31} > 0$ if and only if $\bar{T} < K_1$ and $f_3'(\bar{T}) > \kappa_3$.

□

One can remark that $f_3'(\bar{T}_1) > \kappa_3$ and $f_3'(\bar{T}_2) < \kappa_3$. Then, if $\mu_3 < \mu_{3\max}$ and $0 < \bar{T}_1 < \bar{T}_2 < \min\{K_1, \frac{\kappa_2}{\alpha_2}\}$, the equilibrium $E_{31} := (\bar{T}_1, 0, \bar{C}_1)$ is locally asymptotically stable and the equilibrium $E_{32} := (\bar{T}_2, 0, \bar{C}_2)$ is unstable (see Figure 2).

The analysis of the existence and local asymptotic stability of the equilibrium $E^*$ is much more technical than for other equilibrium points. Before addressing it, we give below the biological meanings of the four equilibrium points.

* $E_0 = (0, 0, 0)$ corresponds to the extinction of all populations. It is not biologically relevant to our study.

* $E_1 = (K_1, 0, 0)$ characterizes treatment failure as it describes the persistence of melanoma cells (at a high level, $K_1$), with the disappearance of CAR T-cells and TAMs population.

* $E_2 = (K_1, K_2 \left(1 - \frac{\kappa_2}{\alpha_2 K_1}\right), 0)$ also characterizes treatment failure with a high level of melanoma cell, the disappearance of CAR T-cells and the persistence of TAMs population. This point intersects $E_1$ at the value $\alpha_2 = \alpha_{2\min}$.

* $E_{3i} = (\bar{T}_i, 0, \bar{C}_i)$, $i = 1, 2$, represent the disappearance of TAMs population and the persistence of melanoma cells and CAR T-cells. These equilibrium points represent a more effective treatment outcome compared to the points $E_1$ and $E_2$, as they are characterized by $\bar{T}_i < K_1$.

In previous sections, we incorporated the term $\mu_3 CT$ to model potential immunosuppressive effects on CAR T-cell dynamics. However, due to the algebraic complexity introduced by this term in the analysis of the coexistence equilibrium $E^*$, we introduce the assumption $\mu_3 = 0$ to simplify the existence conditions and the analysis of the asymptotic stability of this equilibrium. This assumption allows us to focus
on biologically meaningful interactions without compromising the insights derived from the stability analysis of $E^*$. Furthermore, the assumption $\mu_3 = 0$ remains consistent with previous sections where immunosuppressive effects were taken into account in the analysis of the simplest equilibrium points. Therefore, removing this term for the coexistence equilibrium analysis does not affect the general validity of the results, as the immunosuppressive effects can be considered negligible in the specific context of the coexistence equilibrium. Nevertheless, the impact of $\mu_3 > 0$ will be studied in detail in Section 5, where we will numerically explore its influence on the dynamics of the system and the stability of the equilibrium points.

In the remaining part of this section, we assume that $\mu_3 = 0$.

**Proposition 5** *A coexistence equilibrium,*

$$E^* = \left( \frac{\kappa_3 K_3}{\alpha_3 - \kappa_3}, K_2 \left( 1 - \frac{\kappa_2(\alpha_3 - \kappa_3)}{\kappa_3 \alpha_2 K_3} \right), \frac{\alpha_1}{\mu_1} \left( 1 - \frac{T^*}{K_1} \right) (1 + \beta_1 M^*) \right),$$

where the third coordinate is expressed in terms of the first two coordinates

$$T^* = \frac{\kappa_3 K_3}{\alpha_3 - \kappa_3} \quad \text{and} \quad M^* = K_2 \left( 1 - \frac{\kappa_2(\alpha_3 - \kappa_3)}{\kappa_3 \alpha_2 K_3} \right),$$

exists if and only if $\alpha_{3\min} < \alpha_3 < \alpha_{3\max}(\alpha_2)$, where $\alpha_{3\min}$ and $\alpha_{3\max}(\alpha_2)$ are given by (5).

*Proof* Solving the system (4), with $\mu_3 = 0$, to find a coexistence equilibrium is straightforward and is not presented here. It is also not difficult to see that $E^*$ exists if and only if $\alpha_{3\min} < \alpha_3 < \alpha_{3\max}(\alpha_2)$. $\square$

Our objective in this paragraph is to determine some properties of the coexistence equilibrium $E^*$. We show that the last two components of $E^*$ can be written as functions of the first component $T^*$.

**Proposition 6** *The coexistence equilibrium $E^*$ can be written as a function*

$$T^* \in \Omega \mapsto E^*(T^*) := (T^*, M^*(T^*), C^*(T^*)),$$

with the interval $\Omega := (\frac{\kappa_2}{\alpha_2}, K_1)$ and

$$M^*(T^*) = K_2 \left( 1 - \frac{\kappa_2}{\alpha_2 T^*} \right) \text{ and } C^*(T^*) = \frac{\alpha_1}{\mu_1} \left( 1 - \frac{T^*}{K_1} \right) (1 + \beta_1 M^*(T^*)) .$$

Furthermore, $T^* \mapsto M^*(T^*)$ is an increasing function from the interval $\Omega$ to $(0, K_2(1 - \frac{\kappa_2}{\alpha_2 K_1}))$, and there exist two points

$$\bar{T}^* := \sqrt{\frac{\beta_1 K_1 K_2 \kappa_2}{\alpha_2(1 + \beta_1 K_2)}} \quad \text{and} \quad \bar{K}_1 := \frac{\kappa_2(1 + \beta_1 K_2)}{\alpha_2 \beta_1 K_2},$$

such that

(a) if $K_1 > \bar{K}_1$, then $T^* \mapsto C^*(T^*)$ is increasing on $(\kappa_2/\alpha_2, \bar{T}^*)$ and decreasing on $(\bar{T}^*, K_1)$;
(b) if $K_1 \le \bar{K}_1$, then $T^* \mapsto C^*(T^*)$ is decreasing on $\Omega$.

*Proof* The first component of $E^*$ is given by

$$T^* = \frac{\kappa_3 K_3}{\alpha_3 - \kappa_3}.$$

As a function of $\alpha_3 \in (\alpha_{3\min}, \alpha_{3\max}(\alpha_2))$, we obtain a decreasing function $\alpha_3 \mapsto T^*(\alpha_3)$ such that,

$$T^*(\alpha_{3\max}(\alpha_2)) = \frac{\kappa_2}{\alpha_2} < T^*(\alpha_{3\min}) = K_1.$$

This gives the interval $\Omega := (\kappa_2/\alpha_2, K_1)$ to which $T^*$ belongs. On the other hand, a simple calculation shows that

$$M^*(T^*) = K_2 \left(1 - \frac{\kappa_2}{\alpha_2 T^*}\right) \text{ and } C^*(T^*) = \frac{\alpha_1}{\mu_1} \left(1 - \frac{T^*}{K_1}\right) (1 + \beta_1 M^*(T^*)).$$

It is clear that the function $T^* \mapsto M^*(T^*)$ is increasing from the interval $\Omega$ to $(0, K_2(1 - \frac{\kappa_2}{\alpha_2 K_1}))$.

To study the function $T^* \mapsto C^*(T^*)$, we need to calculate its derivative on the interval $\Omega$. This derivative is given by

$$C^{*'}(T^*) = \frac{\alpha_1}{\mu_1} \left[ -\frac{1}{K_1} \left( 1 + \beta_1 K_2 \left( 1 - \frac{\kappa_2}{\alpha_2 T^*} \right) \right) + \left( 1 - \frac{T^*}{K_1} \right) \beta_1 K_2 \left( \frac{\kappa_2}{\alpha_2 T^{*2}} \right) \right].$$

Then, $C^{*'}(T^*) < 0$ if and only if $T^* > \bar{T}^* := \sqrt{\frac{\beta_1 K_1 K_2 \kappa_2}{\alpha_2(1 + \beta_1 K_2)}}$. We have to see if $\bar{T}^*$ belongs to the interval $\Omega := (\frac{\kappa_2}{\alpha_2}, K_1)$. Due to the inequalities $\frac{\beta_1 K_2 \kappa_2}{\alpha_2(1 + \beta_1 K_2)} < \frac{\kappa_2}{\alpha_2} < K_1$, we have $\bar{T}^* < K_1$. But we do not know the position of $\bar{T}^*$ in comparison with $\frac{\kappa_2}{\alpha_2}$. Two cases are possible:

(a) If $K_1 > \bar{K}_1 := \frac{\kappa_2(1 + \beta_1 K_2)}{\alpha_2 \beta_1 K_2}$, then $\bar{T}^* \in \Omega$. This implies that $T^* \mapsto C^*(T^*)$ is increasing on $(\kappa_2/\alpha_2, \bar{T}^*)$ and decreasing on $(\bar{T}^*, K_1)$.

(b) If $K_1 \le \bar{K}_1$, then $\bar{T}^* \le \kappa_2/\alpha_2$. This implies that $T^* \mapsto C^*(T^*)$ is always decreasing on $\Omega$.

The cases analyzed above are shown in Figures 3 and 4.
<figure><img src="image_5.png" /><figcaption>**Fig. 3** These figures represent the first case in which $K_1 > \bar{K}_1$. We use the parameters shown in Table 1 and choose $K_1 = 2 \cdot 10^7$ cells as the value greater than $\bar{K}_1$. Figure 3(a) shows that $T^* \mapsto M^*(T^*)$ is an increasing function on the interval $\Omega := (\kappa_2/\alpha_2, K_1)$. In Figure 3(b), we can see that $T^* \mapsto C^*(T^*)$ is increasing in $(\kappa_2/\alpha_2, \bar{T}^*)$ (see the zoom), the green star represents the point $\bar{T}^*$, and $T^* \mapsto C^*(T^*)$ is decreasing on $(\bar{T}^*, K_1)$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 4** These figures represent the second case, $K_1 \le \bar{K}_1$. We use the parameters given in the Table 1. Figure 4(a) shows that $T^* \mapsto M^*(T^*)$ is an increasing function on $\Omega$. In Figure 4(b) we can observe that $T^* \mapsto C^*(T^*)$ is a decreasing function on $\Omega$.</figcaption></figure>

We now focus on the local asymptotic stability of the coexistence equilibrium $E^*$. We start with the case $\beta_1 = 0$, where tumor proliferation is independent of the action of the TAMs population. We then increase the value of $\beta_1 > 0$ to see how this equilibrium behaves.

**Proposition 7** We assume that $\beta_1 = 0$. When the coexistence equilibrium $E^*$ exists, that is when $T^* \in \Omega$, it is always locally asymptotically stable.
*Proof* The Jacobian matrix for $\beta_1 = 0$, evaluated at $E^*$ is given by

$$J^* = \begin{pmatrix} -\frac{\alpha_1}{K_1} T^* & 0 & -\mu_1 T^* \\ \alpha_2 M^* \left(1 - \frac{M^*}{K_2}\right) & -\frac{\alpha_2}{K_2} T^* M^* & 0 \\ \frac{\alpha_3 C^* K_3}{(K_3 + T^*)^2} & 0 & 0 \end{pmatrix}.$$

So, we have the real eigenvalue

$$\lambda_0 = -\frac{\alpha_2}{K_2} T^* M^* < 0.$$

The two other eigenvalues $\lambda_1$ and $\lambda_2$ are solutions of

$$\lambda^2 + \left(\frac{\alpha_1}{K_1} T^*\right) \lambda + \frac{\mu_1 \alpha_3 C^* K_3 T^*}{(K_3 + T^*)^2} = 0.$$

Then, $\lambda_1 + \lambda_2 = -\frac{\alpha_1}{K_1} T^* < 0$ and $\lambda_1\lambda_2 = \frac{\mu_1\alpha_3 C^* K_3 T^*}{(K_3 + T^*)^2} > 0$. So, $\lambda_1$ and $\lambda_2$ have a negative real part. We conclude that for $\beta_1 = 0$ when the equilibrium $E^*$ exists, it is always locally asymptotically stable. In terms of parameters, the existence of $E^*$ is equivalent to $\alpha_3_{\min} < \alpha_3 < \alpha_3_{\max}(\alpha_2)$, where $\alpha_3_{\min}$ and $\alpha_3_{\max}(\alpha_2)$ are given by (5). $\square$

In the next subsection, we explore various aspects related to the coexistence equilibrium point $E^*$, covering issues such as the emergence of oscillatory solutions around $E^*$ and the Hopf bifurcation phenomenon.

## 3.2 Oscillating solutions and Hopf bifurcation

As in the previous analysis of the equilibrium point $E^*$, we continue to assume in this subsection that $\mu_3 = 0$. This simplification, introduced to manage the algebraic complexity of the system, remains valid throughout our study of the coexistence equilibrium $E^*$ and its surrounding dynamics. The exclusion of the term $\mu_3 CT$ does not affect the qualitative behavior of the system, in particular concerning the existence of oscillations and Hopf bifurcations.

We are now concentrating on the appearance of oscillations around $E^*$. Even in the case where $\beta_1 = 0$, these oscillations, although damped, exist. By increasing the parameter $\beta_1 > 0$, these oscillations can be transformed into periodic solutions through a Hopf bifurcation. Let's first examine the existence of oscillations for $\beta_1 = 0$. We introduce the following notation

$$\tilde{\alpha}_1 := \frac{4\alpha_2^2 \alpha_3 K_1 K_3 (\alpha_2 K_1 - \kappa_2)}{\kappa_2 (\alpha_2 K_3 + \kappa_2)^2}. \qquad (7)$$

**Proposition 8** We assume that $\beta_1 = 0$ and we consider two situations:

* Suppose that $\alpha_1 \ge \tilde{\alpha}_1$. Then, there is no oscillating solution ($E^*$ is a sink).

* Suppose that $\alpha_1 < \tilde{\alpha}_1$. Then, there exists $\tilde{T}^* \in (\frac{\kappa_2}{\alpha_2}, K_1)$, such that

- if $T^* \in (\frac{\kappa_2}{\alpha_2}, \tilde{T}^*)$, then there are damped oscillating solutions ($E^*$ is a spiral),
- if $T^* \in (\tilde{T}^*, K_1)$, then, there is no oscillating solution ($E^*$ is a sink).

*Proof* We compute the discriminant $\Delta$ of the polynomial introduced in the proof of Proposition 7,

$$ \lambda^2 + \left(\frac{\alpha_1}{K_1} T^*\right) \lambda + \frac{\mu_1 \alpha_3 C^* K_3 T^*}{(K_3 + T^*)^2} = 0, \quad T^* \in \Omega := \left(\frac{\kappa_2}{\alpha_2}, K_1\right), $$

with

$$ C^* = \frac{\alpha_1}{\mu_1} \left(1 - \frac{T^*}{K_1}\right). $$

Then, for $T^* \in \Omega$, we get

$$ \Delta = \left(\frac{\alpha_1}{K_1} T^*\right)^2 - 4 \frac{\mu_1 \alpha_3 C^* K_3 T^*}{(K_3 + T^*)^2} = \left[\frac{\alpha_1}{K_1^2} T^* - 4 \frac{\alpha_3 K_3}{(K_3 + T^*)^2} \left(1 - \frac{T^*}{K_1}\right)\right] \alpha_1 T^*. $$

We have to study the sign of the function

$$ \bar{\Delta}(T^*) = \frac{\alpha_1}{K_1^2} T^* - 4 \frac{\alpha_3 K_3}{(K_3 + T^*)^2} \left(1 - \frac{T^*}{K_1}\right), \quad T^* \in \Omega. $$

It is equivalent to study, for $T^* \in \Omega$, the sign of the polynomial function

$$ \tilde{\Delta}(T^*) = T^{*3} + 2K_3 T^{*2} + \left(K_3 + 4 \frac{\alpha_3}{\alpha_1} K_1\right) K_3 T^* - 4 \frac{\alpha_3}{\alpha_1} K_1^2 K_3. $$

Let's first consider the function $\tilde{\Delta}$ defined on the interval $[0, +\infty)$. It is clear that $\tilde{\Delta}$ is increasing on $[0, +\infty)$ with $\tilde{\Delta}(0) < 0$ and $\lim_{T^* \to +\infty} \tilde{\Delta}(T^*) = +\infty$. Furthermore,

$$ \tilde{\Delta}(K_1) = K_1^3 + 2K_3 K_1^2 + K_3^2 K_1 > 0. $$

Then, there exists a unique $\tilde{T}^* \in (0, K_1)$ such that $\tilde{\Delta}(\tilde{T}^*) = 0$. We need to check the position of $\tilde{T}^*$ with respect to $\kappa_2/\alpha_2$. This is given by the sign of

$$ \tilde{\Delta} \left( \frac{\kappa_2}{\alpha_2} \right) = \frac{\kappa_2^3}{\alpha_2^3} + 2K_3 \frac{\kappa_2^2}{\alpha_2^2} + K_3^2 \frac{\kappa_2}{\alpha_2} - 4K_1 K_3 \frac{\alpha_3}{\alpha_1} \left( K_1 - \frac{\kappa_2}{\alpha_2} \right). $$

More precisely, we have

$$ \bar{\Delta} \left( \frac{\kappa_2}{\alpha_2} \right) \ge 0 \quad \text{if and only if} \quad \alpha_1 \ge \tilde{\alpha}_1, $$

where $\tilde{\alpha}_1$ is given by (7). We conclude that, if $\alpha_1 \ge \tilde{\alpha}_1$, then $\tilde{\Delta}(T^*) \ge 0$, for all $T^* \in \Omega$. So, the equilibrium $E^*$ is a sink. Suppose now that

$$ \alpha_1 < \tilde{\alpha}_1. $$

We have $\tilde{T}^* \in (\frac{\kappa_2}{\alpha_2}, K_1)$. So, if $T^* \in (\frac{\kappa_2}{\alpha_2}, \tilde{T}^*)$, then the equilibrium $E^*$ is a spiral, and if $T^* \in (\tilde{T}^*, K_1)$, the equilibrium $E^*$ is a sink. $\square$

We suppose now that $\beta_1 > 0$. The Jacobian matrix evaluated on $E^*$ is given by

$$ J^* = \begin{pmatrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & 0 \\ a_{31} & 0 & 0 \end{pmatrix}, \qquad (8) $$
where

$$a_{11} = -\frac{\alpha_1}{K_1} T^*(1 + \beta_1 M^*) < 0,$$

$$a_{12} = \alpha_1 T^* \left(1 - \frac{T^*}{K_1}\right) \beta_1 > 0,$$

$$a_{13} = -\mu_1 T^* < 0,$$

$$a_{21} = \alpha_2 M^* \left(1 - \frac{M^*}{K_2}\right) > 0,$$

$$a_{22} = -\kappa_2 + \alpha_2 T^* \left(1 - \frac{2M^*}{K_2}\right) < 0,$$

$$a_{31} = \frac{\alpha_3 C^* K_3}{(K_3 + T^*)^2} > 0.$$

Then, the characteristic equation associated to the matrix (8) is given by

$$\lambda^3 + q_1\lambda^2 + q_2\lambda + q_3 = 0,$$

with $q_1 := -(a_{11} + a_{22})$, $q_2 := a_{11}a_{22} - a_{12}a_{21} - a_{13}a_{31}$ and $q_3 := a_{13}a_{22}a_{31}$. We can use the Routh-Hurwitz criterion to evaluate the stability of the equilibrium $E^*$. It is clear that $q_1 > 0$, $q_2 > 0$ and $q_3 > 0$. The only missing condition is the positivity of the expression

$$Q := q_1 q_2 - q_3 = -a_{11}^2 a_{22} + a_{11} a_{12} a_{21} + a_{11} a_{13} a_{31} - a_{11} a_{22}^2 + a_{12} a_{21} a_{22}.$$

Let's consider $Q := Q(\beta_1)$ as a function of $\beta_1 \ge 0$. Only the coefficients $a_{11}$, $a_{12}$ and $a_{31}$ depend on the parameter $\beta_1$, through the expression

$$C^* = \frac{\alpha_1}{\mu_1} \left(1 - \frac{T^*}{K_1}\right) (1 + \beta_1 M^*),$$

because $T^*$ and $M^*$ are independent on $\beta_1$. As $Q(0) > 0$, thanks to the continuity there exists $\beta_1^* > 0$, such that $Q(\beta_1) > 0$ for all $\beta_1 \in [0, \beta_1^*)$. This threshold $\beta_1^*$ can be a bifurcation parameter that destabilizes the equilibrium $E^*$, with the appearance of periodic solutions. If this threshold exists, it will be such that $Q(\beta_1^*) = 0$. Due to the large number of parameters, it is difficult to prove analytically the existence of $\beta_1^*$.

We numerically analyze the situations described above. We consider the two cases $\alpha_1 \ge \tilde{\alpha}_1$ and $\alpha_1 < \tilde{\alpha}_1$, with $\tilde{\alpha}_1$ given by (7). We want to vary $T^* \in \Omega := (\kappa_2/\alpha_2, K_1)$, which also means varying $\kappa_3$ in a certain interval, the parameter that models the rate of apoptosis of CAR T-cells and the effect of immunosuppression. The condition of existence of $E^*$, given by

$$\alpha_{3 \min} := \kappa_3 \left(1 + \frac{K_3}{K_1}\right) < \alpha_3 < \alpha_{3 \max}(\alpha_2) := \kappa_3 \left(1 + \frac{\alpha_2}{\kappa_2} K_3\right),$$

can also be written as follows

$$\frac{\alpha_3 \kappa_2}{\kappa_2 + \alpha_2 K_3} < \kappa_3 < \frac{\alpha_3 K_1}{K_1 + K_3}.$$

Then, by considering $T^*: \kappa_3 \mapsto T^*(\kappa_3) = \frac{\kappa_3 K_3}{\alpha_3 - \kappa_3}$ as a function of $\kappa_3$ on the interval $(\alpha_3\kappa_2/(\kappa_2 + \alpha_2 K_3), \alpha_3 K_1/(K_1 + K_3))$, we obtain an increasing function (see Figure 5). In accordance with the set of parameters presented in Table 1 and taking $\alpha_1 = 0.69 \text{ day}^{-1}$, we determine the specific value for which $\tilde{\Delta}(\tilde{T}^*) = 0$. This value is given by $\tilde{T}^* = 3.86 \cdot 10^8$. We also obtain from (7) that $\tilde{\alpha}_1 = 3.85 \cdot 10^4$. For the values of $\alpha_1$
shown in Table 1, only the condition $\alpha_1 < \tilde{\alpha}_1$ is satisfied, since $\alpha_1 \in [0.69, 0.97]$ to be biologically acceptable [17]. Then, considering the value $\alpha_1 = 0.69 < \tilde{\alpha}_1$, we vary

$$ \kappa_3 \in \left( \frac{\alpha_3 \kappa_2}{\kappa_2 + \alpha_2 K_3}, \frac{\alpha_3 K_1}{K_1 + K_3} \right), $$

such that $T^*(\kappa_3) = \frac{\kappa_3 K_3}{\alpha_3 - \kappa_3}$ crosses the threshold $\tilde{T}^*$ (Figure 5).

<figure><img src="image_7.png" /><figcaption>Fig. 5 The gray area, $T^* \in (\kappa_2/\alpha_2, \tilde{T}^*)$, is the region where the equilibrium $E^*$ is a stable spiral, and the part above the green dashed line, $T^* \in (\tilde{T}^*, K_1)$, is the region where this equilibrium is a stable sink.</figcaption></figure>

With the choice of $\alpha_1 = 0.69 < \tilde{\alpha}_1$ and following Proposition 8, two possibilities for the equilibrium point $E^*$ arise. Employing the parameters from Table 1, with $\beta_1 = 2.3 \cdot 10^{-10}$, and selecting two values for the parameter $\kappa_3$, i.e. $\kappa_3 = 0.5$ and $\kappa_3 = 0.86$, to satisfy the two situations presented in Figure 5, we obtained the following results, as shown in Figure 6.
<figure><img src="image_8.png" /><figcaption>**Fig. 6** In (a), we can observe that when $\kappa_3 = 0.5 \text{ day}^{-1}$, the condition $T^* \in (\kappa_2/\alpha_2, \tilde{T}^*)$ is satisfied, and the equilibrium $E^*$ is a stable spiral. In (b), we observe that when $\kappa_3 = 0.86 \text{ day}^{-1}$, the condition $T^* \in (\tilde{T}^*, K_1)$ is satisfied, and we have that the equilibrium $E^*$ is a stable sink (without oscillations around $T^*$).</figcaption></figure>

We run simulations to assess the influence of parameter $\beta_1 \ge 0$ on the appearance of periodic solutions around equilibrium $E^*$ (Hopf bifurcation). In this scenario, we focus on the case where the equilibrium $E^*$ is a stable spiral for $\beta_1 = 0$, as shown in the grey region of Figure 5. We then increase the value of $\beta_1 > 0$ until the spiral around equilibrium $E^*$ becomes a center (Hopf bifurcation). We use the parameters of Table 1, with $\alpha_1 = 0.69$, $\kappa_3 = 0.35 \text{ day}^{-1}$, and the initial condition $C(0) = 10^7$ cells. In this case, the condition $\alpha_1 < \tilde{\alpha}_1$ is satisfied. We stop the simulations at $t = 1000$ days. Our results (only the component $M(t)$ of the solution) are shown in Figure 7.

<figure><img src="image_9.png" /><figcaption>**Fig. 7** In this figure, we show only the component $M(t)$ of the solution. When $\beta_1 = 0$ in panel (a), the solutions oscillate and converge to the equilibrium $E^*$ (red star). An increase in the parameter $\beta_1$ from 0 leads to the appearance of periodic oscillations through a Hopf bifurcation at $\beta_1 = 5.9 \cdot 10^{-11} \text{ cells}^{-1}$ (panel (b)).</figcaption></figure>

Let's now explain the conditions that give rise to a Hopf bifurcation. To do this, we again consider the characteristic equation associated with the equilibrium $E^*$, in
the case $\beta_1 > 0$,

$$\lambda^3 + q_1\lambda^2 + q_2\lambda + q_3 = 0,$$

with $q_1 := -(a_{11} + a_{22})$, $q_2 := a_{11}a_{22} - a_{12}a_{21} - a_{13}a_{31}$ and $q_3 := a_{13}a_{22}a_{31}$. As mentioned before, we plot the function

$$Q(\beta_1) := q_1q_2 - q_3 = -a_{11}^2a_{22} + a_{11}a_{12}a_{21} + a_{11}a_{13}a_{31} - a_{11}a_{22}^2 + a_{12}a_{21}a_{22}.$$

Figure 8 shows the results obtained.

<figure><img src="image_10.png" /><figcaption>Fig. 8 The function $Q(\beta_1) := q_1q_2 - q_3$ is positive for $\beta_1$ less than $\beta_1^* \approx 5.9 \cdot 10^{-11}$ cells$^{-1}$ and negative when $\beta_1$ is greater than this threshold. The parameters used in this figure are shown in Table 1, with $\alpha_1 = 0.69$ and $\kappa_3 = 0.35$ day$^{-1}$.</figcaption></figure>

If we increase the value of parameter $\beta_1$ from 0, we have that the function $Q(\beta_1) := q_1q_2 - q_3$ is positive until a certain threshold $\beta_1 = \beta_1^*$, which is approximately $5.9 \cdot 10^{-11}$ cells$^{-1}$. The function $Q(\beta_1)$ becomes negative beyond this value. We have $Q(\beta_1^*) = 0$ and $\beta_1^*$ is responsible for the appearance of periodic solutions (Hopf bifurcation) as observed in Figure 7.

Recent experiments have revealed that periodic oscillations in the cell population can occur during melanoma treatment [39]. Our study does not take place under the same conditions as these experiments. However, this study indicates that periodic oscillations in melanoma treatment may exist.

In the next section, we perform a sensitivity analysis to assess the impact of key parameters, including $\mu_3 > 0$, on the behaviour of the model.

## 4 Sensitivity analysis

According to Saltelli, Tarantola, and Chan [40], sensitivity analysis (SA) aims to determine how a model depends on its associated parameters, thus contributing to the understanding of behavior, consistency between the model and the phenomenon being modeled, and how the different parts of the model interact. In this way, it is possible to identify and classify the parameters that have the greatest impact on model output.

The first-order sensitivity index, $S1$, indicates the individual contribution of each input variable to the variation in output. It captures the isolated effect of each input
variable on total model variation and is relevant for discerning which variables have
a significant independent impact [40, 41].

The total sensitivity index, *ST*, is a measure that encompasses both the individual effects of each variable and higher-order interactions. Contrary to the first-order sensitivity index *S1*, the total sensitivity index *ST* encompasses the contribution of a variable across all possible orders of interaction. This index is important for determining the ability to fix a parameter at any point in its range of variability without causing a significant impact on the outcome, *i.e.*, the model output. For instance, if *ST* ≈ 0 for a specific parameter, it is possible to assign it any value within its uncertainty range without substantially affecting the variation in the model output. In this context, total sensitivity indices are particularly relevant for factor-fixing configurations, as they provide valuable insights into the relative influence of variables and their interactions on the model outcome. Parameters whose first-order and total sensitivity indices are greater than 0.05 are considered significant for the model output [41–43].

**4.1 Sensitivity analysis of the mathematical model**

CAR T-cell therapy is understood as a consequence of the interaction of various
phenomena unfolding in diverse biological populations. However, the comprehension
of some side effects, such as cytokine generation and interaction with endogenous
immune cells, such as macrophages, remains limited, with many mechanisms still
unknown.

The proposed mathematical model, as presented in the system (1), incorporates
specific parameters that are presumed to contain uncertainties. In this context, sen-
sitivity analysis aims to examine the most important parameters of the model, i.e.,
those capable of modifying the dynamics of the system by varying their values. Sen-
sitivity analysis was conducted using the Python programming language with the
SALib library, which offers implementations of the eFAST method, [40, 41, 44, 45].

In our analyses, we considered the parameters of the system (1) introducing a
perturbation of 20% to the reference values presented in Table 1. This perturbation
allowed us to observe how these parameters influence the dynamics of tumor cells,
TAMs, and CAR T-cells over a simulated period of one hundred days. Particularly,
the initial conditions considered in the analyses were T(0) = 2 · 10⁶, M(0) = 10⁶,
and C(0) = 10⁷ cells. Additionally, we set α₁ = 0.69 day⁻¹, β₁ = 2.3 · 10⁻¹⁰ cells⁻¹,
κ₃ = 0.35 day⁻¹ and µ₃ = 6.752 · 10⁻¹⁰ as reference values. The results obtained
are presented in Figure 9.
<figure><img src="image_11.png" /><figcaption>**Fig. 9** Sensitivity analysis (SA) of the parameters from Table 1 was conducted using the Python programming language and the SALib package. Blue bars represent first-order sensitivity indices (S1), while orange bars represent total sensitivity indices (ST).</figcaption></figure>

Table 2 presents the values for first-order and total sensitivity indices, along with the 95% confidence interval used in the analyses.

<table><caption>**Table 2** First-order sensitivity indices (S1) and total-order sensitivity indices (ST) for the populations of tumor cells, TAMs, and CAR T-cells with a 95% confidence interval.</caption><thead><tr><th rowspan="2">Parameters</th><th colspan="4">Tumor</th><th colspan="4">TAMs</th><th colspan="4">CAR T-cells</th></tr><tr><th>S1</th><th>ST</th><th>S1 IC 95%</th><th>ST IC 95%</th><th>S1</th><th>ST</th><th>S1 IC 95%</th><th>ST IC 95%</th><th>S1</th><th>ST</th><th>S1 IC 95%</th><th>ST IC 95%</th></tr></thead><tbody><tr><td>α<sub>1</sub></td><td>0.001805</td><td>0.142015</td><td>0.001433</td><td>0.014867</td><td>0.006222</td><td>0.759713</td><td>0.001532</td><td>0.013085</td><td>0.032984</td><td>0.869870</td><td>0.001559</td><td>0.011892</td></tr><tr><td>α<sub>2</sub></td><td>0.000058</td><td>0.071421</td><td>0.001625</td><td>0.012891</td><td>0.008926</td><td>0.047743</td><td>0.001460</td><td>0.011368</td><td>0.000036</td><td>0.049753</td><td>0.001394</td><td>0.013937</td></tr><tr><td>α<sub>3</sub></td><td>0.186196</td><td>0.755235</td><td>0.001754</td><td>0.012913</td><td>0.057875</td><td>0.736822</td><td>0.001509</td><td>0.014630</td><td>0.010143</td><td>0.843995</td><td>0.001271</td><td>0.013829</td></tr><tr><td>β<sub>1</sub></td><td>0.000495</td><td>0.091149</td><td>0.001609</td><td>0.012478</td><td>0.005323</td><td>0.192882</td><td>0.001618</td><td>0.013421</td><td>0.000269</td><td>0.226524</td><td>0.001478</td><td>0.011161</td></tr><tr><td>β<sub>2</sub></td><td>0.000746</td><td>0.113515</td><td>0.001339</td><td>0.011183</td><td>0.002702</td><td>0.250715</td><td>0.001463</td><td>0.013534</td><td>0.020130</td><td>0.340476</td><td>0.001423</td><td>0.013967</td></tr><tr><td>μ<sub>1</sub></td><td>0.000037</td><td>0.079978</td><td>0.001530</td><td>0.013890</td><td>0.016737</td><td>0.153686</td><td>0.001585</td><td>0.012526</td><td>0.000174</td><td>0.152075</td><td>0.001590</td><td>0.013609</td></tr><tr><td>μ<sub>2</sub></td><td>0.010079</td><td>0.231176</td><td>0.001549</td><td>0.013283</td><td>0.000540</td><td>0.068447</td><td>0.001662</td><td>0.013426</td><td>0.000829</td><td>0.092897</td><td>0.001487</td><td>0.012306</td></tr><tr><td>μ<sub>3</sub></td><td>0.114070</td><td>0.674413</td><td>0.001417</td><td>0.012463</td><td>0.084677</td><td>0.430056</td><td>0.001413</td><td>0.013504</td><td>0.004196</td><td>0.383694</td><td>0.001392</td><td>0.012953</td></tr><tr><td>K<sub>1</sub></td><td>0.010250</td><td>0.207585</td><td>0.001420</td><td>0.010248</td><td>0.008481</td><td>0.167951</td><td>0.001449</td><td>0.013162</td><td>0.000441</td><td>0.203998</td><td>0.001784</td><td>0.012409</td></tr><tr><td>K<sub>2</sub></td><td>0.000576</td><td>0.054239</td><td>0.001469</td><td>0.011729</td><td>0.019584</td><td>0.210984</td><td>0.001548</td><td>0.013075</td><td>0.000046</td><td>0.228057</td><td>0.001568</td><td>0.014163</td></tr><tr><td>K<sub>3</sub></td><td>0.008598</td><td>0.200026</td><td>0.001458</td><td>0.014531</td><td>0.008137</td><td>0.079146</td><td>0.001557</td><td>0.014381</td><td>0.000093</td><td>0.141011</td><td>0.001557</td><td>0.012265</td></tr></tbody></table>

Figure 9(a) presents the first-order sensitivity indices (S1) and total sensitivity indices (ST) in the population of tumor cells. The results underscore the importance of parameters specific to tumor dynamics. The two most influential parameters are the proliferation rate of CAR T-cells in response to tumor antigen (α3) and the resistance of CAR T-cells to apoptosis (κ3), emphasizing their significance for tumor
dynamics. The high influence of these parameters suggests that the treatment terms are a key determinant for control growth. The strategies that enhance the proliferation of CAR T-cells in response to the antigen and promote the resistance of these cells may enhance treatment persistence and efficacy.

Figure 9(b) shows the results related to the dynamics of TAMs. The TAMs population is positively influenced by the increase in $\alpha_1$, suggesting that macrophage proliferation is associated with tumor activity. Therapeutically, modulating the proliferative response of macrophages can be explored to impact tumor dynamics. The analysis also highlights the influence of the parameter $\alpha_3$, both through its individual contribution (S1) and total contribution (ST). This indicates an indirect influence of the CAR T-cell response to the tumor antigen on TAMs. Finally, the apoptosis rate of CAR T-cells ($\kappa_3$) also showed a high sensitivity index. These results indicate a complex interaction between the tumor and TAMs, emphasizing the need for strategies that modify these interactions to reduce the pro-tumoral support provided by TAMs cells and enhance the efficacy of CAR-T immunotherapy.

The results related to CAR T-cells shown in Figure 9(c) underscore the crucial importance of the specific antigen proliferation rate ($\alpha_3$) in antitumoral efficacy. The positive influence of $\alpha_3$ suggests that therapeutic strategies promoting a robust and targeted response to the antigen may be crucial for the success of CAR T-therapy. Additionally, the natural proliferation rate of tumor cells ($\alpha_1$) also influences the dynamics, indicating that strategies to improve antigen recognition for CAR T-cell expansion can be explored to optimize therapy.

The results of the sensitivity analysis reveal crucial information for understanding the dynamic interactions among cellular populations in a therapeutic context, highlighting the influence of various parameters on the behaviors of tumor cells, TAMs, and CAR T-cells. The intrinsic tumor proliferation rate ($\alpha_1$) emerges as a significant factor in the overall dynamics, positively affecting all three cellular populations. The antigen-specific response, represented by the increase in $\alpha_3$, is critical for the effectiveness of CAR T-cells against the tumor. Additionally, the pro-tumoral interactions, represented by $\beta_1$, underscore the complexity of the relationships between tumor cells and TAMs, highlighting potential therapeutic targets. The cytotoxic effect of CAR T-cells ($\mu_1$) is identified as a promising strategy for reducing the tumor population. These findings provide valuable insights for the development of personalized therapeutic approaches, directing the intricate interactions among tumor cells, TAMs, and CAR T-cells to optimize clinical outcomes in cancer treatment.

# 5 Computational simulations

Several models involving CAR T-cell kinetics have been developed to present its different phases, namely distribution, expansion, contraction, and persistence. The distribution phase occurs after administration of the immunotherapeutic dose, during which CAR T-cells are disseminated in various tissues of the body, resulting in an initial decrease in the population. After this decline, the cells enter the expansion phase due to interaction with the tumor antigen, reaching a peak of expansion representing the maximum population of CAR T-cells. Finally, CAR T-cells undergo processes such as antigen loss, exhaustion, and apoptosis, which contribute to a subsequent decrease in the population, which enters the contraction and persistence phase. After the contraction and persistence phases, effector CAR T-cells lose their ability to expand [46].
In this section, we proposed several scenarios to simulate the efficacy of CAR T-cell immunotherapy against melanoma and to assess the impact of immunosuppression on treatment. In addition, we examined the characteristic phases of CAR T-cells throughout the therapeutic process.

## 5.1 Scenario 1: immunosuppressive effects have an impact on CAR T-cell immunotherapy

In this scenario, we evaluate the impact of varying the parameter $\mu_3$ in the context of CAR T-cell immunotherapy and identify a bifurcation point that significantly alters the population dynamics.

Using the parameters from Table 1, specifically $\alpha_1 = 0.69 \text{ day}^{-1}$, $\beta_1 = 2.3 \cdot 10^{-10}$ cells$^{-1}$, $\kappa_3 = 0.35 \text{ day}^{-1}$, we set the initial CAR T-cell population at $C(0) = 5 \cdot 10^5$ cells and simulate the system across a range of $\mu_3$ over the interval $[5 \cdot 10^{-10}, 5 \cdot 10^{-9}]$ cells·day$^{-1}$, to observe how the final tumor cell population $T(t)$, and, consequently, the equilibrium points, respond to different levels of immunosuppression. Figure 10 shows the behavior of the tumor population $T(t)$ for different values of $\mu_3$.

The analysis shows a bifurcation in the system's behavior as $\mu_3$ increases. Initially, for lower values of $\mu_3$, the system tends to stabilize at the coexistence equilibrium point $E^*$ (stable spiral). However, as $\mu_3$ approaches a critical value (approximately $\mu_3 \approx 6.752 \cdot 10^{-10}$ cell·day$^{-1}$), the system undergoes a bifurcation, after which the tumor cell population $T(t)$ no longer converges to the coexistence equilibrium $E^*$.

<figure><img src="image_12.png" /><figcaption>**Fig. 10** Dynamics of tumor cell population $T(t)$ under varying 100 values of $\mu_3$ with $t = 60$ days. (a) The time evolution of the tumor population $T(t)$ for different values of $\mu_3$. The figure shows that as $\mu_3$ increases, the tumor cells either stabilize or exhibit oscillatory behavior before reaching equilibrium. (b) The bifurcation diagram shows the final tumor population $T(t)$ as a function of $\mu_3$. A bifurcation occurs around $\mu_3 \approx 6.752 \cdot 10^{-10}$ cells·day$^{-1}$, where the system switches from a stable tumor population to tumor persistence, indicated by the vertical red dashed line.</figcaption></figure>

Beyond this bifurcation point, the system moves into the equilibrium $E_2$, where CAR T-cells are no longer able to develop effectively due to the high immunosuppressive effects, and the tumor population dominates. This is evident in Figure 10(b), where, for higher values of $\mu_3$, the tumor population stabilizes in the carrying capacity $K_1$, indicating the dominance of the tumor and the failure of treatment. The bifurcation therefore represents the threshold at which equilibrium passes from the
coexistence equilibrium point, $E^*$, to treatment failure and tumor dominance, $E_2$, under the effect of increasing immunosuppression.

Using four specific values for $\mu_3$, we can observe the impact on the dynamics of all cell populations with $t = 30$ days. Figure 11 shows the results obtained.

$$ 
\begin{aligned}
\mu_3 &= 10^{-9} (\text{cells} \cdot \text{day})^{-1} \\
\mu_3 &= 5 \cdot 10^{-10} (\text{cells} \cdot \text{day})^{-1}
\end{aligned}
 $$

$$ 
\begin{aligned}
\mu_3 &= 5 \cdot 10^{-9} (\text{cells} \cdot \text{day})^{-1} \\
\mu_3 &= 10^{-11} (\text{cells} \cdot \text{day})^{-1}
\end{aligned}
 $$

<figure><img src="image_13.png" /><figcaption>Fig. 11 For low values of $\mu_3$, the initial dose $C(0) = 5 \cdot 10^5$ cells is able to expand and kill the tumor population.</figcaption></figure>

Given that the parameter $\mu_3$ models the effects of immunosuppression of TAMs, we can note that for high values of $\mu_3$, the CAR T-cell population takes some time to reach its maximum expansion, at which it begins to reduce the number of tumor cells, and this reveals the model's ability to demonstrate phases of expansion, contraction, and persistence (represented by the red lines), suggesting its effectiveness in describing CAR T-cell dynamics. When immunosuppression is high, as with $\mu_3 = 10^{-9}$ and $\mu_3 = 5 \cdot 10^{-9}$, CAR T-cells are unable to expand, leading to treatment failure.

These observations related by immunosuppression may represent, for example, the interaction of PD-L1 ligand present on tumor cells and TAMs with the PD-1 receptor, present on CAR T-cells [28, 34, 47, 48]. Furthermore, these results corroborate data in the literature regarding the presence of TAMs in the tumor microenvironment and its negative influence when using CAR T-cell immunotherapy [49, 50]. Immunosuppressive factors such as interleukin-10 (IL-10) produced by TAMs are also responsible for the failure of CAR T-cells against solid tumors such as melanoma [51].

## 5.2 Scenario 2: increasing the dose of immunotherapy may improve tumor control

In this scenario, we simulate an immunosuppressive tumor, as in the case of cutaneous melanoma. Therefore, we perform simulations by varying the initial doses of CAR T-cells to verify the dose that allows longer control of tumor progression. We consider the parameter values and initial conditions of Table 1, with $\alpha_1 = 0.69 \text{ day}^{-1}$, $\beta_1 = 2.3 \cdot 10^{-10} \text{ cells}^{-1}$, $\kappa_3 = 0.35 \text{ day}^{-1}$ and $\mu_3 = 6.752 \cdot 10^{-10}$ (characterizing the equilibrium point $E_2$). We vary the initial doses of CAR T-cells between $C(0) = 10^6$ cells and $C(0) = 9 \cdot 10^7$ cells. The results are shown in Figures 12.
$$ 
\begin{aligned}
C(0) &= 10^6 \text{ cells} \\
C(0) &= 7 \cdot 10^7 \text{ cells}
\end{aligned}
 $$

$$ 
\begin{aligned}
C(0) &= 10^7 \text{ cells} \\
C(0) &= 9 \cdot 10^7 \text{ cells}
\end{aligned}
 $$

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

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

<figure><img src="image_16.png" /><figcaption>Fig. 12 Time evolution of tumor cells, TAMs, and CAR T-cells for different initial immunotherapeutic doses. The shaded area highlights the range where the tumor population is suppressed below $10^2$ cells.</figcaption></figure>

The results show that, under conditions of high tumor immunosuppression, in particular an immunosuppressive tumor microenvironment, resulting from TAMs-related effects, and CAR T-cell death due to depletion, and to improve treatment efficiency, initial immunotherapeutic doses need to be adjusted. Following the results presented by Barros et al. [27], we adopted a threshold of undetectable tumor cells set at $10^2$ cells.

In Figure 12(a), we observe that from high initial doses of CAR T-cells, tumor cells decrease rapidly, allowing tumor control for 12.4 days for $C(0) = 7 \cdot 10^7$ cells and 21.7 days for $C(0) = 9 \cdot 10^7$ cells, after which the tumor grows to carrying capacity, which usually represents patient death. Figure 12(c) shows that the CAR T-cell population is reduced due to the decrease in the tumor cell population. The accelerated decrease in tumor cells also implies a decrease in the TAMs population, which is to be expected since both populations depend on each other according to the mathematical modeling employed.

Our *in silico* simulations modeled the use of high doses of CAR T-cells, and the results obtained showed that this could overcome the immunosuppressive mechanisms associated with tumor cells and the TAMs population. However, the use of high doses of CAR T-cells is technically not possible due to the complexity of producing them in large quantities [6, 26, 28]. Moreover, in the context of solid tumors such as melanoma, even with the administration of a substantial quantity of CAR T-cells to the patient, only a limited number of cells manage to cross various physical and chemical barriers to localize and combat tumor cells. Our results also showed that, despite the injection of high immunotherapeutic doses, within 30 to 40 days the tumor population resumed its growth, reaching carrying capacity, indicating treatment failure. Anticancer therapies aim to increase patient survival time as much as possible. One strategy for increasing the efficacy of CAR T-cells, and consequently patient survival time, would be to combine them with Immune-Checkpoint Immunotherapy (ICI) such as PD-1/PD-L1 axis blockers [28, 47].
## 5.3 Scenario 3: evaluate the impact of TAMs-associated factors such as $β_1$ on tumor growth in the context of CAR T-cell immunotherapy

In the last scenario, our objective is to simulate the impact of tumor growth factors produced by TAMs, such as growth factors, pro-angiogenic molecules, and immuno-suppressive factors, all represented by the parameter $β_1$, in the context of CAR T-cell immunotherapy. To this end, we carried out simulations by varying the initial doses of CAR T-cells and increasing the parameter $β_1$. To measure the interference of these factors in treatment, we calculated the time of tumor control (in days) for two selected values of $β_1$, i.e. $10^{-7}$ and $2 \cdot 10^{-6}$ cells$^{-1}$ (characterizing the equilibrium point $E_2$). We considered the parameter values and initial conditions presented in Table 1, with $α_1 = 0.69$ day$^{-1}$, $κ_3 = 0.35$ day$^{-1}$ and $μ_3 = 6.752 \cdot 10^{-10}$. We varied the initial doses of CAR T-cells between $C(0) = 5 \cdot 10^6$ cells and $C(0) = 9 \cdot 10^7$ cells. The results are shown in Figure 13.

<figure><img src="image_17.png" /><figcaption>**Fig. 13** Tumor growth, due to the pro-tumor action of TAMs, can reduce the efficacy of immunotherapy, even with the use of high doses.</figcaption></figure>

In the top panel of Figure 13, we can observe that by using $β_1 = 10^{-7}$, tumor cells remained below the threshold considered as control, only when exposed to doses of $7 \cdot 10^7$ and $9 \cdot 10^7$ cells. Under these conditions, the tumor population remained below $10^2$ cells for 12.1 and 21.4 days, respectively. Examining the bottom panel of Figure 13, we see once again that the same immunotherapeutic doses kept tumor cells below $10^2$. However, for the dose of $7 \cdot 10^7$ cells, the control period was 2.1 days,
and for the dose of $9 \cdot 10^7$ cells, the period extended to 15.2 days. Our results in *silico* support information in the existing literature concerning the impact of tumor-associated macrophages (TAMs) on CAR T-cell immunotherapy and, consequently, on tumor promotion [9, 10, 13, 14].

## 6 Conclusions

With this work, the scope of existing models of tumor growth with CAR T-cell immunotherapy has been extended by incorporating tumor-associated macrophages and a specific cancer variant, melanoma. After constructing a system of ordinary differential equations to study the efficacy of CAR T-cell immunotherapy against cutaneous melanoma, taking into account the TAMs population, the system's equilibrium points and their asymptotic local stability were examined using analytical and numerical methods. Local stability conditions were established, with particular emphasis on the coexistence equilibrium. A set of results has been established to ensure local asymptotic stability when $\beta_1 = 0$ and $\mu_3 = 0$. We have shown that around the coexistence equilibrium, oscillations, including periodic oscillations (Hopf bifurcation), occur depending on certain parameters of the proposed model.

In our study, we investigated the impact of CAR T-cell immunotherapy in a model including melanoma and the TAMs population, to establish therapeutic approaches for this type of cancer. Using *in silico* simulations, we observed the expansion, contraction, and persistence phases of CAR T-cells, as well as the cytotoxic effects of immunotherapy against cutaneous melanoma. The proposed model enabled us to observe immunosuppression caused by tumor cells and the TAMs population during immunotherapeutic interventions. These results are in agreement with the findings of the literature, highlighting the immunosuppression induced by tumor cells and the TAMs population in the tumor microenvironment. Immunosuppressive factors produced by the TAMs population appear to contribute to treatment inefficiency. Furthermore, the analyses indicate that, due to tumor- and TAMs-induced immunosuppression, tumor control is possible for a specific period, provided the initial dose of CAR T-cells is adjusted. In high-dose situations, an initial rapid decrease in the tumor and CAR T-cell population was observed, followed by a resumption of tumor growth and attainment of the tumor's carrying capacity.

**Acknowledgments.** This study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior — Brazil (CAPES) — Finance Code 001. MA thanks the Inria international associate team MoCoVec and STIC AmSud BIO-CIVIP.

## Declarations

**Conflict of interest/Competing interests.** The authors declare that there is no conflict/competing of interest.

**Authors' contributions.** GR, MA, and PFAM conceived the mathematical model. GR, JGS, MA, and PFAM did the mathematical analysis and numerical simulations. GR and PFAM wrote the manuscript and arranged the figures and table. GR, JGS, MA, and PFAM contributed to reviewing and editing the paper. All authors read and approved the final manuscript.
References

[1] Arnold, M., Singh, D., Laversanne, M., Vignat, J., Vaccarella, S., Meheus, F., Cust, A.E., de Vries, E., Whiteman, D.C., Bray, F.: Global burden of cutaneous melanoma in 2020 and projections to 2040. JAMA Dermatol. **158**(5), 495–503 (2022). https://doi.org/10.1001/jamadermatol.2022.0160

[2] American Cancer Society: Key Statistics for Melanoma Skin Cancer. Online. Accessed February 21, 2024. (2024). https://www.cancer.org/cancer/types/melanoma-skin-cancer/about/key-statistics.html

[3] Gray-Schopfer, V., Wellbrock, C., Marais, R.: Melanoma biology and new targeted therapy. Nature **445**(7130), 851–857 (2007). https://doi.org/10.1038/nature05661

[4] Rohaan, M.W., Wilgenhof, S., Haanen, J.B.: Adoptive cellular therapies: the current landscape. Virchows Arch **474**, 449–561 (2019). https://doi.org/10.1007/s00428-018-2484-0

[5] Zhao, Z., Chen, Y., Francisco, N.M., Zhang, Y., Wu, M.: The application of CAR T-cell therapy in hematological malignancies: advantages and challenges. Acta Pharm. Sin. B. **8**(4), 539–551 (2018). https://doi.org/10.1016/j.apsb.2018.03.001

[6] Simon, B., Uslu, U.: CAR T-cell therapy in melanoma: A future success story? Exp. Dermatol. **27**(12), 1315–1321 (2018). https://doi.org/10.1111/exd.13792

[7] Hanahan, D., Weinberg, R.A.: The hallmarks of cancer. Cell **100**(1), 57–70 (2000). https://doi.org/10.1016/S0092-8674(00)81683-9

[8] Mantovani, A., Marchesi, F., Malesci, A., Laghi, L., Allavena, P.: Tumour-associated macrophages as treatment targets in oncology. Nat. Rev. Clin. Oncol. **14**(7), 399–416 (2017). https://doi.org/10.1038/nrclinonc.2016.217

[9] Vitale, I., Manic, G., Coussens, L.M., Kroemer, G., Galluzzi, L.: Macrophages and metabolism in the tumor microenvironment. Cell Metab. **30**(1), 36–50 (2019). https://doi.org/10.1016/j.cmet.2019.06.001

[10] Yang, Q., Guo, N., Zhou, Y., Chen, J., Wei, Q., Han, M.: The role of tumor-associated macrophages (TAMs) in tumor progression and relevant advance in targeted therapy. Acta Pharm. Sin. B. **10**(11), 2156–2170 (2020). https://doi.org/10.1016/j.apsb.2020.04.004

[11] Mojsilovic, S.S., Mojsilovic, S., Villar, V.H., Santibanez, J.F.: The
metabolic features of tumor-associated macrophages: Opportunities for
immunotherapy? Anal. Cell. Pathol. **2021**, *12* (2021). https://doi.org/10.
1155/2021/5523055

[12] Shen, M., Du, Y., Ye, Y.: Tumor-associated macrophages, dendritic cells, and neutrophils: biological roles, crosstalk, and therapeutic relevance. Med. Rev. **1**(2), 222–243 (2021). https://doi.org/10.1515/mr-2021-0014

[13] Mäkitie, T., Summanen, P., Tarkkanen, A., Kivelä, T.: Tumor-infiltrating macrophages (CD68+ cells) and prognosis in malignant uveal melanoma. Invest. Ophth. Vis. Sci. **42**(7), 1414–1421 (2001) https://iovs.arvojournals.org/article.aspx?articleid=2199967

[14] Chen, P., Huang, Y., Bong, R., Ding, Y., Song, N., Wang, X., Song, X., Luo, Y.: Tumor-associated macrophages promote angiogenesis and melanoma growth via adrenomedullin in a paracrine and autocrine manner: role of macrophage-derived adrenomedullin in melanoma. Clin. Cancer Res. **17**(23), 7230–7239 (2011). https://doi.org/10.1158/1078-0432.CCR-11-1354

[15] Liu, Z., Zhou, Z., Dang, Q., Xu, H., Lv, J., Li, H., Han, X.: Immuno-suppression in tumor immune microenvironment and its optimization from CAR-T cell therapy. Theranostics **12**(14), 6273–6290 (2022). https://doi.org/10.7150/thno.76854

[16] Kuznetsov, V.A., Makalkin, I.A., Taylor, M.A., Perelson, A.S.: Nonlinear dynamics of immunogenic tumors: parameter estimation and global bifurcation analysis. B. Math. Biol. **56**(2), 295–321 (1994). https://doi.org/10.1016/S0092-8240(05)80260-5

[17] Eftimie, R., Hamam, H.: Modelling and investigation of the CD4+ T cells – Macrophages paradox in melanoma immunotherapies. J. Theor. Biol. **420**, 82–104 (2017). https://doi.org/10.1016/j.jtbi.2017.02.022

[18] de Pillis, L.G., Radunskaya, A.E., Wiseman, C.L.: A validated mathematical model of cell-mediated immune response to tumor growth. Cancer Res. **65**(17), 7950–7958 (2005). https://doi.org/10.1158/0008-5472.CAN-05-0564

[19] de Pillis, L.G., Gu, W., Radunskaya, A.E.: Mixed immunotherapy and chemotherapy of tumors: modeling, applications and biological interpretations. J. of Theor. Biol. **238**(4), 841–862 (2006). https://doi.org/10.1016/j.jtbi.2005.06.037

[20] de Pillis, L., Eladdadi, A., Radunskaya, A.: Modeling cancer-immune responses to therapy. J. Pharmacokinet. Phar. **41**(5), 461–478 (2014). https://doi.org/10.1007/s10928-014-9386-9
[21] Eftimie, R., Bramson, J.L., Earn, D.J.D.: Modeling anti-tumor Th1 and Th2 immunity in the rejection of melanoma. J. Theor. Biol. **265**(3), 467–480 (2010). https://doi.org/10.1016/j.jtbi.2010.04.030

[22] Eftimie, R.: Investigation into the role of macrophages heterogeneity on solid tumour aggregations. Math. Biosci. **322**, 108325 (2020). https://doi.org/10.1016/j.mbs.2020.108325

[23] Eftimie, R., Barelle, C.: Mathematical investigation of innate immune responses to lung cancer: The role of macrophages with mixed phenotypes. J. Theor. Biol. **524**, 110739 (2021). https://doi.org/10.1016/j.jtbi.2021.110739

[24] Shu, Y., Huang, J., Dong, Y., Takeuchi, Y.: Mathematical modeling and bifurcation analysis of pro- and anti-tumor macrophages. Appl. Math. Model. **88**, 758–773 (2020). https://doi.org/10.1016/j.apm.2020.06.042

[25] Sahoo, P., Yang, X., Abler, D., Maestrini, D., Adhikarla, V., Frankhouser, D., Cho, H., Machuca, V., Wang, D., Barish, M., Gutova, M., Branci-amore, S., Brown, C.E., Rockne, R.C.: Mathematical deconvolution of CAR T-cell proliferation and exhaustion from real-time killing assay data. J. R. Soc. Interface **17**(162), 20190734 (2020). https://doi.org/10.1098/rsif.2019.0734

[26] Barros, L.R.C., de Jesus Rodrigues, B., Almeida, R.C.: Car-t cell goes on a mathematical model. Journal of Cellular Immunology **2**(1), 31–37 (2020). https://doi.org/10.33696/immunology.2.016

[27] Barros, L.R.C., Paixão, E.A., Valli, A.M.P., Naozuka, G.T., Fassoni, A.C., Almeida, R.C.: CARTmath–A mathematical model of CAR-T immunotherapy in preclinical studies of hematological cancers. Cancers **13**(12), 2941 (2021). https://doi.org/10.3390/cancers13122941

[28] León-Triana, O., Pérez-Martínez, A., Ramírez-Orellana, M., Pérez-García, V.M.: Dual-target CAR-Ts with on-and off-tumour activity may override immune suppression in solid cancers: A mathematical proof of concept. Cancers **13**(4), 703 (2021). https://doi.org/10.3390/cancers13040703

[29] León-Triana, O., Sabir, S., Calvo, G.F., Belmonte-Beitia, J., Chulián, S., Martínez-Rubio, Á., Rosa, M., Pérez-Martínez, A., Ramirez-Orellana, M., Pérez-García, V.M.: Car t cell therapy in b-cell acute lymphoblastic leukaemia: Insights from mathematical models. Communications in Nonlinear Science and Numerical Simulation **94**, 105570 (2021). https://doi.org/10.1016/j.cnsns.2020.105570

[30] Santurio, D.S., Paixão, E.A., Barros, L.R., Almeida, R.C., Fassoni, A.C.:
Mechanisms of resistance to car-t cell immunotherapy: Insights from a mathematical model. Applied Mathematical Modelling **125**, 1-15 (2024).
https://doi.org/10.1016/j.apm.2023.08.029

[31] Mostolizadeh, R., Afsharnezhad, Z., Marciniak-Czochra, A.: Mathematical model of chimeric anti-gene receptor (car) t cell therapy with presence of cytokine. Numerical Algebra, Control & Optimization **8**(1) (2018).
https://doi.org/10.3934/naco.2018004

[32] Pérez-García, V.M., León-Triana, O., Rosa, M., Perez-Martinez, A.: CAR T-cells for T-cell leukemias: Insights from mathematical models. Communications in Nonlinear Science and Numerical Simulation **96**, 105684 (2021).
https://doi.org/10.1016/j.cnsns.2020.105684

[33] Sica, A., Larghi, P., Mancino, A., Rubino, L., Porta, C., Totaro, M.G., Rimoldi, M., Biswas, S.K., Allavena, P., Mantovani, A.: Macrophage polarization in tumour progression. Semin. Cancer Biol. **18**(5), 349-355 (2008).
https://doi.org/10.1016/j.semcancer.2008.03.004

[34] Soltantoyeh, T., Akbari, B., Karimi, A., Chalbatani, G.M., Ghahri-Saremi, N., Hadjati, J., Hamblin, M.R., Mirzaei, H.R.: Chimeric antigen receptor (CAR) T cell therapy for metastatic melanoma: Challenges and road ahead. Cells **10**(6), 1450 (2021).
https://doi.org/10.3390/cells10061450

[35] Cheever, A., Townsend, M., O'Neill, K.: Tumor Microenvironment Immunosuppression: A Roadblock to CAR T-cell Advancement in Solid Tumors. Cells **11**(22), 3626 (2022).
https://doi.org/10.3390/cells11223626

[36] Tie, Y., Tang, F., Wei, Y.-q., Wei, X.-w.: Immunosuppressive cells in cancer: mechanisms and potential therapeutic targets. Journal of Hematology & Oncology **15**(1), 61 (2022).
https://doi.org/10.1186/s13045-022-01282-8

[37] Lee, S.H., Starkey, P.M., Gordon, S.: Quantitative analysis of total macrophage content in adult mouse tissues. Immunochemical studies with monoclonal antibody F4/80. J. Exp. Med. **161**(3), 475-489 (1985).
https://doi.org/10.1084/jem.161.3.475

[38] Kuo, W.-T., Chang, J.M., Chen, C.C., Tsao, N., Chang, C.-P.: Autophagy drives plasticity and functional polarization of tumor-associated macrophages. Iubmb Life **74**(2), 157-169 (2022).
https://doi.org/10.1002/iub.2543

[39] Gavagnin, E., Vittadello, S.T., Gunasingh, G., Haass, N.K., Simpson, M.J., Rogers, T., Yates, C.A.: Synchronized oscillations in growing cell
populations are explained by demographic noise. Biophysical Journal **120**(8), 1314–1322 (2021). https://doi.org/10.1016/j.bpj.2021.02.017

[40] Saltelli, A., Tarantola, S., Chan, K.P.-S.: A quantitative model-independent method for global sensitivity analysis of model output. Technometrics **41**(1), 39–56 (1999). https://doi.org/10.1080/00401706.1999.10485594

[41] Saltelli, A., Ratto, M., Andres, T., Campolongo, F., Cariboni, J., Gatelli, D., Saisana, M., Tarantola, S.: Global Sensitivity Analysis: the Primer. John Wiley & Sons, Chichester (2008)

[42] Zhang, X.-Y., Trame, M.N., Lesko, L.J., Schmidt, S.: Sobol sensitivity analysis: a tool to guide the development and evaluation of systems pharmacology models. CPT Pharmacometrics Syst. Pharmacol. **4**(2), 69–79 (2015). https://doi.org/10.1002/psp4.6

[43] Ambahera, S.: Mathematical analysis of ordinary differential equations (ode) and fractional order differential equation (fode) models of CD70 CAR T-cell therapy for gliomas. PhD thesis, Middle Tennessee State University (2022)

[44] Herman, J., Usher, W.: SALib: An open-source python library for sensitivity analysis. J. Open Source Softw. **2**(9), 2 (2017). https://doi.org/10.21105/joss.00097

[45] Iwanaga, T., Usher, W., Herman, J.: Toward SALib 2.0: Advancing the accessibility and interpretability of global sensitivity analyses. Socio-Environ. Sys. Model. **4**, 18155 (2022). https://doi.org/10.18174/sesmo.18155

[46] Qi, T., McGrath, K., Ranganathan, R., Dotti, G., Cao, Y.: Cellular kinetics: A clinical and computational review of CAR-T cell pharmacology. Adv. Drug Deliver. Rev., 114421 (2022). https://doi.org/10.1016/j.addr.2022.114421

[47] McGowan, E., Lin, Q., Ma, G., Yin, H., Chen, S., Lin, Y.: PD-1 disrupted CAR-T cells in the treatment of solid tumors: Promises and challenges. Biomed. Pharmacother. **121**, 109625 (2020). https://doi.org/10.1016/j.biopha.2019.109625

[48] Chow, A., Perica, K., Klebanoff, C.A., Wolchok, J.D.: Clinical implications of t cell exhaustion for cancer immunotherapy. Nat. Rev. Clin. Oncol. **19**, 775–790 (2022). https://doi.org/10.1038/s41571-022-00689-z

[49] Mehrabadi, A.Z., Ranjbar, R., Farzanehpour, M., Shahriary, A., Dorostkar, R., Hamidinejad, M.A., Ghaleh, H.E.G.: Therapeutic potential
of CAR T cell in malignancies: A scoping review. Biomed. Pharmacother. **146**, 112512 (2022). https://doi.org/10.1016/j.biopha.2021.112512

[50] Zhang, Z.-Z., Wang, T., Wang, X.-F., Zhang, Y.-G., Song, S.-X., Ma, C.-G.: Improving the ability of CAR-T cells to hit solid tumors: Challenges and strategies. Pharmacol. Res. **175**, 106036 (2022). https://doi.org/10.1016/j.phrs.2021.106036

[51] Rodriguez-Garcia, A., Lynn, R.C., Poussin, M., Eiva, M.A., Shaw, L.C., O'Connor, R.S., Minutolo, N.G., Casado-Medrano, V., Lopez, G., Matsuyama, T., Jr., D.J.P.: CAR-T cell-mediated depletion of immunosuppressive tumor-associated macrophages promotes endogenous antitumor immunity and augments adoptive immunotherapy. Nat. Commun. **12**(1), 1–17 (2021). https://doi.org/10.1038/s41467-021-20893-2