%% Antes de processar este arquivo LaTeX (LaTeX2e) deve
%% verificar que o arquivo TCAM.cls esta na mesma
%% pasta. O arquivo TCAM.cls pode ser obtido do
%% endereco http://tema.sbmac.org.br.

\documentclass{TCAM}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Atenção: 
% O presente arquivo usa a codificação UTF-8. 
%---------------------
% Para usar Overleaf
%--------------------
% Habilitar a opção:     \usepackage[utf8]{inputenc}     % para acentuação em Português UTF-8  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\usepackage[english]{babel}    % para texto em Inglês
%\usepackage[brazil]{babel}      % para texto em Português

\usepackage[utf8]{inputenc}    % para acentuação em Português UTF-8

\usepackage{epsfig}
\usepackage{multirow}
\usepackage{array}
\usepackage{url}
\usepackage{graphics}
\usepackage{subfigure}
\usepackage{psfrag}
\usepackage{dashrule}


\newcommand{\B}{{\tt\symbol{92}}}
\newcommand{\til}{{\tt\symbol{126}}}
\newcommand{\chap}{{\tt\symbol{94}}}
\newcommand{\agud}{{\tt\symbol{13}}}
\newcommand{\crav}{{\tt\symbol{18}}}

%********************************************************
%     LISTE OS PACKAGES NECESSÁRIOS AO ARTIGO AQUI
%********************************************************
\usepackage{framed}
\usepackage{tikz}
\usepackage{amsmath}
\usepackage{float}
\usepackage{multirow}
\usepackage{ragged2e}
%********************************************************
%                INFORMAÇÕES DO ARTIGO
%********************************************************

\title{A mathematical model to describe the melanoma dynamics under effects of macrophage inhibition and CAR T-cell therapy}

\runningtitle{Melanoma dynamics via a mathematical model}

