Introduction

Cancer results from the malignant transformation of cells due to genetic changes or damage that causes cells to grow and spread in an abnormal and uncontrolled way. Mutated cells can form solid tumor tissue at specific sites and spread to distant regions of the body through metastasis. The growth, development, and response to drugs and therapies of tumor masses are highly dependent on biophysical environmental conditions. These involve a not fully understood crosstalk between cancer cells and the surrounding stroma, which consists of cancer-associated fibroblasts, vessels, and immune cells embedded in the extracellular matrix. Indeed, the tumor-stroma ratio turns out to be an independent prognostic factor, with a large proportion of stroma leading to a worse prognosis [1]. In this paper, we focus on modeling tumorstroma interactions in lung cancer which is one of the most common and deadly cancers, with more than 2.2 million new cases diagnosed and 1.8 million deaths worldwide in 2020. Non-small cell lung carcinoma (NSCLC) accounts for the vast majority of lung cancers (85%) [2]. The majority of cells that make up the stroma are fibroblasts and macrophages. Both types are affected and reprogrammed by cancer cells and have been shown to act as tumor promoters and tumor suppressors, depending on the stage of tumor progression and of its mutational status [3, 4]. As cancer progresses, tumor cells often transform surrounding healthy fibroblasts into cancer-associated fibroblasts (CAFs) which, similarly to a wound healing context, produce higher levels of extracellular matrix (ECM) as well as growth factors and cytokines that affect the recruitment of immune or vascular cells and the growth of cancer cells [5].

Experimental observations and evidence of the mutual interactions between cancer cells, the immune response and cancer-associated fibroblasts suggest that a multi-physical model of the tumor microenvironment (TME), including molecular and spatial dynamics, could be highly beneficial in understanding and predicting the complex evolution of tumor growth and progression. Thus, such analysis mayrepresent a critical step in improving T-cell-based therapies. Quantitative models already exist for biochemical interactions between fibroblasts and cancer cells [6, 7, 8], or for the T-cell recruitment in immunotherapy [9], but models for multispecies are less common. A recent paper by Mukherjee et al. analyzed the spatial dynamics of infiltrating splenocytes in an aggregate of cancerous melanocytes using experimental data obtained in vitro [10]. Their conclusion is that a strong persistent random walk and contact energies are important for the ability of T-cells to infiltrate the tumor. However, to the best of our knowledge, there are no approaches in the current literature to model the dynamics behind human lung tumor proliferation due to the role of CAFs in excluding T-cells from the nest. Here, we propose such a model for NSCLC tumors to describe the dynamic interplay between cancer cells, T-cells and fibroblasts in their non-activated and activated states, incorporating the most relevant interactions within the human tumor microenvironment.

To this end, we briefly review the literature on the interactions between different tumor components, including quantitative information on their composition. In particular, we base our model on data extracted from immuno-histochemistry (IHC) staining of the tumor microenvironment performed on a large cohort of human NSCLC tumor samples. This allows us to calibrate our model so that the densities and localization of the various species that the model yields are consistent with the actual data [11] (Section I). We then build the theoretical description in two steps. In a first step, we ignore any spatial organization and model the proliferation and activation phenomena alone, while introducing a pressure-like term to avoid excessive proliferation. This leads to a nonlinear dynamical system for the cell concentrations (Section II). Careful analysis of this dynamical system yields a stable steady state in which T-cells become inefficient against tumor survival for a given set of parameters. More precisely, we demonstrate that two mechanisms control tumor growth: the competition with other species for resources and space, and the efficiency of the T-cell response. In a second step, we enrich the model by considering the spatio-temporal evolution of tumor growth through a continuum approach for cell mixtures based on Onsager’s variational principle. Indeed, it has been shown that the spatial distribution of T-cells inside the tumor is an indicator of the patient survival [12]. Finite element-based simulations in two-dimensional space demonstrate the ambiguous activity of cancer-associated fibroblasts (CAFs) in regulating cancer cell proliferation and invasion. Finally, we discuss our results obtained with this minimal model and suggest directions for future works, such as considering more precisely the geometry of the tumor micro-environment, and introducing anisotropic friction that would depend on the nematic orientation of the fibers (Section III). In the Supporting Information, further details about the results of numerical simulations are provided to show how the proposed dynamical system is able to capture different evolutionary and tumor fates depending on the role played by the species, their initial conditions, as well as the shapes and distributions of the regions occupied by cells. Appendix A is joined to this article with details for the calculation of the dynamical system section, and Appendix B for the anisotropic friction.

I. Lung Cancer Microenvironment

Modeling the tumor environment requires a deep understanding of the interactions between the various components, such as cells of different types, soft materials and fibers as well as fluids and diffusing molecules. These interactions are summarized in section I.1. Next, in section I.2, we present data on the density of the different cell types in the TME. Indeed, our continuous approach incorporates in vivo biological data that vary between patients and tumor sites, and also evolve over time within the tumor. Therefore, after a review of the available values found in the literature, we performed direct measurements from samples stained by multiplex IHC to study the composition of fibroblasts and distribution of T-cells in human

