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 a process known as 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 In this paper, we focus on modeling tumor-stroma 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%) [1]. 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 [2, 3]. 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 [4].

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 may represent a critical step in improving T-cell-based therapies. Quantitative models already exist for biochemical interactions between fibroblasts and cancer cells [5, 6, 7], or for the T-cell recruitment in immunotherapy [8], but models for multi-species 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 [9]. 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 resting 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 immunohistochemistry (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 [10] (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 and jamming contrary to the actual evidence. 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-cells response. In a second step, we enrich the model by considering the spatio-temporal evolution of tumor growth through a continuum mathematical variational approach for cell mixtures based on Onsager’s variational principle. Finite element-based simulations in two-dimensional space demonstrate the ambiguous activity of cancer-associated fibroblasts (CAFs) in regulating cancer cell proliferation and invasion (Section III). In the Supplementary 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.

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 NSCLC tumors in situ [10].

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 ECM and growth factors (cytokine) production and upregulation of activation markers such as FAP [11]. 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 [12]. 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 [13, 14, 15, 16]. 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 [17, 10].

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 [18, 19]. 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 [20, 21]. In addition, biochemical factors expressed by CAFs also help to modify the phenotype of T-cells or inactivate their cytotoxic capacity [22, 23]. 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 [24]. In particular, in the family of inflammatory proteins, also called cytokines, Interlukin-6 (IL6) seems to stimulate the infiltration of CD8⁺ [25]. The origins of CD8⁺ T-cells that infiltrate the tumor are diverse, as they differentiate from both circulating and tissue-resident precursors [26]. 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 [27]. The absence of such chemicals leads to the formation of the so-called immune-desert tumor (Fig. 1A).

I.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 [39], fibroblasts ∼ 16 µm [40], and T-cells between 5-10 µm [41], resulting in volumes and projected areas that can be 10 times smaller for T-cells compared to fibroblasts and cancer cells. 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 [28, 29]. 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 [10]. 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⁺ [10]. 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 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. 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 unpredictable and likely 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.

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 S = C + T + F_NA + F_A = 1 − N, where the mass fraction of cancer cells is represented by C, T-cells by T, quiescent or non-activated fibroblasts NAFs by F_NA, activated fibroblasts CAFs by F_A and normal healthy cells associated with the extracellular medium by N . Note that N is not an active component and therefore does not appear in the following equations. Also, we do not do not model 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: dX/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 δ_TF . With these assumptions, the dynamics of the cancer cell population can be read:

Although proliferation of cytotoxic T-cells has been observed, we do not consider proliferation in our study as we focus on their ability to infiltrate the tumor. 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 α_NA whose role is to maintain an adequate supply given by α_NA/δ_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 δ_NA. NAFs are attracted to the tumor nest at a rate α_NA,CC 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 = δ_NA > 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_NA, F_A. More details on the calculations can be found in the supplementary material.

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_NA = 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 {C₊, T₊} = {ϵ/((a₀ −1)δ), ϵ/δ} where a₀ > 1 (see supplementary material, section 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 δ_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 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 + δ_TF F_A)⁻¹ and Δ_F = δF_A. There are four solution pairs (see supplementary material, section 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 + δ_TF F_A)⁻¹ ∼ ϵ, when activated fibroblasts are abundant F_A = δ⁻¹f_a. This results in δ_TF = d₀δϵ⁻¹.

The only stable equilibrium solution is (C₊, T₊) = (f₀δ⁻¹ + 𝒪 (ϵ), 𝒪 (ϵ)), where the notation 𝒪 (ϵ) denotes a quantity whose order of magnitude is ϵ. Thus and not-surprisingly,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_NA is the equilibrium solution of the equation dF_NA/dt = 0 = α_NA − δF_NA². The parameter α_NA represents the net influx of fibroblasts into the tissue. Assuming the fraction of non-activated fibroblasts to be residual ∼ ϵ², we deduce α_NA ∼ δϵ⁴.

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_NA 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_NA = f_nδ⁻¹ where f_n ∼ 1 is a constant.

Replacing δ_CT = δϵ⁻¹, α_TC = a₀ϵ, δ_TF = 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 supplementary material, section 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_NA, F_A} = δ⁻¹{c₀, t₀ϵ, f_NA, f_A}. The order of magnitude of the attraction parameter is then α_NA,C ∼ 1 (see supplementary material, section 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, given 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 in the phenomena studied. This category includes δ⁻¹, the inverse of the free tumor cell density, and α_NA, the NAF attraction parameter to healthy tissue. The other parameters can be studied as control parameters.

Table 1:

Table 2:

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 α_NA,C + α_{T C} .

T-cells efficiency can be quantified through the parameter α_TCδ_CT /(α_NA,CKδ_TF) (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. To conclude, it 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[42, 43, 44, 45, 46, 47]. Similarly, the localization of immune cells in the tumor environment [9], 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 [48, 49, 50].

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 heterogeneity 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 [46, 51], 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 [52]. 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_NA introduced in the previous sections and ϕ₀ = 1 −∑_i≠0 ϕ_i represents N . Γ_i is the growth rate of each component, which is now space and time dependent (see Eqs. (1)-(4)). 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 [53, 54]. 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:

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 Derivation

Since the dynamics is very slow and completely controlled by dissipation, the Onsager variational principle of least dissipation [55] yields the set of partial differential equations coupling the densities to the velocities in the mixture [51]. This principle, introduced by Lord Rayleigh and further developed by Onsager [56, 57, 58, 59], is widely used in soft matter [60], in the biophysical context [51, 46, 61, 62, 55] as well as other areas of physics [63]. First we define a free energy ℱ from a free energy density F, the associated chemical potentials µ_i, and a dissipation function 𝒲:

where the ξ_ij 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 W and the rate of change of the free energy function F. 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_ij is the friction matrix and µ_i are the chemical potentials. Defining: the friction matrix reads:

Importantly, in our model we assume that all frictions ξ_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: ξ_ij = ξ₀, ξ_iCAF = ξ₀ + ξ₁ϕ_CAF .

After deriving the momentum equations, we construct the free energy density. Following [51] 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 [64, 65]:

The terms D_iϕ_i log(ϕ_i) control both the diffusion and volume exclusion, and the quadratic expansion −α_ijϕ_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_NA, 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_NA) by the chemicals c. We keep the growth rate of cancer cells Γ_C given by Eq. 1 with δ_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, their source rate is driven by an incoming flux due to the chemical c. For the NAFs, the volume term Γ_NA|_Ω does not include the source term α_NA,C (see Eq. 3) which is provided by the boundaries, and α_NA is neglected:

On the other hand, the production of CAFs is not changed and Γ_A defined in Eq. 4 is still relevant.

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 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, κ [51]. 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_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 δ_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 [66, 67, 68]. 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 [69].

Table 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 𝒩 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 mutually interact) their mutual interactions. 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 Fig. 5A). 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. 5C). 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 (see Supplementary information).