\author{G. RODRIGUES%
	\aff{Graduate Programme in Biometrics, Institute of Bioscience, UNESP, 18618-689, Botucatu, SP, Brazil -- E-mail: g.rodrigues2@unesp.br -- ORCID: https://orcid.org/0000-0001-7113-8937}, % %
	J. G. SILVA%
	\aff{Federal Institute of Education, Science and Technology of Mato Grosso, Campus Barra do Garças, 78600-000, Barra do Garças, MT, Brazil -- E-mail:  jairo.gomes@ifmt.edu.br -- ORCID: https://orcid.org/0000-0002-5285-5591}, %
	D. S. SANTURIO%
	\aff{National Laboratory of Scientific Computing, LNCC, 25651-075, Petrópolis, RJ, Brazil -- E-mail:  danielassanturio@gmail.com -- ORCID: https://orcid.org/0000-0002-6041-5748}\\
	\&
    P. F. A. MANCERA%
	\aff{Department of Biodiversity and Biostatistics, Institute of Bioscience, UNESP, 18618-689, Botucatu, SP, Brazil -- E-mail:  paulo.mancera@unesp.br -- ORCID: https://orcid.org/0000-0002-2080-8053}}

\abstracttcam{\justifying Melanoma is considered one of the most aggressive types of cancer due to its high capacity to cause metastases, reducing the chances of survival if discovered late. Moreover, melanoma has very immunogenic characteristics that make its treatment difficult, increasing the need to develop more effective techniques of therapy. Applied to oncology, mathematical modeling enables the analysis and distinction of the various mechanisms involved in tumor progression, a fact that allows the analysis of many scenarios, which would be impractical experimentally. The present work gathers as main objectives to develop a mathematical model to describe melanoma dynamics in the presence of tumor-associated macrophages (TAM) and Chimeric Antigen Receptor (CAR) T-cell immunotherapy in order to evaluate why this therapy fails in fighting solid tumors like melanoma and to understand the role of TAM for the failure. We performed stability analysis of the equilibrium points of the model as well as sensitivity analysis of the parameters and some numerical solutions. Our results showed that immunosuppression caused by TAM has a negative impact on the immunotherapeutic dose and varying the cytotoxicity of CAR T-cells together with dose adjustments can help to increase tumour control. In addition, tumour size may also influence treatment efficacy.}

\keywords{Immunotherapy, sensitivity analysis, tumour-associated macrophages.}

\acknowledgments{This study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior – Brasil (CAPES) – Finance Code 001.} 

%********************************************************

\begin{document}
	
	\maketitle
	
	%********************************************************
	\section{Introduction}
	
	One of the most aggressive types of cancer is the cutaneous melanoma. It is estimated for the year 2023 that there will be about 8,980 new cases of cutaneous melanoma in Brazil \cite{inca}. Global projections for new cases of this disease and number of deaths indicate that there will be more than 500,000 new cases and approximately 100,000 deaths by 2040 \cite{Arnold:2022}.	One of the striking features of melanoma is its high resistance to chemotherapy, a fact that compromises the long-term success of treatment, even if it is well targeted. If melanoma is discovered late, the chance of recovery and survival rate decreases dramatically due to its potential to cause metastases \cite{gray-schopfer:2007}. 

    Adoptive immunotherapy with CAR T-cells has been used as a new strategy to improve cancer treatment, especially for haematological malignancies, with promising results \cite{zhao:2018}. In addition, several researches involving the use of this therapy against melanoma are ongoing \cite{simon-cart:2018, soltantoyeh:2021}. In this treatment, the patient its own T-lymphocytes are genetically reprogrammed to recognize antigens expressed on the surface of tumour cells, and upon contact with the malignant cells, the modified T-lymphocyte cells (CAR T-cells) are activated and lead an antitumor immune response \cite{waldman:2020}. There are several challenges for applicability of CAR T-cell immunotherapy against melanoma, and the immunosuppressive tumour microenvironment (TME) is a classical example of challenge to be overcome \cite{soltantoyeh:2021}. 

    Among the mechanisms that induce the rapid progression of melanoma is the recruitment of immune system cells \cite{hallmarks:2000}. Of all the immune system cells recruited to tumour sites, macrophages constitute the majority and collaborate to the pathway of tumour progression through a variety of pathways. Macrophages that infiltrate solid tumours are called as tumour-associated macrophages and participate in tumour progression by secreting growth factors, pro-angiogenic molecules and immunosuppressive factors such as interleukin-10 (IL-10) \cite{mojsilovic:2021, shen:2021, yang:2020, vitale:2019, mantovani:2017}. Some studies indicate that TAM density in the tumour microenvironment correlates with a poor prognosis of melanoma and acts to prevent the cytotoxic action of T-lymphocytes against tumour cells \cite{makitie:2001, chen:2011, vitale:2019}.
    
    There are many papers published involving mathematical models, cancer and immune cells \cite{kuznetsov:1994, eftimie:2010, eftimie:2020investigation, paixao:2022, swanson:2022}. Recently, some papers were published related to mathematical models to CAR T-cells. For example, \cite{eftimie:2017} investigated the possible mechanisms that could explain the elimination of B16 mouse melanoma cells by immune system cells, such as T lymphocytes and M2 macrophages\footnote{Tumour-associated macrophages are also called M2-type macrophages.} considering, among other hypotheses, the properties of immune cells and their role in tumour progression. A mathematical model based on prey-predator dynamics to explore the death kinetics of CAR T-cells in a solid tumour (glioma) was developed by \cite{sahoo:2020} and they observed that the death rate of CAR T-cells has a direct relationship with the proliferation rate and exhaustion, suggesting that lower doses of immunotherapy treatment are more effective, although these cells become more exhausted compared to higher doses of treatment. Mathematical model to investigate the interaction between tumour cells, effector and memory CAR T-cells was present in \cite{barros:2021}. The aim, among others, was to investigate the performance of memory CAR T-cells against a solid tumour over time.
  
    Considering the challenges of using CAR T-cell immunotherapy against cutaneous melanoma, the present work aims to develop a mathematical model based on Ordinary Differential Equations (ODE) to describe the tumour dynamics of advanced melanoma considering its immunosuppressive TME in the presence of TAM and answer the question: Why does CAR T-cell therapy fail to fight melanoma and what role does TAM play in this failure?

    This paper is structured as follows: In Section \ref{sec:model}, we present the proposed mathematical model based on ODE and describe its variables and parameters. Furthermore, we present the equilibrium points and the linear stability analysis. In Section \ref{sec:results} we present the sensitivity analysis and simulations considering some scenarios. Finally, in Section \ref{sec:conclusions}, we end this paper with conclusions.
    
    \section{Mathematical Model}
    \label{sec:model}
    Based on the works by \cite{eftimie:2017}, \cite{sahoo:2020} and \cite{barros:2021}, we formulate an ODE model describing the interactions between melanoma cells, TAM and CAR T-cells as a function of time $t$. Denoting by $T$ the melanoma cell population, by $M$ the TAM population, and by $C$ the number of CAR T-cells. Table \ref{tab1} shows descriptively the biological assumptions embedded in the mathematical expressions.

\begin{table}[H]
\caption{\label{tab1} Mathematical expressions and their meanings.}
\centering
\resizebox{\linewidth}{!}{%
\begin{tabular}{ccl} 
\hline
\multicolumn{1}{c}{\textbf{Equation}} & \textbf{Term}                                 & \multicolumn{1}{c}{\textbf{Description}}                              \\ 
\hline
\\
\multirow{3}{*}{$\dfrac{dT}{dt}$}    & $\alpha_{1}T\left(1-\dfrac{T}{K_{1}}\right)$   & Melanoma cells grow logistically \cite{eftimie:2017, sica:2008}.                         \\
\\
                                     & $\alpha_{1}T(1+\beta_{1}M)$                    & TAM induces melanoma growth \cite{eftimie:2017}.                              \\
\\
                                     & $-\mu_1CT$                                     & Melanoma death due to cytotoxic action of CAR T-cells \cite{barros:2021}.       \\
\\
\multirow{2}{*}{$\dfrac{dM}{dt}$}    & $\alpha_{2}MT\left(1-\dfrac{M}{K_{2}} \right)$ & Logistic growth of TAM due to recruitment by melanoma \cite{eftimie:2017}.  \\
\\
                                     & $-\kappa_{2}M$                                 & Natural death of TAM \cite{eftimie:2017}.                                               \\
\\
\multirow{4}{*}{$\dfrac{dC}{dt}$}    & $\alpha_{3}CT$                                 & Growth of CAR T-cells due to contact with antigen \cite{barros:2021}.                \\
\\
                                     & $-\gamma_{3}CT$                                & Exhaustion of CAR T-cells due to the tumour.                         \\
\\
                                     & $-\theta_{3}MC$                                & Immunosuppression of CAR T-cells caused by TAM.                  \\
\\
                                     & $-\kappa_{3}C$                                 & Natural death of CAR T-cells \cite{barros:2021}.                                     \\
\hline
\end{tabular}
}
\end{table}

The following biological assumptions are incorporated into the model:
\begin{itemize}
    \item[1.] TAM favour tumour progression by secreting growth factors, pro-angiogenic molecules and immunosuppressive factors such as IL-10 \cite{vitale:2019, yamaguchi:2016}.
    
    \item[2.] CAR T-cells kill melanoma cells due to recognition of antigens such as CD126, MCSP, VEGFR2, GD2, among others \cite{soltantoyeh:2021, razavi:2021}.
    
    \item[3.] TAM infiltrate the TME due to cytokines and chemokines produced by the \cite{eftimie:2017} melanoma cells. In addition, they have their population growth described by logistics and compete for space with melanoma cells up to a carrying capacity.
    
    \item[4.] We assume that the expansion of CAR T-cells occurs due to contact with the tumour antigen and, due to immunosuppressive mechanisms present in the TME related to both melanoma cells and TAM, they are inhibited and suffer exhaustion. Moreover, they die naturally. 
\end{itemize}

In particular, the main interest in this modelling proposal is to understand why the immunotherapy with CAR T-cells proves ineffective in the treatment against melanoma and what is the role played by TAM for such ineffectiveness. To do this, we have incorporated the term $-\theta_{3}MC$ to simulate \textit{in silico} the impact of TAM in CAR T-cells population. Furthermore, the term $-\gamma_{3}CT$ models the immunosuppressive TME of melanoma. Then, we propose the following mathematical model:

\begin{equation}
\label{model}
    \left\{\begin{array}{lll}
         \dfrac{dT}{dt}=\alpha_{1}T\left(1-\dfrac{T}{K_{1}}\right)(1+\beta_{1}M)-\mu_1CT,\\
         \\
         \dfrac{dM}{dt}=\alpha_{2}MT\left(1-\dfrac{M}{K_{2}} \right)-\kappa_{2}M,\\
         \\
         \dfrac{dC}{dt}=\alpha_{3}CT-\gamma_{3}CT-\theta_{3}MC-\kappa_{3}C.
    \end{array}\right.
\end{equation}

Table \ref{tab2} shows the parameters with, their respective values and units, and the initial conditions.

\begin{table}[H]
\begin{center}
\begin{minipage}{\textwidth}
\caption{\label{tab2} Descriptions, values and units of the parameters.}
\resizebox{\linewidth}{!}{%
\begin{tabular}{clccc} 
\hline
\multicolumn{1}{c}{\textbf{Parameter}} & \multicolumn{1}{c}{\textbf{Description}}     &  \multicolumn{1}{c}{\textbf{Value}}   & \multicolumn{1}{c}{\textbf{Unit}} & \multicolumn{1}{c}{\textbf{Reference}}  \\ 
\hline
    $\alpha_{1}$    & Natural rate of tumour cells proliferation &  0.69    &     day$^{-1}$       &    \cite{eftimie:2017}                             \\
    $\alpha_{2}$   & Rate of macrophage proliferation promoted by the tumor              &    $10^{-7}$     &  (cells$\cdot$day)$^{-1}$          &      \cite{eftimie:2017}                        \\
    $\alpha_{3}$   & Proliferation rate of CAR T-cells due to contact with tumor antigen        &    $5\cdot10^{-8}$       &     (cells$\cdot$day)$^{-1}$           &   Assumed         \\
    $\beta_{1}$   & Tumor proliferation rate due to pro-tumor action of TAM        &     $2.3\cdot10^{-10}$             &     cells$^{-1}$           &   \cite{eftimie:2017}  \\
    $K_{1}$   & Tumor cell carrying capacity                 &    $10^{9}$                  &     cells       & \cite{eftimie:2017}    \\
    $K_{2}$   & Macrophage carrying capacity                       &    $10^{9}$             &   cells         &    \cite{eftimie:2017} \\
    
    $\mu_{1}$   & Tumor cell death rate due to cytotoxic action of CAR T-cells   &      $[10^{-7}, 10^{-6}]$   &       (cells$\cdot$day)$^{-1}$           &    Assumed        \\
    $\kappa_{2}$   & Macrophage apoptosis rate                        &     0.34      &       day$^{-1}$              & \cite{eftimie:2017}      \\
    $\gamma_3$  &   CAR T-cells immunosuppression rate caused by melanoma cells & $4\cdot10^{-8}$ & (cells$\cdot$day)$^{-1}$ & Assumed \\
    $\theta_3$  &    CAR T-cells immunosuppression rate caused by TAM &[$0, 5\cdot10^{-10}$] & (cells$\cdot$day)$^{-1}$ & Assumed \\
    $\kappa_{3}$   & CAR T-cells apoptosis rate                  &  0.13     &    day$^{-1}$      &  Assumed                                \\
    $T(0)$      & Initial condition of the tumour cells         &   $[2.5\cdot10^{5},5\cdot10^6]$          & cells     &   Assumed\\
    $M(0)$      & Initial condition of the TAM         &   $10^{6}$          & cells     &   \cite{vitale:2019}\\
    $C(0)$      & Initial condition of the CAR T-cells         &   $[2.5\cdot10^{5}, 5\cdot10^6]$          & cells     &   Assumed\\

\hline
\end{tabular}
}
\end{minipage}
\end{center}
\end{table}

\subsection{Equilibrium points and stability}
\label{subsec:stability}

Without loss of generality, let us consider $\alpha_{3}-\gamma_{3}=\beta_{3}>0$. The equilibrium points can be found from the solution of the following system of equations:

\begin{eqnarray}
    \left\{\begin{array}{lll}
        \left(\alpha_{1}\left(1-\dfrac{T^{*}}{K_{1}}\right)(1+\beta_{1}M^{*})-\mu_1C^{*}\right)T^{*}=0,\\
        \\
        \left(\alpha_{2}T^{*}\left(1-\dfrac{M^{*}}{K_{2}} \right)-\kappa_{2}\right)M^{*}=0,\\
        \\
        \left(\beta_{3}T^{*}-\theta_{3}M^{*}-\kappa_{3}\right)C^{*}=0.
    \end{array}\right.
\end{eqnarray}

The equilibrium points and their respective biological meanings are:

\begin{itemize}
    \item $E_0=(0,0,0)$. The trivial equilibrium point, which represents the extinction of all cell populations, which is not biologically relevant.
    
    \item $E_1=(K_1,0,0)$. This equilibrium point represents the failure of immunotherapy treatment and describes the persistence of melanoma cells up to the carrying capacity $K_1$. In addition, disappearance of CAR T-cells and TAM population occurs.
    
    \item $E_2=\left(K_{1},K_{2}\left(1-\dfrac{\kappa_2}{\alpha_2K_1} \right),0\right)$. It exists if, and only if $\alpha_2>\dfrac{\kappa_2}{K_1}$. Like the equilibrium point $E_1$, this also represents the failure of immunotherapeutic treatment, in which the disappearance of CAR T-cells and the persistence of the macrophage population occurs. Furthermore, this equilibrium point intersects $E_1$ when $\alpha_2=\dfrac{\kappa_2}{K_1}$.
    
    \item $E_3=\left(\dfrac{\kappa_3}{\beta_3}, 0, \dfrac{\alpha_1}{\mu_1}\left(1-\dfrac{\kappa_3}{\beta_3K_1} \right)\right)$ which exists if, and only if, $\beta_3>\dfrac{\kappa_3}{K_1}$. This point represents the disappearance of the TAM population and the persistence of melanoma cells and CAR T-cells, and equals $E_1$ when $\beta_3=\dfrac{\kappa_3}{K_1}$.
    \item The state of coexistence satisfies:
    \begin{eqnarray}
    \label{sys}
        \left\{\begin{array}{lll}
            \alpha_{1}\left(1-\dfrac{T^{*}}{K_{1}}\right)(1+\beta_{1}M^{*})-\mu_1C^{*}=0,\\
            \\
            \alpha_{2}T^{*}\left(1-\dfrac{M^{*}}{K_{2}} \right)-\kappa_{2}=0,\\
            \\
            \beta_{3}T^{*}-\theta_{3}M^{*}-\kappa_{3}=0.
        \end{array}\right.
    \end{eqnarray}
    Let $M$ be expressed in terms of $T^{*}$, such that $M^{*}(T^{*})=\dfrac{\beta_3T^{*}-\kappa_3}{\theta_3}$. Then, from the second equation of \eqref{sys} we obtain
    \[T_{1}^*=\dfrac{\alpha_2(K_2\theta_3+\kappa_3)+\sqrt{\delta}}{2\alpha_2\beta_3}\quad\text{and}\quad T_{2}^*=\dfrac{\alpha_2(K_2\theta_3+\kappa_3)-\sqrt{\delta}}{2\alpha_2\beta_3},\]
    where
    \[\delta:=\alpha_2((K_2\theta_3+\kappa_3)^2\alpha_2-4K_2\beta_3\kappa_2\theta_3).\]
    The discriminant will be positive if, and only if,
    \[\beta_3\leq\dfrac{(K_2\theta_3+\kappa_3)^2\alpha_2}{4\theta_3K_2\kappa_2}:=\bar{\delta}.\]
    Then, $T_{1}^{*}\geq0$ since all parameters are positive. Substituting $T_{1}^{*}$ in $M^*(T^*)$, we have
    \[M_{1}^{*}=\dfrac{\alpha_2(K_2\theta_3-\kappa_3)+\sqrt{\delta}}{2\alpha_2\theta_3},\]
    that exists if, and only if, $\beta_3\leq\dfrac{\alpha_2\kappa_3}{\kappa_2}$.
    Therefore, the coexistence equilibrium point $E_{1}^{*}$ is given by:
    \[E_1^*=\left(T_1^*, M_1^*, \dfrac{\alpha_1}{\mu_1}\left(1-\dfrac{T_1^{*}}{K_{1}}\right)(1+\beta_{1}M_1^{*})\right).\]
    Let us now consider the coordinate $T_{2}^{*}$. Substituting $T_{2}^{*}$ in $M^*(T^*)$, we have that
    \[M_2^{*}=\dfrac{\alpha_2(K_2\theta_3-\kappa_3)-\sqrt{\delta}}{2\alpha_2\theta_3},\]
    that exists if, and only if, $\beta_3\geq\dfrac{\alpha_2\kappa_3}{\kappa_2}$. 
    In particular, since the parameters of the model are all positive, it implies that $\beta_3\leq\dfrac{\alpha_2\kappa_3}{\kappa_2}$ and $M_2^{*}<0$. Therefore, $E_2^{*}$ has no biological relevance.  
\end{itemize}

Now we present the local stability of the equilibrium points. Let $\Bar{E}:=(\Bar{T}, \Bar{M}, \Bar{C})$ be any steady-state of the system. We denote by $F_{1}$, $F_{2}$ and $F_{3}$ the three components of the nonlinear function of the system \eqref{model}. Then, the partial derivatives at the equilibrium point $\Bar{E}:=(\Bar{T}, \Bar{M}, \Bar{C})$ are given by
\[\dfrac{\partial F_{1}}{\partial T}=\dfrac{-2(\bar{M}\beta_1+1)(\bar{T}-\frac{K_1}{2})\alpha_1-\mu_1\bar{C}K_1}{K_1},\quad\dfrac{\partial F_{1}}{\partial M}=\dfrac{\alpha_1\bar{T}(K_1-\bar{T})\beta_1}{K_1},\]
\[\dfrac{\partial F_{1}}{\partial C}=-\mu_1\bar{T},\quad\dfrac{\partial F_{2}}{\partial T}=\dfrac{\alpha_2\bar{M}(K_2-\bar{M})}{K_2},\quad\dfrac{\partial F_{2}}{\partial M}=\dfrac{(\alpha_2\bar{T}-\kappa_2)K_2-2\alpha_2\bar{M}\bar{T}}{K_2},\]
\[\dfrac{\partial F_{2}}{\partial C}=0,\quad\dfrac{\partial F_{3}}{\partial T}=\bar{C}\beta_3,\quad\dfrac{\partial F_{3}}{\partial M}=-\bar{C}\theta_3,\quad \dfrac{\partial F_{3}}{\partial C}=-\bar{M}\theta_3+\bar{T}\beta_3-\kappa_3.\]

Thus, analysing the equilibrium point $E_0:=(0,0,0)$, the Jacobian matrix $J$ is given as

\begin{equation}
\label{E0}
    J_0=\begin{pmatrix}
    \alpha_{1} & 0 & 0\\
    \\
    0 & -\kappa_{2} & 0\\
    \\
    0 & 0 & -\kappa_{3}
    \end{pmatrix}.
\end{equation}

The eigenvalues are $\lambda_1=\alpha_1$, $\lambda_2=-\kappa_{2}$ and $\lambda_3=-\kappa_{3}$ and then the equilibrium point $E_0$ is unstable.

Analysing the stability of the equilibrium point $E_{1}:=(K_{1},0,0)$, the Jacobian matrix $J$ is
\begin{equation}
\label{E1}
    J_1=\begin{pmatrix}
    -\alpha_{1} & 0 & -\mu_1K_{1}\\
    \\
    0 & K_{1}\alpha_{2}-\kappa_{2} & 0\\
    \\
    0 & 0 & K_{1}\beta_3-\kappa_{3}
    \end{pmatrix}.
\end{equation}

The eigenvalues associated to the matrix \eqref{E1} are
 \[\lambda_{1}=-\alpha_{1},\quad\lambda_{2}=K_{1}\alpha_{2}-\kappa_{2},\quad\text{and}\quad\lambda_{3}=K_{1}\beta_3-\kappa_{3}.\]

We can conclude that the equilibrium point $E_1$ is linearly asymptotically stable if, and only if,
\[ \alpha_2<\dfrac{\kappa_2}{K_1}\quad\text{and}\quad \beta_3<\dfrac{\kappa_3}{K_1}.\]

Let us now consider the equilibrium point $E_2=\left(K_{1},K_{2}\left(1-\dfrac{\kappa_2}{\alpha_2K_1},0 \right)\right)$, $\alpha_2>\dfrac{\kappa_2}{K_1}$. The Jacobian matrix is given by
\begin{equation}
\label{J2}
    J_2=\begin{pmatrix}
    a_{11} & 0  & a_{13} \\
    \\
    a_{21} & a_{22} & 0\\
    \\
    0 & 0  & a_{33}
    \end{pmatrix},
\end{equation}
where,
\[a_{11}=-\alpha_1\left(K_2\beta_1\left(1-\dfrac{\kappa_2}{K_1\alpha_2}\right)+1 \right),\quad a_{13}=-\mu_1K_1,\quad a_{21}=\dfrac{K_2\kappa_2(K_1\alpha_2-\kappa_2)}{K_{1}^{2}\alpha_2},\]
\[a_{22}=-\alpha_2K_1\left(1-\dfrac{\kappa_2}{K_1\alpha_2}\right),\quad a_{33}=-K_2\theta_3\left(1-\dfrac{\kappa_2}{K_1\alpha_2} \right)+K_1\beta_3-\kappa_3.\]

The characteristic equation associated to \eqref{J2} is given by
\begin{equation}
\label{ecj2}
    (\lambda-a_{11})(\lambda-a_{22})(\lambda-a_{33})=0.
\end{equation}
According to the terms of \eqref{ecj2}, the eigenvalues are given by $\lambda_1:=a_{11}$, $\lambda_2:=a_{22}$ and $\lambda_3:=a_{33}$. We know that $\alpha_2>\dfrac{\kappa_2}{K_1}$ is the existence condition of the point $E_3$, implying that $\lambda_1<0$ and $\lambda_2<0$. Considering the third term of \eqref{ecj2}, we have that $\lambda_{3}:=a_{33}<0$ if, and only if, $\beta_3<\dfrac{\kappa_3}{K_1}$.

Considering the equilibrium point $E_3=\left(\dfrac{\kappa_3}{\beta_3}, 0, \dfrac{\alpha_1}{\mu_1}\left(1-\dfrac{\kappa_3}{\beta_3K_1} \right)\right)$, $\beta_3>\dfrac{\kappa_3}{K_1}$, then the Jacobian matrix is
\begin{equation}
\label{J3}
    J_3=\begin{pmatrix}
    b_{11} & b_{12}  & b_{13} \\
    \\
    0 & b_{22} & 0\\
    \\
    b_{31} & b_{32}  & 0
    \end{pmatrix},
\end{equation}
where,
\[b_{11}=-\dfrac{\alpha_1\kappa_3}{\beta_3 K_1},\quad b_{12}=\dfrac{\alpha_1\kappa_3\beta_1(K_1\beta_3-\kappa_3)}{\beta_{3}^{2}K_1},\quad b_{13}=-\dfrac{\mu_1\kappa_3}{\beta_3},\]
\[b_{22}=\dfrac{\alpha_2\kappa_3}{\beta_3}-\kappa_2,\quad b_{31}=\dfrac{\alpha_1(K_1\beta_3-\kappa_3)}{\mu_1K_1},\quad b_{32}=-\dfrac{\alpha_1\theta_3(K_1\beta_3-\kappa_3)}{\mu_1K_1\beta_3}.\]

The characteristic equation associated to the matrix \eqref{J3} is
\begin{equation}
\label{ecJ3}
    (-b_{22}+\lambda)(\lambda^{2}-\lambda b_{11}-b_{13}b_{31})=0.
\end{equation}

From the first term of \eqref{ecJ3}, we have that $\lambda_1=b_{22}$. It is negative if, and only if, $\beta_3>\dfrac{\alpha_2\kappa_3}{\kappa_2}.$ For the second term of \eqref{ecJ3}, we use the Routh-Hurwitz criterion:
\begin{itemize}
    \item $p_1:=-b_{11}>0$, which is always satisfied.
    \item $p_2:=-b_{13}b_{31}>0$ if, and only if, $\beta_3>\dfrac{\kappa_3}{K_1}$.
\end{itemize}

Therefore, if the above conditions are satisfied, then $E_3$ is linearly asymptotically stable.

Let us analyse the coexistence equilibrium point
\[E_1^*=\left(T_1^*, M_1^*, \dfrac{\alpha_1}{\mu_1}\left(1-\dfrac{T_1^{*}}{K_{1}}\right)(1+\beta_{1}M_1^{*})\right),\]
where $T_1^{*}:=\dfrac{\alpha_2(K_2\theta_3+\kappa_3)+\sqrt{\delta}}{2\alpha_2\beta_3}$, and $M_1^{*}:=\dfrac{\alpha_2(K_2\theta_3-\kappa_3)+\sqrt{\delta}}{2\alpha_2\theta_3}$, $\beta_3\leq\dfrac{\alpha_2\kappa_3}{\kappa_2}.$

The Jacobian matrix is expressed by
\begin{equation}
\label{JE1*}
    J^*=\begin{pmatrix}
    c_{11} & c_{12}  & c_{13} \\
    \\
    c_{21} & c_{22} & 0\\
    \\
    c_{31} & c_{32}  & 0
    \end{pmatrix},
\end{equation}
where

\begin{eqnarray*}
    c_{11}&=&-\dfrac{\alpha_1(K_2\alpha_2\theta_3+\alpha_2\kappa_3+\sqrt{\delta})(K_2\alpha_2\beta_1\theta_3-\alpha_2\beta_1\kappa_3+\sqrt{\delta}\beta_1+2\alpha_2\theta_3)}{4\alpha_2^{2}\beta_3 K_1\theta_3},\\
    c_{12}&=&-\dfrac{\alpha_1(K_2\alpha_2\theta_3+\alpha_2\kappa_3+\sqrt{\delta})(-2\alpha_2\beta_3K_1+K_2\alpha_2\theta_3+\alpha_2\kappa_3+\sqrt{\delta})\beta_1}{4\alpha_2^{2}{\beta_3}^{2}K_1},\\
    c_{13}&=&-\dfrac{\mu_1(K_2\alpha_2\theta_3+\alpha_2\kappa_3+\sqrt{\delta})}{2\alpha_2\beta_3},\\
    c_{21}&=&-\dfrac{(K_2\alpha_2\theta_3-\alpha_2\kappa_3+\sqrt{\delta})(-K_2\alpha_2\theta_3-\alpha_2\kappa_3+\sqrt{\delta})}{4\theta_3^{2}K_2\alpha_2},\\
    c_{22}&=&\dfrac{K_2\alpha_2^{2}\kappa_3\theta_3-2\kappa_2\alpha_2\theta_3\beta_{3}K_2-\sqrt{\delta}K_2\alpha_2\theta_3+\alpha_2^2\kappa_3^2-\delta}{2\alpha_2\theta_3\beta_3K_2},\\
    c_{31}&=&-\dfrac{(\sqrt{\delta}\beta_1+\alpha_2(K_2\beta_1\theta_3-\beta_1\kappa_3+2\theta_3))\alpha_1(\sqrt{\delta}+\alpha_2(-2K_1\beta_3+K_2\theta_3+\kappa_3))}{4\mu_1\alpha_2^2K_1\theta_3},\\
    c_{32}&=&\dfrac{(\sqrt{\delta}\beta_1+\alpha_2(K_2\beta_1\theta_3-\beta_1\kappa_3+2\theta_3))\alpha_1(\sqrt{\delta}+\alpha_2(-2K_1\beta_3+K_2\theta_3+\kappa_3))}{4\mu_1\alpha_2^2K_1\beta_3}.
\end{eqnarray*}

Because of the algebraic complexity characterized by the equilibrium point $E_1^{*}$ evaluated in the Jacobian matrix, the stability will be studied numerically using Python and the values for the parameters present in Table \ref{tab2}, where we assumed $\mu_1=10^{-7}$ (cells$\cdot$day)$^{-1}$. For different values of $\theta_3$ we obtain the following equilibrium point and its eigenvalues associated:

\begin{itemize}
    \item $\theta_3=5\cdot10^{-11}$ (cells$\cdot$day)$^{-1}$ $\Rightarrow E_1^{*}=(1.7\cdot10^{7}, 8.0\cdot10^8, 8.03\cdot10^6)$. The eigenvalues are $\lambda_1=-1.383$, $\lambda_2=0.0048+0.355i$ and $\lambda_3=0.0048-0.355i$. Thus, the point is unstable.
    \item $\theta_3=10^{-10}$ (cells$\cdot$day)$^{-1}$ $ \Rightarrow E_1^{*}=(2.14\cdot10^{7}, 8.41\cdot10^8, 8.05\cdot10^6)$. The eigenvalues are $\lambda_1=-1.8179$, $\lambda_2=-0.0004+0.397i$ and $\lambda_3=-0.0004-0.397i$. Thus, the point is a stable spiral.
    \item $\theta_3=5\cdot10^{-10}$ (cells$\cdot$day)$^{-1}$ $ \Rightarrow E_1^{*}=(6.017\cdot10^{7}, 9.434\cdot10^8, 7.892\cdot10^6)$. The eigenvalues are $\lambda_1=-5.681$, $\lambda_2=-0.023+0.672i$ and $\lambda_3=-0.023-0.672i$. Thus, the point is a stable spiral. 
\end{itemize}    

\section{Results and Discussions}
\label{sec:results}

\subsection{Sensitivity Analysis}
\label{subsec:sensitivity}

    The dynamics underlying CAR T-cell therapy are recognized to be the result of an interplay of many phenomena that occur at different biological populations; however, some indirect effects, such as cytokine production, the interaction of endogenous immune cells (macrophages, NK cells and so on) are not well understood, and many mechanisms remain unanswered. Under this perspective, defining the focus of the investigation, as well as acquiring and integrating data into the model,  are sources of uncertainty but are critical steps toward developing a predictive mathematical model. As a result, how uncertainties are quantified and propagated in the model during  the prediction process determines a model's predictability.
    
    The ODE system presented in Equations \eqref{model} has 10 parameters that have been considered from previous research. In models with several uncertain inputs, it is important to assess which parameters have a significant impact on the outputs. To this end, we perform a sensitivity analysis (SA) using the elementary effects (EE) technique, by measuring the effect of small parameter changes (a perturbation of $\pm$ 25\% from the values established in Table \ref{tab2}) on different quantities of interest (QoI) of the model using an eighty-day simulation of the system's evolution. 
    
     On the EE method, the selected output of interest can be represented by a function $Y(\textbf{X})$, where $\textbf{X} = (X1, . . . , X_d)$ is a vector of $d$ independent input variables (parameters) defined within the range of a continuous interval. The elementary effect of the $i$-th input factor for a given value of $\textbf{X}$ is defined as:
\begin{equation}
    \label{EEf1}
    EE_i = \dfrac{Y(X_1,...,X_i+\Delta,...,X_d) - Y(\textbf{X})}{\Delta} = \dfrac{Y(\textbf{X} + \textbf{e}_i\Delta) - Y(\textbf{X})}{\Delta},
\end{equation}
where $\Delta$ is a predetermined perturbation factor of $X_i$ and $Y(\textbf{X} + \textbf{e}_i\Delta)$ is the model output corresponding to a $\Delta$ change in $X_i$. The distribution of the elementary effects ($\rho_i$) for each input parameter is obtained by constructing $r$-trajectories in a grid of $d$-dimensions and $p$-levels, where each input is randomly sampled using a one-at-a-time (OAT) sampling strategy. %Morris \cite{Morris:1991} proposed this economical design, which generates a cost of $r(d+1)$ model evaluations.

In order to rank the parameters, we chose to estimate the mean $\mu_i$ of the distribution ($\rho_i$). This measure evaluates the general impact of the $i$-th parameter on the output (QoI) and the average of absolute elementary effects is:
\begin{equation}
    \label{EEf12}
    {\mu_i}^*= \dfrac{1}{r} \sum\limits_{j=1}^{r}|{EE_i}^j|.
\end{equation}

The results of this parameter sensitivity analysis for the mathematical model are shown in Figure \ref{sens}. In the case where the QoI is the tumour load, Figure \ref{sens}(a) shows that the system is found to be most sensitive to the balance between the CAR-T proliferation and inhibition rates ($\alpha_3-\gamma_3$), as well as to the tumour growth parameter and its carrying capacity $K_1$. This suggests that, in addition to the tumour's aggressiveness, as represented by growth parameter $\alpha_1$, even minor differences in the absolute number of CAR-T cells in the system during the dynamics can affect clinical outcome. It is possible to notice that, during tumour escape, all parameters associated with changes in CAR-T cell concentration appear to be relevant, including $\theta_3$ and $\kappa_3$. This would indicate that any treatment which might enhance CAR T-cell expansion should be pursued. By contrast, the size of the tumour after eighty days is not very sensitive to the TAM participation in tumour progression expressed by $\beta_1$. According to this model, then, the pro-tumour action of TAM alone is not a determining factor in the eventual size of the tumour, and the cytolytic activity should be considered as the conjunction of number of CAR T-cells binding with tumour cells as well as its killing efficacy.

\begin{figure}
    \centering
    \subfigure[]{\includegraphics[scale=0.30]{SENS ADAPTED-TUMOR.png}}
    \subfigure[]{\includegraphics[scale=0.30]{SENS ADAPTED-MAC.png}}
    \subfigure[]{\includegraphics[scale=0.30]{SENS ADAPTED-CAR.png}}
    \caption{\label{sens}{Sensitivity analysis was performed using the EE method for different quantities of interest (QoI) with $r = 40$ and $p = 4$. The average of absolute elementary effects is shown in (a) Tumour/Melanoma Cells, (b) TAM population and (c) CAR T-cells. The parameters from Table \ref{tab2} and $\theta_3=10^{-10}$ (cells$\cdot$day)$^{-1}$ were perturbed by $\pm$ 25\% and the system's evolution was simulated for eighty days.}}
\end{figure}

In the scenario where the output is the concentration of macrophages, Figure \ref{sens}(b) shows that during an initial expansion phase, this population is highly sensitive to the TAM’s intrinsic characteristics ($\alpha_2$ and $\kappa_2$) together with tumour proliferation ($\alpha_1$). Then, as the macrophage population is being depleted and similar to what was described previously for the tumour, the system is also sensitive to the balance between the CAR T-cell proliferation and inhibition rates ($\alpha_3-\gamma_3$). This reveals that, while CAR T-cells do not directly interact with TAM, their action in reducing tumour population will have an indirect effect on the interplay between tumour cells and macrophages. Matter of fact, during tumour progression, TAM are again recruited to the tumour site, and CAR T-cells attempt to control these populations at first, but then these populations escape. %The sensitivity analysis shows that between 40-60 days, $\kappa_3$ and $\theta_3$ play a significant role in macrophage re-expansion; however, once this population has reached its carrying capacity, $K_2$ is the most influential parameter.

At last, when QoI is the CAR T-cell density, during an initial phase, the model is most sensitive to the CAR T expansion strength, given by antigen recognition ($\alpha_3-\gamma_3$) and its natural apoptotic rate $\kappa_3$. After the peak, during the contraction phase, the inhibition of CAR T-cells by macrophages ($\theta_3$), and the absolute number of TAM in the microenvironment might affect the residual number of CAR T in the system. This finding, indicates that investigating CAR T-cell interactions with the host immune system could shed light on the mechanisms underlying (non)durable responses. The sensitivity analysis shows that, during the final phase, antigen binding and CAR T proliferation alone are not determining factors in CAR T-cell persistence, and its efficiency in killing tumour cells must also be considered. %, since tumour control may imply coexistence. However, care must be taken on interpreting this state, once the model does not yet account memory formation and the CAR-T exhaustion.
\subsection{Numerical Simulations}
\label{subsec:simulations}

The simulations were performed by implementing the fourth order Runge-Kutta method adapted and applied to the \eqref{model} system in Python programming language. Several models involving CAR T-cell kinetics were developed to show the distinct phases that comprise it, being the phases of distribution, expansion, contraction and persistence \cite{qi:2022}. As shown in the Subsection \ref{subsec:stability}, some equilibrium points are stable spirals, indicating oscillatory dynamics in the model solutions, including in the CAR T-cell dynamics, which is not biologically possible. To adapt the simulations for a realistic scenario, such as when CAR T-cells are exhausted or the tumour loses antigen, we used a numerical procedure that the CAR T-cells population will not re-expand even if the tumour escapes. Furthermore, the algorithm was projected to hold the period of time that the tumour population $T$ remains below $10^2$ cells \cite{barros:2021}.

\subsubsection{TAM inhibition affects tumour control}
The immunosuppressive TME represents a major barrier to effective tumour-specific T-cell responses to cancer \cite{rodriguez:2021}. As TAM are major constituents of the TME, and are comprised of a vastly heterogeneous population, we first studied the impact caused by an increasingly immunosuppressive TME in a patient under CAR T-cell treatment. To do so, we performed several simulations varying the parameter $\theta_3$ and the results are shown in Figure \ref{sim1}.    

\begin{figure}[H]
    \centering
    $\theta_3 = 0$ (cells$\cdot$day)$^{-1}$  \hspace{0.10cm} \textcolor{black}{\tikz[baseline]{\draw[thick] (0,.5ex)--++(.5,0) ;}}\hspace{1.50cm}
    $\theta_3$ = 5$\cdot$10$^{-11}$ (cells$\cdot$day)$^{-1}$ \hspace{0.10cm} \textcolor{black}{\tikz[baseline]{\draw[dashed] (0,.5ex)--++(.5,0) ;}}
    
    $\theta_3$ = 10$^{-10}$ (cells$\cdot$day)$^{-1}$ \hspace{0.10cm} \textcolor{black}{\tikz[baseline]{\draw[dotted] (0,.5ex)--++(.5,0) ;}}\hspace{1.50cm}
    $\theta_3$ = 5$\cdot$10$^{-10}$ (cells$\cdot$day)$^{-1}$ \hspace{0.10cm} \textcolor{black}{\tikz[baseline]{\draw[dash dot] (0,.5ex)--++(.5,0) ;}}

    \subfigure[]{\includegraphics[scale=0.27]{TUM2 - theta3.png}}
    \subfigure[]{\includegraphics[scale=0.27]{MAC2 - theta3.png}}
    \subfigure[]{\includegraphics[scale=0.27]{CAR - theta3.png}}
    \caption{\label{sim1}{An immunosuppressive TME rules tumour control \textit{in-silico}. Dynamics of the number of tumour cells (blue curves), TAM cells (black curves) and CAR T-cells (red curves) ruled by Equation \eqref{model} in different immunosuppression scenarios. Initial data used in the simulations $T(0)=2\cdot10^6$ cells, $M(0)=10^6$ cells and $C(0)=5\cdot10^5$ cells and parameter values presented in Table \ref{tab2}. The grey shaded area represent the non-detection limit of $10^2$ cells.}}
\end{figure}

TAM-induced immunosuppression has a negative impact on treatment, as illustrated in Figure \ref{sim1}(a). It is possible to see that in a highly immunosuppressive TME, the tumour is not controlled at all; however, in an environment where the immunosuppression is gradually decreasing, $\theta_3=5\cdot10^{-11}$, $\theta_3=10^{-10}$, $\theta_3=5\cdot10^{-10}$ (cells$\cdot$day)$^{-1}$, the tumour is controlled during 6, 17, and 25 days, respectively. The simulations agree with clinical observations that TAM's action correlates with poor survival in many types of solid cancer. These observations also suggests that therapeutics strategies targeting  TAM subsets should be explored, even if it means a longer window of opportunity for the application of other therapies.

Figure \ref{sim1}(b) shows that TAM dynamics occurs similarly as in the melanoma cells. This scenario is expected considering that TAM are recruited to the tumour site, promoting its progression. Therefore, depletion of tumour cells by the action of immunotherapy also impacts the TAM population. This means that even if TAM do not express the target-antigen of CAR T-cells, the treatment will have an indirect effect on this population dynamics.

From Figure \ref{sim1}(c), we can observe that macrophages play a pivotal role in the development of effective immunotherapy, the simulations show that increasing immunosuppression caused by TAM causes the CAR T-cell population to decrease more quickly, affecting treatment. Moreover, from the parameter set used, the model was able to capture the phases of expansion, contraction and persistence, indicating that the functional response of the model is adequate to describe the dynamics of CAR T-cells \cite{qi:2022}. Furthermore, our results corroborate data from the literature about the immunosuppression caused by TAM when CAR T-cell therapies are used \cite{mehrabadi:2022, zhang:2022}. Immunosuppressive factors such as IL-10 produced by TAM are among the possible causes for the failure of CAR T-cell immunotherapy to combat solid tumours such as melanoma \cite{rodriguez:2021}.

\subsubsection{Killing efficiency affects CAR-T dynamics}
Next, we investigated the effect of the initial number of CAR T-cells injected on the system's dynamics in the case where CAR T-cells have different killing efficiency when targeting the tumour antigen. To do so, we performed simulations of equations \eqref{model} and we found a dependence of the dynamics on the number of injected CAR T-cells. Results shown in Figure \ref{sim2} present some examples for numbers of cells initially injected, ranging from $2.5\cdot10^{5}$ to $5\cdot10^{6}$ cells.
\begin{figure}[H] 
    \centering
    $C(0) = 2.5\cdot10^{5}$ cells  \hspace{0.10cm} \textcolor{black}{\tikz[baseline]{\draw[thick] (0,.5ex)--++(.5,0) ;}}\hspace{0.7cm}
    $C(0) = 5\cdot10^{5}$ cells \hspace{0.10cm} \textcolor{black}{\tikz[baseline]{\draw[dashed] (0,.5ex)--++(.5,0) ;}}   
    
    $C(0) = 2.5\cdot10^{6}$ cells \textcolor{black}{\tikz[baseline]{\draw[dash dot] (0,.5ex)--++(.5,0) ;}}\hspace{0.7cm}
    $C(0) = 5\cdot10^{6}$ cells \hspace{0.10cm} \textcolor{black}{\tikz[baseline]{\draw[dotted] (0,.5ex)--++(.5,0) ;}}
    
\hspace{-1cm}\rotatebox[origin=l]{90}{\small{$\mu_1 = 10^{-6}$ (cells$\cdot$day)$^{-1}$}}
    \subfigure[]{\includegraphics[scale=0.27]{TUM2 - CICAR 1e-6.png}}
    \subfigure[]{\includegraphics[scale=0.27]{MAC2 - CICAR 1e-6.png}}
    \subfigure[]{\includegraphics[scale=0.27]{CAR - CICAR 1e-6.png}} 
    
 \hspace{-1cm}\rotatebox[origin=l]{90}{\small{$\mu_1 = 10^{-7}$ (cells$\cdot$day)$^{-1}$}}
    \subfigure[]{\includegraphics[scale=0.27]{TUM2 - CICAR 1e-7.png}}
    \subfigure[]{\includegraphics[scale=0.27]{MAC2 - CICAR 1e-7.png}}
    \subfigure[]{\includegraphics[scale=0.27]{CAR - CICAR 1e-7.png}}
    \caption{\label{sim2}{Killing efficiency and CAR T dose change the system dynamics. System evolution  of the total number of tumour cells (blue curves), TAM cells (black curves) and CAR T-cells (red curves). The curves correspond to different values of CAR T-cells injected in a patient with initial tumour burden of $T(0)=2\cdot10^6$ cells, $M(0)=10^6$ cells. The parameter values used in the simulations are in Table \ref{tab2} and $\theta_3=5\cdot10^{-11}$ (cells$\cdot$day)$^{-1}$.}}
\end{figure}

The top panel of Figure \ref{sim2} shows there were two different types of dynamics in the CAR T population depending on the initial choice of $C(0)$, see Figure \ref{sim2}(c). For a highly cytotoxic CAR T-cell ($\mu_1 = 10^{-6}$ (cell$\cdot$day)$^{-1}$) and the parameters listed in the Table \ref{tab2}, there was a qualitative change of the dynamics around $C(0) = 2.5\cdot10^{6}$ cells. Thus, small doses of CAR T-cells led to an early tumour progression followed by a small reduction in the tumour load (Figure \ref{sim2}(a)), while for doses larger than the mentioned threshold, the therapy was able to control tumour growth for a longer period, a behaviour also demonstrated by TAM (Figure \ref{sim2}(b)). Nevertheless, the threshold is related to this particular set parameter and would change under different conditions, which may provide a durable response.

The bottom panel of Figure \ref{sim2} shows a different dynamics when the cytotoxic activity against tumour is low ($\mu_1 = 10^{-7}$ (cell$\cdot$day)$^{-1}$). The change in the magnitude of the initial CAR T-cell dose, had a significant effect on the maximum expansion achieved (Figure \ref{sim2}(c)). We simulated \textit{in silico}, and saw that a reduction in the efficiency of stimulation in the CAR T-cells led to a slower growth of this population. Since CAR T expansion relies on the contact with tumour cells, Figure \ref{sim2}(a) shows that a higher CAR T infusion, quickly decrease the tumour population which impairs CAR T levels on the system.

\subsubsection{Initial tumour burden and CAR T interaction with antigen changes the overall dynamics}
The initial tumour load on has a dual role: a high tumour load may favour the initial expansion of CAR T-cells on one hand, but it may enhance tumour immune suppression capabilities on the other. As a result, the question of how to use CAR T-cell treatments in combination with other neoadjuvant therapies arises. To shed some light on the question, we computationally investigated the idea of using CAR T-cells after a partial surgical resection or a primary regimen of radiotherapy, which is a common scenario in the context of melanoma.
\begin{figure}[H] 
    \centering
    $T(0) = 2.5\cdot10^{5}$ cells  \hspace{0.10cm} \textcolor{black}{\tikz[baseline]{\draw[thick] (0,.5ex)--++(.5,0) ;}}\hspace{0.7cm}
    $T(0) = 5\cdot10^{5}$ cells \hspace{0.10cm} \textcolor{black}{\tikz[baseline]{\draw[dashed] (0,.5ex)--++(.5,0) ;}}   
    
    $T(0) = 2.5\cdot10^{6}$ cells \textcolor{black}{\tikz[baseline]{\draw[dash dot] (0,.5ex)--++(.5,0) ;}}\hspace{0.7cm}
    $T(0) = 5\cdot10^{6}$ cells \hspace{0.10cm} \textcolor{black}{\tikz[baseline]{\draw[dotted] (0,.5ex)--++(.5,0) ;}}
    
    \subfigure[]{\includegraphics[scale=0.29]{TUM2 - CITUMOR 1e-8.png}}
    \subfigure[]{\includegraphics[scale=0.29]{MAC2 - CITUMOR 1e-8.png}}
    \subfigure[]{\includegraphics[scale=0.29]{CAR - CITUMOR 1e-8.png}}

    \caption{\label{sim3}{\textit{In silico} simulations under different initial tumour loads. System evolution  of the total number of tumour cells (blue curves), TAM cells (black curves) and CAR T-cells (red curves). The curves correspond to different values of tumour burden in a patient with initial CAR T injection of $C(0)=5\cdot10^5$ cells, $M(0)=10^6$ cells. The parameter values used in the simulations are in Table \ref{tab2} and $\theta_3=5\cdot10^{-11}$ (cells$\cdot$day)$^{-1}$.}}
\end{figure}

In Figure \ref{sim3} we notice that during the first weeks of treatment, the reduction in initial tumour load resulted in an initial tumour progression, which delays the expansion phase of CAR T-cells and gradually recruits TAM to the tumour site. A higher tumour burden, on the other hand, shortens the duration of the CAR T distribution phase, keeping more CAR T-cells in the system and slowing the initial tumour proliferation. Since a lower number of tumour cells implies a lower density of CAR T-cells, the CAR-peak levels in this case scenario are smaller. However, following the CAR T contraction phase, the greater the tumour burden, the greater the concentration of CAR T-cells during the persistence phase.

It is important to notice that this behaviour is most probably due to the tumour aggressiveness and to the TAM support to tumour growth. These findings suggest that a new CAR structure, designed to reprogram the TME, might indirectly limit tumour progression. However, when making broad statements about specific responses to treatments, caution is necessary, and any quantitative data must be interpreted as one possibility, with further investigation required for a different set of parameters.

\section{Conclusions}
\label{sec:conclusions}
From a proposed ODE model, it was possible to capture the expansion, contraction and persistence phases of CAR T-cells and their cytotoxic effects against cutaneous melanoma. The sensitivity analysis of the parameters showed that the pro-tumour action of TAM alone does not represent a determining factor for tumour growth. Furthermore, sensitivity analysis revealed that during the contraction phase of CAR T-cells the absolute number of TAM can affect the residual number of CAR T-cells in the system, revealing that interactions between CAR T-cells and immune system cells could shed light on the mechanisms that lead treatment to fail. 

Analytical results and numerical simulations revealed the challenges for the applicability of CAR T-cell immunotherapy against melanoma. In the first scenario we investigated the inhibitory effects of TAM on CAR T-cells and the results showed that increased suppression inherent to TAM has a negative impact on treatment and, furthermore, our results corroborate clinical observations about poor survival of CAR T-cells in fighting several types of solid tumours, including melanoma. In the second scenario, we investigated the cytotoxic efficiency of CAR T-cells against melanoma and obtained as results that increasing the cytotoxicity of CAR T-cells in conjunction with initial immunotherapeutic dose adjustments may contribute to greater tumour control. Finally, in the third scenario we varied the initial tumour burden and found that this variation affects CAR T-cell dynamics in general, suggesting that the development of new CAR structures with the ability to reprogram the TME may be a strategy to limit tumour progression and increase immunotherapeutic efficacy. Indeed, the reasons why CAR T-cell therapy fails in fighting melanoma are many and our \textit{in silico} analyses evidenced that TAM has an important role for this failure, opening the way for further investigations and improvements in therapeutic techniques.


	\bibliographystyle{ieeetr}
	\bibliography{TCAM_modelo}
	
\end{document}
