Almost all tissues contain a population of fibroblasts that provide the tissue architecture, and serve as sentinels for tissue dysfunctions. When a solid tumor grows, quiescent or progenitor lung fibroblasts activate an initial wound-healing response with increased production of ECM and growth factor (cytokine) and upregulation of activation markers such as FAP (Salmon et al., 2019). Due to mechanotransduction and/or biochemical signals from tumor or immune cells, CAFs often increase their level of contractility, which affects the maintenance of the tumor stroma (Dvorak, 1986). CAFs are characterized by their increased mobility, proliferation, and ECM remodeling and, unlike wound-associated fibroblasts, they seem to undergo poorly reversible activation in absence of an appropriate therapy (Nakamura et al., 2004; Wang et al., 2009; Peng et al., 2013; Kalluri, 2016). The diversity of this population, resulting from phenotypic modifications, explains their diverse functions and localization in the stroma, close or distant from the tumor nest (Barrett and Puré, 2020; Grout et al., 2022).

Activated fibroblasts produce fibers that act as a mechanical barrier around the tumor, impeding the movement of immune cells and limiting the interaction between cytotoxic T-cells and cancer cells (Salmon et al., 2012; Yaegashi et al., 2021). However, by creating such a barrier around tumors, CAFs may also prevent the spread of cancer, as mechanical stress can reduce cell spreading and promote cell apoptosis (Cheng et al., 2009; Barbazan et al., 2023). In addition, biochemical factors expressed by CAFs also help to modify the phenotype of T-cells or inactivate their cytotoxic capacity (Nazareth et al., 2007; Lakins et al., 2018; Koppensteiner et al., 2022). Thus, CAFs have an ambiguous role as tumor promoters, inhibiting T-cell invasion into the nest, and as tumor suppressors, limiting the cancer growth and giving rise to an immune-excluded tumor (Figure 2B), which is less proliferative than a free-growing tumor (Figure 2A). Failure to confine the tumor induces high levels of cytotoxic T-cell infiltration, creating a hot tumor (Figure 2C). Therefore, one aim of this study is to predict susceptibility of a tumor to T-cell infiltration according to the scheme given in Figure 2.

Figure 2

Download asset Open asset

In the context of non-small cell lung carcinoma (NSCLC), the recruitment of CD8 ⁺ T-cells seems to be modulated by a specific tumor-associated antigen present on the surface of the cancer cells (Hiraoka et al., 2006). In particular, in the family of inflammatory proteins, also called cytokines, Interlukin-6 (IL6) and (IL8) seem, among others to stimulate the infiltration of CD8⁺ (Liu et al., 2022; McKeown et al., 2004; Brenner et al., 2017). The origins of CD8⁺ T-cells that infiltrate the tumor are diverse, as they differentiate from both circulating and tissue-resident precursors (Gueguen et al., 2021). However, infiltration alone does not imply that the immune response is efficient. In fact, an immunotherapy such as anti-PD1/PDL1 antibodies is often needed to boost the response of these T-cells. Chemical factors also play an important role in attracting T-cells via chemotaxis, such as the chemokines produced by dendritic cells or cancer cells (Brown et al., 2007). The absence of such chemicals leads to the formation of the so-called immune-desert tumor (Figure 2A).

For a quantitative modeling, it is critical to use biologically well-identified data on both tumor and stroma, in terms of the density and the localization of the different species. Therefore, we survey the literature to build a quantitative view for the composition of the TME. We recapitulate these values in Table 1. However, human data show a great diversity depending on many factors such as tumor edge, patient age, non-cancer health status, as well as method of analysis. The fibroblast population in the stroma, which is very heterogeneous and not easy to identify, is the best example. In particular, although fibroblasts are often considered to be a major component of the stroma, they are only found in small proportions in scRNA-seq with 10 X single-cell systems, which may be rather surprising. This situation is well known and may be due to tissue digestion processes, to the fact that part of the stroma is not extracted with the tumor nodule, and to a lower efficiency of 10 X single-cell systems for fibroblasts. Another caveat is that identifying the surface fraction of cells by staining can be very different from counting individual cells. In fact, cell types can vary greatly in size. For example, lung cancer cells have a diameter between 13–18 μm (Hosokawa et al., 2013), fibroblasts ∼16 μm (Mitsui and Schneider, 1976), and T-cells between 5–10 μm (Tasnim et al., 2018), resulting in volumes and projected areas that can be 10 times smaller for T-cells compared to fibroblasts and cancer cells for the same densities. Moreover, in some studies, authors focus on well-characterized zones of enrichment for the different species (tumor nest, fibrotic areas), and do not consider each cell type per se (Sieren et al., 2010; Sieren et al., 2011). Despite these remarks, cancer cells, T-cells, and α TME. The limitations of using data coming from the literature led us to perform our own measurements. In this study, we analyzed data from 13 patients with lung squamous cell carcinoma (LUSC) and 22 patients with lung adenocarcinoma (LUAD), based on a recent publication by some of the present authors (Grout et al., 2022). In this previous study, different fibroblast types were identified using multiplex IHC imaging. Tumor islets were stained with keratin, T-cells with CD3, and fibroblasts with αSMA, FAP, and ADH1B. The coverage was evaluated for each fibroblast type and the total fibroblast coverage corresponds to the sum of these different coverages. Here, in the case of LUSC, we found that fibroblasts (composed of 6% ADH1B, 37% FAP⁺ αSMA, 25% FAP⁺ αSMA⁺, 31% FAP⁺ αSMA⁺) occupy 35% of the stroma with a tumor stroma ratio (TSR) of 1.29. In LUAD, these proportions were partly modified: fibroblasts (composed of 48% ADH1B, 18% FAP⁺ αSMA⁻, 7% FAP⁺⁺ αSMA⁺, 26% FAP⁻ αSMA⁺) occupy 31% of the stroma, with a TSR of 1.07. The fibroblasts most responsible for T-cell exclusion are shown to be FAP⁺ αSMA⁺(Grout et al., 2022). In Figure 3, we compare two microenvironments. In the top one, FAP⁺ αSMA⁺ are lining the tumor (Figure 3a, c, d), and T-cells are excluded from it (Figure 3b). On the contrary, in cases where no FAP⁺ αSMA⁺ are observed, fibroblasts are homogeneously distributed in the stroma (Figure 3g, e, h), and T-cells infiltrate the tumor (Figure 3f). Although no data were directly obtained for the surface fraction of T-cells, the estimated number of CD3⁺ T-cells measured in both tumor islands and stroma suggests a low surface fraction compared to cancer cells and fibroblasts (see Appendix C, Raw data for Table 1).