Figure 5:

IV. Discussion

Our work provides a physical model to quantify the role of the immune system in human lung tumors. 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 proposed 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 microenvironment 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 established the spatial organization of the different cell types inside and outside the tumor core. The model was 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 considered them through the friction increase. An extension of this modeling would be to include anisotropic friction due to a nematic field, as done in Ref. [70], where the nematic fields, the fibroblasts and the fibronectin cause anisotropic friction.

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.1 we used both this latter strategy and the spatial structure of the mixing. In addition, the different combinations of the parameters’ 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.1 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 predicting the outcome of therapy [54]. Indeed, both pharmacodynamics and pharmacokinetics at the tumor site are highly dependent on the spatial structure and composition of the tumor. 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 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 to extend our model to follow tumor evolution under treatment and the related internal dynamics starting from the early alveolar stage of the lung by including the actual three-dimensional environment would be a major step forward as well.

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.

Figure. Sup. 1:

Figure. Sup. 2:

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

Supplementary Table 1:

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

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.

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.

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.

Declaration of interests

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

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). C.B acknowledges the Italian Ministry of Universities and Research (MUR) for the project DICAM-EXC, University of Trento, Departments of Excellence 2023-2027 (grant DM 230/2022). We thank Philippe Benaroch, Layla Mathieson, Xudong Xing for useful discussions.