I.1 Interactions Between TME Components

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 [13]. 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 [14]. 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 [15, 16, 17, 18]. 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 [19, 11].

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 [20, 21]. 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 [22, 23]. In addition, biochemical factors expressed by CAFs also help to modify the phenotype of T-cells or inactivate their cytotoxic capacity [24, 25, 26]. 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 (Fig. 1B), which is less proliferative than a free-growing tumor (Fig. 1A). Failure to confine the tumor induces high levels of cytotoxic T-cell infiltration, creating a hot tumor (Fig. 1C). Therefore, one aim of this study is to predict susceptibility of a tumor to T-cell infiltration according to the scheme given in Fig. 1.

Figure 1:

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 [27]. 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⁺ [28, 29, 30]. The origins of CD8⁺ T-cells that infiltrate the tumor are diverse, as they differentiate from both circulating and tissue-resident precursors [31]. 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 [32]. The absence of such chemicals leads to the formation of the so-called immune-desert tumor (Fig. 1A).

1.2. Lung TME Composition

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 10X 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 10X 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 can be between 13-18 µ m [33], fibroblasts ∼ 16 µm m [34], and T-cells between 5-10 µm [35], resulting in volumes and projected areas that can be 10 times smaller for T-cells compared to fibroblasts and can-cer 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 [36, 37]. Despite these remarks, cancer cells, T-cells, and macrophages (fibroblasts are underestimated as mentioned above) appear to be a significant part of the 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 [11]. 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⁺ [11]. In Fig.2, we compare two LUSC microenvironments. In the top one, FAP⁺αSMA⁺ are lining the tumor (Fig.2(a),(c),(d)), and T-cells are excluded from it (Fig.2(b)). On the contrary, in cases where no FAP⁺αSMA⁺ are observed, fibroblasts are homogeneously distributed in the stroma (Fig.2(g),(e),(h)), and T-cells infiltrate the tumor (Fig. 2(f)).

Figure 2:

Although no data for T-cells were obtained from the last measurements, it turns out that since the surface is occupied at 80% by cancer cells and fibroblasts, and that the rest of the TME occupies less than 20%, the surface fraction of T-cells can be considered as small compared to cancer cells and fibroblasts.

II. Dynamic Modeling in the Lung Cancer TME

Our theoretical and numerical analyses consist of two steps in the spirit of Ref. [47] for the study of cancer stem cells. In a first analysis, we examine only the dynamics of an interactive ecological system in order to evaluate the physical parameters that quantify these interactions and how the dynamics depend on them. Spatial constraints are represented by a pressure term avoiding overcrowding. We present this approach step by step in order to set the parameters one by one, highlighting the physical importance of each choice through stability analyzes of the system. In addition, this first step allows to explore efficiently the parameter space without time consuming numerical resolution, and to already produce a first classification of the global outcomes. The second step is the spatial description of the tumor growth in Section III. In both approaches, the case where fibroblasts, cancer cells and T-cells are present is intended to correspond to the patient data presented in the previous section.