Figure 3

Download asset Open asset

As we saw in the previous section, the complexity of the microenvironment makes the role of the immune system hardly predictable and highly dependent on the tumor being studied. In the case of lung tumors, the immune system is triggered as the carcinoma expands, but T-cells may be excluded from the tumor nest by activated fibroblasts. Therefore, the goal of this work is to physically and quantitatively understand the process of T-cell exclusion from the tumor mass in the simplest way possible, and to explore different possible scenarios, since the dynamical system can be seen as the spatial integration of the different processes over the domain under study.

We focus on the interaction between different cell types in the case of the NSCLC. In particular, we consider a closed system including the cancer cells, T-cells, non-activated fibroblasts (NAFs), cancer-associated fibroblasts (CAFs) and healthy cells with the extracellular medium. Diffusive signaling molecules are not explicitly introduced: their production by one cell type and their effect on another cell type is modeled as a direct interaction between the two. For example, the attraction of T-cells to cancer cells by chemotaxis is introduced in the mathematical model as a source term proportional to the product of the two concentrations T and C in the T-cell equation (see below Equation 2). We also hypothesize that the main difference between NAFs and CAFs is the fiber production of the latter, which prevents T-cells from infiltrating the tumor. Furthermore, our model does not consider transformation of NAFs into CAFs as a reversible process. All these cells have the same mass density and the sum of their mass fraction satisfies the relationship ${\displaystyle \mathcal{S}=C+T+F_{NA}+F_A=1-N}$, where ${\displaystyle N}$ is a healthy non active component such as healthy cells and interstitial fluid, for example. The mass fraction of cancer cells is represented by ${\displaystyle C}$, T-cells by ${\displaystyle T}$, quiescent or non-activated fibroblasts NAFs by ${\displaystyle F_{NA}}$, activated fibroblasts CAFs by ${\displaystyle F_A}$. Note that ${\displaystyle N}$ which is not an active component does not appear in the following equations. Also, we do not consider here the recruitment of macrophages to better highlight the competing mechanisms related to the sole role of T-cells and CAFs in the tumor mass development. With this in mind, we write an evolution equation for each component of the system: ${\displaystyle dX/dt={\mathrm{\Gamma }}_X}$, where ${\displaystyle {\mathrm{\Gamma }}_X}$ corresponds to a source term modeling the proliferation, death, differentiation, or fluxes into/out of the system under study. The source terms for each species are described in detail below.

The dynamics of the cancer cells are driven by their proliferation, controlled by a growth rate coefficient ${\displaystyle {\alpha }_C}$, and limited by a death rate ${\displaystyle {\tilde{\delta }}_C}$. It takes into account the population pressure caused by self-inhibition as well as by the inhibition of the other species. In the following, we will choose ${\displaystyle {\alpha }_C^{-1}}$ as the unit of time, and all coefficients introduced in the following will be pure constants without unit, so that ${\displaystyle {\tilde{\delta }}_C={\alpha }_C{\delta }_C}$. In addition, cytotoxic T-cells eliminate cancer cells if their anti-tumor activity is not inhibited by the activated fibroblasts (although they do not remove T-cells from the mixture). This process is quantified by the cytotoxic coefficient ${\displaystyle {\delta }_{CT}}$, and by the coefficient of T-cell inhibition by CAFs, ${\displaystyle {\delta }_{TF}}$. With these assumptions, the dynamics of the cancer cell mass fraction can be read:

Although proliferation of cytotoxic T-cells has been observed, we do not consider explicitly proliferation in our study as we focus on their ability to infiltrate the tumor. Rather, we consider that T-cells proliferate outside the domain boundaries, so that this proliferation is included in the boundary source contributions. Therefore, their only source is their attraction to cancer cells, which occurs at a recruitment rate of ${\displaystyle {\alpha }_{TC}C}$, while the crowding limits this recruitment:

Fibroblasts are the key regulators of tumor immunity and progression. Their dynamics involve a recruitment of non-activated fibroblasts ${\displaystyle {\alpha }_{NA}}$ whose role is to maintain an adequate supply given by ${\displaystyle {\alpha }_{NA}/{\delta }_{NA}}$ in healthy tissue, i.e., in the absence of cancer cells, and a death rate due to the pressure exerted by the cells and controlled by ${\displaystyle {\delta }_{NA}}$. NAFs are attracted to the tumor nest at a rate ${\displaystyle {\alpha }_{NA,C}C}$ by the cancer cells that activate them to become CAFs. This process is accounted by introducing a specific transformation rate controlled by the plasticity coefficient ${\displaystyle {\text{K}}_A}$, so that the dynamic equation for NAFs reads:

Thus, the dynamics of the CAFs is purely determined by the transformation of the NAFs and by the pressure through ${\displaystyle {\delta }_A}$, which leads to:

Therefore, the role of fibroblasts in human lung carcinoma can be investigated by studying the interplay between cancer cells and the TME. The interactions discussed above are described by 10 parameters and lead to a system of four coupled nonlinear differential equations concerning four unknowns. Within this framework, we proceed to the estimation of the parameters and the steady states of the system.

When all the parameters vary, the steady states or the equilibrium points of the dynamics can be quite impractical, so we start by assuming that all the coefficients at the origin of a pressure are equivalent: ${\displaystyle \delta ={\delta }_C={\delta }_T={\delta }_A={\delta }_{NA}>0}$. This reduces the number of independent parameters to 7. Fixed points are obtained by setting the time derivatives in Equations 1–4 to 0 which gives the long-term behavior when the system reaches equilibrium.

In order to evaluate the values of the parameters in Equations 1–4, we study different simplified situations that can be reproduced in experiments in vitro or in controlled experiments in vivo. We also use the conclusions obtained in Section Lung TME Composition for the cell densities. We will start by analyzing the case of cancer cells alone and we will successively add all the other cell types with fractions ${\displaystyle C,T,F_{NA},F_A}$. More details on the calculations can be found in Appendix 1, Section A, Dynamical system.

Numerical study of cell population dynamics

The dynamics of each cellular component of the mixture can be systematically studied according to the full set of parameters summarized in Table 2, with the corresponding orders of magnitude. Some of them can be considered as fixed in the system, i.e., they do not vary significantly between the different tumors. This category includes ${\displaystyle {\delta }^{-1}}$, the inverse of the free tumor cell density, and ${\displaystyle {\alpha }_{NA}}$, the NAF attraction parameter to healthy tissue. The other parameters can be studied as control parameters.

Table 2

	${\displaystyle {\delta }_{CT}}$	${\displaystyle {\alpha }_{TC}}$	${\displaystyle {\delta }_{TF}}$	${\displaystyle {\alpha }_{NA,C}}$	${\displaystyle K_A}$
Control of:	Killing C-c by T-c	Attraction of T-c by C-c	Inhibition of T-c by F-c	Attraction of F-c to C-c	Activation of F-c
Efficient T-cells (T-c) but no attraction by (C-c)	${\displaystyle \delta {ϵ}^{-1}}$ ${\displaystyle \sim 10}$	${\displaystyle {\displaystyle {ϵ}^2}}$ ${\displaystyle \sim 0.01}$	No role	No role	No role
Efficient T-cells but inhibited by (F-c)	${\displaystyle \delta {ϵ}^{-1}}$ ${\displaystyle \sim 10}$	${\displaystyle ϵ}$ ${\displaystyle \sim 0.1}$	${\displaystyle \delta {ϵ}^{-1}}$ ${\displaystyle \sim 10}$	1∼1	${\displaystyle \delta }$ ${\displaystyle \sim 1}$
Inefficient T-cellsno need of fibroblasts	${\displaystyle \delta }$ ${\displaystyle \le 1}$	${\displaystyle ϵ}$ ${\displaystyle \sim 0.1}$	No role	No role	No role
Efficient T-cells not inhibited by fibroblasts	${\displaystyle \delta {ϵ}^{-1}}$ ${\displaystyle \sim 10}$	${\displaystyle ϵ}$ ${\displaystyle \sim 0.1}$	${\displaystyle \delta }$ ${\displaystyle \sim 1}$	1∼1	${\displaystyle \delta }$ ${\displaystyle \sim 1}$
Efficient T-cells and fibroblasts not attracted	${\displaystyle \delta {ϵ}^{-1}}$ ${\displaystyle \sim 10}$	${\displaystyle ϵ}$ ${\displaystyle \sim 0.1}$	${\displaystyle \delta {ϵ}^{-1}}$ ${\displaystyle \sim 10}$	${\displaystyle ϵ}$ ${\displaystyle \sim 0.1}$	${\displaystyle \delta }$ ${\displaystyle \sim 1}$
Efficient T-cells fibroblasts not activated	${\displaystyle \delta {ϵ}^{-1}}$ ${\displaystyle \sim 10}$	${\displaystyle ϵ}$ ${\displaystyle \sim 0.1}$	${\displaystyle \delta {ϵ}^{-1}}$ ${\displaystyle \sim 10}$	${\displaystyle ϵ}$ ${\displaystyle \sim 1}$	${\displaystyle \delta ϵ}$ ${\displaystyle \sim 0.1}$

Figure 4

Download asset Open asset

The simulated time-dependent densities of each cell type are displayed in Figure 4 for different sets of parameters in the system of equations in Equations 1–4. At time ${\displaystyle t=0}$, we assume small mass fractions of cancer cells, T-cells, activated and quiescent fibroblasts. Over time, the particular choice of a quadratic model allows the dynamics to reach a plateau for each cell type, confirming the stability of the fixed points found in Section Model Parameters and Fixed Point Analysis.

In the case of an immune-desert tumor, i.e., when T-cells are not attracted to the tumor nest or are unable to penetrate it (Figure 4a–d), cancer cell growth is not limited by the immune response. However, this growth saturates when it reaches a steady state due to the competition for space and resources between the different species and controlled by ${\displaystyle \delta }$ and ${\displaystyle {\alpha }_{NA,C}+{\alpha }_{TC}}$.

T-cells efficiency can be quantified through the parameter ${\displaystyle {\alpha }_{TC}{\delta }_{CT}/({\alpha }_{NA,C}K{\delta }_{TF})}$ (Figure 4i). If the T-cells are efficient, the tumor is said to be immune-inflamed and several scenarios can be discussed. Their response is triggered by the proliferation of cancer cells. Thus, in the absence of CAF inhibition, T-cells inhibit cancer growth (Figure 4f–h). In contrast, in the immune-excluded tumor, CAFs impede the T-cell response and thus indirectly promote cancer cell growth, as observed in Figure 4e.

In conclusion, our dynamical system, limited to four cell types, recapitulates the different possible scenarios that could evolve according to the interaction between cancer cell, T-cells, and fibroblasts in lung adenocarcinoma. With respect to cancer cells, three main outcomes can be identified: the tumor can invade the tissue (Figure 4a), the cancer cell population can be limited without disappearing (Figure 4b–e), or it can be eradicated (Figure 4f–h). Each outcome can have multiple origins, that are all related to either competition for space and resources or T-cell efficiency, see Figure 4i. Each origin has a particular signature with respect to the different species densities, so that the TME properties can be inferred from the measured densities. However, the dynamical system does not provide information about the morphology of the tumor, and its local composition. Indeed, this analysis poorly represents the structure of the tumor, which is divided into a nest, a cancer-associated stroma, and a healthy stroma.

For these reasons, in what follows we extend our model by incorporating the effects of cell species diffusion in order to obtain a more faithful spatio-dynamical description of the system, thus re-analyzing the cases above mentioned in this more general framework.

Since the dynamics are very slow and completely controlled by dissipation, the Onsager variational principle of least dissipation (Ackermann and Ben Amar, 2023) yields the set of partial differential equations coupling the densities to the velocities in the mixture (Ackermann et al., 2021). This principle, introduced by Lord Rayleigh and further developed by Onsager (Strutt, 1871; Onsager, 1931b; Onsager, 1931a; Onsager and Machlup, 1953), is widely used in soft matter (Doi, 2011), in the biophysical context (Ackermann et al., 2021; Balois and Ben Amar, 2014; Wang et al., 2021; Borja da Rocha et al., 2022; Ackermann and Ben Amar, 2023) as well as other areas of physics (Minguzzi, 2015). First we define a free energy functional ${\displaystyle \mathcal{F}}$ from a free energy density ${\displaystyle F}$ integrated over the volume ${\displaystyle V}$, the associated chemical potentials ${\displaystyle {\mu }_i}$, and a dissipation function ${\displaystyle \mathcal{W}}$:

where the ${\displaystyle {\xi }_{ij}}$ are the relative friction coefficients between components ${\displaystyle i}$ and ${\displaystyle j}$, and ${\displaystyle {\mathbf{v}}_{\mathbf{i}}}$ and ${\displaystyle {\mathbf{v}}_{\mathbf{j}}}$ are the velocities of the different cell types. The Rayleighian is then defined as the sum of the dissipation function ${\displaystyle \mathcal{W}}$ and the rate of change of the free energy function ${\displaystyle \mathcal{F}}$. Within this framework, the final equations for the local velocities are obtained by minimizing the Rayleighian with respect to each velocity ${\displaystyle {\mathbf{v}}_{\mathbf{i}}}$:

where ${\displaystyle A_{ij}}$ is the friction matrix and ${\displaystyle {\mu }_i}$ are the chemical potentials. Defining: ${\displaystyle {\tilde{\varphi }}_i={\varphi }_i{\varphi }_0^{-1}}$,the friction matrix reads:

Importantly, in our model, we assume that all frictions ${\displaystyle {\xi }_{ij}}$ are equal, except for the friction with the CAFs. In fact, the mass fraction of CAFs is assumed to reflect the amount of matrix produced by the fibroblasts, resulting in a very high friction in the medium. Therefore, we write: ${\displaystyle {\xi }_{ij}={\xi }_0,{\xi }_{iCAF}={\xi }_0+{\xi }_1{\varphi }_{CAF}}$.

After deriving the momentum equations, we construct the free energy density. Building a free-energy density for a biological material is justified, because, although biological materials are out of equilibrium, their behavior often resembles that dictated by thermodynamics. It is, therefore, useful to write a free energy in terms of state variables. Following (Ackermann et al., 2021) and inspired by the Cahn-Hilliard approach, we define a free energy density ${\displaystyle F}$ as a sum of the interaction potential ${\displaystyle f}$ and of the cost induced by the mass fraction gradients ${\displaystyle \kappa (\mathrm{\nabla }{\varphi }_i{)}^2/2}$. We also assume for ${\displaystyle f}$, the Flory-Huggins free energy density of mixing that depends on the local mass fractions (or equivalently the volume fractions) of each component ${\displaystyle {\varphi }_i}$ (Flory, 1942; Huggins, 1941):