II.1. Dynamical System for Immune and Cancer Cells in Interaction

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 Eq. 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 , where N is a healthy non active compo-nent such as healthy cells and interstitial fluid, for example. The mass fraction of cancer cells is represented by C, T-cells by T, quiescent or non-activated fibroblasts NAFs by F_{N A}, activated fibroblasts CAFs by F_A. Note that N is not an active component and therefore 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: d X /dt = Γ_X, where Γ_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 is driven by their proliferation, controlled by a growth rate coefficient α_C, and limited by a death rate . 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 as the unit of time, and all coefficients introduced in the following will be pure constants without unit, so that . 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 δ_CT, and by the coefficient of T-cell inhibition by CAFs, δ_{T F}. 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 α_TCC, 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 α_{N A} whose role is to maintain an adequate supply given by α_{N A}/δ_{N A} 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 δ_{N A}. NAFs are attracted to the tumor nest at a rate α_{N A,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 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 δ_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 11 parameters and lead to a system of 4 coupled nonlinear differential equations concerning 4 unknowns. Within this framework, we proceed to the estimation of the parameters and the steady states of the system.

II.2 Model Parameters and Fixed Point Analysis

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: δ = δ_C = δ_T = δ_A = δ_{N A} > 0. This reduces the number of independent parameters to 7. Fixed points are obtained by setting the time derivatives in Eqs. (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 Eqs. (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. I.2 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 C, T, F_{N A}, F_A. More details on the calculations can be found in Appendix A.

II.2.1. Cancer cells and T-cells

The dynamic evolution of cancer cells alone is limited to: dC /dt = C − δC ². This equation leads to two equilibrium fixed points: C = 0 and C = δ⁻¹ and only the second one C = δ⁻¹ is stable. When cancer cells are isolated from other active cellular components, they are expected to invade the system or to be its major component, leading to a cancer death rate of δ ∼ 1.

We now examine the interaction between T-cells and cancer cells in the nest in the absence of fibroblasts (so F_{N A} = F_A = 0) and the relevant parameters scaled by δ are δ_CT and α_TC. Then, we study the equilibrium regime:

There are three equilibrium solution pairs {C, T}, including the trivial one: {0, 0}. To analyze whether the solutions found are physically relevant (0 ≤ C, T ≤ 1) and dynamically stable, we estimate possible scaling for the two parameters, focusing on an effective immune response against cancer. For efficient elimination of cancer cells by T-cells, the killing rate δ_CT must be much larger than the natural death rate δ, so we introduce a small parameter 0 < ϵ ≪ 1, so that δ_CT = δϵ⁻¹. We also assume that the T-cell recruitment is slow compared to cancer cells, which means that α_TC = a₀ϵ, where a₀ being of order one. So the only stable solution is where a₀ > 1 (see Appendix A.1). Thus, even if the inflammation level is low, resulting in a small number of T-cells, the immune action on the cancer cells remains efficient.

II.2.2. Role of activated fibroblasts on T-cells

Fibroblasts play an active role in the exclusion of the T-cells from the tumor nest. We isolate a subsystem composed of cancer cells, T-cells and active fibroblasts to determine the inhibition rate δ_{T F}, responsible for the marginalization of the T-cells and subsequently for the increase of cancer cells at a fixed volume fraction of fibroblasts. We consider the fixed points corresponding to an ensemble {C, T, 0, F_A}. For simplicity, we first assume that F_A is constant, which leads to the dynamical system:

where f₀ = 1 − δF_A, Δ_CT = δ_CT (1 + δ_{T F} F_A)⁻¹ and Δ_F = δF_A. There are four solution pairs (see Appendix A.2).

We use the scaling of both parameters already established in the previous paragraph: δ_CT = δϵ⁻¹, α_TC = a₀ϵ. We also assume that the inhibition of the T-cells may counteract their cytotoxic effect on cancer cells, i.e. (1+δ_{T F} F_A)⁻¹ ∼ ϵ, when activated fibroblasts are abundant F_A = δ⁻¹ f_a. This results in δ_{T F} = d₀δϵ⁻¹.

The only stable equilibrium solution is , where the notation denotes a quantity whose order of magnitude is ϵ. Thus, when fibroblasts inactivate T-cells, they promote the cancer cell proliferation.

II.2.3. Residual Fibroblasts in Healthy Tissue

In healthy tissue, most of the fibroblasts are in a quiescent state and are not activated in the absence of pathologies such as wounds, allergic reactions or cancer cells. Therefore, in such tissues, for a quiescent fibroblast population, the density of F_{N A} is the equilibrium solution of the equation d F_{N A}/dt = 0 = α_{N A} − δF_{N A}². The parameter α_{N A} represents the net influx of fibroblasts into the tissue. Assuming the fraction of non-activated fibroblasts to be residual ∼ ϵ², we deduce α_{N A} ∼ δϵ⁴.

II.2.4. Fibroblast Plasticity

Fibroblast plasticity is the phenotypic change responsible for T-cell inhibition and for more active fiber production. Cancer cells drive this transformation of the current population F_{N A} into F_A. This process is quantified by the constant K_A. We first estimate that the fibroblast population is comparable to the cancer cell population when they are alone, which leads to: F_{N A} = f_nδ⁻¹ where f_n ∼ 1 is a constant.

Replacing δ_CT = δϵ⁻¹, α_TC = a₀ϵ, δ_{T F} = d₀δϵ⁻¹, leads to:

Because of the effect of the fibroblasts on them, the T-cells have little effect on the cancer cells, so cancer cell proliferation is little affected. At equilibrium, this leads to C ∼ δ⁻¹, T ∼ ϵδ⁻¹. The equation for T gives F_A ∼ a₀δ⁻¹/d₀. This means that the activated fibroblasts must be relatively abundant to fully inhibit the activity of the T-cells. It follows that the plasticity parameter is of low order in ϵ, since the non-activated fibroblasts are maintained at a high density: K_A ∼ δ (see Appendix A.3).

II.2.5. Tumor Fibroblast Attraction

In this last step, we want to determine the parameter that controls the attraction of the fibroblasts to the tumor. To do this, we write K_A = kδ in the system of equations Eqs. 1→4, which we rewrite according to the previous findings:

Then, we look for solutions of the type: {C, T, F_{N A}, F_A} = δ⁻¹{c₀, T₀ϵ, f_{N A}, f _A}. The order of magnitude of the attraction parameter is then α_{N A,C} ∼ 1 (see Appendix A.4).

In the next section, we summarize all scaling laws possible, according to the obtained results and present the different scenarios related to the state of the tumor, shown in Fig.1.

II.2.6. 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 δ⁻¹, the inverse of the free tumor cell density, and α_{N A}, the NAF attraction parameter to healthy tissue. The other parameters can be studied as control parameters.

Table 1:

Table 2:

Table 3:

The simulated time-dependent densities of each cell type are displayed in Fig. 3 for different sets of parameters in the system of equations in Eqs. (1)-(4). At time 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 II.2.

Figure 3:

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 (Fig. 3(a)-3(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 δ and α_{N A,C} +α_TC.

T-cells efficiency can be quantified through the parameter α_TC δ_CT /(α_{N A,C} K δ_{T F}) (Fig. 3(i)). 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 (Fig. 3(f)-3(h)). In contrast, in the immune-excluded tumor, CAFs impede the T-cell response and thus indirectly promote cancer cell growth, as observed in Fig. 3(e).

In conclusion, our dynamical system, limited to 4 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 (Fig. 3(a)), the cancer cell population can be limited without disappearing (Fig. 3(b)-3(e)), or it can be eradicated (Fig. 3(f)-3(h)). Each outcome can have multiple origins, that are all related to either competition for space and resources or T-cell efficiency, see Fig. 3(i). 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-analysing the cases above mentioned in this more general framework.

III. Spatio-Temporal Behavior of Tumor Growth

Spatial study is an important tool in cancer diagnosis since the tumor shape reveals the aggressiveness of cancer cells and the role of their microenvironment[48, 49, 50, 51, 52, 53, 54, 55]. Similarly, the localization of immune cells in the tumor environment [10], as well as proliferating and pre-metastatic cells often found in niches, are active areas of research in oncology and a valuable support for clinical prognosis [56, 57, 58].

In the previous section, we described our cell mixing through the global time evolution given by a dynamical system. Here, we aim to complete the modeling by considering the spatial het-erogeneity of the interplay between fibroblasts, T-cells and cancer cells and its consequences for tumor cell localization and proliferation. We first present the mixture model [52, 59], which is able to incorporate the spatial distribution and evolution of active cells leading to a set of partial differential equations that we solve with the finite element (FEM) software COMSOL Multiphysics (COMSOL Multiphysics® v. 6.2. www.comsol.com. COMSOL AB, Stockholm, Sweden.). Last, we propose different mechanisms as sources of anisotropy in the tumor nest, including the introduction of a nematic tensor for the orientation of fibroblasts.

As in the previous section, we consider 4 different cell types namely cancer cells, T-cells, activated and non-activated fibroblasts, with different attraction properties. A component that does not play an active role in the mixture is also added via an inactive N fraction. It concerns the intercellular fluid, the healthy cells and the dead cells. This last component is also a source of material, for the proliferation of cancer cells. Each component is described by a local mass fraction φ_i, a velocity v_i, and a proliferation/death term Γ_i, with φ and Γ denoting generic quantities that are now space and time dependent. More specifically, φ_i corresponds to the local value of C, T, F_A, F_{N A} introduced in the previous sections and φ₀ = 1 − ∑_i≠0 φ_i repre-sents N. Γ_i is the growth rate of each component, which is now space and time dependent (see Eqs. (1)-(4)) and may be positive only for cancer cells, since we consider that fibroblasts and T-cells are recruited from the surrounding environment. CAFs produce a significant amount of fibers, resulting in a higher friction between different species, chosen to be proportional to the local amount of CAFs. This friction, which is also a space-time dependent quantity, will impact the dynamics of the mixture.

The sample also contains diffusive signaling molecules that are at the origin of the immune activity. In the previous section, chemicals were not introduced, although they were associated with some coefficients of the dynamical system. Here, the chemicals that determine the cell behavior, are represented by a concentration c _j [47, 60]. Note that these chemicals have no mass and diffuse through the mixture with the diffusion coefficient D_c. For simplicity, we restrict ourselves to a single chemical of concentration c, that mediates both chemotaxis and activation of fibroblasts. The balance between its production and degradation rate writes , since it is produced by cancer cells and naturally degraded. With these considerations, we study the case where the tumor is well supplied with nutrients, which are therefore not explicitly mentioned, and we write a set of conservation equations for each component of the mixture:

When integrated over the whole space and when the only boundary conditions for the velocity are considered, the first equation in (8) recovers the dynamical system. Besides, we consider that chemicals equilibrate much faster than the tissue dynamics, so that: λ_c∇²_c+α_cC φ_C −c = 0, with the penetration length and α_cC = (δ_cτ_cC)⁻¹. We now present a derivation for the average local velocity v_i of each species by evaluating its momentum equation and the expressions for the various source terms.

III.1 Momentum and Free Energy Density Derivation

Since the dynamics is very slow and completely controlled by dissipation, the Onsager variational principle of least dissipation [61] yields the set of partial differential equations coupling the densities to the velocities in the mixture [59]. This principle, introduced by Lord Rayleigh and further developed by Onsager [62, 63, 64, 65], is widely used in soft matter [66], in the biophysical context [59, 52, 67, 68, 61] as well as other areas of physics [69]. First we define a free energy functional ℱ from a free energy density F integrated over the volume V, the associated chemical potentials µ_i, and a dissipation function :

where the ξ_{i j} are the relative friction coefficients between components i and j, and v_i and v_j are the velocities of the different cell types. The Rayleighian is then defined as the sum of the dissipation function and the rate of change of the free energy function ℱ. Within this framework, the final equations for the local velocities are obtained by minimizing the Rayleighian with respect to each velocity v_i:

where A_{i j} is the friction matrix and µ_i are the chemical potentials. Defining: friction matrix reads:

Importantly, in our model we assume that all frictions ξ_{i j} 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: ξ_ij = ξ₀,ξ_iCAF = ξ₀+ξ₁ϕ_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 [59] and inspired by the Cahn-Hilliard approach, we define a free energy density F as a sum of the interaction potential f and of the cost induced by the mass fraction gradients κ(∇φ_i)²/2. We also assume for f, the Flory-Huggins free energy density of mixing that depends on the local mass fractions (or equivalently the volume fractions) of each component φ_i [70, 71]:

The terms D_i φ_i log(φ_i) control both the diffusion and volume exclusion, and the quadratic expansion −α_{i j} φ_i φ_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.

III.2. Effective Free Energy Density and Source Terms

We now apply the general formalism presented in the previous section to the cancer cell mixture using the same notation as Section II for C, T, F_{N A}, F_A. However, whereas 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 f defined in Eq. 9.

The attraction between cancer cells is represented by −λ_CCC ², between cancer cells and T-cells by −λ_CT CT, and between all types of fibroblasts and cancer cells by −λ_CF C (F_NA+F_A), where the λ 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 Eq. 8. Note that we must distinguish the cells produced in situ, such as the cancer cells C and the activated fibroblasts F_A, from the cells attracted to the tumor nest (T and F_{N A}) by the chemical c. We keep the growth rate of cancer cells Γ_C given by Eq. 1 with δ_{T F} = 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 S, their source rate is driven by an incoming flux due to the chemical c. For the NAFs, the volume term Γ_{N A}|_V does not include the source term α_{N A,C} (see Eq. 3) which is provided by the boundaries, and α_{N A} is neglected:

On the other hand, the production of CAFs is not changed and Γ_A defined in Eq. 4 is still relevant, provided that the plasticity contribution is now proportional to the chemical concentration 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 ar-ticle. For example, the tumor nest is a dense phase of cancer cells with mass fraction C = 0.8 that an interface of size d∼ 10 µ m separates from the dilute phase C ∼ 0. This choice, together with the typical value for the energy-density of 10⁶ kg.µm⁻¹.d⁻² found in literature imposes the values of D_C, D₀, λ_CC, κ [59]. Similarly, we assume a weak infiltration of T-cells and fibroblasts into the nest, resulting in low values for the attraction parameters λ_CT, λ_CF when compared to λ_CC. Since D_{N AF}, D_{C AF}, 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 δ_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 τ_T and τ_F, respectively. The NAFs transform into CAFs in a time controlled by K. We use the letter T₀ 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 one day and the same time is needed for the eradication of tumor cells [72, 73, 74]. 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 [75]. 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 30 µm. At the same time, the ratio α between the chemical production and degradation parameters is set to α_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.

III.3. Ambiguous Role of Fibroblasts in Tumorigenesis

To illustrate the ambiguous role of fibroblasts in tumor progression, we explore different scenarios through a 2-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 Supplementary Information, we present other scenarios of free growth: a TME in the presence of only NAFs and inefficient T-cells (Fig. Sup. 1B), in presence of NAFs (Fig. Sup. 1C), and a TME in which only cancer cells are present (Fig. Sup. 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 t = 0 (Fig. 4A).

Figure 4:

As shown in Fig. 4E,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 Eq. 1 and leads to a growth that would be less important than in a case without stroma or composed only of NAFs (Fig. Sup. 1A→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 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 t = 35 T₀, 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 (Fig. 4B,C), with a mass fraction around the nest reaching 20-30%. As explained earlier (in Section III.1), 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 35 T₀, with a fibrotic area reaching 70%, while the average cancer cell decreases to 17% (Fig. 4E,F,G cyan curve). The fibroblast population is composed of 30% NAFs and 70% CAFs, as the average fibroblast mass fraction is 26% (Fig. 4I,K cyan curve).

At the same time, even in a situation where T-cells are efficient, the barrier precludes T-cells from the tumor (Fig. 4C). 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% (Fig. 4E,F,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 (Fig. 4I-red curve), the tumor dies as cancer cells are eliminated and the compact nest disappears (Fig. 4E-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 λ_CT) and allows tumor reduction (Fig. 4D). Therefore, their invasion is efficient only in the absence of fibers. (Fig. 4G-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 Fig. 4B, 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 immuneinflamed 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 Fig. 5C). When T-cell activity is marginalized by the presence of CAFs, the two nests slowly coalesce to form a single solid tumor nest (see Fig. 5B). 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 Fig. 5A). In the latter cases, coalescence leads to anisotropic shapes of the tumor. Al-though 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 5:

III.4. Toward anisotropy

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 Fig. 2, where blood vessels are disk-like patterns in the tissue slice (Fig. 6A), 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 (Fig. 2). 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.

Figure 6:

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 Q is introduced which characterizes the average orientation of the fibroblasts and fibers and their degree of order (see Appendix B). Indeed, it has been shown that the density of the extracellular matrix is not always sufficient to induce T-cell exclusion [12]. We therefore provide the CAFs with an orientation, that is the result of there coupling with the matrix they deposit [76, 77]. This long range order has been shown to induce an anisotropic friction [78] that may be related to immune cell exclusion [79]. 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 Fig. 6B,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 (σ = −ζQ with ζ 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 [61], and a precise description of the nematic order may require to describe separately the fibroblasts and the matrix [76, 77, 78].

IV. Discussion

Our work provides a physical model to quantify the role of the immune system in human lung tumors. This is in line with different similar mathematical models, that study through this lens the inhibition/activation of the immune system by cancer cells either by the mean of nonlinear compartment models resembling our dynamical system, for instance regarding macrophage enrolment and cytokine signaling [80, 81], or mixture models [82]. We thus con-nect the two approaches in order to rigorously derive the parameters of the model and derive insights from both. In our article, the focus is on early tumor growth prior to angiogenesis and the development of its stroma rich in T-cells and fibroblasts, the activation of non-activated fibroblasts into cancer-associated fibroblasts and the marginalization of T-cells from the tumor nest. Other immune cells, such as macrophages, could complete this study in a further investigation.

After analyzing the data from the literature and from patients with LUSC or LUAD pathologies, we propose a physical model, both theoretical and numerical, for the interactions between cancer cells, T-cells and fibroblasts during tumor progression and their different roles. In particular, fibroblasts are introduced in the model in their inactive and active states (NAFs and CAFs). The modeling involves two steps: firstly, a dynamic study aimed at elucidating the key factors that can characterize the properties of the micro-environment on which the ability of T-cells to inhibit tumor growth depends. Scenarios scaled by only two parameters control the dynamics and evaluate the aggressiveness of the tumor. In particular, we have established the spatial organization of the different cell types inside and outside the tumor core. The model is based on the continuous mixture model derived from Onsager’s variational principle. Simulations were performed with the FEM software COMSOL Multiphysics in a two-dimensional framework. This provided a complete description of the tumor morphology and composition. In fact, we explored different scenarios to fully appreciate the role of all cell types involved. The results confirmed the experimental evidence that CAFs play a dual role in promoting and suppressing cancer cell proliferation. First, they alter the tumor microenvironment with a fiber barrier that reduces the motility and activity of T-cells, thereby promoting cancer cell growth. At the same time, this growth is limited and results to be lower than in the absence of active fibroblasts. Moreover, we have also shown that although the different parameters can take a wide range of values and continuously change the results of the numerical study, different scenarios can be drawn by looking only at the orders of magnitude.

In this study, we did not consider the fibers of the extracellular matrix as a component per se, but we consider them through the friction increase. We also provided the basics and started to describe the consequences of an anisotropic friction in relation to a nematic fiber orientation (see Appendix B.2). The same formalism could also integrate non trivial active stresses as done in Ref. [61].

Our physical model aims to limit the number of parameters as much as possible. This may seem far-fetched given the complexity of the biological system, especially in vivo. However, there are several reasons for this limitation. First, determining the range of values for these parameters is a complex task in the absence of direct measurement. The strategy we used in Section II is to isolate each process corresponding to each parameter and thus construct the model step by step. In Section III we used both this latter strategy and the spatial structure of the mixing. In addition, the different combinations of the parameter values provide a wide variety of scenarios and it would be sufficient to modify these values to model new chemical entering the system. Limiting the complexity of the different processes allows us to focus on other aspects of the problem, such as the time evolution and the spatial structure, as well as a more precise quantitative study. In the same spirit, we numerically report in Section III a 2-dimensional system without substrate, which yields results that we can consider similar as a 3-dimensional system with an invariance in the third dimension.

Identifying the different scenarios that can occur in the TME and their likelihood is critical to predict the outcome of therapy [60]. Indeed, both pharmacodynamics and pharmacokinetics are highly dependent on the spatial structure and composition at the tumor site. In particular, the dynamic of arrival of the T-cell population is over simplified here via chemotaxis and discard auto-chemotaxis at the origin of a swarming displacement of an increased number of cytotoxic T-cells eventually, CAR-T cells [83]. This mechanism, rare in the solid-tumor case, is sequential in time leading to a competition between killing efficiency and cancerous cell repair [84, 85] which will depend on the cancer type considered. Classification and quantitative characterization of the various outcomes can also help to monitor the treatment in real time according to the evolution of the TME.

In conclusion, this numerical investigation may help to understand the limitations of the immune system in the face of solid tumor growth. Understanding quantitatively how immune cells are excluded from the tumor nest may be helpful in the drug design of T-cell-based therapies. Conversely, the numerical reproduction of various processes that can be observed in vivo is an important step in the context of personalized medicine. A next step in the theoretical and numerical investigation would therefore be to introduce drug molecules into the framework we have presented. Indeed extending our model to follow the tumor evolution under treatment and the associated internal dynamics from the early alveolar stage of the lung by incorporating the actual three-dimensional environment would also be a major step forward.

Supporting information

In this appendix, we show cases of free growth but in simple stroma: the case of a stroma composed by NAFs and inefficient T-cells (Fig. Sup. 1B), a stroma composed by NAFs only (Fig. Sup. 1C), and the case without stroma (Fig. Sup. 1D), the tumor nest surface fraction in the simulation window for the different scenarios (Fig. Sup. 1E), and the average mass fraction of the cancer cells and the NAFs (Fig. Sup. 1F,G. The values of the parameters depending on the scenario are shown in Table. Sup 1.

Table Sup 1:

Besides, we show the evolution of the tumor in the case with efficient T-cells and NAFs without CAFs (Fig. Sup. 1H).

We also obtain different configurations by assuming deformed domains in order to study the anisotropic response of the tumor stroma (Fig. Sup. 2A,B).

Figure Sup 1:

In particular, the results obtained differ from those presented above by introducing T-cells and NAFs only from the upper side of the boundaries, this assuming the possibility that T-cells are able to reach the nest from different parts. Finally, different regions have been highlighted: in light blue the cancer cells nest, in red the barrier of activated fibroblasts and in green the healthy stroma composed of NAFs and T-cells.

Figure Sup 2:

Data and materials availability

All data are available in the main text, in the figures and in the tables. Comsol Multiphysics files and raw data used in this work are available upon request.

Acknowledgements

M.BA and H.S acknowledge the financial support from ITMO Cancer of Aviesan within the framework of the 2021-2030 Cancer Control Strategy, on funds administrated by Inserm (PCSI 2021,MCMP 2022). J.A acknowledges the financial support from ANR COLLAMOEBOID (ANR-20-CE13-0031). M.F. and C.B. acknowledge the financial support under the National Recovery and Resilience Plan (NRRP), Mission 4, Component 2, Investment 1.1, Call for tender No. 104 published on 2.2.2022 by the Italian Ministry of University and Research (MUR), funded by the European Union – NextGenerationEU– Project Title 2022ATZCJN AMPHYBIA – CUP E53D23003040006 - Grant Assignment Decree No. 961 adopted on 30.06.2023 by the Italian Ministry of Ministry of University and Research (MUR). We thank Philippe Benaroch, Layla Mathieson, Xudong Xing for useful discussions.

Additional information

Author contributions

M.BA and J.A conceived and designed the project, built the physical model and performed the theoretical analysis, assisted by M.F and C.B. H.S and P.S provided the biological data and images. C.B and J.A performed the numerical simulations. M.BA, J.A, and C.B wrote the manuscript with contributions from P.S, H.S, and M.F.

Appendix A. Dynamical System

In this section, we detail the calculations referred to in the dynamical system part of the article (Section II). More precisely, we demonstrate the scaling of the various parameters in the case where activated fibroblasts prevent T cells from killing the cancer cells. With the found scaling, we show that the stationary solutions of this dynamical system are the only stable solutions.

First, we recall the complete system of equations describing the dynamics of the tumor microenvironment. C, T, F_A, F_{N A} represent the cancer cells, T cells, active fibroblasts, and inactive fibroblasts, respectively.

where δ_CT is the killing rate of cancer cells by T cells in the absence of active fibroblasts, δ_{T F} is the inhibition constant of T cells by the active fibroblasts, α_TC is the attraction rate of T cells to the tumor, α_{N A} is the attraction rate of inactive fibroblasts in the healthy tissue, α_{N A,C} is the attraction rate of non-active fibroblasts in the presence of cancer cells, K_A is the plasticity constant between active and non-active fibroblasts, and δ_C, δ_T, δ_{N A}, δ_A are the death rates due to space and resource constraints of the different species. For simplicity, in this article we consider δ_C = δ_T = δ_{N A} = δ_A = δ. We also define the small parameter ϵ ≪ 1.

Appendix A.1. Interaction between T-cells and cancer cells

We examine the interaction between T cells and cancer cells in the nest in the absence of fibroblasts F_{N A} = F_A = 0 as described in Section II.2.2. This case would lead to the tumor eradication. We look for parameters that lead to significant growth inhibition such that C ∼ ϵ at equilibrium.

The equilibrium equations read:

Three pairs of solutions exist:

with Δ = (α_TC δ+α_TC δ_CT +δ)² − 4α_TC δδ_CT.

Assuming δ_CT = δϵ⁻¹, α_TC = a₀ϵ, with ϵ ≪ 1, a₀ ∼ 1, the solutions, eq.(B.3) simplify:

The only solution leading to the tumor eradication is , and a₀ > 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 [86]. The solution found is also the only stable solution provided that a₀ > 1. Here,the different eigenvalues of the Jacobian matrix are for the different solutions:

The second solution with a₀ > 1 is the only stable solution.

Appendix A.2. Role of activated fibroblasts on T-cells

We isolate a subsystem composed of cancer cells, T-cells and active fibroblasts, as in Section II.2.2, in order to determine the inhibition rate δ_{T F}, 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 F_A is constant. The equilibrium equations read:

where f₀ = 1 −δF_A, Δ_CT = δ_CT (1 +δ_{T F} F_A)⁻¹ and Δ_F = δF_A.

The different pairs of solutions read:

We assume δ_CT = δϵ⁻¹, α_TC = a₀ϵ, δ_{T F} = d₀δϵ⁻¹, as justified in the main text (Section II.2.2). We obtain the following solution pairs:

The only physical solution is the third pair: .

The eigenvalues of the Jacobian matrix read:

and the third solution is the only stable one.

Appendix A.3. Fibroblast plasticity

We consider the system described in Section II.2.4, where the non-active fibroblasts population is maintained at a fraction F_{N A} = f_nδ⁻¹ where f_n < 1 is a constant.

We write the equilibrium equations, where we replace δ_CT = δϵ⁻¹, α_TC = a₀ϵ, δ_{T F} = d₀δϵ⁻¹:

where .

Equating and leads to F_A = K_A f_nT /(δa₀ϵ). The system rewrites:

Since the active fibroblasts should prevent T-cells to kill cancer cells, but without reducing the T-cell population, we write C ∼ δ⁻¹ and T ∼ ϵδ⁻¹. In addition, the term related to the killing of cancer-cells by T-cells should be of order . This leads to: K_a ∼ δ.

We now look for the solutions of the dynamical system, as expansions in ϵ. Obvious solutions are: C = 0, T = 0 and C = 0, T = −f_nϵ/(δϵ+K_a f_n/a₀). We look for the other solutions in the form: C = c₀ +c₁ϵ, T = t₀ + t₁ϵ. The order 1/ϵ in the dynamical system lead to t₀ = 0. The orders 1 and ϵ 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 is the only stable solution.

Appendix A.4. Non active fibroblasts attraction to tumor

Last, we demonstrate the results of Section II.2.5. In order to analyse non active fibroblasts to tumor, we now consider the full dynamical system:

We look for solutions of the type: for the case of fibroblasts preventing T-cells getting access to the tumor.

The dynamical system becomes to the first order in ϵ:

In order for f _A, f_{N A} to be of order 0 in ϵ, α_{N A,C} must also be of order 0.

At the lowest order in ϵ we obtain:

The solutions of the full system write:

Only the second solution is physical.

We now check the stability of those solutions. The linear stability analysis yields the eigenvalues for the stability rates:

The only stable solution is the second solution where all cell fractions are above 0 and below 1.

Appendix B. Model for an anisotropic friction with the matrix

Appendix B.1. Nematic order

We introduce a refined model for the fiber role in the tumor microenvironment. We first add a nematic order to the CAF population since its structure involves a more or less disordered nets of fibers. This nematic order is described by a traceless nematic tensor Q. We recall that the nematic tensor may be written as:

where S is the amplitude of the anisotropy, and θ the angle for the nematic orientation. Besides, Q is determined by the minimization of the following free-energy:

where α = −F_A, and Φ indicates the orientation of the tumor nest boundary: . The orientation of the fibers with regards to the tumor nest is controlled by the parameter γ. In the spirit of our work, we assume that fibers orient fast compared to the growth times, so that we consider the following equation:

which allows to evaluate Q.

Appendix B.2. Anisotropic friction

We consider an anisotropic friction between the different components and the CAFs in the dissipation function.

where is the velocity component along the unit vector of the nematic order orientation of the matrix while X represents one component of the mixture . In the same way, . after some algebra, this leads to the nematic friction tensor that depends only on the tensor Q. More justifications can be found in Ref. [78].