The terms ${\displaystyle D_i{\varphi }_i\log ({\varphi }_i)}$ control both the diffusion and volume exclusion, and the quadratic expansion ${\displaystyle -{\alpha }_{ij}{\varphi }_i{\varphi }_j}$ controls the attraction or repulsion between species.

In the following, we explain in more detail the effective free energy density of our system, as well as the expression of the different proliferation rates, and the values of the parameters, focusing on the attraction/repulsion and the chemotactic terms.

We now apply the general formalism presented in the previous section to the cancer cell mixture using the same notation as Section Dynamic Modeling in the Lung Cancer TME for ${\displaystyle C,T,F_{NA},F_A}$. However, while in the first section, these quantities represented averaged mass fractions over the whole tumor, here they are local mass fractions averaged over a small volume of the tumor but large compared to the cell size. We now detail the effective free energy, focusing first on the interaction potential ${\displaystyle f}$ defined in Equation 9.

The attraction between cancer cells is represented by ${\displaystyle -{\lambda }_{CC}{C}^2}$, between cancer cells and T-cells by ${\displaystyle -{\lambda }_{CT}CT}$, and between all types of fibroblasts and cancer cells by ${\displaystyle -{\lambda }_{CF}C(F_{NA}+F_A)}$, where the ${\displaystyle \lambda }$ coefficients are positive. These interactions are accounted for the free energy density that reads:

Finally, we focus on the proliferation rates to complete the system Equation 8. Note that we must distinguish the cells produced in situ, such as the cancer cells ${\displaystyle C}$ and the activated fibroblasts ${\displaystyle F_A}$, from the cells attracted to the tumor nest (${\displaystyle T}$ and ${\displaystyle F_{NA}}$) by the chemical ${\displaystyle c}$. We keep the growth rate of cancer cells ${\displaystyle {\mathrm{\Gamma }}_C}$ given by Eq. 1 with ${\displaystyle {\delta }_{TF}=0}$ because the inhibition of T-cells by fibroblasts is now treated by increasing the friction created by the fiber barrier. T-cells and NAFs do not proliferate and are not derived from precursors at the tumor site, but are instead generated far away from the tumor and are attracted to it by chemotaxis. Thus, at the boundaries ${\displaystyle S}$, their source rate is driven by an incoming flux due to the chemical ${\displaystyle c}$. For the NAFs, the volume term ${\displaystyle {\mathrm{\Gamma }}_{NA}{|}_V}$ does not include the source term ${\displaystyle {\alpha }_{NA,C}}$ (see Equation 3) which is provided by the boundaries, and ${\displaystyle {\alpha }_{NA}}$ is neglected:

On the other hand, the production of CAFs is not changed and ${\displaystyle {\mathrm{\Gamma }}_A}$ defined in Equation 4 is still relevant, provided that the plasticity contribution is now proportional to the chemical concentration ${\displaystyle c}$ and not the cancer cells fraction.

We scale the various physical parameters in Table 3. Their values are derived from the literature and from the spatial structure of the TME, as well as the expected outcomes for the different scenarios, while the growth parameters are related to those in the first part of the article. For example, the tumor nest is a dense phase of cancer cells with mass fraction ${\displaystyle C=0.8}$ that an interface of size ${\displaystyle d\sim 10\mu m}$ separates from the dilute phase ${\displaystyle C\sim 0}$. This choice, together with the typical value for the energy-density of ${\displaystyle {10}^{6}kg.\mu m^{-1}.d^{-2}}$ found in literature imposes the values of ${\displaystyle D_C,D_0,{\lambda }_{CC},\kappa }$ (Ackermann et al., 2021). Similarly, we assume a weak infiltration of T-cells and fibroblasts into the nest, resulting in low values for the attraction parameters ${\displaystyle {\lambda }_{CT},{\lambda }_{CF}}$ when compared to ${\displaystyle {\lambda }_{CC}}$. Since ${\displaystyle D_{NAF},D_{CAF},D_T}$ control the diffusion of the fibroblasts and T-cells, whose source is located at the boundaries, their values are determined by the density gradients between the boundaries and the tumor nest. At the same time, the times required for T-cells to kill the tumor nest is given by ${\displaystyle {\delta }_{CT}}$, and the times for the arrival of the T-cells and NAFs into the domain, as well as their average mass fractions, are governed by ${\displaystyle {\tau }_T}$ and ${\displaystyle {\tau }_F}$, respectively. The NAFs transform into CAFs in a time controlled by ${\displaystyle K}$. We use the letter ${\displaystyle T_0}$ as the time unit which corresponds to the typical division time for tumor development. For experiments in vitro in perfect conditions of nutrient access, the relevant time is the largest one, corresponding to the slowest process, so the typical division time of cancer cells. It is about 1 d and the same time is needed for the eradication of tumor cells (Cellosaurus, 2023; Xie et al., 2010; Kloss et al., 2013). But in vivo, it is obviously a very different time. It is more related to the time scale necessary for the tumor to grow in its natural environment and is of the order of months (Winer-Muram et al., 2002). The plasticity of fibroblasts is controlled by the chemical penetration length. Since we assume that only the fibroblasts lining the tumor nest are activated, we choose a low penetration length of ${\displaystyle 30\mu m}$. At the same time, the ratio ${\displaystyle \alpha }$ between the chemical production and degradation parameters is set to ${\displaystyle {\alpha }_{cC}=3}$.

In the following, we will show that, once activated by tumor cells, fibroblasts can inhibit tumor growth through a confinement effect but also limit the cytotoxic role of T-cells or simply prevent their infiltration into the tumor.

To illustrate the ambiguous role of fibroblasts in tumor progression, we explore different scenarios through a two-dimensional numerical study. This allows us to compare a fibrotic and a non-fibrotic tumor in the presence or absence of an immune response. Thus, the different growth cases we present are: a tumor without CAFs and T-cells, a tumor with both CAFs and T-cells, a tumor with T-cells but with a low level of CAFs, and a case with CAFs but no efficient T-cells. In the appendix, we present other scenarios of free growth: a TME in the presence of only NAFs and inefficient T-cells (Figure 5—figure supplement 1B), in presence of NAFs (Figure 5—figure supplement 1C), and a TME in which only cancer cells are present (Figure 5—figure supplement 1D). We first consider a single tumor nest before analyzing the case of two adjacent nests. Therefore, there is only one tumor nest at time ${\displaystyle t=0}$ (Figure 5A).

Figure 5 with 1 supplement see all

Download asset Open asset

As shown in Figure 5E and H (blue curve), in the case of a free-growing tumor, i.e., in the presence of NAFs and inactive T-cells, cancer cell growth is not hindered by any obstacles, so this is the most severe situation. However, the tumor provokes the formation of a stroma composed of the NAFs and inefficient T-cells. The presence of a stroma plays a role in the cancer cell growth described in Equation 1 and leads to a growth that would be less important than in a case without stroma or composed only of NAFs (Figure 5—figure supplement 1A and G). In fact, the pressure exerted by the stroma increases the mass fraction of cancer cells in the core of the nest, as well as the total mass fraction ${\displaystyle \mathcal{N}}$ at the stroma-nest interface. This phenomenon corresponds to the competition for space and resources between the different species already described in the previous section. In the case of a free cancer growth without stroma, the tumor nest rapidly invades the environment and its surface fraction in the simulation window is almost 40% after reaching ${\displaystyle t=35T_0}$, with a trend that is not yet saturated and an average mass fraction of cancer cells of 25%.

When fibroblasts are activated, they make the environment around the tumor nest fibrotic (Figure 5B and C), with a mass fraction around the nest reaching 20-30%. As explained earlier (in Section Momentum and Free Energy Density Derivation), the fibers are introduced in our model by an increase of the friction coefficient between the different species and the activated fibroblasts, which is directly related to the fiber concentration. The tumor cells are then trapped behind a barrier with a very high friction which prevents the nest from expanding, and the surface area of the tumor nest decreases to 25% after ${\displaystyle 35T_0}$, with a fibrotic area reaching 70%, while the average cancer cell decreases to 17% (Figure 5E, F and G cyan curve). The fibroblast population is composed of 30% NAFs and 70% CAFs, as the average fibroblast mass fraction is 26% (Figure 5I and K cyan curve).

At the same time, even in a situation where T-cells are efficient, the barrier precludes T-cells from the tumor (Figure 5C). In the latter case, tumor integrity is maintained and CAFs play a tumor-promoting role, inhibiting the immune response and stabilizing the tumor nest at 7.5% of the domain, with a fibrotic zone of 50% (Figure 5E, F and H, green curve). In this scenario, the average mass fractions are 5% for cancer cells, 10% for CAFs, 10% for NAFs (so 20% for the fibroblast population), and 20% for T-cells. It is interesting to note that this scenario quickly leads to a stable steady state. In contrast, the other simulations take longer times to reach a steady state. This could be due to the fact that the stroma builds up quickly compared to cell death when the cell population is not renewed, and that the steady state corresponds to a small nest that is reached after a short period of growth.

When T-cells are introduced without transformation from NAFs to CAFs (Figure 5I-red curve), the tumor dies as cancer cells are eliminated and the compact nest disappears (Figure 5E- red curve). As the stroma initially builds up due to the presence of a tumor, the fibroblast and T-cell populations slowly relax. However, the tumor stroma takes a long time to retract. This suggests that even if the cancer is cured, the effects on the tissue can be long-lasting. In this case of non-activation of NAFs, T-cells infiltrate the core of the nest thanks to the attraction of cancer cells (through the parameter ${\displaystyle {\lambda }_{CT}}$) and allows tumor reduction (Figure 5D). Therefore, their invasion is efficient only in the absence of fibers. (Figure 5G-red curve).

There are many mechanisms able to inhibit the immune system. In addition to the exclusion of T-cells from the tumor and the absence of active feedback, low attractiveness can also reduce the immune response. Indeed, when T-cells are less attracted, for example, because chemotaxis is not efficient enough, the tumor is free to expand. In Figure 5B, the cytotoxicity of T-cells is impaired, but their chemotaxis from the boundaries of the domains is not. This leads to an accumulation of inefficient T-cells around the tumor nest.

Next, we analyze the case where two tumor nests are nucleated and interact. We assume that they have the same initial size and mass fraction of cancer cells. In the absence of CAFs, an immune response is triggered and T-cell infiltration can occur. As expected, in the immune-inflamed tumor, the activity of T-cells reduces the ability of the two nests to coalesce by also reducing their mass fraction, leading to very low values of cancer cells (see Figure 6C). When T-cell activity is marginalized by the presence of CAFs, the two nests slowly coalesce to form a single solid tumor nest (see Figure 6B). In this case, growth is actually limited by the tumor-promoting function of CAFs. Although fibroblasts surround the tumor, its growth is enhanced due to the lack of immune cells inside. On the contrary, in the absence of T-cell chemotaxis, the tumor growth is unrestricted and the nests create a larger tumor by also increasing their mass fraction (see Figure 6A). In the latter cases, coalescence leads to anisotropic shapes of the tumor. Although relaxation in our model eventually leads to a round shape, this relaxation is slower for large tumors and even more for fibrotic tumors where friction significantly slows down this process. Patterns associated with coalescence may thus provide insights into the interpretation of anisotropic patterns in tumors. These are also related to the particular geometry of the system, such as the shape of the organ and the location of various blood vessels. In the next paragraph, we provide different directions to better describe this anisotropy.

Figure 6

Download asset Open asset

To gain insight into the mechanisms leading to anisotropy, we numerically study different cases where the blood vessels are not spatially homogeneously distributed in the vicinity of the tumor nest. Thus, we introduce in our systems localized sources of T-cells and NAF, in the sense of Figure 3, where blood vessels are disk-like patterns in the tissue slice (Figure 7A), or sources corresponding to only a part of the boundaries, as it would be the case for a blood vessel located in the plane of the tissue slice (Figure 7—figure supplement 1). We note that the anisotropy in the sources can lead to non-trivial shapes of the tumor tissue that are not due to instability but rather to the geometry of the system. These vessels are considered static, as dynamic vessels would require a more accurate modeling of angiogenesis. In the case of localized blood vessels modeled as disks, since the interface between the vessel and the system is smaller than in the case of boundary sources, the parameters driving T-cell population growth and T-cell efficiency in killing cancer cells must be increased in order to eradicate the tumor in the absence of cancer-associated fibroblasts. This highlights a peculiar phenomenon: increasing the density of blood vessels has ambiguous consequences. On the one hand, it allows a more efficient immune response and activate the arrival of T-cells. However, on the other hand, it also provides more nutrients to cancer cells and could ultimately favors metastasis. Therefore, critical processes such as the presence of a fibrous barrier around the tumor may have dual effects on and during tumor evolution. Furthermore, we introduce a more complete description of the fibrous stroma by adding a nematic order to the matrix and an anisotropic friction with respect to this matrix. A tensor order parameter ${\displaystyle \mathit{Q}}$ is introduced which characterizes the average orientation of the fibroblasts and fibers and their degree of order (see Appendix 1 Section B, Model for an anisotropic friction with the matrix). Indeed, it has been shown that the density of the extracellular matrix is not always sufficient to induce T-cell exclusion (Carstens et al., 2017). We, therefore, provide the CAFs with an orientation, that is the result of there coupling with the matrix they deposit (Li et al., 2017; Bell et al., 2025). This long range order has been shown to induce an anisotropic friction (Jacques et al., 2023) that may be related to immune cell exclusion (Sun et al., 2021). At the same time, matrix orientation strongly influences metastasis and tumor growth by providing directional cues. Thus, there is a strong difference in the outcome of a tumor with fibers normal to the tumor surface and one with fibers along the tumor, both at the level of immune response and cancer escape. In Figure 7B and C we show the nematic order of the matrix and the flux field of the cancer cells in the case of orthonormal and normal orientation with respect to the tumor nest boundary, and due to the influx of fibroblasts from the right side of the system. The flux of cancer cells is much higher in the case of normal orientation of the fibers than in the case of parallel orientation. In conclusion, the description of T-cell exclusion by the sole density would be adapted to a case with fibers orientation aligned with the cancer nest surface, but not in the case with different orientations. Further investigations may, therefore, introduce these two elements to our model: a nematic order in the CAF layer, whose orientation is determined by coupling with the orientation of the interfaces, and an anisotropic friction, which is higher in the direction normal to the nematic matrix order than in the orientation of the matrix. Note that active stresses (${\displaystyle \mathit{\sigma }=-\zeta \mathbf{Q}}$ with ${\displaystyle \zeta }$ the activity) and couplings with the proliferation can develop in interaction with this nematic order, as explained and modeled in a previous article by some of us (Ackermann and Ben Amar, 2023), and a precise description of the nematic order may require to describe separately the fibroblasts and the matrix (Li et al., 2017; Bell et al., 2025; Jacques et al., 2023).

Figure 7 with 1 supplement see all

Download asset Open asset

We examine the interaction between T-cells and cancer cells in the nest in the absence of fibroblasts ${\displaystyle F_{NA}=F_A=0}$ as described in the main text. This case would lead to the tumor eradication. We look for parameters that lead to significant growth inhibition such that ${\displaystyle C\sim ϵ}$ at equilibrium.

The equilibrium equations read:

Three pairs of solutions exist:

with ${\displaystyle \mathrm{\Delta }=({\alpha }_{TC}\delta +{\alpha }_{TC}{\delta }_{CT}+\delta {)}^2-4{\alpha }_{TC}\delta {\delta }_{CT}}$.

Assuming ${\displaystyle {\delta }_{CT}=\delta {ϵ}^{-1},{\alpha }_{TC}=a_0ϵ}$, with ${\displaystyle ϵ\ll 1,a_0\sim 1}$, the solutions, Equation A4 simplify:

The only solution leading to the tumor eradication is ${\displaystyle \{C,T\}=\{\frac{ϵ}{(a_0-1)\delta },\frac{ϵ}{\delta }\}}$, and ${\displaystyle a_0>1}$.

We now check the linear stability of the different solutions by evaluating the Jacobian eigenvalues. A solution is stable if all its eigenvalues have a negative real part (Strogatz, 2018). The solution found is also the only stable solution provided that ${\displaystyle a_0>1}$. Here, the different eigenvalues of the Jacobian matrix are for the different solutions:

The second solution with ${\displaystyle a_0>1}$ is the only stable solution.

We isolate a subsystem composed of cancer cells, T-cells, and active fibroblasts, as in Lung TME Composition, in order to determine the inhibition rate ${\displaystyle {\delta }_{TF}}$, responsible for the marginalization of the T-cells and subsequently for the increase of cancer cells at a fixed volume fraction of fibroblasts. We first assume that ${\displaystyle F_A}$ is constant. The equilibrium equations read:

where ${\displaystyle {\mathrm{\Delta }}_F=\delta F_A}$, ${\displaystyle {\mathrm{\Delta }}_{CT}={\delta }_{CT}(1+{\delta }_{TF}{F}_A{)}^{-1}\mathrm{a}\mathrm{n}\mathrm{d}{\mathrm{\Delta }}_F=\delta F_A}$.

The different pairs of solutions read:

We assume ${\displaystyle {\delta }_{CT}=\delta {ϵ}^{-1},{\alpha }_{TC}=a_0ϵ,{\delta }_{TF}=d_0\delta {ϵ}^{-1}}$, as justified in the main text. We obtain the following solution pairs:

The only physical solution is the third pair: ${\displaystyle \{C,T\}=\{\frac{f_0}{\delta }+\mathcal{O}(ϵ),\frac{a_0{f}_0}{\delta (f_0+f_a)}ϵ\}}$.

The eigenvalues of the Jacobian matrix read:

${\displaystyle \{{\omega }_1,{\omega }_2\}=\left\{\begin{array}{ccc}-f_a& ,& f_0\\ f_a& ,& (\frac{1}{d_0}+f_0+f_a)\\ -(f_0+f_a)& ,& -f_0\\ \frac{f_a-i\sqrt{f_a}\sqrt{4f_0(d_0{f}_0+2d_0{f}_a+1)+f_a(4d_0{f}_a+3)}}{2}& ,& \frac{f_a+i\sqrt{f_a}\sqrt{4f_0(d_0{f}_0+2d_0{f}_a+1)+f_a(4d_0{f}_a+3)}}{2}\end{array}\right\},}$(A9) and the third solution is the only stable one.

We consider the system described in the main text in the subsection Fibroblast plasticity, where the non-active fibroblasts population is maintained at a fraction ${\displaystyle F_{NA}=f_n{\delta }^{-1}}$ where ${\displaystyle f_n<1}$ is a constant.

We write the equilibrium equations, where we replace ${\displaystyle {\delta }_{CT}=\delta {ϵ}^{-1},{\alpha }_{TC}=a_0ϵ}$, ${\displaystyle {\delta }_{TF}=d_0\delta {ϵ}^{-1}}$:

where ${\displaystyle {\displaystyle \mathcal{S}=C+T+F_A+F_{NA}}}$.

Equating ${\displaystyle C=\delta T\mathcal{S}/(a_0ϵ)}$ and ${\displaystyle C={\delta }^2{F}_A\mathcal{S}/(K_A{f}_n)}$ leads to ${\displaystyle F_A=K_A{f}_{n}T/(\delta a_0ϵ)}$. The system rewrites:

Since the active fibroblasts should prevent T-cells to kill cancer cells, but without reducing the T-cell population, we write ${\displaystyle C\sim {\delta }^{-1}}$ and ${\displaystyle T\sim ϵ{\delta }^{-1}}$. In addition, the term related to the killing of cancer-cells by T-cells should be of order ε. This leads to: ${\displaystyle K_a\sim \delta }$.

We now look for the solutions of the dynamical system, as expansions in ${\displaystyle ϵ}$. Obvious solutions are: ${\displaystyle C=0,T=0}$ and ${\displaystyle C=0,T=-f_nϵ/(\delta ϵ+K_a{f}_n/a_0)}$. We look for the other solutions in the form: ${\displaystyle C=c_0+c_1ϵ,T=t_0+t_1ϵ}$. The order ${\displaystyle 1/ϵ}$ in the dynamical system lead to ${\displaystyle t_0=0}$. The orders 1 and ${\displaystyle ϵ}$ in the dynamical system leads to:

Finally the different solutions are:

with only one physical solution apart from the disappearance of both T-cells and cancer cells.

We now check the stability of the different solutions. The eigenvalues for the growth rates are:

and the solution ${\displaystyle C=\frac{1-f_n}{K_a{f}_n+\delta },T=\frac{a_0ϵ(1-f_n)}{K_a{f}_n+\delta },F_A=\frac{K_A{f}_n(1-f_n)}{\delta (K_a{f}_n+\delta )}}$ is the only stable solution.