# Mathematical modelling of the contribution of senescent fibroblasts to basement membrane digestion during carcinoma invasion

Luís Almeida, Alexandre Poulain, Albin Dr Pourtier, Chiara Villa

► To cite this version:

Luís Almeida, Alexandre Poulain, Albin Dr Pourtier, Chiara Villa. Mathematical modelling of the contribution of senescent fibroblasts to basement membrane digestion during carcinoma invasion. 2024. [hal-04574340](https://hal.science/hal-04574340)

HAL Id: hal-04574340

https://hal.science/hal-04574340

Preprint submitted on 14 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
# Mathematical modelling of the contribution of senescent fibroblasts to basement membrane digestion during carcinoma invasion

Luís Almeida, Alexandre Poulain, Albin Pourtier, Chiara Villa

► To cite this version:

Luís Almeida, Alexandre Poulain, Albin Pourtier, Chiara Villa. Mathematical modelling of the contribution of senescent fibroblasts to basement membrane digestion during carcinoma invasion. 2025. hal-04574340v4

HAL Id: hal-04574340

https://hal.science/hal-04574340v4

Preprint submitted on 23 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
# Mathematical modelling of the contribution of senescent fibroblasts to basement membrane digestion during carcinoma invasion

Luís Almeida¹, Alexandre Poulain²,*, Albin Pourtier³, Chiara Villa⁴,*,

1 Sorbonne Université, CNRS, Université de Paris, Laboratoire Jacques-Louis Lions UMR 7598, 75005 Paris, France.

2 Université de Lille, CNRS, UMR 8524 Laboratoire Paul Painlevé, F-59000 Lille, France. Email: alexandre.poulain@univ-lille.fr.

3 Univ. Lille, CNRS, Inserm, CHU Lille, Institut Pasteur Lille, UMR9020-U1277 -CANTHER -Cancer Heterogeneity Plasticity and Resistance to Therapies, F-59000 Lille, France.

4 Sorbonne Université, CNRS, Université de Paris, Inria, Laboratoire Jacques-Louis Lions UMR 7598, 75005 Paris, France. Email: chiara.villa.1@sorbonne-universite.fr.

* Corresponding authors.

## Abstract

Senescent cells have been recognized to play major roles in tumor progression and are nowadays included in the hallmarks of cancer [46]. Our work aims to develop a mathematical model capable of capturing a pro-invasion effect of senescent fibroblasts located in the conjunctive tissue. We focus in the present article on the first moments of the invasion cascade. Considering a localized epithelial tumor, we model the digestion of the collagen fibers of the basement membrane by the proteolytic enzyme MMP-2. The activation of MMP-2 is modelled in detail, as MT1-MMPs bound to the surface of tumor cells interact with proMMPs and TIMPs, proteins enriched in the secretome of senescent Cancer-Associated Fibroblasts, along with its inhibition by TIMPs. Using numerical simulations of the model, calibrated via an extensive literature search, reproducing biologically relevant scenarios, we test the model's suitability to investigate the effect on basement membrane digestion of fibroblasts presenting a senescence-associated secretory phenotype. Via model reduction, steady state and global sensitivity analyses, we identify the most influential parameters in view of their calibration with empirical data. We conclude the paper discussing mathematical and interdisciplinary perspectives.

# 1 Introduction

## 1.1 Biological context

### 1.1.1 Carcinoma progression

Cancer is a disease characterized by a disorder of the regulatory processes, i.e. a disease of cellular homeostasis. Despite the tremendous progress made in the past decades to understand tumor progression and find efficient treatment strategies, cancer persists as the first cause of death in many countries. This is partly due to the extreme complexity of this disease, which involves numerous actors and interlaced biological and physical processes. In particular, other cells (e.g. endothelial cells, fibroblasts, immune cells) in the tumor micro-environment (TME) can act alternatively in favour and against tumor advancement, depending on their activated or altered status along cancer etiology or progression [47, 48].

Carcinomas are tumors that originate in an epithelium (such as lung, breast, colorectum, prostate, stomach, liver, thyroid, cervix, bladder, ...) prior to their progression toward the invasion of their cellular and matrix micro-environment. Carcinomas account for the major part of cancer incidence worldwide (about 20 millions cases) with a mean mortality rate of about 50% worldwide, although the latter varies a lot depending on the country/region [37].

The risk of being diagnosed for a carcinoma is related to ageing [12, 110]. Among the various carcinomas, non-melanoma skin cancers such as basal cell carcinomas or squamous cell carcinomas are amongst those whose incidence will be the most dynamic in the next decades [56]. Hence, although the rate of mortality is moderate in most of these cancers, their massive incidence is already a critical issue.

Typical in solid tumors, carcinoma evolution can be broken down into a series of steps common to rather most types. Tumor initiation occurs from an accumulation of mutations in a single cell, or a few cells, resulting in abnormal proliferation. After a commonly avascular growth, tumors generally stabilise at a diameter of 0.1–0.2cm [53], from lack of nutrients and/or space. Regulation by the environment may
<figure><img src="image_2.png" /></figure>

Peer-reviewed &
recommended by

PCI

Math &
Comp Biol

# Mathematical modelling of the contribution of senescent fibroblasts to basement membrane digestion during carcinoma invasion

Luís Almeida¹, Alexandre Poulain²,¶, Albin Pourtier³,*, Chiara Villa⁴,⁵,¶,*,

1 Sorbonne Université, CNRS, Université de Paris, Laboratoire Jacques-Louis Lions UMR 7598, 75005 Paris, France.

2 Univ. Lille, CNRS, UMR 8524 - Laboratoire Paul Painlevé, F-59000 Lille, France.

3 Univ. Lille, CNRS, Inserm, CHU Lille, Institut Pasteur Lille, UMR9020-U1277 -CANTHER -Cancer Heterogeneity Plasticity and Resistance to Therapies, F-59000 Lille, France.

4 Sorbonne Université, CNRS, Université de Paris, Inria, Laboratoire Jacques-Louis Lions UMR 7598, 75005 Paris, France.

5 Current affiliation: Université Paris-Saclay, Inria, Centre Inria de Saclay, 91120, Palaiseau, France.

* Corresponding authors. Emails: albin.pourtier@univ-lille.fr, chiara.villa@inria.fr

¶ These authors contributed equally to this work.

## Abstract

Senescent cells have been recognized to play major roles in tumor progression and are nowadays included in the hallmarks of cancer [26]. Our work aims to develop a mathematical model capable of capturing a pro-invasion effect of senescent fibroblasts located in the conjunctive tissue. We focus in the present article on the first moments of the invasion cascade. Considering a localized epithelial tumor, we model the digestion of the collagen fibers of the basement membrane by the proteolytic enzyme MMP-2. The activation of MMP-2 is modelled in detail, as MT1-MMPs bound to the surface of tumor cells interact with proMMPs and TIMPs, proteins enriched in the secretome of senescent Cancer-Associated Fibroblasts, along with its inhibition by TIMPs. Using numerical simulations of the model, calibrated via an extensive literature search, reproducing biologically relevant scenarios, we test the models suitability to investigate the effect on basement membrane digestion of fibroblasts presenting a senescence-associated secretory phenotype. Via model reduction, steady state and global sensitivity analyses, we identify the most influential parameters in view of their calibration with empirical data. We conclude the paper discussing mathematical and interdisciplinary perspectives.

## Key points

* We develop a mathematical model to test the hypothesized role of senescent fibroblasts located in the conjunctive tissue on the digestion of the basement membrane by tumor cells.

* Numerical simulations suggest a faster degradation of the basement membrane closer to the senescent fibroblasts, revealing the importance of senescent fibroblasts' proximity to the tumor for the timeline of invasion.

* Sensitivity analysis reveals the most influential parameters for basement membrane digestion, i.e.: the removal or decay rates of MT1-MMP dimers and TIMP-2, the activation rate of MMP-2, the rate of MT1-MMP monomer expression by cancer cells, and the tumor cell density.

* Model reduction results in an analytically tractable model with fewer parameters, which is nonetheless more sensitive to parameter uncertainties – as revealed by sensitivity analysis – and will require extra care during experimental calibration.

## 1 Introduction

**Biological context.** Carcinomas are tumors that originate in an epithelium (such as lung, breast, colorectum, prostate, stomach, liver, thyroid, cervix, bladder, ...) prior to their progression toward the invasion of their cellular and matrix micro-environment. Carcinomas account for the major part of cancer incidence worldwide (about 20 millions cases) with a mean mortality rate of about 50% worldwide,
also occur. At this stage, the carcinoma remains *in situ*, meaning that the epithelial cells are still in their original compartment, touching their underlying basal lamina, physical and anchoring border with the neighbouring conjunctive tissue [93, 88]. Carcinoma progression presents then the basal lamina rupture as a critical point prior to invasion, *i.e.* progression within the TME. The invasion of cancer cells within the surrounding tissue is associated with the development of malignancy, the acquisition of aggressive traits by the cells, and metastasis. Although cancer risk is commonly attributed to intrinsic or inherited mutations, the involvement of age-related perturbations within the TME in the increase of the cancer risk or its progression is clear [71, 110]. In the present article, we are interested in the contributions of senescent stromal cell (fibroblasts) secretions, and particularly MMPs, in the progression of carcinomas in the initial stages of cancer invasion via the digestion of the basement membrane.

### 1.1.2 Basement membrane digestion

During invasion, tumor cells need to migrate within the extracellular matrix (ECM) – *i.e.* the network consisting of extracellular macromolecules providing cells with structural support – in the conjunctive tissue (CT), after they have successfully degraded the basement membrane (BM) – *i.e.* a thin layer of dense ECM separating the tumor from the CT. The ECM is made of a network of fibers, with collagen being the most abundant structural protein (amongst others such as elastin, laminin, fibronectin, proteoglycans, glycoproteins, ...), and is crucial for the homeostasis and the regulation of the cell migration, growth and differentiation. Its cellular and molecular composition varies along ageing [39, 70].

Cancer cells can break through dense barriers of ECM using proteolytic enzymes that remodel the ECM’s fibers network and expose new routes of invasion. The matrix metalloproteinases (MMPs) class is a large family of proteins, secreted in nonactive forms (proMMPs) and become proteolytically active in the extracellular medium under the action of other proteinases such as Membrane-Type Metalloproteinases (MT-MMPs), but also serine proteinases, furin, or other MMPs. In fact, MMPs also have cross-activating effects on diffusing or soluble MMPs: MMP-2, for instance, may be activated both by MT1-MMP and interstitial MMP-1 [72]. MMP activation occurs after successive proteolytic cleavages within the pro-domain and a bait domain leading to a cystein-zinc ion switch that will drive the intra-molecular remodeling of the MMP [108]. The various forms of the MMPs are sensitive to their binding to proteins known as Tissue-Inhibitors of MetalloProteinases (TIMPs). The main function of TIMPs is to act as MMP inhibitors, playing a crucial role in regulating proteolytic degradation in a variety of physiological processes, hence ensuring homeostasis [106]. In cancer, however, the physiological ratio between TIMPs and proMMPs may be disrupted, favoring tumor invasion. TIMPs are commonly found in the TME [27], and produced by stromal cells such as fibroblasts (see e.g. [14]).

Being interested in BM digestion, in this paper we focus on MMP-2, which is known to digest gelatin and type IV collagen, this latter being one of the main components of the BM [104]. MMP-2 may be overexpressed both by tumor and stromal cells, such as fibroblasts, and it is generally produced in its inactive form. The activation depends on the activity of MT1-MMP that relies on the stabilizing effect of TIMP-2. Indeed, after a stabilized complex composed of two MT1-MMP and one TIMP is formed, pro-MMP-2 is presented to this latter by integrin, resulting in the cleavage of a part of it and, ultimately, the activation of the MMP-2. This effect is therefore dependent on the amount of TIMP-2 in a dose-dependent manner, as an excess of this inhibitory enzyme will disorganize the dimeric MT1-MMP complex and inhibit the interaction with the proMMP-2 protein and this latter activation [54, 72, 85]. Membrane bound MT1-MMP is enriched at tumor invadipodia, *i.e.* specific membrane protrusions where tumor cells invade their TME. MT1-MMP membrane expression shows a turn-over which has been indicated as a necessary function to control proMMP-2 activation and increase collagen digestion [109]. Clear gelatinolytic activity was described in the stromal compartment proximal to damaged BMs in human skin biopsies, from aged and young donors, in association with the detection of cells displaying markers of cancer initiation [81]. The basal membrane in aged tissues is thinner in association with profound modifications of molecular and cellular components in the underlying microenvironment contributing to epithelial dysfunctions with potential promotion of cancer [39]. Among the stromal cells in the TME, senescent fibroblasts (SFs) have been observed to secrete pro- and active MMP1 and MMP-2 [81].

### 1.1.3 The role of senescent fibroblasts

Senescence is a cellular state characterized by an arrest of proliferation, resistant to apoptosis, resulting from a panel of replicative, oncogenic, oxydative, DNA-impacting, metabolic stresses [11, 76]. Senescent cells accumulate with age in various tissues as in the skin [28]. Among many others, cancer is an age-related disease [12, 110]. In association with the cell-cycle arrest, inducing senescence of cancer cells have
although the latter varies a lot depending on the country/region [23]. Typical in solid tumors, carcinoma evolution can be broken down into a series of steps common to rather most types. During the first stages of the disease, carcinomas are denoted as *in situ*, meaning that the epithelial cells are still in their original compartment, touching their underlying basal lamina, physical and anchoring border with the neighboring conjunctive tissue (CT) [55, 50]. Carcinoma progression presents then the basal lamina rupture as a critical point prior to invasion, *i.e.* progression within the tumor microenvironment (TME). This step consists in the degradation by cancer cells of the basement membrane (BM), the thin layer separating the tumor from the CT made of dense extra-cellular matrix (ECM), i.e. of a dense network of fibers with collagen being the most abundant. Cancer cells can break through dense barriers of the ECM using proteolytic enzymes that remodel the ECM's fibers network and expose new routes of invasion.

The matrix metalloproteinases (MMPs) class is a large family of proteins, secreted in nonactive forms (proMMPs) in the TME. These latter become proteolytically active in the extracellular medium under the action of other proteinases such as Membrane-Type Metalloproteinases (MT-MMPs), but also serine proteinases, furin, or other MMPs. In particular, MMP-2 is known to digest gelatin and type IV collagen, this latter being one of the main components of the BM [63], and the activation of proMMP-2 is known to be regulated by membrane bound MT1-MMP, typically expressed on the surface of invading cancer cells [66]. The various forms of the MMPs are sensitive to their binding to proteins known as Tissue-Inhibitors of MetalloProteinases (TIMPs), which mainly act as MMP inhibitors but are also known to be involved in the chemical reactions responsible for the activation of certain MMPs [65]. TIMPs are commonly found in the TME [15], and produced by stromal cells such as fibroblasts (see *e.g.* [5]). Indeed, depending on their activated or altered status, various cell types in the TME may account for cancer progression [27]. We here focus on that of senescent cells – which are known to accumulate with age in various tissues as in the skin [16] – given that cancer, among many others, is an age-related disease [4]. As highlighted in Hanahan's latest presentation of the hallmarks of cancer [26], senescent cells are nowadays recognized as important contributors to cancer progression via the release a huge amount of proteins – *e.g.* growth factors, chemokines, MMPs – in the TME [22, 37, 45]. Specifically, among the stromal cells in the TME, senescent fibroblasts (SFs) have been observed to secrete pro- and active MMP-1 and MMP-2 [45]. Despite accumulating evidence on the impact of SFs in carcinoma progression, much more needs to be achieved in order to identify the net of cellular contributions and adapt potential treatment strategies targeting SFs to prevent, or slow down, cancer invasion. We seek to develop a theoretical platform, initially fed by data from literature and biological data from patient samples or biological models, in order to provide a tool to highlight and manipulate critical actors of the first steps of the age-related cancer risk progression, *i.e.* the digestion of the BM as initial mark of invasion.

**Mathematical state of the art.** Enormous modelling efforts have been made to describe ECM digestion and cancer invasion of the CT, *e.g.* see [20, 59] and references therein. Various works proposed mathematical models to study the pro-invasive effect of stromal cells on progressing cancers, with some PDE [39, 38] or agent-based [28] models capturing the crosstalk between cancer cells and fibroblasts. Among these, the model proposed by Kim *et al.* [38] captured the secretion of proMMPs by fibroblasts, their diffusion across a thin layer of ECM and activation once in contact with cancer cells. While this paved the way to subsequent works including BM digestion by active MMPs [25, 37] with interesting results reproducing the transition from the *in situ* stage to the invasive phase, none of these works capture the complex biochemical interactions that regulate MMP/TIMP activation and inhibition. Attempts have been made to model these latter kinetics alone, see for instance [36, 58, 18, 30], discarding spatial dynamics. Works integrating simplified modeling of the MMP/TIMP activation/inhibition and spatial dynamics can be found in the literature for the invasion of cancer cells in the CT [9, 8, 51, 35], but we here focus on the preceding step of BM rupture. Unlike the ECM in the CT, the BM is a thin barrier that cells need to cross and may be better represented as a $d-1$ surface separating tissues of dimension $d$. A multi-dimensional coupling has been employed in previous studies [56, 25, 44, 38, 39] allowing an efficient representation of the BM without having the need to mesh it.

**Objectives of our work.** The short review we presented showcases high levels of biological complexity involved in SF-mediated carcinoma invasion and a vast mathematical modelling literature capturing many elements of the biological system we seek to model – which we expand on in more detail in Supplementary Material 1 for the interested reader. Recent empirical observations on the role of SF secretome on carcinoma invasion, which includes both proMMP-2 and TIMP-2, calls for further development of these modelling frameworks.

In the present work, we propose a bio-physically consistent mathematical model of BM digestion by
been aimed and reported to regulate tumor growth. Cancer progression should however be considered within its microenvironment network of influences. Indeed, while normal mammary fibroblasts can regulate cancerous features of breast epithelial cells in collagen 3D matrices via their secretions and ECM contraction level [78], senescence of stromal cells contribute to tumor progression via a specific secretome the Senescence-Associated Secretory Phenotype, SASP. Therefore, as highlighted in Hanahan's latest presentation of the hallmarks of cancer [46], senescent cells are nowadays recognized as important contributors to cancer progression via the release a huge amount of proteins – e.g. growth factors, chemokines, MMPs – in the TME [36, 23, 24, 49, 64, 69, 81]

Senescent fibroblasts as Cancer-Associated Fibroblasts (CAFs) share common features as they display different subsets/impacts depending on their molecular or activation state. Activated CAFs and late senescent fibroblasts display mostly similar secretomes. Hence, the acquisition of a senescence-associated secretory phenotype by SFs turns them equivalent to a subtype of Cancer-Associated Fibroblasts (CAFs) [5] that similarly secrete proMMP-2 on top of many growth factors, cytokines such as the interleukin IL-8 [107], ECM components and other enzymes that affect the TME similarly to the senescent fibroblasts SASP. In line with the antagonistic pleiotropy hypothesis [97, 35], targeting senescent cells is under active therapeutic trials. The removal of senescent cells in vivo was shown to reduce age-related diseases in mice, among which cancer risk [6, 45]. Thus, SFs may promote carcinogenesis and carcinoma progression in aged organisms. The SASP have numerous possible impacts leading to different therapeutic concepts [43, 76]. Early steps of cancer invasions rely on basal membrane disruption. The fibroblast's SASP displays strong MMP levels and activities, and fibroblasts as CAFs contribute to cancer invasion via their MMP-related matrix digestion or via the remodelling of the latter (stiff tracks) [5, 111]. Despite accumulating evidence on the impact of SFs in carcinoma progression, much more needs to be achieved in order to identify the net of cellular contributions in order to adapt potential treatment strategies targeting SFs to prevent or slow down cancer invasion. We seek to develop a theoretical platform, initially fed by data from literature and biological data from patient samples or biological models, in order to provide a tool to highlight and manipulate critical actors of the first steps of the age-related cancer risk progression, i.e. the digestion of the BM as initial mark of invasion.

## 1.2 Mathematical modelling

Mathematical and computational models are nowadays commonly used to study biological systems. They allow to test hypotheses, uncover hidden mechanisms and explore scenarios in silico, with limited cost and time engagement. In cancer research, this ultimately helps to improve our understanding of tumor progression and propose treatment optimisation strategies.

The field of mathematical oncology has undergone enormous growth in the past thirty years (see for instance the review articles [1, 3, 8, 10, 13, 95] and references therein). The integration of mathematical and computational models into precision oncology is only emerging in recent times, particularly through the development of digital twins [52, 112]. This includes, for instance, their application to radiation treatment of non-melanoma skin cancers [1].

Enormous modelling efforts have been made to describe ECM digestion and cancer invasion of the CT, e.g. see [33, 99] and references therein. Models comprising systems of partial differential equations (PDEs), which allow for a mechanistic description of the invasion process at the tissue-scale, comprise a particularly prolific branch of this field following the first model proposed by Gatenby and Gawlinski [42]. However, while many models have been proposed for different types of movements in the CT, fewer are the works focussing on the role of stroma cells in boosting cancer invasion, on the biochemical reactions responsible for activation of proteolytic enzymes and on the digestion of the BM.

### 1.2.1 Models of cancer invasion promoted by stroma cells

Among the many stroma cells that play a role in cancer progression, particular attention has been given to macrophages and fibroblasts, whose prominent role in tumor spread has lead to dedicated studies on Tumor-Associated Macrophages (TAMs) and CAFs. The focus on these stroma cells and cancer-associated phenotypes has indeed spilled over to the mathematical modelling literature.

Various works proposed PDE and computational models of tumor spread under the cytotoxic action of macrophages, or TAMs, e.g. see [9, 31, 68, 79, 87, 103]. In these works, the crosstalk between tumor cells and macrophages results in the emergence of spatial patterns in the tumor spread. Most models predict these features whilst ignoring ECM degradation, with the emergent behaviour resulting from macrophage infiltration in the tumor, enhanced cell aggregation and cooperative migration. Eftimie and
MMP-2, which is assumed to be produced in its inactive form together with TIMP-2 by SFs in the stroma, and activated by TIMP-2 and MT1-MMPs expressed on the surface of cancer cells. In summary, our model, comprising a system of PDEs and ODEs, combines the following features:

i. secretome spatiotemporal dynamics in the CT, including secretion by SFs modelled as local sources in the stroma;

ii. transmission of SF secretome from the stroma to the BM;

iii. detailed biochemical reactions regulating MMP-2 activation at the BM;

iv. BM digestion by active MMP-2.

Compared to previous works, our model is new on several aspects. First, it comprises a detailed description of the bio-chemical reactions leading to activation and inhibition of MMP-2. Second, it integrates the spatial component by representing the diffusion of the enzymes in the TME. And third, we model the BM as a $d-1$ dimensional manifold (where $d$ is the dimension of the domain modelling the stroma), and let the biochemical reactions regulating MMP activation occur in this domain, aided by cancer cells in the strict vicinity of the BM. It is important to consider that the amount of soluble molecules coming from the CT contributing to the MMP activation dynamics may be lowered as these molecules pass from the CT to the BM, the latter being more dense and thus less permeable. Moreover, the utilization of these molecules in the activation cascade may affect their overall dynamics, e.g. by resulting in a sink at the boundary corresponding to the BM. To ensure these biophysical effects are not neglected, we formally derive the transmission conditions of SF secretome at the BM by performing a dimensionality reduction procedure on the extended model in which the BM is modelled as a $d$-dimensional domain. These phenomena, which were not captured in [10, 25], represent an important element of novelty introduced in our work, together with the comprehensive framework including all of the dynamics listed in points i-iv above. To the authors' knowledge, a mathematical model that captures accurately these dynamics while representing the evolution of the BM density in both time and space does not exist in the literature.

Our model aims to capture accurately the described biological phenomena and to be tractable. To improve analytical and computational tractability, also in view of model calibration with experimental data and model extensions, we derive a reduced model under simplifying assumptions and study this latter analytically. Via numerical simulations and global sensitivity analyses, we explore biological scenarios and compare the sensitivity of the two models, to ensure key elements of the model's predictive capability are not lost during model reduction. Using relevant parameter values carefully drawn from the literature (that rely on dosages and experimental measurements), the numerical simulations of our model illustrate the model's ability to accurately describe the rupture of the BM, and the potential effect of SFs on tumor progression from the *in situ* stage to the invasive one.

In summary, our objectives are to construct a new biologically accurate mathematical model to qualitatively study the role of SFs during the first steps of the invasion of the CT by carcinoma cells, and to understand the benefits – and the disadvantages – of reducing a lengthy but detailed mathematical models into a smaller system that may be preferable for the calibration of the model with data.

**Outline of the paper.** The paper is organized as follows: the model assumptions and equations are presented in Section 2; the reduced model is summarized in Section 2.3 together with selected key analytical results; results of the numerical simulations and global sensitivity analyses are outlined in Section 3; a discussion of the results and research perspectives is given in Section 4.2, particularly in view of the model's extension to a comprehensive framework investigating tumor invasion aided by stroma cells. A list of frequently used acronyms can be consulted at the end of the manuscript (Table 4).

# 2 Methods

## 2.1 The mathematical model

### 2.1.1 Modelling assumptions

We consider a potentially invasive tumour in a healthy non-proliferating tissue, focusing on the stage at which the tumour is isolated and separated from the CT by a dense BM. The BM permeability depends on the local density of ECM, and cancer cells may invade the CT only if the ECM density is sufficiently lowered by proteolytic degradation. In order to study the effect of the position of the SFs and the rate at which they secrete TIMPs and proMMPs on the digestion of the BM, we consider a reduced model of SFs' secretome dynamics in the CT, and ECM degradation at the BM mediated by MMP-2 activation. In particular, the model is based on the following assumptions and simplifications:
coworkers [9, 103] additionally showed that tumor cell motility can significantly increase in the presence of TAMs, which may secrete MMPs involved in the degradation of the ECM.

Other works focused on modelling, using PDEs [66, 65] or agent-based [50] models, the dynamics stemming from the molecular crosstalk between cancer cells and fibroblasts, leading to fibroblast differentiation into myofibroblasts [66, 65] or CAFs [50], which then enhanced cancer cell proliferation and motility. Kim *et al.* [66], in particular, reproduced *in silico* a Tumor Chamber Invasion Assay, in which cancer cells and fibroblasts are plated in two chambers separated by a semi-permeable membrane which allows soluble molecules to pass through, but not cells, and its extension including a layer of ECM on the membrane in [65]. In the latter work, myofibroblasts secrete proMMPs which are assumed to become active and degrade the ECM in combination with cancer cells, an assumption which indirectly models the MMP-activating role of MT1-MMPs on the membrane of cancer cells. This simplified approach to modelling MMP activation was also employed by Gallinato *et al.* [41], who proposed a PDE model for the transition of a tumor from the *in situ* stage to the invasive phase, in which fibroblasts in the CT secreted proMMPs. The proteolytic degradation of the BM by active MMPs led to an increase in the BM permeability over time, ultimately allowing cells to pass through.

The invasion-boosting effect of SFs was captured by the hybrid model presented by Kim *et al.* [64], which reproduced experimental observations on melanoma invasion. In their work, matrix-degrading enzymes are assumed to be secreted by SFs already in their activated form. Nonetheless, as detailed in Section 1.1.2, MMP activation may involve complex interactions and ultimately TIMPs may play an important role in regulating ECM degradation given its dual role of MMP activator and inhibitor [77].

### 1.2.2 Models of TIMPs-regulated proteolytic degradation of the ECM

The first model of MMP-mediated ECM degradation accounting for the complex biochemical reactions responsible for MMP activation was proposed by Karagiannis and Popel [62]. Ignoring spatial dynamics or sources of proMMPs and TIMPs, the model comprises a system of ordinary differential equations (ODEs) and was employed to investigate how varying the initial concentrations of proMMP-2 and TIMP-2 affects the final percentage of collagen digested. The authors found that a high initial TIMP-2 concentration saturates ECM degradation, as expected from its inhibitory function, and that the highest levels of collagen digestions were obtained for intermediate concentrations of TIMP-2, resulting from the key role this protein has in the activation of MMP-2. Analogous results were obtained in [98], who employed a similar model to explore the efficacy of different MMP inhibitors, motivated by the large variety of interactions between complexes formed during MMP-2 activation. The model proposed in [62] was extended by Donze *et al.* [30] to analyse the existence and stability of oscillatory proteolytic dynamics assuming a periodic switch in the role of TIMP-2 as activator or inhibitor of MMP-2. Another extension was proposed by Hoshino *et al.* [55] who included MT1-MMP turnover on the membrane of cancer cells, following experimental results highlighting the importance of MT1-MMP expression on invadopodia for ECM digestion. They found that low turnover rates hindered ECM degradation also due to the key role of MT1-MMP in the MMP-2 activation process, providing a systematic explanation for the experimental observations.

Deakin and Chaplain [17] implemented simplified MMP-2 activation dynamics, mediated by TIMP-2 and MT1-MMP, in a spatially explicit model of cancer invasion. The model, comprising a system of PDEs, built on the previous work of Chaplain and Lolas [18] which considered ECM digestion by cancer cells to be regulated by the plasminogen activation system. The model in [17], which could already predict complex invasion patterns, assumed proMMP-2 to be already attached to the MT1-MMP-TIMP-2 complex and TIMP-2 to be only secreted by cancer cells, and did not include MMP-2 inhibition by TIMP-2. Conversely, in a model of glioma invasion [89] the authors considered the inhibitory effect of TIMP-2 on MMP-2, respectively secreted by healthy and cancer cells in the presence of healthy tissue, although MMP-2 was assumed to be secreted already in its active form. The aim of this work was to investigate how the production of TIMP-2 and that of MMP-2 may balance each other to contain invasion. A similar investigation was conducted by Joshi *et al.* [61], who considered a well-mixed system of stem cells to reproduce *in vitro* experiments to elucidate the mechanisms underlying stem cell-based transplantation therapies.

Despite the variety of models that have been proposed which consider the role of TIMPs in the proteolytic degradation of the ECM, the few spatially explicit ones focus on the invasion of the CT by cancer cells, assuming the digestion of the BM has already occurred. Indeed, we highlight that in [41] the inhibitory effect of TIMPs is considered using a constant decaying rate in the active MMPs concentration equation. Hence, the dynamics of TIMPs concentration and the dual role (activation/inhibition) that TIMP plays is not represented.
* We do not model explicitly the healthy cells but we take into account the ECM remodelling effect of healthy cells near the BM.

* We assume cancer cell proliferation has reached saturation near the BM due to the physical constraints imposed by it. Therefore, we neglect the spatio-temporal dynamics of cancer cells, but do take into account the fact that a higher cell density near the BM results in higher concentrations of molecules involved in MMP-2 activation.

* We neglect the proteolytic action of MT1-MMPs in order to focus on their role in MMP-2 activation. We also neglect the secretion of endogenous proMMP-2 by tumor cells.

* We do not model dynamics of specific integrins (e.g. avβ3 or β1) and their important roles in the activation cascade of MMP-2, but we assume that this effect is captured within specific parameters (i.e. the reaction rates of our system).

We propose a model which translates into mathematical terms the following biological assumptions:

A1. The collagen fibers in the BM can only be decomposed by active MMPs.

A2. The reverse reaction decomposing the active MMP-ECM protein complex into the original enzyme and substrate occur at fast rates, i.e. ECM degradation follows Michaelis-Menten kinetics.

A3. The ECM can be remodelled (i.e. produced) by healthy cells, provided space is available.

A4. MT1-MMPs, which are present on the cancer cell surfaces in monomeric form due to protein expression [7] and recycling of internalised MT1-MMPs [66], can activate proMMP-2 only when in dimeric form.

A5. MT1-MMP dimerization, proMMP-2 activation and TIMPs-mediated MMP-2 inhibition occur according to the biochemical reactions depicted in Figure 1.

A6. MMP and TIMP molecules decay over time (naturally and due to uptake/elimination by other cells).

A7. MT1-MMP and complexes involving MT1-MMP on the surface of cancer cells decay over time due to the MT1-MMP internalization by tumour cells involved in the rapid turnover of this enzyme [30].

A8. SFs secrete both proMMP-2 and TIMPs, increasing the physiological levels of these molecules to pathological levels.

A9. While MT1-MMPs and associated complexes are bound to the cell surfaces, soluble MMPs and TIMPs secreted by SFs are diffusive.

Note that the biochemical reactions summarized in Figure 1 allow us to capture the role TIMP-2 has in the proMMP-2 activation process and the fact that they may hinder the proteolytic action of MMP-2 in three possible ways [43]:

(i) by excessive binding to MT1-MMP, reducing the number of free MT1-MMP molecules available for proMMP-2 fission;

(ii) by directly binding to proMMP-2, thus sequestering it from the activation complex;

(iii) by directly binding to active MMP-2, blocking its proteolytic action.

### 2.1.2 Model equations

Let $t \in \mathbb{R}_{>0}$ denote time and $x = (x_1, x_2) \in \Omega$ denote the spatial position in the two-dimensional domain $\Omega \subset \mathbb{R}^2$, modelling a 2D cross-section of the CT. The boundary of the domain $\Omega$ is denoted $\partial\Omega$. Namely, we assume the domain to be the square $\Omega = [0, L] \times [0, L]$ and that the BM is located at the top boundary $\Gamma = \{x = (x_1, L)\}$ with $x_1 \in [0, L]$, as illustrated in Figure 2.

In this article, we use the following units: seconds (s) for time, decimeter (dm) for space, and nanomolar (nM=nmol/dm³) for molecular concentrations and BM density. Moreover, we use the following notation (see also Figure 1):

* The density of tumor cells in a sufficiently close neighbourhood of the BM is denoted by $\rho^0$ (in number of cell/dm³) and is assumed to be constant for $t \in [0, T]$ where $T > 0$ is a finite time horizon.
### 1.2.3 Models of cancer invasion mediated by the digestion of the basal membrane

Ribba et al. [94] proposed a modelling framework that represents the digestion of the BM and the healthy tissue to investigate the role of therapies targeting MMPs, modelling both the BM and healthy tissue as *d*-dimensional domains. Nonetheless, the BM is very thin compared to the rest of the tissue, measuring between 1 and 1.5 µm in width [80], and modelling frameworks considering the BM to be a *d*−1 dimensional manifold portray more accurately the dimensionality of the invasion process. In the works of Kim and coworkers introduced in Section 1.2.1 [65, 66], such multi-dimensional coupling was employed to model the semi-permeable membrane in the Invasion Chamber Assay they were reproducing, while keeping the ECM in [66] as a *d*-dimensional domain. The thin membrane modelled was permeable only to soluble molecules, which satisfied the transmission conditions derived by Friedman et al. [40] via homogenization, and not cells, which satisfied no-flux boundary conditions.

Such multi-dimensional coupling has been employed in the PDE model of BM invasion by tumor cells proposed by Gallinato et al. [41]. In their model, the authors apply nonlinear Kedem-Katchalsky transmission conditions [63, 67] across the BM, which describe the discontinuity of the pressure across the thin weakly porous membrane as a linear function of the flux. This discontinuity depends on permeability of the membrane, which in [41] is assumed to increase over time under the proteolytic action of MMPs, as previously described in Section 1.2.1. Chaplain et al. [19] presented a formal derivation of biophysically-consistent transmission conditions which, compared to those in [41], allow the cell volume fraction to also be discontinuous across the membrane, while ensuring mass conservation. These conditions, which were later extended to the case of multiple populations in [44] and rigorously derived in [21], can be regarded as a nonlinear generalisation of the classical Kedem-Katchalsky interface conditions [63]. In [19], in particular, they capture BM digestion leading to cancer invasion by letting the effective mobility coefficient depend on the local concentration of MMPs, released by cancer cells in the active form.

The selected literature review we presented showcases the very rich research in tumor invasion modelling, with many elements of the biological system we seek to model. However, recent empirical observations on the role of SF secretome on carcinoma invasion, which includes both proMMP-2 and TIMP-2, calls for further development of these modelling frameworks.

## 1.3 Scope and outline of the paper

In the present work, we propose a biophysically consistent mathematical model of BM digestion by MMP-2, which are produced in their inactive form together with TIMP-2 by SFs in the stroma, and activated by TIMP-2 and MT1-MMPs expressed on the surface of cancer cells. To the authors' knowledge a mathematical model that captures accurately this dynamics while representing the evolution of the BM density in both time and space does not exist in the literature.

The model, comprising a system of PDEs and ODEs, combines the following features:

i. secretome spatiotemporal dynamics in the CT, including secretion by SFs modelled as local sources in the stroma;

ii. transmission of SF secretome from the stroma to the BM;

iii. detailed biochemical reactions regulating MMP-2 activation at the BM;

iv. BM digestion by active MMP-2.

We model the BM as a *d*−1 dimensional manifold (where *d* is the dimension of the domain modelling the stroma), and let the biochemical reactions regulating MMP activation occur in this domain, aided by cancer cells in the strict vicinity of the BM. It is important to consider that the amount of soluble molecules coming from the CT contributing to the MMP activation dynamics may be lowered as these molecules pass from the CT to the BM, the latter being more dense and thus less permeable. Moreover, the utilization of these molecules in the activation cascade may affect their overall dynamics, e.g. by resulting in a sink at the boundary corresponding to the BM. To ensure these biophysical effects are not neglected, we formally derive the transmission conditions of SF secretome at the BM by performing a dimensionality reduction procedure on the extended model in which the BM is modelled as a *d*-dimensional domain. These phenomena, which were not captured in [19, 41], represent an important element of novelty introduced in our work, together with the comprehensive framework including all of the dynamics listed in points i-iv above.

The model aims to capture accurately the biological phenomena described above and to be tractable. To improve analytical and computational tractability, also in view of model calibration with experimental data and model extensions, we derive a reduced model under simplifying assumptions and study this latter analytically. Via numerical simulations and global sensitivity analyses, we explore biological scenarios and compare the sensitivity of the two models, to ensure key elements of the model's predictive capability
<figure><img src="image_3.png" /><figcaption>Figure 1: **Biochemical reactions included in the model.** Visual summary of biochemical reactions between MT1-MMP monomers, MT1-MMP dimers, proMMP-2s, TIMP-2s and active MMP-2s included in the model. MT1-MMPs, and all derived complexes, reside on the membrane of cancer cells, and may be subject to turnover (via expression and internalisation), a simplified version of which is included in our model. The association and dissociation rate $k_i$ depend on specific integrin availability (e.g. beta 3 integrin is responsible for recruiting proMMP-2 [14]) and may *de facto* be subject to temporal fluctuations, but are assumed to be constant in our model for simplicity.</figcaption></figure>

• The density of ECM in the BM, i.e. at $x \in \Gamma$, is denoted by $M(x,t)$ (in nM).

• The concentration of MT1-MMPs in monomeric and dimeric form at $x \in \Gamma$ are denoted, respectively, by $c_m(x,t)$ and $c_d(x,t)$ (in nM).

• The concentration of unbound proMMP-2, TIMP-2 and active MMP-2 at $x \in \Gamma$ are denoted, respectively, by $c_p(x,t)$, $c_t(x,t)$ and $c_a(x,t)$ (in nM).

• The concentration of unbound soluble proMMP-2 and TIMP-2 secreted by SFs in the CT, i.e. $x \in \Omega$, are denoted, respectively, by $\bar{c}_p(x,t)$ and $\bar{c}_t(x,t)$ (in nM).

• The concentration of the different MT1-MMP dimer (‘2MT1’) complexes observed at $x \in \Gamma$ during MMP-2 activation by $c_i(x,t)$ with $i = 1, ..., 3$, i.e. 2MT1-TIMP ($i=1$), 2MT1-TIMP-proMMP ($i=2$) and 2MT1-TIMP-TIMP ($i=3$) (in nM).

<figure><img src="image_4.png" /><figcaption>Figure 2: **Schematic illustration of the epithelium, CT, and the BM.** The red vertical plan represents our two-dimensional domain $\Omega$ that is a cross section of the CT. We also represented that the CT contains a SF (denoted by the star-shaped cell).</figcaption></figure>
are not lost during model reduction. Using relevant parameter values carefully drawn from the literature (that rely on dosages and experimental measurements), the numerical simulations of our model illustrate the model's ability to accurately describe the rupture of the BM, and the potential effect of SFs on tumor progression from the *in situ* stage to the invasive one.

The paper is organized as follows: the model assumptions and equations are presented in Section 2; the reduced model is summarized in Section 3 together with selected key analytical results; results of the numerical simulations and global sensitivity analyses are outlined in Section 4; a discussion of the results and research perspectives is given in Section 5.2, particularly in view of the model's extension to a comprehensive framework investigating tumor invasion aided by stroma cells. A list of frequently used acronyms can be consulted at the end of the manuscript (Table 4).

# 2 The mathematical model

## 2.1 Model assumptions

We consider a potentially invasive tumour in a healthy non-proliferating tissue, focusing on the stage at which the tumour is isolated and separated from the CT by a dense BM. The BM permeability depends on the local density of ECM, and cancer cells may invade the CT only if the ECM density is sufficiently lowered by proteolytic degradation. In order to study the effect of the position of the SFs and the rate at which they secrete TIMPs and proMMPs on the digestion of the BM, we consider a reduced model of SFs' secretome dynamics in the CT, and ECM degradation at the BM mediated by MMP-2 activation. In particular, the model reduction is based on the following assumptions and simplifications:

* We do not model explicitly the healthy cells but we take into account the ECM remodelling effect of healthy cells near the BM.

* We assume cancer cell proliferation has reached saturation near the BM due to the physical constraints imposed by it. Therefore, we neglect the spatio-temporal dynamics of cancer cells, but do take into account the fact that a higher cell density near the BM results in higher concentrations of molecules involved in MMP-2 activation.

* We neglect the proteolytic action of MT1-MMPs in order to focus on their role in MMP-2 activation. We also neglect the secretion of endogenous proMMP-2 by tumor cells.

* We do not model dynamics of specific integrins (e.g. αvβ3 or β1) and their important roles in the activation cascade of MMP-2, but we assume that this effect is captured within specific parameters (i.e. the reaction rates of our system).

We propose a model which translates into mathematical terms the following biological assumptions:

A1. The collagen fibers in the BM can only be decomposed by active MMPs.

A2. The reverse reaction decomposing the active MMP-ECM protein complex into the original enzyme and substrate occur at fast rates, i.e. ECM degradation follows Michaelis-Menten kinetics.

A3. The ECM can be remodelled (i.e. produced) by healthy cells, provided space is available.

A4. MT1-MMPs, which are present on the cancer cell surfaces in monomeric form due to protein expression [16] and recycling of internalised MT1-MMPs [109], can activate proMMP-2 only when in dimeric form.

A5. MT1-MMP dimerization, proMMP-2 activation and TIMPs-mediated MMP-2 inhibition occur according to the biochemical reactions depicted in Figure 1.

A6. MMP and TIMP molecules decay over time (naturally and due to uptake/elimination by other cells).

A7. MT1-MMP and complexes involving MT1-MMP on the surface of cancer cells decay over time due to the MT1-MMP internalization by tumour cells involved in the rapid turnover of this enzyme [55].

A8. SFs secrete both proMMP-2 and TIMPs, increasing the physiological levels of these molecules to pathological levels.

A9. While MT1-MMPs and associated complexes are bound to the cell surfaces, soluble MMPs and TIMPs secreted by SFs are diffusive.
• The concentration of TIMP-bound proMMP-2 and MMP-2 at $x \in \Gamma$ are denoted, respectively, by $c_{tp}(x, t)$ and $c_{ta}(x, t)$ (in nM).

Our model is given by the following ordinary differential equation system at $x \in \Gamma$ for $t \in (0, T)$

$$ \frac{dM}{dt} = -\gamma c_a \frac{M}{K_m + M} + r_M \left(1 - \frac{M}{M_{\max}}\right)_+ , \quad (2.1) $$

$$ \frac{dc_m}{dt} = \alpha_m \rho^0 - k_0 c_m^2 + k_{-0} c_d - \beta_m c_m, \quad (2.2) $$

$$ \frac{dc_d}{dt} = k_0 c_m^2 - k_{-0} c_d - k_1 c_d c_t + k_{-1} c_1 - \beta_d c_d, \quad (2.3) $$

$$ \frac{dc_p}{dt} = I_p(t) - k_2 c_p (c_t + c_1) + k_{-2} (c_{tp} + c_2) - \beta_p c_p, \quad (2.4) $$

$$ \frac{dc_t}{dt} = I_t(t) - c_t[k_1(c_d + c_1) + k_2(c_p + c_a)] + k_{-1}(c_1 + c_3) + k_{-2}(c_{tp} + c_{ta}) - \beta_t c_t, \quad (2.5) $$

$$ \frac{dc_a}{dt} = k_3 c_2 - k_2 c_t c_a + k_{-2} c_{ta} - \beta_a c_a, \quad (2.6) $$

$$ \frac{dc_1}{dt} = k_1 c_d c_t - c_1[k_1 c_t + k_2 c_p] + k_{-1}(c_3 - c_1) + (k_{-2} + k_3)c_2 - \beta_1 c_1, \quad (2.7) $$

$$ \frac{dc_2}{dt} = k_2 c_1 c_p - (k_{-2} + k_3)c_2 - \beta_2 c_2, \quad (2.8) $$

$$ \frac{dc_3}{dt} = k_1 c_1 c_t - k_{-1} c_3 - \beta_3 c_3, \quad (2.9) $$

$$ \frac{dc_{tp}}{dt} = k_2 c_t c_p - k_{-2} c_{tp} - \beta_{tp} c_{tp}, \quad (2.10) $$

$$ \frac{dc_{ta}}{dt} = k_2 c_t c_a - k_{-2} c_{ta} - \beta_{ta} c_{ta}. \quad (2.11) $$

modelling the digestion of the ECM in the BM, coupled – though the terms $I_p(t)$ and $I_t(t)$ in (2.4) and (2.5), defined in (2.14), the initial conditions and boundary conditions (2.18)-(2.19) – with the following PDE system for $(x, t) \in \Omega \times (0, T)$,

$$ \partial_t \bar{c}_p - \bar{D}_p \Delta \bar{c}_p = \bar{S}_p(x) - \beta_p \bar{c}_p, \quad (2.12) $$

$$ \partial_t \bar{c}_t - \bar{D}_t \Delta \bar{c}_t = \bar{S}_t(x) - \beta_t \bar{c}_t, \quad (2.13) $$

modelling the spatio-temporal dynamics of the concentration fields of TIMP-2 and proMMP-2 in the CT.

### 2.1.3 Explanation of the terms appearing the model

**ECM degradation and remodeling in the BM.** The first term on the right-hand-side of equation (2.1) models the proteolytic degradation of ECM proteins by active MMP-2s (A1), modelled according to Michaelis-Menten kinetics (A2), with turnover number $\gamma > 0$ and MichaelisMenten constant $K_m > 0$. Following the modelling strategy adopted in [17], the second term on the right-hand-side of equation (2.1) models ECM remodeling (A3), which occurs at a maximum rate $r_M > 0$ and only if space is available, where $M_{\max} > 0$ is the ECM density carrying capacity at the BM and where we have used the notation $(\cdot)_+ = \max(\cdot, 0)$.

**MMP-2 activation near the BM.** The first term on the right-hand-side of equation (2.2) models the expression of MT1-MMP in monomeric form on the surface of cancer cells (A4), occurring at a rate $\alpha_m > 0$. The terms highlighted in *teal* in equations (2.2)-(2.11) model the biochemical reactions summarised in Figure 1 (A5) according to the law of mass action, with reaction rates $k_i > 0$ ($i = 0, -0, 1, -1, 2, -2, 3$). The last terms on the right-hand-side of equations (2.4)-(2.6), (2.10) and (2.11)
Note that the biochemical reactions summarized in Figure 1 allow us to capture the role TIMP-2 has in the proMMP-2 activation process and the fact that they may hinder the proteolytic action of MMP-2 in three possible ways [77]:

(i) by excessive binding to MT1-MMP, reducing the number of free MT1-MMP molecules available for proMMP-2 fission;

(ii) by directly binding to proMMP-2, thus sequestering it from the activation complex;

(iii) by directly binding to active MMP-2, blocking its proteolytic action.

<figure><img src="image_2.png" /><figcaption>Figure 1: **Biochemical reactions included in the model.** Visual summary of biochemical reactions between MT1-MMP monomers, MT1-MMP dimers, proMMP-2s, TIMP-2s and active MMP-2s included in the model. MT1-MMPs, and all derived complexes, reside on the membrane of cancer cells, and may be subject to turnover (via expression and internalisation), a simplified version of which is included in our model. The association and dissociation rate $k_i$ depend on specific integrin availability (e.g. beta 3 integrin is responsible for recruiting proMMP-2 [26]) and may *de facto* be subject to temporal fluctuations, but are assumed to be constant in our model for simplicity.</figcaption></figure>

## 2.2 Model equations

Let $t \in \mathbb{R}_{\ge 0}$ denote time and $x = (x_1, x_2) \in \Omega$ denote the spatial position in the two-dimensional domain $\Omega \subset \mathbb{R}^2$, modelling a 2D cross-section of the CT. The boundary of the domain $\Omega$ is denoted $\partial\Omega$. Namely, we assume the domain to be the square $\Omega = [0, L] \times [0, L]$ and that the BM is located at the top boundary $\Gamma = \{x = (x_1, L)\}$ with $x_1 \in [0, L]$, as illustrated in Figure 2.

In this article, we use the following units: seconds (s) for time, decimeter (dm) for space, and nanomolar (nM=nmol/dm³) for molecular concentrations and BM density. Moreover, we use the following notation (see also Figure 1):

* The density of tumor cells in a sufficiently close neighbourhood of the BM is denoted by $\rho^0$ (in number of cell/dm³) and is assumed to be constant for $t \in [0, T]$ where $T > 0$ is a finite time horizon.

* The density of ECM in the BM, i.e. at $x \in \Gamma$, is denoted by $M(x, t)$ (in nM);

* The concentration of MT1-MMPs in monomeric and dimeric form at $x \in \Gamma$ are denoted, respectively, by $c_m(x, t)$ and $c_d(x, t)$ (in nM).

* The concentration of unbound proMMP-2, TIMP-2 and active MMP-2 at $x \in \Gamma$ are denoted, respectively, by $c_p(x, t)$, $c_t(x, t)$ and $c_a(x, t)$ (in nM).

* The concentration of unbound soluble proMMP-2 and TIMP-2 secreted by SFs in the CT, i.e. $x \in \Omega$, are denoted, respectively, by $\bar{c}_p(x, t)$ and $\bar{c}_t(x, t)$ (in nM).
model natural decay of $c_i(t, x)$ at a rate $\beta_i > 0$, for $i = \{p, t, a, tp, ta\}$ (A6). The last terms on the right-hand-side of equations (2.2)-(2.3) and (2.7)-(2.9) model internalization of MT1-MMP (A7) at rates $\beta_i > 0$, for $i = \{m, d, 1, 2, 3\}$ (assumed to be equal). Finally, the first terms on the right-hand-side of equations (2.4) and (2.5) model the inflow/outflow of proMMP-2 and TIMP-2 due to environmental sources (A8), and are given by

$$I_p(t) = \hat{\kappa}_p (\bar{c}_p|_{x=0} - c_p) \quad \text{and} \quad I_t(t) = \hat{\kappa}_t (\bar{c}_t|_{x=0} - c_t), \qquad (2.14)$$

modelling the transmission of proMMP-2 and TIMP-2 from the CT to the BM at the right boundary $\Gamma$, with transmission rates $\hat{\kappa}_i > 0$, as derived in the formal model reduction procedure inspired by [41] and detailed in Appendix A.

**SFs' secretome dynamics in the CT.** The second terms on the left-hand-side of equations (2.12) and (2.13) model diffusion of $\bar{c}_p(t, x)$ and $\bar{c}_t(t, x)$ at a rate $D_p > 0$ and $D_t > 0$, respectively (A9). The first terms on the right-hand-side of equations (2.12) and (2.13), i.e. $\bar{S}_p(x), \bar{S}_t(x)$, respectively model the sources of proMMP-2 and TIMP-2 in the CT. They are of the form

$$\bar{S}_p(x) = s_p^{ph} + f_p(x) \quad \text{and} \quad \bar{S}_t(x) = s_t^{ph} + f_t(x), \qquad (2.15)$$

where $s_i^{ph}$ ($i = \{p, t\}$) is a constant source term modelling the physiological production of proMMP-2 and TIMP-2 by cells of the CT, and $f_i(x)$ ($i = \{p, t\}$) model the constant secretion of proMMP-2 and TIMP-2 by SFs. In particular, we take

$$f_p(x) = \sum_{i=1}^{N} r_p \delta_{x_i, \sigma}(x) \quad \text{and} \quad f_t(x) = \sum_{i=1}^{N} r_t \delta_{x_i, \sigma}(x)$$

where we model proMMP-2 and TIMP-2 production at maximum rates $r_p > 0$ and $r_t > 0$, respectively, as sharp gaussians centered at the positions $x_i \in \Omega$ ($i = 1, ..., N$) of the N SFs present in the domain, i.e. we define

$$\delta_{x_i, \sigma}(x) := \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left(-\frac{(x-x_i)^2}{2\sigma^2}\right), \qquad (2.16)$$

for some small variance $0 < \sigma^2 \ll 1$ (we take $\sigma^2 = 0.04 \text{ dm}^2$). Finally, as in equations (2.4) and (2.5), the last terms in (2.12) and (2.13) model natural decay (A6).

**Boundary conditions.** We assume that soluble proMMPs and TIMPs produced by SFs cannot cross the boundaries of the domain, and complement equations (2.12) and (2.13) with the Neumann boundary conditions

$$\bar{D}_p \nabla \bar{c}_p \cdot \nu = \bar{D}_t \nabla \bar{c}_t \cdot \nu = 0 \quad \text{on } \partial\Omega \setminus \Gamma \times (0, T), \qquad (2.17)$$

with $\nu$ the outward normal vector to the boundary $\partial\Omega$. While *in vivo* these soluble molecules may escape the tissue through epithelial membranes, the physical walls of experimental devices justify this assumption *in vitro*. The transmission of proMMP-2 and TIMP-2 from the CT to the BM is captured by Robin boundary conditions at $\Gamma$, as derived in the model reduction procedure detailed in Appendix A, i.e. equations (2.12) and (2.13) are also complemented with (which are compatible with (2.14))

$$-\bar{D}_p \nabla \bar{c}_p \cdot \nu = \kappa_p (\bar{c}_p - c_p), \quad \text{on } \Gamma \times (0, T), \qquad (2.18)$$

$$-\bar{D}_t \nabla \bar{c}_t \cdot \nu = \kappa_t (\bar{c}_t - c_t), \quad \text{on } \Gamma \times (0, T), \qquad (2.19)$$

where $\kappa_p > 0$ and $\kappa_t > 0$ describe the diffusive permeability of proMMP-2 and TIMP-2, respectively, in the BM. While these molecules can cross the BM, we assume that they will immediately interact with cancer cells in the close proximity of the BM. Such dynamics can be adequately captured by the chosen boundary conditions at $\Gamma$.

**Initial conditions.** We consider the initial conditions to reflect the scenario in which no ECM degradation, biochemical reaction, or secretion by SFs has taken place. Thus, for simplicity, we consider the initial ECM density to be at carrying capacity, i.e.

$$M(x, 0) = M_{\max} \quad \text{on } \Gamma. \qquad (2.20)$$
<figure><img src="image_3.png" /><figcaption>Figure 2: Schematic illustration of the epithelium, CT, and the BM. The red vertical plan represents our two-dimensional domain $\Omega$ that is a cross section of the CT. We also represented that the CT contains a SF (denoted by the star-shaped cell).</figcaption></figure>

* The concentration of the different MT1-MMP dimer ('2MT1') complexes observed at $x \in \Gamma$ during MMP-2 activation by $c_i(x, t)$ with $i = 1, ..., 3$, i.e. 2MT1-TIMP ($i = 1$), 2MT1-TIMP-proMMP ($i = 2$) and 2MT1-TIMP-TIMP ($i = 3$) (in nM).

* The concentration of TIMP-bound proMMP-2 and MMP-2 at $x \in \Gamma$ are denoted, respectively, by $c_{tp}(x, t)$ and $c_{ta}(x, t)$ (in nM).

Then, the model is given by the following ordinary differential equation system at $x \in \Gamma$ for $t \in (0, T)$

$$ \frac{dM}{dt} = -\gamma c_a \frac{M}{K_m + M} + r_M \left(1 - \frac{M}{M_{\max}}\right)_+ \quad (2.1) $$

$$ \frac{dc_m}{dt} = \alpha_m \rho^0 - k_0 c_m^2 + k_{-0} c_d - \beta_m c_m \quad (2.2) $$

$$ \frac{dc_d}{dt} = k_0 c_m^2 - k_{-0} c_d - k_1 c_d c_t + k_{-1} c_1 - \beta_d c_d \quad (2.3) $$

$$ \frac{dc_p}{dt} = I_p(t) - k_2 c_p(c_t + c_1) + k_{-2}(c_{tp} + c_2) - \beta_p c_p \quad (2.4) $$

$$ \frac{dc_t}{dt} = I_t(t) - c_t[k_1(c_d + c_1) + k_2(c_p + c_a)] + k_{-1}(c_1 + c_3) + k_{-2}(c_{tp} + c_{ta}) - \beta_t c_t \quad (2.5) $$

$$ \frac{dc_a}{dt} = k_3 c_2 - k_2 c_t c_a + k_{-2} c_{ta} - \beta_a c_a \quad (2.6) $$

$$ \frac{dc_1}{dt} = k_1 c_d c_t - c_1[k_1 c_t + k_2 c_p] + k_{-1}(c_3 - c_1) + (k_{-2} + k_3)c_2 - \beta_1 c_1 \quad (2.7) $$

$$ \frac{dc_2}{dt} = k_2 c_1 c_p - (k_{-2} + k_3)c_2 - \beta_2 c_2 \quad (2.8) $$

$$ \frac{dc_3}{dt} = k_1 c_1 c_t - k_{-1} c_3 - \beta_3 c_3 \quad (2.9) $$

$$ \frac{dc_{tp}}{dt} = k_2 c_t c_p - k_{-2} c_{tp} - \beta_{tp} c_{tp} \quad (2.10) $$

$$ \frac{dc_{ta}}{dt} = k_2 c_t c_a - k_{-2} c_{ta} - \beta_{ta} c_{ta} \quad (2.11) $$

modelling the digestion of the ECM in the BM, coupled – though the terms $I_p(t)$ and $I_t(t)$ in (2.4) and (2.5), defined in (2.14), the initial conditions and boundary conditions (2.18)-(2.19) – with the following PDE
All other soluble molecules, in the CT and at the BM, are taken to be at physiological levels at day 0, i.e.

$$c_i(x, 0) = c_i^{\text{ph}} \quad \text{for} \quad i = \{m, d, a, 1, 2, 3, tp, ta\} \quad \text{on } \Gamma, \qquad (2.21)$$

$$c_p(x, 0) = c_p^{\text{ph}} \quad \text{and} \quad c_p(x, 0) = c_t^{\text{ph}} \quad \text{on } \Gamma, \qquad (2.22)$$

$$\bar{c}_p(x, 0) = c_p^{\text{ph}} \quad \text{and} \quad \bar{c}_t(x, 0) = c_t^{\text{ph}} \quad \text{a.e. in } \bar{\Omega}. \qquad (2.23)$$

## 2.2 Parameter values

The baseline values used for the model parameters are summarized in Table 1. We combine results from experiments and computational works to characterize the maximum number of parameters from the existing literature (see Reference publications in Table 1). Some parameters are estimated here, consistently with our model equations (marked with a * next to the Reference publications in Table 1). In particular:

* the BM carrying capacity is estimated from empirical measurements of collagen IV, the main component of the BM, and its molecular weight;

* internalisation rates of MT1-MMP dimers and complexes are assumed to be the same as those of MT1-MMP monomers;

* decay rates of proMMP-2, TIMP-2 and MMP-2 are computed from half-life values, those of TIMP-2/proMMP-2 and TIMP-2/MMP-2 complexes are estimated from physiological levels assuming they correspond to equilibrium concentrations;

* physiological production rates of proMMP-2 and TIMP-2 are estimated from physiological levels assuming they correspond to equilibrium concentrations;

* physiological levels of MT1-MMPs, proMMP-2, TIMP-2, MMP-2 and TIMP-2/MMP-2 complex are estimated from empirical measurements of serum levels of these proteins and their molecular weight;

* values of the diffusive permeability and transmission rates of proMMP-2 and TIMP-2 in the BM are estimated consistently with the physical considerations and definitions introduced in the formal BM dimension reduction procedure, detailed in appendix A.

Details about the references used and the computations made to obtain the above-listed parameter values can be found in Appendix B.

**Remark 1.** It is worth mentioning that we use measurements in blood and transpose the values to tissue concentrations. This may be untrue and most probably the values depend on the tissue. However, as we were unable to find measurement in tissue samples, we use these values in this work. This can easily be updated whenever better data is available.

Few parameters lack empirical measurements or estimates in the literature (marked with only a * in the last column of Table 1). For simplicity, we set the physiological levels of the remaining molecular complexes to zero. Relevant parameter values for the tumor cell density in contact with BM and the production rates of proMMP-2 and TIMP-2 by SFs are found using numerical testing. In particular, baseline values of TIMP-2 and proMMP-2 production rates by SFs are chosen to add 1% and 200%, respectively, to the respective physiological production rate at the SF location.

## 2.3 Reduced system from simplifying assumptions

To improve analytical and computational tractability, in view of future work comprising model calibration with experimental data and model extensions, we derived a reduced system under simplifying assumptions.

In summary, we assume the inhibitory action of TIMP-2 is irreversible, that the evolution of the complexes' concentrations $c_1, c_2, c_3$ and MT1-MMP turnover are very fast compared to other dynamics in the system, and that natural decay in the BM is negligible in the timescale considered (cf. SA1-SA4 in Appendix C). In the reduced, simplified model we have that $\forall x \in \Gamma$ and $t \in (0, T)$

$$\frac{dM}{dt} = -\gamma c_a \frac{M}{K_M + M} + r_M \left(1 - \frac{M}{M_{\text{max}}}\right), \qquad (2.24)$$

$$\frac{dc_t}{dt} = \hat{\kappa}_t \bar{c}_t - k_2 c_t (c_p + c_a) - \hat{\kappa}_t c_t, \qquad (2.25)$$

$$\frac{dc_p}{dt} = \hat{\kappa}_p \bar{c}_p - k_2 c_p c_t \left(1 + \frac{k_1 k_3}{k_{-1}(k_{-2} + k_3)} \bar{c}_d\right) - \hat{\kappa}_p c_p, \qquad (2.26)$$

$$\frac{dc_a}{dt} = \frac{k_1 k_2 k_3}{k_{-1}(k_{-2} + k_3)} \bar{c}_d c_t c_p - k_2 c_t c_a, \qquad (2.27)$$
system for $(x, t) \in \Omega \times (0, T)$,

$$ \partial_t \bar{c}_p - \bar{D}_p \Delta \bar{c}_p = \bar{S}_p(x) - \beta_p \bar{c}_p, \qquad (2.12) $$

$$ \partial_t \bar{c}_t - \bar{D}_t \Delta \bar{c}_t = \bar{S}_t(x) - \beta_t \bar{c}_t, \qquad (2.13) $$

modelling the spatio-temporal dynamics of the concentration fields of TIMP-2 and proMMP-2 in the CT.

**ECM degradation and remodeling in the BM.** The first term on the right-hand-side of equation (2.1) models the proteolytic degradation of ECM proteins by active MMP-2s (A1), modelled according to Michaelis-Menten kinetics (A2), with turnover number $\gamma > 0$ and Michaelis-Menten constant $K_m > 0$. Following the modelling strategy adopted in [29], the second term on the right-hand-side of equation (2.1) models ECM remodeling (A3), which occurs at a maximum rate $r_M > 0$ and only if space is available, where $M_{max} > 0$ is the ECM density carrying capacity at the BM and where we have used the notation $(\cdot)_+ = \max(\cdot, 0)$.

**MMP-2 activation near the BM.** The first term on the right-hand-side of equation (2.2) models the expression of MT1-MMP in monomeric form on the surface of cancer cells (A4), occurring at a rate $\alpha_m > 0$. The terms highlighted in *teal* in equations (2.2)-(2.11) model the biochemical reactions summarised in Figure 1 (A5) according to the law of mass action, with reaction rates $k_i > 0$ ($i = 0, -0, 1, -1, 2, -2, 3$). The last terms on the right-hand-side of equations (2.4)-(2.6), (2.10) and (2.11) model natural decay of $c_i(t, x)$ at a rate $\beta_i > 0$, for $i = \{p, t, a, tp, ta\}$ (A6). The last terms on the right-hand-side of equations (2.2)-(2.3) and (2.7)-(2.9) model internalization of MT1-MMP (A7) at rates $\beta_i > 0$, for $i = \{m, d, 1, 2, 3\}$ (assumed to be equal, see Table 1). Finally, the first terms on the right-hand-side of equations (2.4) and (2.5) model the inflow/outflow of proMMP-2 and TIMP-2 due to environmental sources (A8), and are given by

$$ I_p(t) = \hat{\kappa}_p (\bar{c}_p|_{x=0} - c_p) \quad \text{and} \quad I_t(t) = \hat{\kappa}_t (\bar{c}_t|_{x=0} - c_t), \qquad (2.14) $$

modelling the transmission of proMMP-2 and TIMP-2 from the CT to the BM at the right boundary $\Gamma$, with transmission rates $\hat{\kappa}_i > 0$, as derived in the formal model reduction procedure inspired by [74] and detailed in Appendix A.

**SFs' secretome dynamics in the CT.** The second terms on the left-hand-side of equations (2.12) and (2.13) model diffusion of $\bar{c}_p(t, x)$ and $\bar{c}_t(t, x)$ at a rate $D_p > 0$ and $D_t > 0$, respectively (A9). The first terms on the right-hand-side of equations (2.12) and (2.13), i.e. $\bar{S}_p(x), \bar{S}_t(x)$, respectively model the sources of proMMP-2 and TIMP-2 in the CT. They are of the form

$$ \bar{S}_p(x) = s_p^{ph} + f_p(x) \quad \text{and} \quad \bar{S}_t(x) = s_t^{ph} + f_t(x), \qquad (2.15) $$

where $s_i^{ph}$ ($i = \{p, t\}$) is a constant source term modelling the physiological production of proMMP-2 and TIMP-2 by cells of the CT, and $f_i(x)$ ($i = \{p, t\}$) model the constant secretion of proMMP-2 and TIMP-2 by SFs. In particular, we take

$$ f_p(x) = \sum_{i=1}^{N} r_p \delta_{x_i, \sigma}(x) \quad \text{and} \quad f_t(x) = \sum_{i=1}^{N} r_t \delta_{x_i, \sigma}(x) $$

where we model proMMP-2 and TIMP-2 production at maximum rates $r_p > 0$ and $r_t > 0$, respectively, as sharp gaussians centered at the positions $x_i \in \Omega$ ($i = 1, ..., N$) of the N SFs present in the domain, i.e. we define

$$ \delta_{x_i, \sigma}(x) := \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left(-\frac{(x-x_i)^2}{2\sigma^2}\right), \qquad (2.16) $$

for some small variance $0 < \sigma^2 \ll 1$ (we take $\sigma^2 = 0.04 \text{ dm}^2$). Finally, as in equations (2.4) and (2.5), the last terms in (2.12) and (2.13) model natural decay (A6).

**Boundary conditions.** We assume that soluble proMMPs and TIMPs produced by SFs cannot cross the boundaries of the domain, and complement equations (2.12) and (2.13) with the Neumann boundary conditions

$$ \bar{D}_p \nabla \bar{c}_p \cdot \nu = \bar{D}_t \nabla \bar{c}_t \cdot \nu = 0 \quad \text{on } \partial\Omega \setminus \Gamma \times (0, T), \qquad (2.17) $$
<table><caption>Table 1: List of parameters found in equations and initial conditions (2.1)-(2.23) and baseline parameter values used in this paper, drawn from the literature. Parameters which are not explicitly found in the literature and were calculated here are marked with a * in the last column. Details about the calculations used to estimate parameters from the referenced papers can be found in Appendix B.</caption><thead><tr><th>Parameter</th><th>Biological meaning</th><th>Value &amp; Unit</th><th>Reference</th></tr></thead><tbody><tr><td>γ</td><td>Maximum ECM degradation rate</td><td>0.236 s<sup>-1</sup></td><td>[35]</td></tr><tr><td>K<sub>m</sub></td><td>Michaelis-Menten constant of ECM proteolytic degradation</td><td>0.1357 nM</td><td>[35]</td></tr><tr><td>r<sub>M</sub></td><td>Maximum ECM remodeling rate</td><td>6.18 × 10<sup>-4</sup> nM/s</td><td>[35]</td></tr><tr><td>M<sub>max</sub></td><td>ECM carrying capacity</td><td>62.5 × 10<sup>3</sup> nM</td><td>* [1]</td></tr><tr><td>α<sub>m</sub></td><td>Production rate of monomeric MT1-MMP by cancer cells</td><td>5 × 10<sup>-4</sup> nmol/s/cell</td><td>[8]</td></tr><tr><td>ρ<sup>0</sup></td><td>Number of tumor cells adjacent to the BM</td><td>1 × 10<sup>2</sup> cells/dm<sup>2</sup></td><td>*</td></tr><tr><td>k<sub>0</sub></td><td>Maximum MT1-MMP dimerization rate</td><td>2 1/nM/s</td><td>[30]</td></tr><tr><td>k<sub>-0</sub></td><td>Dissociation rate of MT1-MMP dimers</td><td>1 × 10<sup>-2</sup> s<sup>-1</sup></td><td>[30]</td></tr><tr><td>k<sub>1</sub></td><td>Association rate of TIMP-2 and MT1-MMP</td><td>2.71 1/nM/s</td><td>[64]</td></tr><tr><td>k<sub>-1</sub></td><td>Dissociation rate of bound MT1-MMP and TIMP-2</td><td>1 × 10<sup>-4</sup> s<sup>-1</sup></td><td>[64]</td></tr><tr><td>k<sub>2</sub></td><td>Association rate of proMMP-2 and TIMP-2</td><td>0.14 1/nM/s</td><td>[49]</td></tr><tr><td>k<sub>-2</sub></td><td>Dissociation rate of bound proMMP-2 and TIMP-2</td><td>1 × 10<sup>-4</sup> s<sup>-1</sup></td><td>[49]</td></tr><tr><td>k<sub>3</sub></td><td>MMP-2 activation rate</td><td>0.02 s<sup>-1</sup></td><td>[36]</td></tr><tr><td>β<sub>m</sub></td><td>MT1-MMP monomer internalization rate</td><td>3.85 × 10<sup>-2</sup></td><td>[30]</td></tr><tr><td>β<sub>d</sub></td><td>MT1-MMP dimer internalization rate</td><td>3.85 × 10<sup>-2</sup></td><td>* [30]</td></tr><tr><td>β<sub>1</sub>, β<sub>2</sub>, β<sub>3</sub></td><td>Internalization rate of MT1-MMP complexes</td><td>3.85 × 10<sup>-2</sup></td><td>* [30]</td></tr><tr><td>β<sub>p</sub></td><td>Decay rate of proMMP-2</td><td>6.08 × 10<sup>-4</sup> s<sup>-1</sup></td><td>* [47]</td></tr><tr><td>β<sub>t</sub></td><td>Decay rate of TIMP-2</td><td>4.45 × 10<sup>-5</sup> s<sup>-1</sup></td><td>* [6]</td></tr><tr><td>β<sub>a</sub></td><td>Decay rate of active MMP-2</td><td>1.22 × 10<sup>-3</sup> s<sup>-1</sup></td><td>* [47]</td></tr><tr><td>β<sub>tp</sub></td><td>Decay rate of TIMP-2/proMMP-2 complex</td><td>3.60 s<sup>-1</sup></td><td>*</td></tr><tr><td>β<sub>ta</sub></td><td>Decay rate of TIMP-2/MMP-2 complex</td><td>0.035 s<sup>-1</sup></td><td>*</td></tr><tr><td>D<sub>p</sub>, D<sub>t</sub></td><td>Diffusion coefficients of proMMP-2, TIMP-2</td><td>1.29 × 10<sup>-6</sup> dm<sup>2</sup>/s</td><td>[8, 12]</td></tr><tr><td>s<sup>ph</sup><sub>p</sub></td><td>Physiological production rate of proMMP-2</td><td>1.09 × 10<sup>-2</sup> nM/s</td><td>* [52, 47]</td></tr><tr><td>s<sup>ph</sup><sub>t</sub></td><td>Physiological production rate of TIMP-2</td><td>4.28 × 10<sup>-4</sup> nM/s</td><td>* [6, 52]</td></tr><tr><td>r<sub>p</sub></td><td>Production rate of proMMP-2 by senescent fibroblasts</td><td>8.34 × 10<sup>-4</sup> nM/s</td><td>*</td></tr><tr><td>r<sub>t</sub></td><td>Production rate of TIMP-2 by senescent fibroblasts</td><td>2.22 × 10<sup>-6</sup> nM/s</td><td>*</td></tr><tr><td>κ<sub>p</sub>, κ<sub>t</sub></td><td>Diffusive permeability of proMMP-2 and TIMP-2 in the BM</td><td>1.935 × 10<sup>-6</sup> dm/s</td><td>* [8, 12, 57]</td></tr><tr><td>κ̃<sub>p</sub>, κ̃<sub>t</sub></td><td>Transmission rate of proMMP-2 and TIMP-2 across the BM</td><td>0.9675 s<sup>-1</sup></td><td>* [8, 12, 57]</td></tr><tr><td>c<sup>ph</sup><sub>m</sub></td><td>Physiological levels of monomeric MT1-MMP</td><td>0.14 nM</td><td>* [40]</td></tr><tr><td>c<sup>ph</sup><sub>d</sub></td><td>Physiological levels of dimeric MT1-MMP</td><td>0.14 nM</td><td>* [11]</td></tr><tr><td>c<sup>ph</sup><sub>p</sub></td><td>Physiological levels of proMMP-2</td><td>18 nM</td><td>* [52, 61]</td></tr><tr><td>c<sup>ph</sup><sub>t</sub></td><td>Physiological levels of TIMP-2</td><td>9.6 nM</td><td>* [52, 62]</td></tr><tr><td>c<sup>ph</sup><sub>a</sub></td><td>Physiological levels of active MMP-2</td><td>0.44 nM</td><td>* [52, 61]</td></tr><tr><td>c<sup>ph</sup><sub>ta</sub></td><td>Physiological levels of TIMP-2/MMP-2 complex</td><td>17 nM</td><td>* [52]</td></tr><tr><td>c<sup>ph</sup><sub>tp</sub></td><td>Physiological levels of TIMP-2/proMMP-2 complex</td><td>6.72 nM</td><td>* [2]</td></tr><tr><td>c<sup>ph</sup><sub>1</sub>, c<sup>ph</sup><sub>2</sub>, c<sup>ph</sup><sub>3</sub></td><td>Physiological levels of MMP-2 activation complexes</td><td>0 nM</td><td>*</td></tr></tbody></table>

in which $\bar{c}_d$ is a given parameter representing the constant concentration of available dimeric MT1-MMP and $\bar{c}_i \equiv \bar{c}_i(x, t)|_\Gamma (i = \{t, p\})$, i.e. the concentrations of TIMP-2 and proMMP-2 at the BM coming from the CT. System (2.24)-(2.27) is coupled to System (2.12)-(2.13).

**Equilibrium concentrations.** We analysed the steady state of the reduced system (2.25)-(2.27), cf. Appendix C.1.1. The corresponding equilibrium concentrations are plotted in Figure 3 for varying values of the concentrations of TIMP-2 and proMMP-2 reaching the BM from the CT. We see that the equilibrium concentration of active MMP-2: increases (linearly) with the amount of proMMP-2 coming from the CT; (linearly) decreases with the amount of TIMP-2 coming from the CT; is zero in the absence of TIMP-2 coming from the CT; increases with the amount of MT1-MMP dimers (cf. Supplementary Figure S1). These results are consistent with the dual role of TIMP-2 in MMP-2 activation.

**BM rupture and cancer invasion.** We define BM rupture as the decrease of the BM density below an arbitrary value that corresponds to the critical density below which tumor cells can invade. We derived a condition for this long-term sustained BM rupture, allowing cells to begin invading the CT, cf. Appendix C.1.2. From this condition, we have that BM rupture and cancer invasion is more likely to occur for: higher values of active MMP-2 at equilibrium, which depends on the amount of proMMP-2 and TIMP-2 coming from the CT and MT1-MMP2 dimers; higher BM degradation rates; lower BM
with $\nu$ the outward normal vector to the boundary $\partial\Omega$. While in vivo these soluble molecules may escape the tissue through epithelial membranes, the physical walls of experimental devices justify this assumption in vitro. The transmission of proMMP-2 and TIMP-2 from the CT to the BM is captured by Robin boundary conditions at $\Gamma$, as derived in the model reduction procedure detailed in Appendix A, i.e. equations (2.12) and (2.13) are also complemented with (which are compatible with (2.14))

$$-\bar{D}_p \nabla \bar{c}_p \cdot \nu = \kappa_p (\bar{c}_p - c_p), \quad \text{on } \Gamma \times (0, T), \quad (2.18)$$

$$-\bar{D}_t \nabla \bar{c}_t \cdot \nu = \kappa_t (\bar{c}_t - c_t), \quad \text{on } \Gamma \times (0, T), \quad (2.19)$$

where $\kappa_p > 0$ and $\kappa_t > 0$ describe the diffusive permeability of proMMP-2 and TIMP-2, respectively, in the BM. While these molecules can cross the BM, we assume that they will immediately interact with cancer cells in the close proximity of the BM. Such dynamics can be adequately captured by the chosen boundary conditions at $\Gamma$.

**Initial conditions.** We consider the initial conditions to reflect the scenario in which no ECM degradation, biochemical reaction, or secretion by SFs has taken place. Thus, for simplicity, we consider the initial ECM density to be at carrying capacity, i.e.

$$M(x, 0) = M_{\max} \quad \text{on } \Gamma. \quad (2.20)$$

All other soluble molecules, in the CT and at the BM, are taken to be at physiological levels at day 0, i.e.

$$c_i(x, 0) = c_i^{\text{ph}} \quad \text{for} \quad i = \{m, d, a, 1, 2, 3, tp, ta\} \quad \text{on } \Gamma, \quad (2.21)$$

$$c_p(x, 0) = c_p^{\text{ph}} \quad \text{and} \quad c_p(x, 0) = c_t^{\text{ph}} \quad \text{on } \Gamma, \quad (2.22)$$

$$\bar{c}_p(x, 0) = c_p^{\text{ph}} \quad \text{and} \quad \bar{c}_t(x, 0) = c_t^{\text{ph}} \quad \text{a.e. in } \bar{\Omega}. \quad (2.23)$$

## 2.3 Parameter values

The baseline values used for the model parameters are summarized in Table 1. We combine results from experiments and computational works to characterize the maximum number of parameters from the existing literature (see Reference publications in Table 1). Some parameters are estimated here, consistently with our model equations (marked with a * next to the Reference publications in Table 1). In particular:

* the BM carrying capacity is estimated from empirical measurements of collagen IV, the main component of the BM, and its molecular weight;

* internalisation rates of MT1-MMP dimers and complexes are assumed to be the same as those of MT1-MMP monomers;

* decay rates of proMMP-2, TIMP-2 and MMP-2 are computed from half-life values, those of TIMP-2/proMMP-2 and TIMP-2/MMP-2 complexes are estimated from physiological levels assuming they correspond to equilibrium concentrations;

* physiological production rates of proMMP-2 and TIMP-2 are estimated from physiological levels assuming they correspond to equilibrium concentrations;

* physiological levels of MT1-MMPs, proMMP-2, TIMP-2, MMP-2 and TIMP-2/MMP-2 complex are estimated from empirical measurements of serum levels of these proteins and their molecular weight;

* values of the diffusive permeability and transmission rates of proMMP-2 and TIMP-2 in the BM are estimated consistently with the physical considerations and definitions introduced in the formal BM dimension reduction procedure, detailed in appendix A.

Details about the references used and the computations made to obtain the above-listed parameter values can be found in Appendix B.

**Remark 1.** It is worth mentioning that we use measurements in blood and transpose the values to tissue concentrations. This may be untrue and most probably the values depend on the tissue. However, as we were unable to find measurement in tissue samples, we use these values in this work. This can easily be updated whenever better data is available.

Few parameters lack empirical measurements or estimates in the literature (marked with only a * in the last column of Table 1). For simplicity, we set the physiological levels of the remaining molecular complexes to zero. Relevant parameter values for the tumor cell density in contact with BM and the production rates of proMMP-2 and TIMP-2 by SFs are found using numerical testing. In particular, baseline values of TIMP-2 and proMMP-2 production rates by SFs are chosen to add 1% and 200%, respectively, to the respective physiological production rate at the SF location.
<figure><img src="image_5.png" /><figcaption>Figure 3: The analytical equilibrium of the reduced system demonstrates that the dual role of TIMP-2, as activator and inhibitor of MMP-2, is still captured in the reduced model. Concentration of TIMP-2 (left), proMMP-1 (center) and active MMP-2 (right) at steady state of the reduced ODE system (2.25)–(2.27), for varying fixed concentrations of proMMP-2 ($\overline{c}_p^*$) and TIMP-2 ($\overline{c}_t^*$) coming from the CT. The steady state is given by (C.3) for $\overline{c}_t^* > 0$, and by (C.4) otherwise, under the parameter values listed in Table 1 and $C_a = 0$. The latter corresponds to the value obtained at equilibrium under initial conditions $c_t = c_p = c_a = 0$ at $t = 0$. We consider TIMP-2 and proMMP-2 concentrations from the CT in the ranges $\overline{c}_t^* \in [0, c_t^{ph} + 1\%c_t^{ph}]$ and $\overline{c}_p^* \in [0, c_p^{ph} + 200\%c_p^{ph}]$, in line with the parameter choices reported in Section 2.2.</figcaption></figure>

remodeling rates. These results are consistent with the biological meaning of these quantities.

## 2.4 Numerical method and set up

To test the model's ability to predict biologically relevant scenarios, we conduct numerical simulations of the full System (2.1)–(2.13), complemented with definitions (2.14)–(2.16), boundary conditions (2.17)–(2.19), and initial conditions (2.20)–(2.23), under the baseline parameter values stated in Table 1. The results are discussed in Section 3.1, while section 3.2 is dedicated to a global sensitivity analysis of both the full and reduced models, using the SAFE (Sensitivity analysis For Everyone) toolbox [54].

In Section 3.1 we set $dim(\Omega) = 2$ and $dim(\Gamma) = 1$ (“2D” case), while in Section 3.2 we set $dim(\Omega) = 1$ and $dim(\Gamma) = 0$ (“1D” case). In particular, in 2D we assume that the domain $\Omega = [0, L] \times [0, L]$ is a square of side length $L$ dm (or, in 1D, a line of length $L$ dm). We consider a uniform discretization of the domain $\Omega$ using $N_x \times N_x$ nodes (or, in 1D, using $N_x$ nodes), and using $N_x$ nodes for $\Gamma$. We construct a cell-centered finite volume spatial discretization of our system and integrate it in time using an explicit Runge-Kutta integration method, up to a final time $T$. For details about our numerical scheme, we refer to Appendix D, and the values of $T$, $L$ and $N_x$ used are reported in Table 2.

<table><caption>Table 2: Simulation parameters used in Section 3.1 and Section 3.2.</caption><thead><tr><th>Parameter</th><th>Meaning</th><th>Unit</th><th>Value S.3.1</th><th>Value S.3.2</th></tr></thead><tbody><tr><td>dim(Ω)</td><td>Dimension of the CT</td><td></td><td>2</td><td>1</td></tr><tr><td>T</td><td>End time</td><td>days</td><td>7</td><td>1</td></tr><tr><td>L</td><td>Length in each dimension of Ω</td><td>dm</td><td>0.1</td><td>0.1</td></tr><tr><td>N<sub>x</sub></td><td>Number of nodes in the x<sub>i</sub> direction</td><td></td><td>50</td><td>20</td></tr></tbody></table>

Our code (with scripts to replicate the figures of this article) can be downloaded from Zenodo https://doi.org/10.5281/zenodo.12654067 or from GitHub https://github.com/alexandrepoulain/TumInvasion-BM.

# 3 Results

## 3.1 Model's ability to reproduce biologically relevant scenarios

### 3.1.1 Biologically relevant test cases

To test the model's ability to reproduce biologically relevant scenarios, we consider three test cases:
<table><caption>Table 1: List of parameters found in equations and initial conditions (2.1)-(2.23) and baseline parameter values used in this paper, drawn from the literature. Parameters which are not explicitly found in the literature and were calculated here are marked with a * in the last column. Details about the calculations used to estimate parameters from the referenced papers can be found in Appendix B.</caption><thead><tr><th>Parameter</th><th>Biological meaning</th><th>Value &amp; Unit</th><th>Reference</th></tr></thead><tbody><tr><td>γ</td><td>Maximum ECM degradation rate</td><td>0.236 s<sup>-1</sup></td><td>[61]</td></tr><tr><td>K<sub>m</sub></td><td>Michaelis-Menten constant of ECM proteolytic degradation</td><td>0.1357 nM</td><td>[61]</td></tr><tr><td>r<sub>M</sub></td><td>Maximum ECM remodeling rate</td><td>6.18 &times; 10<sup>-4</sup> nM/s</td><td>[61]</td></tr><tr><td>M<sub>max</sub></td><td>ECM carrying capacity</td><td>62.5 &times; 10<sup>3</sup> nM</td><td>* [2]</td></tr><tr><td>α<sub>m</sub></td><td>Production rate of monomeric MT1-MMP by cancer cells</td><td>5 &times; 10<sup>-4</sup> nmol/s/cell</td><td>[17]</td></tr><tr><td>ρ<sup>0</sup></td><td>Number of tumor cells adjacent to the BM</td><td>1 &times; 10<sup>2</sup> cells/dm<sup>2</sup></td><td>*</td></tr><tr><td>k<sub>0</sub></td><td>Maximum MT1-MMP dimerization rate</td><td>2 1/nM/s</td><td>[55]</td></tr><tr><td>k<sub>-0</sub></td><td>Dissociation rate of MT1-MMP dimers</td><td>1 &times; 10<sup>-2</sup> s<sup>-1</sup></td><td>[55]</td></tr><tr><td>k<sub>1</sub></td><td>Association rate of TIMP-2 and MT1-MMP</td><td>2.71 1/nM/s</td><td>[105]</td></tr><tr><td>k<sub>-1</sub></td><td>Dissociation rate of bound MT1-MMP and TIMP-2</td><td>1 &times; 10<sup>-4</sup> s<sup>-1</sup></td><td>[105]</td></tr><tr><td>k<sub>2</sub></td><td>Association rate of proMMP-2 and TIMP-2</td><td>0.14 1/nM/s</td><td>[86]</td></tr><tr><td>k<sub>-2</sub></td><td>Dissociation rate of bound proMMP-2 and TIMP-2</td><td>1 &times; 10<sup>-4</sup> s<sup>-1</sup></td><td>[86]</td></tr><tr><td>k<sub>3</sub></td><td>MMP-2 activation rate</td><td>0.02 s<sup>-1</sup></td><td>[62]</td></tr><tr><td>β<sub>m</sub></td><td>MT1-MMP monomer internalization rate</td><td>3.85 &times; 10<sup>-2</sup></td><td>[55]</td></tr><tr><td>β<sub>d</sub></td><td>MT1-MMP dimer internalization rate</td><td>3.85 &times; 10<sup>-2</sup></td><td>* [55]</td></tr><tr><td>β<sub>1</sub>, β<sub>2</sub>, β<sub>3</sub></td><td>Internalization rate of MT1-MMP complexes</td><td>3.85 &times; 10<sup>-2</sup></td><td>* [55]</td></tr><tr><td>β<sub>p</sub></td><td>Decay rate of proMMP-2</td><td>6.08 &times; 10<sup>-4</sup> s<sup>-1</sup></td><td>* [83]</td></tr><tr><td>β<sub>t</sub></td><td>Decay rate of TIMP-2</td><td>4.45 &times; 10<sup>-5</sup> s<sup>-1</sup></td><td>* [15]</td></tr><tr><td>β<sub>a</sub></td><td>Decay rate of active MMP-2</td><td>1.22 &times; 10<sup>-3</sup> s<sup>-1</sup></td><td>* [83]</td></tr><tr><td>β<sub>tp</sub></td><td>Decay rate of TIMP-2/proMMP-2 complex</td><td>3.60 s<sup>-1</sup></td><td>*</td></tr><tr><td>β<sub>ta</sub></td><td>Decay rate of TIMP-2/MMP-2 complex</td><td>0.035 s<sup>-1</sup></td><td>*</td></tr><tr><td>D<sub>p</sub>, D<sub>t</sub></td><td>Diffusion coefficients of proMMP-2, TIMP-2</td><td>1.29 &times; 10<sup>-6</sup> dm<sup>2</sup>/s</td><td>[17, 22]</td></tr><tr><td>s<sup>ph</sup><sub>p</sub></td><td>Physiological production rate of proMMP-2</td><td>1.09 &times; 10<sup>-2</sup> nM/s</td><td>* [90, 83]</td></tr><tr><td>s<sup>ph</sup><sub>t</sub></td><td>Physiological production rate of TIMP-2</td><td>4.28 &times; 10<sup>-4</sup> nM/s</td><td>* [15, 90]</td></tr><tr><td>r<sub>p</sub></td><td>Production rate of proMMP-2 by senescent fibroblasts</td><td>8.34 &times; 10<sup>-4</sup> nM/s</td><td>*</td></tr><tr><td>r<sub>t</sub></td><td>Production rate of TIMP-2 by senescent fibroblasts</td><td>2.22 &times; 10<sup>-6</sup> nM/s</td><td>*</td></tr><tr><td>κ<sub>p</sub>, κ<sub>t</sub></td><td>Diffusive permeability of proMMP-2 and TIMP-2 in the BM</td><td>1.935 &times; 10<sup>-6</sup> dm/s</td><td>* [17, 22, 96]</td></tr><tr><td>κ̂<sub>p</sub>, κ̂<sub>t</sub></td><td>Transmission rate of proMMP-2 and TIMP-2 across the BM</td><td>0.9675 s<sup>-1</sup></td><td>* [17, 22, 96]</td></tr><tr><td>c<sup>ph</sup><sub>m</sub></td><td>Physiological levels of monomeric MT1-MMP</td><td>0.14 nM</td><td>* [73]</td></tr><tr><td>c<sup>ph</sup><sub>d</sub></td><td>Physiological levels of dimeric MT1-MMP</td><td>0.14 nM</td><td>* [20]</td></tr><tr><td>c<sup>ph</sup><sub>p</sub></td><td>Physiological levels of proMMP-2</td><td>18 nM</td><td>* [90, 101]</td></tr><tr><td>c<sup>ph</sup><sub>t</sub></td><td>Physiological levels of TIMP-2</td><td>9.6 nM</td><td>* [90, 102]</td></tr><tr><td>c<sup>ph</sup><sub>a</sub></td><td>Physiological levels of active MMP-2</td><td>0.44 nM</td><td>* [90, 101]</td></tr><tr><td>c<sup>ph</sup><sub>ta</sub></td><td>Physiological levels of TIMP-2/MMP-2 complex</td><td>17 nM</td><td>* [90]</td></tr><tr><td>c<sup>ph</sup><sub>tp</sub></td><td>Physiological levels of TIMP-2/proMMP-2 complex</td><td>6.72 nM</td><td>* [4]</td></tr><tr><td>c<sup>ph</sup><sub>1</sub>, c<sup>ph</sup><sub>2</sub>, c<sup>ph</sup><sub>3</sub></td><td>Physiological levels of MMP-2 activation complexes</td><td>0 nM</td><td>*</td></tr></tbody></table>

## 3 Reduced system from simplifying assumptions

To improve analytical and computational tractability, in view of future work comprising model calibration with experimental data and model extensions, we derived a reduced system under simplifying assumptions.

In summary, we assume the inhibitory action of TIMP-2 is irreversible, that the evolution of the complexes' concentrations $c_1, c_2, c_3$ and MT1-MMP turnover are very fast compared to other dynamics in the system, and that natural decay in the BM is negligible in the timescale considered (cf. SA1-SA4 in Appendix C). In the reduced, simplified model we have that $\forall x \in \Gamma$ and $t \in (0, T)$

$$ \frac{dM}{dt} = -\gamma c_a \frac{M}{K_M + M} + r_M \left(1 - \frac{M}{M_{\max}}\right), \quad (3.1) $$

$$ \frac{dc_t}{dt} = \hat{\kappa}_t \bar{c}_t - k_2 c_t (c_p + c_a) - \hat{\kappa}_t c_t, \quad (3.2) $$

$$ \frac{dc_p}{dt} = \hat{\kappa}_p \bar{c}_p - k_2 c_p c_t \left( 1 + \frac{k_1 k_3}{k_{-1}(k_{-2} + k_3)} \bar{c}_d \right) - \hat{\kappa}_p c_p, \quad (3.3) $$

$$ \frac{dc_a}{dt} = \frac{k_1 k_2 k_3}{k_{-1}(k_{-2} + k_3)} \bar{c}_d c_t c_p - k_2 c_t c_a, \quad (3.4) $$

in which $\bar{c}_d$ is a given parameter representing the constant concentration of available dimeric MT1-MMP
*   **Healthy test case:** we consider a healthy tissue, free from any tumor cells, i.e. $ρ_0 = 0$, and assume that there are no SFs in the CT.

*   **Tumor test case:** we assume that a tumor is present, but no SFs, and the tumor cell density in contact with the BM is given by $ρ_0 = 10^3 χ_{[L/4,3L/4]}$.

*   **Tumor with a senescent fibroblast test case:** we assume a tumour is present, and again consider $ρ_0 = 10^3 χ_{[L/4,3L/4]}$, as well as a SF in the CT at a position ($x_1^{SF}, x_2^{SF}$).

### 3.1.2 Results of the full model under the baseline parameter set

The BM density at day 7 predicted by the model for each test case is displayed in Figure 4 (left), where we also mark the critical BM density for BM rupture (assumed to be $M_{crit} = M_{max}/2$), and the corresponding concentration of TIMP-2 and proMMP-2 in the CT at day 7 is displayed in Figure 5. We also show the time evolution of the BM density in Figure 4 (right), which is still evolving at day 7 in the presence of a tumour – in fact, its dynamics span several weeks (cf. Supplementary Figure S2). This is not the case for the remaining molecular concentrations in our system, which have all reached a steady state by day 7 – in fact, their dynamics only span several hours (cf. Supplementary Figures S3-S5).

<figure><img src="image_6.png" /><figcaption>Figure 4: The model predicts that the presence of a SF speeds up BM digestion, with the fastest digestion rate occurring at the point closest to the SF. BM density predicted by the full model for the healthy (green), tumour (blue) and tumour+SF (red) test cases. The BM density $M(t, x)$ is plotted over time (up to day 7) at $x_1 = x_1^{SF}$, i.e. at the point in the BM closest to the SF, on the left. On the right we plot the spatial distribution of the BM density at day 7, marking the position $x_1 = x_1^{SF}$ closest to the SF with a dotted line. The dashed horizontal line marks the critical BM density $M_{crit}$ below which we expect cells to be able to invade BM. These results are obtained solving System (2.1)–(2.13), under definitions (2.14)–(2.16), boundary conditions (2.17)–(2.19), and initial conditions (2.20)–(2.23), under the parameter values in Table 1 or as specified in Section 3.1.1 for each test case.</figcaption></figure>

**Healthy test case.** In the healthy test case we observe no BM degradation, as displayed in Figure 4 (green line), and as expected in a healthy epithelial tissue. As we assumed no tumor cells are present in a healthy tissue, there is no production of monomeric MT1-MMP and consequently no MMP-2 activation: even though their initial concentrations were positive, they becomes zero exponentially fast (cf. Supplementary Figure S3) due to internalisation and natural decay. Interestingly, we notice that a gradient of TIMP-2 and proMMP-2 in the CT forms in the $x_2$ direction, i.e. orthogonally to the BM, as displayed in Figure 5 (first column). This is due to the transmission conditions at the BM and the reactions occurring at this location, which result in an increase in TIMP-bound proMMPs (cf. Supplementary Figure S3).

**Tumor test case.** In the tumour test case we observe partial degradation of the BM in the region where the tumour is present, as displayed in Figure 4 (blue line). This localized degradation of the BM in the tumor region is explained by the localization of MT1-MMP produced by tumor cells, which ensure
and $\overline{c}_i \equiv \overline{c}_i(x,t)|_\Gamma (i = \{t,p\})$, i.e. the concentrations of TIMP-2 and proMMP-2 at the BM coming from the CT. System (3.1)–(3.4) is coupled to System (2.12)–(2.13).

**Equilibrium concentrations.** We analysed the steady state of the reduced system (3.2)–(3.4), cf. Appendix C.1.1. The corresponding equilibrium concentrations are plotted in Figure 3 for varying values of the concentrations of TIMP-2 and proMMP-2 reaching the BM from the CT. We see that the equilibrium concentration of active MMP-2: increases (linearly) with the amount of proMMP-2 coming from the CT; (linearly) decreases with the amount of TIMP-2 coming from the CT; is zero in the absence of TIMP-2 coming from the CT; increases with the amount of MT1-MMP dimers (cf. Supplementary Figure S1). These results are consistent with the dual role of TIMP-2 in MMP-2 activation.

<figure><img src="image_4.png" /><figcaption>Figure 3: The analytical equilibrium of the reduced system demonstrates that the dual role of TIMP-2, as activator and inhibitor of MMP-2, is still captured in the reduced model. Concentration of TIMP-2 (left), proMMP-1 (center) and active MMP-2 (right) at steady state of the reduced ODE system (3.2)–(3.4), for varying fixed concentrations of proMMP-2 ($\overline{c}_p^*$) and TIMP-2 ($\overline{c}_t^*$) coming from the CT. The steady state is given by (C.3) for $\overline{c}_t^* > 0$, and by (C.4) otherwise, under the parameter values listed in Table 1 and $C_a = 0$. The latter corresponds to the value obtained at equilibrium under initial conditions $c_t = c_p = c_a = 0$ at $t = 0$. We consider TIMP-2 and proMMP-2 concentrations from the CT in the ranges $\overline{c}_t^* \in [0, c_t^{ph} + 1\%c_t^{ph}]$ and $\overline{c}_p^* \in [0, c_p^{ph} + 200\%c_p^{ph}]$, in line with the parameter choices reported in Section 2.3.</figcaption></figure>

**BM rupture and cancer invasion.** We define BM rupture as the decrease of the BM density below an arbitrary value that corresponds to the critical density below which tumor cells can invade. We derived a condition for this long-term sustained BM rupture, allowing cells to begin invading the CT, cf. Appendix C.1.2. From this condition, we have that BM rupture and cancer invasion is more likely to occur for: higher values of active MMP-2 at equilibrium, which depends on the amount of proMMP-2 and TIMP-2 coming from the CT and MT1-MMP2 dimers; higher BM degradation rates; lower BM remodeling rates. These results are consistent with the biological meaning of these quantities.

# 4 Numerical results

To test the model's ability to predict biologically relevant scenarios, we conduct numerical simulations of the full System (2.1)–(2.13), complemented with definitions (2.14)–(2.16), boundary conditions (2.17)–(2.19), and initial conditions (2.20)–(2.23), under the baseline parameter values stated in Table 1. The results are discussed in Section 4.1, while section 4.2 is dedicated to a global sensitivity analysis of both the full and reduced models, using the SAFE (Sensitivity analysis For Everyone) toolbox [92].

**Numerical method and set up.** In Section 4.1 we set $\dim(\Omega) = 2$ and $\dim(\Gamma) = 1$ ("2D" case), while in Section 4.2 we set $\dim(\Omega) = 1$ and $\dim(\Gamma) = 0$ ("1D" case). In particular, in 2D we assume that the domain $\Omega = [0, L] \times [0, L]$ is a square of side length $L$ dm (or, in 1D, a line of length $L$ dm). We consider a uniform discretization of the domain $\Omega$ using $N_x \times N_x$ nodes (or, in 1D, using $N_x$ nodes), and using $N_x$ nodes for $\Gamma$. We construct a cell-centered finite volume spatial discretization of our system and integrate it in time using an explicit Runge-Kutta integration method, up to a final time $T$. For details about our numerical scheme, we refer to Appendix D, and the values of $T$, $L$ and $N_x$ used are reported in Table 2.
<figure><img src="image_7.png" /><figcaption>Figure 5: The inclusion of a SF in the model results in the prediction of accurate gradients of proMMP-2 in the CT, and faster digestion of the BM. TIMP-2 (top) and proMMP-2 (bottom) concentrations in the CT at day 7 predicted by the full model for the healthy (left), tumour (center) and tumour+SF (right) test cases. On the right plots, the cross marks the position of the SF in the CT, i.e. the center of diffusion of SF secretome. The BM is marked in red at $x_2 = L$ along with its status at day 7: intact (solid line), partially degraded (dashed) or ruptured (dotted), i.e. degraded below a critical value, cf. Figure 4 (left). These results are obtained solving System (2.1)-(2.13), under definitions (2.14)-(2.16), boundary conditions (2.17)-(2.19), and initial conditions (2.20)-(2.23), under the parameter values in Table 1 or as specified in Section 3.1.1 for each test case.</figcaption></figure>

positive levels of active MMP-2 persist locally (cf. Supplementary Figure S3). This supported by the experimental results of [45], in which we observe a localized activation of proMMP-2 close to the BM. Nonetheless, the resulting levels of active MMP-2 do not suffice to predict cancer invasion within one week – this will happen as time advances (cf. Supplementary Figure S2). On the other hand, we observe in Figure 5 (second column) a gradient of TIMP-2 and proMMP-2 concentrations in the CT analogous to that observed in the healthy test case, with slightly lower levels of TIMP-2 recorded in the CT at day 7. This is likely due to the recruitment of these molecules in the MMP-2 activation cascade, as demonstrated by the increase in the levels of molecular complexes (cf. Supplementary Figure S3, bottom row).

**Tumour with a SF test case.** With the addition of a SF close to the BM – specifically, at a distance of 0.5 mm in Figures 4 and 5 – we instead observe BM rupture within 5 days, as displayed in Figure 4 (red line). This is due to the large amount of proMMP-2 secreted by the SF that reached the BM in the tumour region, creating a larger imbalance between activation and inhibition compared to the tumour test case, resulting in a larger amount of active MMP-2 (cf. Supplementary Figure S5) and ultimately faster BM digestion. Interestingly, the digestion of the BM is asymmetric in space, with the fastest digestion rate occurring at $x_1 = x_1^{SF}$, i.e. in correspondence of the SF position, at which point we expect complete BM digestion shortly after 10 days (cf. Supplementary Figure S2). This is due to the radial diffusion of the SF secretome from the SF position in the CT, as displayed in Figure 5 (bottom right plot), which creates spatial heterogeneity in the enzymes concentration fields in the BM mirroring the pattern observed in the BM, with the largest amount of proMMP-2 in the BM occurring at the point closest to the BM (cf. Supplementary Figure S5), i.e. at $x_1 = x_1^{SF}$ in our geometric setting. The spatial pattern of proMMP-2 in the CT is not observed for TIMP-2 in Figure 5 (top right plot) only because production of TIMP-2 by the SF is assumed to be small, and thus the effect of the SF on this concentration field appears negligible.
<table><caption>Table 2: Simulation parameters used in Section 4.1 and Section 4.2.</caption><thead><tr><th>Parameter</th><th>Meaning</th><th>Unit</th><th>Value S.4.1</th><th>Value S.4.2</th></tr></thead><tbody><tr><td>dim(Ω)</td><td>Dimension of the CT</td><td></td><td>2</td><td>1</td></tr><tr><td>T</td><td>End time</td><td>days</td><td>7</td><td>1</td></tr><tr><td>L</td><td>Length in each dimension of Ω</td><td>dm</td><td>0.1</td><td>0.1</td></tr><tr><td>N<sub>x</sub></td><td>Number of nodes in the x<sub>i</sub> direction</td><td></td><td>50</td><td>20</td></tr></tbody></table>

# 4.1 Model's ability to reproduce biologically relevant scenarios

## 4.1.1 Biologically relevant test cases

To test the model's ability to reproduce biologically relevant scenarios, we consider three test cases:

* *Healthy test case:* we consider a healthy tissue, free from any tumor cells, i.e. $\rho_0 = 0$, and assume that there are no SFs in the CT.

* *Tumor test case:* we assume that a tumor is present, but no SFs, and the tumor cell density in contact with the BM is given by $\rho_0 = 10^3 \chi_{[L/4,3L/4]}$.

* *Tumor with a senescent fibroblast test case:* we assume a tumour is present, and again consider $\rho_0 = 10^3 \chi_{[L/4,3L/4]}$, as well as a SF in the CT at a position $(x_1^{SF}, x_2^{SF})$.

## 4.1.2 Results of the full model under the baseline parameter set

The BM density at day 7 predicted by the model for each test case is displayed in Figure 4 (left), where we also mark the critical BM density for BM rupture (assumed to be $M_{crit} = M_{max}/2$), and the corresponding concentration of TIMP-2 and proMMP-2 in the CT at day 7 is displayed in Figure 5. We also show the time evolution of the BM density in Figure 4 (right), which is still evolving at day 7 in the presence of a tumour – in fact, its dynamics span several weeks (cf. Supplementary Figure S2). This is not the case for the remaining molecular concentrations in our system, which have all reached a steady state by day 7 – in fact, their dynamics only span several hours (cf. Supplementary Figures S3-S5).

<figure><img src="image_5.png" /><figcaption>Figure 4: The model predicts that the presence of a SF speeds up BM digestion, with the fastest digestion rate occurring at the point closest to the SF. BM density predicted by the full model for the healthy (green), tumour (blue) and tumour+SF (red) test cases. The BM density $M(t, x)$ is plotted over time (up to day 7) at $x_1 = x_1^{SF}$, i.e. at the point in the BM closest to the SF, on the left. On the right we plot the spatial distribution of the BM density at day 7, marking the position $x_1 = x_1^{SF}$ closest to the SF with a dotted line. The dashed horizontal line marks the critical BM density $M_{crit}$ below which we expect cells to be able to invade BM. These results are obtained solving System (2.1)-(2.13), under definitions (2.14)-(2.16), boundary conditions (2.17)-(2.19), and initial conditions (2.20)-(2.23), under the parameter values in Table 1 or as specified in Section 4.1.1 for each test case.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 6: The model predicts that the presence of a SF in the CT can induce BM rupture according to a delicate balance between its distance from the BM and the rates at which it secretes proMMP-2 and TIMP-2. BM density at the point closest to the SF (i.e. at $x_1 = x_1^{SF}$) predicted by the full model for the tumour+SF test case, varying the distance between the SF and the BM. In particular, we consider a SF located at $x_2^{SF} = 0.005$ dm (left), at $x_2^{SF} = 0.05$ dm (center) and at $x_2^{SF} = 0.095$ dm (right), recalling that the BM is located at $x_2 = 0.1$ dm. In each plot, the BM density predicted by the model at day 7 is shown for varying values of the proMMP-2 secretion rate $r_p \in [0, 4.17 \times 10^{-3}]$ s⁻¹ and the TIMP-2 secretion rate $r_t \in [0, 2.22 \times 10^{-4}]$ s⁻¹ by the SF. The black line marks the critical BM density $M_{crit}$ below which ECM pore size are sufficiently large for cancer cells to pass through, here referred to as BM rupture. Following this criterion, the regions of the parameter space in which the model predicts “RUPTURE” or “NO RUPTURE” by day 7 are marked on each plot. These results are obtained solving System (2.1)-(2.13), under definitions (2.14)-(2.16), boundary conditions (2.17)-(2.19), and initial conditions (2.20)-(2.23), under the parameter values in Table 1 or as specified in Section 3.1.1 for the Tumour+SF test case.</figcaption></figure>

### 3.1.3 Further exploring the role of SF location and secretion rates

Under the baseline parameter set it is clear that the presence of a SF close to the BM may significantly speed up BM digestion and, consequently, cancer invasion due to their secretome. We recall that, due to a lack of empirical measurements, the secretion rates of TIMP-2 and proMMP-2 by SFs were arbitrarily set in the baseline scenario. We thus explored how the BM density at $x_1 = x_1^{SF}$, i.e. the point at the BM most affected by SF secretome, varying the distance of the SF from the BM, and its TIMP-2 and proMMP-2 secretion rates, as displayed in Figure 6. We observe that higher secretion rates of proMMP-2 and lower secretion rates of TIMP-2 by SFs result in a lower BM density at day 7 and thus higher change of observing BM rupture within a week. This is in line with the results of the steady state analysis of the reduced system and the equilibrium concentration of active MMP-2 shown in Figure 3. However, if the SF does not secrete any TIMP-2, BM degradation may still occur, suggesting that physiological levels of TIMP-2 in the CT suffice to initiate the MMP-2 activation cascade, and any TIMP-2 produced by SF will likely play a saturating role. Finally, it is clear from Figure 6 that decreasing the distance of the SF from the BM increases its influence on BM dynamics, and BM rupture may be observed at day 7 for a wider range of TIMP-2 and proMMP-2 secretion rates. We expect that increasing the proportion of SFs in the CT close to the BM would lead analogous results.

## 3.2 Global sensitivity analysis of the full and reduced models

### 3.2.1 Description of the global sensitivity analysis method

Global sensitivity analysis (GSA) is a powerful tool to study how the uncertainty in the model's outputs depends on the uncertainties of the inputs. To better understand the relative importance of the parameters of the model and identify key ones, we conduct GSA using the SAFE toolbox [54, 53] for Matlab. In particular, we employ the Elementary Effects (EEs) Test included in the toolbox, encoding a revised version of the Morris method [48]. This GSA method is based on the calculation of the EEs associated with each model input (parameter values, in our case), and revised versions of the Morris method consider the absolute value of EEs to avoid potential cancellation effects in the summation of EEs with opposite
<figure><img src="image_6.png" /><figcaption>Figure 5: The inclusion of a SF in the model results in the prediction of accurate gradients of proMMP-2 in the CT, and faster digestion of the BM. TIMP-2 (top) and proMMP-2 (bottom) concentrations in the CT at day 7 predicted by the full model for the healthy (left), tumour (center) and tumour+SF (right) test cases. On the right plots, the cross marks the position of the SF in the CT, i.e. the center of diffusion of SF secretome. The BM is marked in red at $x_2 = L$ along with its status at day 7: intact (solid line), partially degraded (dashed) or ruptured (dotted), i.e. degraded below a critical value, cf. Figure 4 (left). These results are obtained solving System (2.1)-(2.13), under definitions (2.14)-(2.16), boundary conditions (2.17)-(2.19), and initial conditions (2.20)-(2.23), under the parameter values in Table 1 or as specified in Section 4.1.1 for each test case.</figcaption></figure>

**Healthy test case.** In the healthy test case we observe no BM degradation, as displayed in Figure 4 (green line), and as expected in a healthy epithelial tissue. As we assumed no tumor cells are present in a healthy tissue, there is no production of monomeric MT1-MMP and consequently no MMP-2 activation: even though their initial concentrations were positive, they becomes zero exponentially fast (cf. Supplementary Figure S3) due to internalisation and natural decay. Interestingly, we notice that a gradient of TIMP-2 and proMMP-2 in the CT forms in the $x_2$ direction, i.e. orthogonally to the BM, as displayed in Figure 5 (first column). This is due to the transmission conditions at the BM and the reactions occurring at this location, which result in an increase in TIMP-bound proMMPs (cf. Supplementary Figure S3).

**Tumor test case.** In the tumour test case we observe partial degradation of the BM in the region where the tumour is present, as displayed in Figure 4 (blue line). This localized degradation of the BM in the tumor region is explained by the localization of MT1-MMP produced by tumor cells, which ensure positive levels of active MMP-2 persist locally (cf. Supplementary Figure S3). This supported by the experimental results of [81], in which we observe a localized activation of proMMP-2 close to the BM. Nonetheless, the resulting levels of active MMP-2 do not suffice to predict cancer invasion within one week – this will happen as time advances (cf. Supplementary Figure S2). On the other hand, we observe in Figure 5 (second column) a gradient of TIMP-2 and proMMP-2 concentrations in the CT analogous to that observed in the healthy test case, with slightly lower levels of TIMP-2 recorded in the CT at day 7. This is likely due to the recruitment of these molecules in the MMP-2 activation cascade, as demonstrated by the increase in the levels of molecular complexes (cf. Supplementary Figure S3, bottom row).

**Tumour with a SF test case.** With the addition of a SF close to the BM – specifically, at a distance of 0.5 mm in Figures 4 and 5 – we instead observe BM rupture within 5 days, as displayed in Figure 4 (red line). This is due to the large amount of proMMP-2 secreted by the SF that reached the BM in the tumour region, creating a larger imbalance between activation and inhibition compared to the tumor test case, resulting in a larger amount of active MMP-2 (cf. Supplementary Figure S5) and ultimately faster BM digestion. Interestingly, the digestion of the BM is asymmetric in space, with the fastest digestion
signs. In particular:

* a large mean of the absolute value of the EEs indicates the parameter has individual important effect and, viceversa, a small mean indicates the parameter has a negligible effect on the output;

* a large standard deviation of the absolute value of the EEs indicates the parameter has a non-linear effect or one dependent on its interaction with other parameters.

The EE Test is the least computationally expensive GSA method, and it is therefore particularly suited for screening and ranking of model inputs in more involved mathematical models, like the one presented here. For an overview of GSA methods, we refer the interested reader to the review article [32].

As GSA methods need to run a lot of simulations of the model, we restrict our analysis to the one dimensional case, i.e. dim($\Omega$) = 1 and dim($\Gamma$) = 0. As the addition of a SF in the CT can particularly speed up BM digestion, as seen in Section 3.1, we set up our GSA to target variations of the BM density $M$. To capture variations in the time evolution of this quantity, we consider $N_t = 100$ equally spaced time points $t_j$ ($j = 1, \dots, 100$) at which we save the BM density $M(t_j)$, and compute

$$Y = \Delta t \sum_{j=1}^{N_t} \frac{|M(t_j)|}{|M^0|}, \qquad (3.1)$$

where $\Delta t = T/N_t$ is the spacing between two saved time points. From now on, the quantity $Y$ will be referred to as the output of our simulations. The spatio-temporal parameters used in each model simulation (i.e. $T$, $L$, $N_x$) for the GSA are specified in Table 2 (last column). The parameter sets considered as inputs during the GSA are specified in following sections for the full and reduced model, respectively. For each parameters set, the GSA relies on 500 points sampled from the parameter space using the Latin Hypercube technique with a radial sampling method. The parameter space is defined considering each parameter has a uniform probability distribution in a range chosen to span several orders of magnitude both below and above that of the baseline parameter value. These ranges are specified in Table 3 for the parameters of the full model, and for $\bar{c}_d$ we consider values in the range $[10^{-7}, 10^{-3}]$.

<table><caption>Table 3: Lower and upper bounds of the parameter values considered in the GSA of the full model, sorted in parameter Set 1 (left) and Set 2 (right). The unit of measurement of each parameter is given in Table 1.</caption><thead><tr><th>Parameter</th><th>Lower bound</th><th>Upper bound</th><th>Parameter</th><th>Lower bound</th><th>Upper bound</th></tr></thead><tbody><tr><td>$k_0$</td><td>0.02</td><td>200</td><td>$\gamma$</td><td>2.36 &times; 10<sup>-3</sup></td><td>23.6</td></tr><tr><td>$k_{-0}$</td><td>1 &times; 10<sup>-4</sup></td><td>1</td><td>$K_M$</td><td>1.357 &times; 10<sup>-3</sup></td><td>13.57</td></tr><tr><td>$k_1$</td><td>0.0271</td><td>271</td><td>$r_M$</td><td>6.18 &times; 10<sup>-6</sup></td><td>6.18 &times; 10<sup>-2</sup></td></tr><tr><td>$k_{-1}$</td><td>1 &times; 10<sup>-6</sup></td><td>0.01</td><td>$\rho_0$</td><td>10</td><td>10000</td></tr><tr><td>$k_2$</td><td>0.0014</td><td>14</td><td>$\alpha_m$</td><td>5 &times; 10<sup>-6</sup></td><td>0.05</td></tr><tr><td>$k_{-2}$</td><td>1 &times; 10<sup>-6</sup></td><td>0.01</td><td>$D_t$</td><td>1.29 &times; 10<sup>-8</sup></td><td>1.29 &times; 10<sup>-4</sup></td></tr><tr><td>$k_3$</td><td>2 &times; 10<sup>-4</sup></td><td>2</td><td>$D_p$</td><td>1.29 &times; 10<sup>-8</sup></td><td>1.29 &times; 10<sup>-4</sup></td></tr><tr><td>$\beta_t$</td><td>4.45 &times; 10<sup>-7</sup></td><td>4.45 &times; 10<sup>-3</sup></td><td>$\kappa_t$</td><td>1.935 &times; 10<sup>-8</sup></td><td>1.935 &times; 10<sup>-4</sup></td></tr><tr><td>$\beta_p$</td><td>6.08 &times; 10<sup>-6</sup></td><td>6.08 &times; 10<sup>-2</sup></td><td>$\kappa_p$</td><td>1.935 &times; 10<sup>-8</sup></td><td>1.935 &times; 10<sup>-4</sup></td></tr><tr><td>$\beta_a$</td><td>1.29 &times; 10<sup>-8</sup></td><td>1.29 &times; 10<sup>-4</sup></td><td>$s_t^{ph}$</td><td>4.272 &times; 10<sup>-6</sup></td><td>4.272 &times; 10<sup>-2</sup></td></tr><tr><td>$\beta_m$</td><td>3.85 &times; 10<sup>-4</sup></td><td>3.85</td><td>$s_p^{ph}$</td><td>1.0944 &times; 10<sup>-4</sup></td><td>1.0944</td></tr><tr><td>$\beta_d$</td><td>3.85 &times; 10<sup>-4</sup></td><td>3.85</td><td>$r_t$</td><td>2.22 &times; 10<sup>-8</sup></td><td>2.22 &times; 10<sup>-4</sup></td></tr><tr><td>$\beta_{tp}$</td><td>0.036</td><td>360</td><td>$r_p$</td><td>8.28 &times; 10<sup>-6</sup></td><td>8.28 &times; 10<sup>-2</sup></td></tr><tr><td>$\beta_{ta}$</td><td>3.5 &times; 10<sup>-4</sup></td><td>3.5</td><td>$x_1^{SF}$</td><td>0</td><td>0.1</td></tr></tbody></table>

### 3.2.2 Sensitivity analysis of the full model

We perform GSA on the full model (2.1)–(2.13), by analysing separately two parameters sets:

* Set 1: The biochemical reaction coefficients ($k_0, k_{-0}, k_1, k_{-1}, k_2, k_{-2}, k_3$) and the molecular decay coefficients ($\beta_t, \beta_t, \beta_p, \beta_m, \beta_d, \beta_{tp}, \beta_{ta}$);

* Set 2: The parameters related to the BM degradation and production ($\gamma, K_M, r_M$), the cancer cells and their MT1-MMP expression ($\rho_0, \alpha_m$), the production and diffusion of TIMP-2 and proMMP-2 in CT ($s_t^{ph}, s_p^{ph}, r_t, r_p, D_p, D_t$) and their transmission to the BM ($\kappa_t, \kappa_p$), the position of the SF ($x_1^{SF}$).
<figure><img src="image_7.png" /><figcaption>Figure 6: The model predicts that the presence of a SF in the CT can induce BM rupture according to a delicate balance between its distance from the BM and the rates at which it secretes proMMP-2 and TIMP-2. BM density at the point closest to the SF (i.e. at $x_1 = x_1^{SF}$) predicted by the full model for the tumour+SF test case, varying the distance between the SF and the BM. In particular, we consider a SF located at $x_2^{SF} = 0.005$ dm (left), at $x_2^{SF} = 0.05$ dm (center) and at $x_2^{SF} = 0.095$ dm (right), recalling that the BM is located at $x_2 = 0.1$ dm. In each plot, the BM density predicted by the model at day 7 is shown for varying values of the proMMP-2 secretion rate $r_p \in [0, 4.17 \times 10^{-3}]$ s⁻¹ and the TIMP-2 secretion rate $r_t \in [0, 2.22 \times 10^{-4}]$ s⁻¹ by the SF. The black line marks the critical BM density $M_{crit}$ below which ECM pore size are sufficiently large for cancer cells to pass through, here referred to as BM rupture. Following this criterion, the regions of the parameter space in which the model predicts “RUPTURE” or “NO RUPTURE” by day 7 are marked on each plot. These results are obtained solving System (2.1)-(2.13), under definitions (2.14)-(2.16), boundary conditions (2.17)-(2.19), and initial conditions (2.20)-(2.23), under the parameter values in Table 1 or as specified in Section 4.1.1 for the Tumour+SF test case.</figcaption></figure>

rate occurring at $x_1 = x_1^{SF}$, i.e. in correspondence of the SF position, at which point we expect complete BM digestion shortly after 10 days (cf. Supplementary Figure S2). This is due to the radial diffusion of the SF secretome from the SF position in the CT, as displayed in Figure 5 (bottom right plot), which creates spatial heterogeneity in the enzymes concentration fields in the BM mirroring the pattern observed in the BM, with the largest amount of proMMP-2 in the BM occurring at the point closest to the BM (cf. Supplementary Figure S5), i.e. at $x_1 = x_1^{SF}$ in our geometric setting. The spatial pattern of proMMP-2 in the CT is not observed for TIMP-2 in Figure 5 (top right plot) only because production of TIMP-2 by the SF is assumed to be small, and thus the effect of the SF on this concentration field appears negligible.

### 4.1.3 Further exploring the role of SF location and secretion rates

Under the baseline parameter set it is clear that the presence of a SF close to the BM may significantly speed up BM digestion and, consequently, cancer invasion due to their secretome. We recall that, due to a lack of empirical measurements, the secretion rates of TIMP-2 and proMMP-2 by SFs were arbitrarily set in the baseline scenario. We thus explored how the BM density at $x_1 = x_1^{SF}$, i.e. the point at the BM most affected by SF secretome, varying the distance of the SF from the BM, and its TIMP-2 and proMMP-2 secretion rates, as displayed in Figure 6. We observe that higher secretion rates of proMMP-2 and lower secretion rates of TIMP-2 by SFs result in a lower BM density at day 7 and thus higher change of observing BM rupture within a week. This is in line with the results of the steady state analysis of the reduced system and the equilibrium concentration of active MMP-2 shown in Figure 3. However, if the SF does not secrete any TIMP-2, BM degradation may still occur, suggesting that physiological levels of TIMP-2 in the CT suffice to initiate the MMP-2 activation cascade, and any TIMP-2 produced by SF will likely play a saturating role. Finally, it is clear from Figure 6 that decreasing the distance of the SF from the BM increases its influence on BM dynamics, and BM rupture may be observed at day 7 for a wider range of TIMP-2 and proMMP-2 secretion rates. We expect that increasing the proportion of SFs in the CT close to the BM would lead analogous results.
**Remark 2.** We emphasize that variations of the parameter $\beta_d$ affect Equation (2.3) and Equations (2.7)-(2.9) as $\beta_d = \beta_1 = \beta_2 = \beta_3$. Similarly, variations of $\kappa_i$ ($i = \{p, t\}$) affect Equations (2.18)-(2.19), but also Equations (2.4)-(2.5) via definition (2.14) as $\hat{\kappa}_i = \kappa_i/\varepsilon$ (see (A.16)).

<figure><img src="image_9.png" /><figcaption>Figure 7: GSA highlights the most influential parameters of the full model (2.1)–(2.13) on the speed of BM digestion. Using Morris' method, we quantify the impact of changes in parameter values on the time evolution of the BM density. The boxplots display mean and standard deviation of the EEs of each parameter in Set 1 (a) and Set 2 (b), as prescribed in Section 3.2.2, on the output $Y$ defined in (3.1) (see Section 3.2.1 for details on the methodology). The most influential parameters are the decay rates of MT1-MMP dimers $\beta_d$ and TIMP-2 $\beta_t$, the activation rate of MMP-2 $k_3$, the rate of MT1-MMP monomer expression by cancer cells $\alpha_m$ and the tumor cell density $\rho_0$. Changes of these parameters lead to large linear and nonlinear effects on the output.</figcaption></figure>

Figure 7 summarizes the results of the GSA of the full model. We observe that, for parameter Set 1 (Figure 7a), the most important parameter is the decay rate $\beta_d$. This is explained by the fact that changes in this value not only affect the concentration of MT1-MMP dimers, which is crucial in the activation process, but also that of all complexes involving this compound (cf. Remark 2), resulting in an amplified effect of this parameter. Changes in the TIMP-2 decay rate $\beta_t$ also lead to large variations in the output, visibly larger than those observed for the decay rates of proMMP-2 and MMP-2. As this parameter directly affects the amount of TIMP-2, this result further highlights the importance of the dual role of TIMP-2 as both activator and inhibitor of MMP-2. The parameter with third largest impact on the model output is $k_3$, the rate of MMP-2 activation at the end of the activation cascade, which directly impacts the amount of active MMP-2. The larger effect of this parameter, compared to those of the association coefficients involved in the activation process $k_0$, $k_1$ and $k_2$, could be attributed to the irreversibility of the reaction it relates to. On the other hand, the dissociation rates $k_{-i}$ ($i = \{0, 1, 2\}$) and decay rates of TIMP-bound active and inactive MMP-2 appear to have a negligible effect. In parameter Set 2 (Figure 7b), we notice that the MT1-MMP expression rate by cancer cells $\alpha_m$ has the largest effect on the outcome, immediately followed by the cancer cell density $\rho_0$. This re-stresses the importance for BM digestion of MT1-MMP production on the surface of cancer cells, indeed one of the corner stones of MMP-2 activation. We note that the Michaelis-Menten constant of BM proteolytic degradation $K_M$, the BM production rate $r_M$ and TIMP-2 production rate by the SF $r_t$ appear to have negligible effect, likely due to the relatively small value of these parameters. Last, but not least, we remark that all most influential parameters induce both linear and non-linear effects. This means that the influence of an individual change of any of these parameters lead to an important variation in the output but also that changing two or more of these influential parameters at the same time results in an even larger variation of the output due to cumulative effects.

### 3.2.3 Global sensitivity analysis of the reduced system

We perform GSA on the reduced system (2.24)–(2.27) with (2.12)-(2.13), considering two sets of parameters:
## 4.2 Global sensitivity analysis of the full and reduced models

### 4.2.1 Description of the global sensitivity analysis method

Global sensitivity analysis (GSA) is a powerful tool to study how the uncertainty in the model's outputs depends on the uncertainties of the inputs. To better understand the relative importance of the parameters of the model and identify key ones, we conduct GSA using the SAFE toolbox [92, 91] for Matlab. In particular, we employ the Elementary Effects (EEs) Test included in the toolbox, encoding a revised version of the Morris method [84]. This GSA method is based on the calculation of the EEs associated with each model input (parameter values, in our case), and revised versions of the Morris method consider the absolute value of EEs to avoid potential cancellation effects in the summation of EEs with opposite signs. In particular:

* a large mean of the absolute value of the EEs indicates the parameter has individual important effect and, viceversa, a small mean indicates the parameter has a negligible effect on the output;

* a large standard deviation of the absolute value of the EEs indicates the parameter has a non-linear effect or one dependent on its interaction with other parameters.

The EE Test is the least computationally expensive GSA method, and it is therefore particularly suited for screening and ranking of model inputs in more involved mathematical models, like the one presented here. For an overview of GSA methods, we refer the interested reader to the review article [58].

As GSA methods need to run a lot of simulations of the model, we restrict our analysis to the one dimensional case, i.e. $dim(\Omega) = 1$ and $dim(\Gamma) = 0$. As the addition of a SF in the CT can particularly speed up BM digestion, as seen in Section 4.1, we set up our GSA to target variations of the BM density $M$. To capture variations in the time evolution of this quantity, we consider $N_t = 100$ equally spaced time points $t_j$ ($j = 1, \dots, 100$) at which we save the BM density $M(t_j)$, and compute

$$Y = \Delta t \sum_{j=1}^{N_t} \frac{|M(t_j)|}{|M^0|}, \qquad (4.1)$$

where $\Delta t = T/N_t$ is the spacing between two saved time points. From now on, the quantity $Y$ will be referred to as the output of our simulations. The spatio-temporal parameters used in each model simulation (i.e. $T$, $L$, $N_x$) for the GSA are specified in Table 2 (last column). The parameter sets considered as inputs during the GSA are specified in following sections for the full and reduced model, respectively. For each parameters set, the GSA relies on 500 points sampled from the parameter space using the Latin Hypercube technique with a radial sampling method. The parameter space is defined considering each parameter has a uniform probability distribution in a range chosen to span several orders of magnitude both below and above that of the baseline parameter value. These ranges are specified in Table 3 for the parameters of the full model, and for $\bar{c}_d$ we consider values in the range $[10^{-7}, 10^{-3}]$.

<table><caption>Table 3: Lower and upper bounds of the parameter values considered in the GSA of the full model, sorted in parameter Set 1 (left) and Set 2 (right). The unit of measurement of each parameter is given in Table 1.</caption><thead><tr><th>Parameter</th><th>Lower bound</th><th>Upper bound</th><th>Parameter</th><th>Lower bound</th><th>Upper bound</th></tr></thead><tbody><tr><td>$k_0$</td><td>0.02</td><td>200</td><td>$\gamma$</td><td>2.36 &times; 10<sup>-3</sup></td><td>23.6</td></tr><tr><td>$k_{-0}$</td><td>1 &times; 10<sup>-4</sup></td><td>1</td><td>$K_M$</td><td>1.357 &times; 10<sup>-3</sup></td><td>13.57</td></tr><tr><td>$k_1$</td><td>0.0271</td><td>271</td><td>$r_M$</td><td>6.18 &times; 10<sup>-6</sup></td><td>6.18 &times; 10<sup>-2</sup></td></tr><tr><td>$k_{-1}$</td><td>1 &times; 10<sup>-6</sup></td><td>0.01</td><td>$\rho_0$</td><td>10</td><td>10000</td></tr><tr><td>$k_2$</td><td>0.0014</td><td>14</td><td>$\alpha_m$</td><td>5 &times; 10<sup>-6</sup></td><td>0.05</td></tr><tr><td>$k_{-2}$</td><td>1 &times; 10<sup>-6</sup></td><td>0.01</td><td>$D_t$</td><td>1.29 &times; 10<sup>-8</sup></td><td>1.29 &times; 10<sup>-4</sup></td></tr><tr><td>$k_3$</td><td>2 &times; 10<sup>-4</sup></td><td>2</td><td>$D_p$</td><td>1.29 &times; 10<sup>-8</sup></td><td>1.29 &times; 10<sup>-4</sup></td></tr><tr><td>$\beta_t$</td><td>4.45 &times; 10<sup>-7</sup></td><td>4.45 &times; 10<sup>-3</sup></td><td>$\kappa_t$</td><td>1.935 &times; 10<sup>-8</sup></td><td>1.935 &times; 10<sup>-4</sup></td></tr><tr><td>$\beta_p$</td><td>6.08 &times; 10<sup>-6</sup></td><td>6.08 &times; 10<sup>-2</sup></td><td>$\kappa_p$</td><td>1.935 &times; 10<sup>-8</sup></td><td>1.935 &times; 10<sup>-4</sup></td></tr><tr><td>$\beta_a$</td><td>1.29 &times; 10<sup>-8</sup></td><td>1.29 &times; 10<sup>-4</sup></td><td>$s_t^{ph}$</td><td>4.272 &times; 10<sup>-6</sup></td><td>4.272 &times; 10<sup>-2</sup></td></tr><tr><td>$\beta_m$</td><td>3.85 &times; 10<sup>-4</sup></td><td>3.85</td><td>$s_p^{ph}$</td><td>1.0944 &times; 10<sup>-4</sup></td><td>1.0944</td></tr><tr><td>$\beta_d$</td><td>3.85 &times; 10<sup>-4</sup></td><td>3.85</td><td>$r_t$</td><td>2.22 &times; 10<sup>-8</sup></td><td>2.22 &times; 10<sup>-4</sup></td></tr><tr><td>$\beta_{tp}$</td><td>0.036</td><td>360</td><td>$r_p$</td><td>8.28 &times; 10<sup>-6</sup></td><td>8.28 &times; 10<sup>-2</sup></td></tr><tr><td>$\beta_{ta}$</td><td>3.5 &times; 10<sup>-4</sup></td><td>3.5</td><td>$x_1^{SF}$</td><td>0</td><td>0.1</td></tr></tbody></table>
Set 1: The biochemical reaction coefficients ($k_1, k_{-1}, k_2, k_{-2}, k_3$) and the fixed concentration of free dimeric MT1-MMPs ($\bar{c}_d$).

Set 2: The proMMP-2 and TIMP-2 decay rates ($\beta_p, \beta_t$), their production and diffusion rates in the CT ($s_t^{ph}, s_p^{ph}, r_t, r_p, D_p, D_t$), their transmission to the BM ($\kappa_t, \kappa_p$), the parameters related to the digestion and the production of the BM ($\gamma, r_M, K_M$), and the position of the SF in the CT ($x_1^{SF}$).

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

Figure 8: **GSA indicates that the reduced model (2.24)-(2.27) with (2.12)-(2.13) is more sensitive to uncertainty in the most influential parameters compared to the full model.** Using Morris' method, we quantify the impact of changes in parameter values on the time evolution of the BM density. The boxplots display mean and standard deviation of the EEs of each parameter in Set 1 (a) and Set 2 (b), as prescribed in Section 3.2.3, on the output Y defined in (3.1) (see Section 3.2.1 for details on the methodology). The most influential parameter in the reduced system is the fixed concentration of free dimeric MT1-MMPs $\bar{c}_d$, followed by the BM degradation rate $\gamma$ and the proMMP-2 physiological production rate $s_p^{ph}$. EEs means and standard deviation of the most influential parameters of the reduced system are larger than those of the full system, indicating increased sensitivity of the reduced model.

Figure 8 summarizes the results of the GSA of the reduced model. The main observation is that the value of $\bar{c}_d$ (Figure 8a), i.e. the fixed concentration of free dimeric MT1-MMPs, is crucial for the output of the reduced system, with EEs mean notably larger than those obtained for the full system. Indeed this parameter encapsulates the effect of MT1-MMP expression and internalisation by cancer cells, and that of MT1-MMP dimerisation and de-dimerisation, i.e. it includes the effect of $\rho_0$, $\alpha_m$ and $\beta_d$ which were among the most influential parameters in the full model. Notably, the MMP-2 activation rate $k_3$ appears to have negligible effects in the reduced system, as opposed to its strong impact in the full model. This might result from the quasi steady-state assumption on the activation complexes imposed during model reduction – mathematically, this is explained by the fact that $k_3$ appears at the nominator and denominator of the same fraction in equations (2.26) and (2.27). The reduced impact of $k_3$ likely favours the striking effect of $\bar{c}_d$. Another parameter that has a remarkably larger EEs mean in the GSA of the reduced model, compared to that found for the full model, is the BM degradation rate $\gamma$ (Figure 8b). Interestingly, we observe that the parameters in Set 2 relating to proMMP-2 dynamics appear to have consistently larger effects on the output than the analogous TIMP-2 related parameters, with the physiological production rate of proMMP-2 $s_p^{ph}$ having gained significant importance. While the model's sensitivity to uncertainty in these inputs has changed upon model reduction, the negligible effect of $K_M, r_M$ and $r_t$ remains consistent across models. Finally, we again have that all most influential parameters induce both linear and non-linear effects, yet with a stronger impact than those identified for the full system, indicating an increased sensitivity upon model reduction.

# 4 Discussion

It is now widely recognised that cells in the stroma may play major roles in tumour progression. Our goal is to develop a mathematical model which can implement empirical data and suggest promising
4.2.2 Sensitivity analysis of the full model

We perform GSA on the full model (2.1)–(2.13), by analysing separately two parameters sets:

Set 1: The biochemical reaction coefficients ($k_0, k_{-0}, k_1, k_{-1}, k_2, k_{-2}, k_3$) and the molecular decay coefficients ($\beta_t, \beta_t, \beta_p, \beta_m, \beta_d, \beta_{tp}, \beta_{ta}$);

Set 2: The parameters related to the BM degradation and production ($\gamma, K_M, r_M$), the cancer cells and their MT1-MMP expression ($\rho_0, \alpha_m$), the production and diffusion of TIMP-2 and proMMP-2 in CT ($s_t^{ph}, s_p^{ph}, r_t, r_p, D_p, D_t$) and their transmission to the BM ($\kappa_t, \kappa_p$), the position of the SF ($x_1^{SF}$).

**Remark 2.** We emphasize that variations of the parameter $\beta_d$ affect Equation (2.3) and Equations (2.7)–(2.9) as $\beta_d = \beta_1 = \beta_2 = \beta_3$. Similarly, variations of $\kappa_i$ ($i = \{p, t\}$) affect Equations (2.18)–(2.19), but also Equations (2.4)–(2.5) via definition (2.14) as $\hat{\kappa}_i = \kappa_i/\epsilon$ (see (A.16)).

<figure><img src="image_8.png" /><figcaption>Figure 7: GSA highlights the most influential parameters of the full model (2.1)–(2.13) on the speed of BM digestion. Using Morris' method, we quantify the impact of changes in parameter values on the time evolution of the BM density. The boxplots display mean and standard deviation of the EEs of each parameter in Set 1 (a) and Set 2 (b), as prescribed in Section 4.2.2, on the output $Y$ defined in (4.1) (see Section 4.2.1 for details on the methodology). The most influential parameters are the decay rates of MT1-MMP dimers $\beta_d$ and TIMP-2 $\beta_t$, the activation rate of MMP-2 $k_3$, the rate of MT1-MMP monomer expression by cancer cells $\alpha_m$ and the tumor cell density $\rho_0$. Changes of these parameters lead to large linear and nonlinear effects on the output.</figcaption></figure>

Figure 7 summarizes the results of the GSA of the full model. We observe that, for parameter Set 1 (Figure 7a), the most important parameter is the decay rate $\beta_d$. This is explained by the fact that changes in this value not only affect the concentration of MT1-MMP dimers, which is crucial in the activation process, but also that of all complexes involving this compound (cf. Remark 2), resulting in an amplified effect of this parameter. Changes in the TIMP-2 decay rate $\beta_t$ also lead to large variations in the output, visibly larger than those observed for the decay rates of proMMP-2 and MMP-2. As this parameter directly affects the amount of TIMP-2, this result further highlights the importance of the dual role of TIMP-2 as both activator and inhibitor of MMP-2. The parameter with third largest impact on the model output is $k_3$, the rate of MMP-2 activation at the end of the activation cascade, which directly impacts the amount of active MMP-2. The larger effect of this parameter, compared to those of the association coefficients involved in the activation process $k_0, k_1$ and $k_2$, could be attributed to the irreversibility of the reaction it relates to. On the other hand, the dissociation rates $k_{-i}$ ($i = \{0, 1, 2\}$) and decay rates of TIMP-bound active and inactive MMP-2 appear to have a negligible effect. In parameter Set 2 (Figure 7b), we notice that the MT1-MMP expression rate by cancer cells $\alpha_m$ has the largest effect on the outcome, immediately followed by the cancer cell density $\rho_0$. This re-stresses the importance for BM digestion of MT1-MMP production on the surface of cancer cells, indeed one of the corner stones of MMP-2 activation. We note that the Michaelis-Menten constant of BM proteolytic degradation $K_M$, the BM production rate $r_M$ and TIMP-2 production rate by the SF $r_t$ appear to have negligible effect, likely due to the relatively small
experimental investigations, in order to reduce the cost and timeline of the empirical work necessary to promote therapeutic strategies targeting SFs against carcinoma progression. We presented in this article the results of our initial modelling efforts, which focus on the first step of the invasion cascade: the rupture of the BM. We proposed a tractable and detailed model of BM digestion by cancer cells involving the contribution of SFs located in the CT, which builds on the current biological knowledge of this phenomenon.

## 4.1 Summary of modelling efforts and key results

**Our model contains several elements of novelty which ensure that various aspects of biological complexity relevant to carcinoma invasion are captured within the same framework.** To the best of our knowledge, the present work is the first to mathematically represent tumor cells mediated rupture of the BM aided by SF secretome. Remarkable modelling progress to capture cancer invasion following BM digestion was achieved by Gallinato et al. [25], and following extensions [10]. However, in these works the BM digestion was modelled as the result of the activity of degrading enzymes secreted only by the tumor. Our model captures the diffusion of the SF secretome inside the CT and its transmission to the BM, at the other end of which these enzymes may interact with MT1-MMP located on the membrane of tumor cells, and activate MMP-2. The inclusion of these dynamics in the mathematical model enables the study of the partnerships that exist between tumor cells and helping cells in the stroma, which are nowadays well known to aid cancer invasion. While this framework, including fibroblast facilitating ECM degradation by cancer cells, was partially considered in [39] for a Tumour Chamber Invasion Assay, our model captures better the dimensionality of the problem by considering the BM to be a $d-1$ dimensional domain, where $d$ is the dimension of the CT, as similarly done in [25]. Altogether, our model allows for a more complete and biophysically consistent representation of the phenomenon, and the results of numerical simulations highlight the importance of considering these pro-tumour cells located in the stroma, as the digestion of the BM appears faster closer to a SF rich CT area.

Compared to the aforementioned works and previous mathematical models of tumor invasion (see e.g. [8, 9, 35, 50]), our model also comprises a comprehensive description of the biochemical reactions regulating the activation and inhibition of MMP-2, capturing the dual role of TIMP-2 as only previously seen in dedicated studies ignoring spatial dynamics [36, 58]. We remark that such dual role of TIMP-2 is still captured in the reduced model that we derived using typical simplifying assumptions, as demonstrated by the results of our steady state analysis, summarised in Figure 3.

Moreover, additional efforts were made in this work to extract parameter values from the literature, deriving several missing ones from steady state assumptions and reference enzyme concentrations reported in empirical studies. Hence, except for 2 of them (the production rates of enzymes by SF, i.e. $r_p, r_t$), most of the parameters are directly found from experimental studies or previous modelling works that fit these parameters using experimental data. The remaining parameters are obtained from simple computations: starting from empirical measurements of physiological enzyme concentrations, we compute the remaining parameters assuming that the physiological conditions correspond to steady-states in the reaction-diffusion equations. This procedure allowed us to find biologically relevant parameter values. We conducted a sensitivity analysis to estimate how errors on these parameter values could influence the output of the model and that of a reduced model obtained under typical simplifying assumptions. This latter part of the work is crucial, as it also highlights changes in the model's sensitivity upon reduction, which is not very often considered for invasion models.

**The results of numerical simulations are in good qualitative agreement with the current biological understanding of the phenomenon.** The flexibility of the model allowed us to explore the behavior of the system on the relevant test cases, reproducible under different parameter sets, via numerical simulations. The model predicts that the presence of an epithelial tumour may result in slow BM digestion at the tumour location, as opposed to the case of a healthy tissue where the BM remains intact, thanks to the presence of MT1-MMPs on the surface of tumour cells which initiate the MMP-2 activation cascade. Moreover, our results clarify that SF secretome may aggravate the imbalance between the activation and inhibition dynamics of MMP-2, boosting cancer invasion. In fact, the inclusion of SFs in the CT results in the emergence of more accurate gradients of proMMP-2 in the CT, extending radially from the position of the SF, and faster digestion of the BM localized at the points closest to the SF. It is at this point that we expect to first observe BM rupture, i.e. the decrease of BM density below a threshold allowing cells to pass through the BM, and the initiation of cancer invasion. Our results also highlight that the possibility of observing BM rupture within a given time frame increases as the distance between the SF and the BM decreases, as the rate of proMMP-2 secretion by the SF increases and as its TIMP-2
value of these parameters. Last, but not least, we remark that all most influential parameters induce both linear and non-linear effects. This means that the influence of an individual change of any of these parameters lead to an important variation in the output but also that changing two or more of these influential parameters at the same time results in an even larger variation of the output due to cumulative effects.

### 4.2.3 Global sensitivity analysis of the reduced system

We perform GSA on the reduced system (3.1)–(3.4) with (2.12)–(2.13), considering two sets of parameters:

Set 1: The biochemical reaction coefficients ($k_1, k_{-1}, k_2, k_{-2}, k_3$) and the fixed concentration of free dimeric MT1-MMPs ($\bar{c}_d$).

Set 2: The proMMP-2 and TIMP-2 decay rates ($\beta_p, \beta_t$), their production and diffusion rates in the CT ($s_t^{ph}, s_p^{ph}, r_t, r_p, D_p, D_t$), their transmission to the BM ($\kappa_t, \kappa_p$), the parameters related to the digestion and the production of the BM ($\gamma, r_M, K_M$), and the position of the SF in the CT ($x_1^{SF}$).

<figure><img src="image_9.png" /><figcaption>Figure 8: GSA indicates that the reduced model (3.1)–(3.4) with (2.12)–(2.13) is more sensitive to uncertainty in the most influential parameters compared to the full model. Using Morris’ method, we quantify the impact of changes in parameter values on the time evolution of the BM density. The boxplots display mean and standard deviation of the EEs of each parameter in Set 1 (a) and Set 2 (b), as prescribed in Section 4.2.3, on the output Y defined in (4.1) (see Section 4.2.1 for details on the methodology). The most influential parameter in the reduced system is the fixed concentration of free dimeric MT1-MMPs $\bar{c}_d$, followed by the BM degradation rate $\gamma$ and the proMMP-2 physiological production rate $s_p^{ph}$. EEs means and standard deviation of the most influential parameters of the reduced system are larger than those of the full system, indicating increased sensitivity of the reduced model.</figcaption></figure>

Figure 8 summarizes the results of the GSA of the reduced model. The main observation is that the value of $\bar{c}_d$ (Figure 8a), i.e. the fixed concentration of free dimeric MT1-MMPs, is crucial for the output of the reduced system, with EEs mean notably larger than those obtained for the full system. Indeed this parameter encapsulates the effect of MT1-MMP expression and internalisation by cancer cells, and that of MT1-MMP dimerisation and de-dimerisation, i.e. it includes the effect of $\rho_0$, $\alpha_m$ and $\beta_d$ which were among the most influential parameters in the full model. Notably, the MMP-2 activation rate $k_3$ appears to have negligible effects in the reduced system, as opposed to its strong impact in the full model. This might result from the quasi steady-state assumption on the activation complexes imposed during model reduction – mathematically, this is explained by the fact that $k_3$ appears at the nominator and denominator of the same fraction in equations (3.3) and (3.4). The reduced impact of $k_3$ likely favours the striking effect of $\bar{c}_d$. Another parameter that has a remarkably larger EEs mean in the GSA of the reduced model, compared to that found for the full model, is the BM degradation rate $\gamma$ (Figure 8b). Interestingly, we observe that the parameters in Set 2 relating to proMMP-2 dynamics appear to have consistently larger effects on the output than the analogous TIMP-2 related parameters, with the physiological production rate of proMMP-2 $s_p^{ph}$ having gained significant importance. While the model’s sensitivity to uncertainty in
secretion rate decreases. These secretion rates are hard to estimate as they may vary across experimental conditions, tissue samples and among SF population.

While these results are in good qualitative agreement with the current biological understanding of the phenomenon, we remark that the amount of time after which the model predicts BM rupture may not accurately convey the empirical timescale of this phenomenon, which may be faster. First of all, the critical BM density for cancer invasion, which we arbitrarily set to half the carrying capacity, may be higher. Secondly, we considered the presence of only one SF in the CT for illustrative purposes – although the model is flexible enough to represent numerous SFs in different locations within the CT – and *in vivo* the number of SFs may be significantly larger, yielding an enhanced effect of combined SFs' secretome. Lastly, despite the great efforts made in identifying appropriate parameter values, these are taken or estimated from a variety of works in the literature, and model calibration with experimental data is necessary to ensure a wholly coherent parameter set.

**Sensitivity analysis reveals that model prediction accuracy depends on knowledge of the amount of MT1-MMPs that epithelial tumor cell express at the invadopodia, especially after model reduction.** Given the large number of parameters in the full model, we derived a reduced version of the model under biologically motivated simplifying assumptions on the relative scales of the dynamics prescribed within the MMP-2 activation cascade. We remark that these assumptions were further justified by the results of numerical simulations, during which concentrations of soluble molecules reached equilibrium much faster than the BM density. Model reduction not only lowers the degrees of freedom of the system, facilitating future experimental calibration of the model, but also allows us to draw analytical conclusions on the equilibrium of the system and necessary conditions for sustained BM rupture. Moreover, it allows for a comparison with previous cancer invasion mathematical models that included some MMP-2 activation or inhibition dynamics under similar assumptions.

In view of this calibration, we performed GSA to investigate the models’ robustness to uncertainty in the parameter values and identify the most influential parameters on the speed of BM digestion by cancer cells in the presence of a SF in the CT. The most influential parameters of both versions of the model displayed both linear and non-linear effects, indicating uncertainty in these inputs may influence the model output also through its interaction with other parameters, as expected given the complexity of the equations and dynamics prescribed in the model. Notably, we found that the most influential parameters of the full model have lower EEs mean and standard deviation than those of the reduced model. This suggests that the predictive capability of the full model is more robust to uncertainty of parameter values than the reduced one, and it is therefore better suited for *in silico* investigations when the parameters are not accurately measured, despite the higher computational cost.

From the GSA we found that uncertainty in the parameters associated with the amount of MT1-MMPs available on the surface of cancer cells – *i.e.* tumour cell number, MT1-MMP expression rate by tumour cells, internalisation rate of MT1-MMP dimers and complexes involving it – has the largest impact on the output of the full model. This was followed by an important role played by the rate of the irreversible reaction at the end of the MMP-2 activation cascade. Interestingly, we found that this latter parameter had negligible effect on the output of the reduced model, which may be an advantage in the absence of accurate estimates for this rate. However, this advantage is counterbalanced by the intensified effect of the parameter modelling the fixed concentration of free dimeric MT1-MMPs, which only appears in the reduced model capturing the collective effect of the MT1-MMP dynamics on the surface of cancer cells. Such simplification is the main drawback of the reduced model, as enormous importance if given to a parameter which is even more difficult to measure empirically.

## 4.2 Modelling and interdisciplinary perspectives

**Mathematical modelling perspectives.** In this work we use mathematical modelling to represent and study the contribution of SF secretome in initial cancerous progression of carcinoma toward BM rupture and invasion. In particular, we investigate the potential impact of this partnership on the digestion of the BM, which is the first step of carcinoma invasion. While our model is able to capture numerous biological effects, there are several planned extensions that would allow the inclusion of additional biologically relevant features.

The transmission coefficient $\kappa_i$ ($i = t, p$) models the rate at which proteins from the CT are transmitted to and across the BM. It is computed using the diffusion coefficient $D_i$ of the protein and the width of the BM. In particular, we assumed that the diffusivity of each protein in the BM is $10^{-6}$ smaller than in the CT (see Appendix B), a coarse assumption explained by the fact that the BM is a dense sheet of collagen fibers. A more accurate description of the protein diffusivity in the BM could be achieved
these inputs has changed upon model reduction, the negligible effect of $K_M$, $r_M$ and $r_t$ remains consistent across models. Finally, we again have that all most influential parameters induce both linear and non-linear effects, yet with a stronger impact than those identified for the full system, indicating an increased sensitivity upon model reduction.

# 5 Discussion

It is now widely recognised that cells in the stroma may play major roles in tumour progression. Our goal is to develop a mathematical model which can implement empirical data and suggest promising experimental investigations, in order to reduce the cost and timeline of the empirical work necessary to promote therapeutic strategies targeting SFs against carcinoma progression. We presented in this article the results of our initial modelling efforts, which focus on the first step of the invasion cascade: the rupture of the BM. We proposed a tractable and detailed model of BM digestion by cancer cells involving the contribution of SFs located in the CT, which builds on the current biological knowledge of this phenomenon.

## 5.1 Summary of modelling efforts and key results

**Our model contains several elements of novelty which ensure that various aspects of biological complexity relevant to carcinoma invasion are captured within the same framework.** To the best of our knowledge, the present work is the first to mathematically represent tumor cells mediated rupture of the BM aided by SF secretome. Remarkable modelling progress to capture cancer invasion following BM digestion was achieved by Gallinato et al. [41], and following extensions [19]. However, in these works the BM digestion was modelled as the result of the activity of degrading enzymes secreted only by the tumor. Our model captures the diffusion of the SF secretome inside the CT and its transmission to the BM, at the other end of which these enzymes may interact with MT1-MMP located on the membrane of tumor cells, and activate MMP-2. The inclusion of these dynamics in the mathematical model enables the study of the partnerships that exist between tumor cells and helping cells in the stroma, which are nowadays well known to aid cancer invasion. While this framework, including fibroblast facilitating ECM degradation by cancer cells, was partially considered in [66] for a Tumour Chamber Invasion Assay, our model captures better the dimensionality of the problem by considering the BM to be a $d-1$ dimensional domain, where $d$ is the dimension of the CT, as similarly done in [41]. Altogether, our model allows for a more complete and biophysically consistent representation of the phenomenon, and the results of numerical simulations highlight the importance of considering these pro-tumour cells located in the stroma, as the digestion of the BM appears faster closer to a SF rich CT area.

Compared to the aforementioned works and previous mathematical models of tumor invasion (see e.g. [17, 18, 61, 88]), our model also comprises a comprehensive description of the biochemical reactions regulating the activation and inhibition of MMP-2, capturing the dual role of TIMP-2 as only previously seen in dedicated studies ignoring spatial dynamics [62, 98]. We remark that such dual role of TIMP-2 is still captured in the reduced model that we derived using typical simplifying assumptions, as demonstrated by the results of our steady state analysis, summarised in Figure 3.

Moreover, additional efforts were made in this work to extract parameter values from the literature, deriving several missing ones from steady state assumptions and reference enzyme concentrations reported in empirical studies. Hence, except for 2 of them (the production rates of enzymes by SF, i.e. $r_p, r_t$), most of the parameters are directly found from experimental studies or previous modelling works that fit these parameters using experimental data. The remaining parameters are obtained from simple computations: starting from empirical measurements of physiological enzyme concentrations, we compute the remaining parameters assuming that the physiological conditions correspond to steady-states in the reaction-diffusion equations. This procedure allowed us to find biologically relevant parameter values. We conducted a sensitivity analysis to estimate how errors on these parameter values could influence the output of the model and that of a reduced model obtained under typical simplifying assumptions. This latter part of the work is crucial, as it also highlights changes in the model's sensitivity upon reduction, which is not very often considered for invasion models.

**The results of numerical simulations are in good qualitative agreement with the current biological understanding of the phenomenon.** The flexibility of the model allowed us to explore the behavior of the system on the relevant test cases, reproducible under different parameter sets, via numerical simulations. The model predicts that the presence of an epithelial tumour may result in slow
integrating information on the physical properties of the BM in this parameter. Specifically, it would be relevant to use a diffusion coefficient dependent on the tortuosity of the BM, which is related to the pore sizes of the ECM network within the BM, estimated to be between 10 and 130 nm in diameter [33]. This is particularly relevant in our modelling framework because as the BM is degraded, the size of the pores increase, ultimately increasing the protein diffusivity within the BM. This effect could be implemented in our model by considering a BM density dependent transmission coefficient $\kappa_i(M) = K_i(M) \frac{D_i}{\varepsilon}$ with $0 \le K_i(M) < 1$ being the partitioning coefficient (as in [24]).

Secondly, in our model, MT1-MMP turnover is modelled as a combination of the production of monomeric MT1-MMP by cancer cells, with rate $\alpha_m$, and the internalization of monomeric, dimeric and complexes involving MT1-MMPs, via decay terms with rates $\beta_i$ ($i = m, d, 1, 2, 3$). While this was sufficiently detailed for the purpose of this first work, it might be an overly simplified representation of MT1-MMP turnover. As MT1-MMP turnover also plays an important role in subsequent steps of invasion in the CT, it might be relevant to adapt the model proposed in [30] to our setting.

In the present work, we only consider the digestion of the BM by MMP-2 enriched in SF secretome and activated by tumor cells. Nonetheless, degradation of the BM may also be mediated by MT1-MMPs or by other MMPs directly produced by tumor cells. This process was neglected for the moment to discriminate the effect of SFs located in the CT in the digestion of the BM. The addition of a source term of proMMP, of the form of $\alpha_p\rho_0$ ($\alpha_p$ corresponding to the rate of production of proMMP by the tumor cells, having unit nmol/cell/s) in the equation for the proMMP-2 in the BM may be sufficient to include this phenomenon.

Furthermore, while we only consider the action of MMP-2 in this work, it may be relevant to extend this framework to include the dynamics of other MMP subtypes. For instance, MMP-9 may play an important role in the digestion of the BM and the invasion process (see [31] and references therein), even though no relevant amount of MMP-9 was measured in SF secretome in [45]. In this natural extension of our work, MMP-9 dynamics may be captured with a modelling strategy analogous to that used for MMP-2.

Finally, as already demonstrated by experimental observations, SFs' secretome plays a role not only in the digestion of the BM, but also in the subsequent invasion of carcinoma cells in the CT [45]. It may, in particular, promote epithelial-to-mesenchymal transition of carcinoma cells, enhancing their motility capabilities. Our next efforts will be devoted to the development of a model of invasion of carcinoma cells inside the CT, and we will focus on capturing the cancer cell motility-boosting effect of SF secretome. Coupling a model for invasion inside the CT and the present model will lead to a complete mathematical representation of the partnership between carcinoma cells and the SFs during invasion.

**Interdisciplinary perspectives.** A limitation of our work was our ability to find precise parameter values that correspond to the specific pathophysiological context under study. Although a thorough calibration from the literature was conducted, some parameters are not found in the current literature or originate from different contexts, not specifically measured using epithelial samples. For example, reference enzyme concentrations used in this work were measured in plasma serum, due to a lack of measurements in tissues. Most probably, concentrations differ from one tissue to another. The importance of this limitation is highlighted by the results of our sensitivity analysis, which demonstrate that uncertainties on limited parameter values could tremendously impact the simulations' outputs. This will be a crucial point for our future work, where we will compare our model's numerical results to the results of biological experiments. We aim to accurately replicate the experimental setup *in silico* to achieve quantitative comparisons.

Therefore our numerical model needs to be fed by data gathered from biological models mimicking skin carcinoma progression, which will allow us to harvest critical information on specific factors involved in early cancer progression. Specifically, organotypic 3D models of skin carcinoma progression are being developed to allow us to gather data on the impact of cellular and matricial densities, proximity of stroma cells and their secretory activity in the senescent context. Moreover, data on protein concentrations and gradients analyzed in human biopsy specimen from the clinic will also be relevant, in order to identify the histological expression patterns of critical molecules associated with different stages of cancer progression and invasion features, or in samples from patients of different ages.

Data from these experiments, and from patients, will nourish the mathematical model with quantitative information that will improve the quality and specificity of parameter values upon empirical calibration. Therefore, the concentrations of critical actors identified in the model from the sensitivity analysis will be measured in the biological model in order to replace those inferred from related literature. This will help improve the predictive capability of the model and allow the simulation of scenarios in which molecules or processes are targeted to interfere with cancer progression. In summary, the mathematical and biological model will mutually feed one another, in a synergistic numerical and experimental approach.
BM digestion at the tumour location, as opposed to the case of a healthy tissue where the BM remains intact, thanks to the presence of MT1-MMPs on the surface of tumour cells which initiate the MMP-2 activation cascade. Moreover, our results clarify that SF secretome may aggravate the imbalance between the activation and inhibition dynamics of MMP-2, boosting cancer invasion. In fact, the inclusion of SFs in the CT results in the emergence of more accurate gradients of proMMP-2 in the CT, extending radially from the position of the SF, and faster digestion of the BM localized at the points closest to the SF. It is at this point that we expect to first observe BM rupture, *i.e.* the decrease of BM density below a threshold allowing cells to pass through the BM, and the initiation of cancer invasion. Our results also highlight that the possibility of observing BM rupture within a given time frame increases as the distance between the SF and the BM decreases, as the rate of proMMP-2 secretion by the SF increases and as its TIMP-2 secretion rate decreases. These secretion rates are hard to estimate as they may vary across experimental conditions, tissue samples and among SF population.

While these results are in good qualitative agreement with the current biological understanding of the phenomenon, we remark that the amount of time after which the model predicts BM rupture may not accurately convey the empirical timescale of this phenomenon, which may be faster. First of all, the critical BM density for cancer invasion, which we arbitrarily set to half the carrying capacity, may be higher. Secondly, we considered the presence of only one SF in the CT for illustrative purposes – although the model is flexible enough to represent numerous SFs in different locations within the CT – and *in vivo* the number of SFs may be significantly larger, yielding an enhanced effect of combined SFs’ secretome. Lastly, despite the great efforts made in identifying appropriate parameter values, these are taken or estimated from a variety of works in the literature, and model calibration with experimental data is necessary to ensure a wholly coherent parameter set.

**Sensitivity analysis reveals that model prediction accuracy depends on knowledge of the amount of MT1-MMPs that epithelial tumor cell express at the invadopodia, especially after model reduction.** Given the large number of parameters in the full model, we derived a reduced version of the model under biologically motivated simplifying assumptions on the relative scales of the dynamics prescribed within the MMP-2 activation cascade. We remark that these assumptions were further justified by the results of numerical simulations, during which concentrations of soluble molecules reached equilibrium much faster than the BM density. Model reduction not only lowers the degrees of freedom of the system, facilitating future experimental calibration of the model, but also allows us to draw analytical conclusions on the equilibrium of the system and necessary conditions for sustained BM rupture. Moreover, it allows for a comparison with previous cancer invasion mathematical models that included some MMP-2 activation or inhibition dynamics under similar assumptions.

In view of this calibration, we performed GSA to investigate the models’ robustness to uncertainty in the parameter values and identify the most influential parameters on the speed of BM digestion by cancer cells in the presence of a SF in the CT. The most influential parameters of both versions of the model displayed both linear and non-linear effects, indicating uncertainty in these inputs may influence the model output also through its interaction with other parameters, as expected given the complexity of the equations and dynamics prescribed in the model. Notably, we found that the most influential parameters of the full model have lower EEs mean and standard deviation than those of the reduced model. This suggests that the predictive capability of the full model is more robust to uncertainty of parameter values than the reduced one, and it is therefore better suited for *in silico* investigations when the parameters are not accurately measured, despite the higher computational cost.

From the GSA we found that uncertainty in the parameters associated with the amount of MT1-MMPs available on the surface of cancer cells – *i.e.* tumour cell number, MT1-MMP expression rate by tumour cells, internalisation rate of MT1-MMP dimers and complexes involving it – has the largest impact on the output of the full model. This was followed by an important role played by the rate of the irreversible reaction at the end of the MMP-2 activation cascade. Interestingly, we found that this latter parameter had negligible effect on the output of the reduced model, which may be an advantage in the absence of accurate estimates for this rate. However, this advantage is counterbalanced by the intensified effect of the parameter modelling the fixed concentration of free dimeric MT1-MMPs, which only appears in the reduced model capturing the collective effect of the MT1-MMP dynamics on the surface of cancer cells. Such simplification is the main drawback of the reduced model, as enormous importance if given to a parameter which is even more difficult to measure empirically.
Carcinoma invasive growth also comprises other factors such as inflammation, collective and individual migration of tumor cells, mechanical effects such as pressure, stiffening of the tissue, and partnerships existing with other cells of the stroma (e.g. macrophages, cancer-associated fibroblasts, ...). The ultimate goal of our work will be to progressively add new cellular contributions to build a complete model of carcinoma invasion that comprises more of the previously mentioned effects. Even though such a model will be highly complex, the present work demonstrates that designing and simulating an already complex subsystem, as well as characterizing a large parameter set using the literature, is manageable.

## 4.3 Conclusion

Our model shows that an enhanced digestion of the BM can occur due to the presence of SFs close to the location of epidermal carcinomas, due to the secretion of TIMP-2 and proMMP-2 by SFs and the expression of MT1-MMP by cancer cells. This finding reveals the potential role of these stromal cells in the initial steps of carcinoma invasion and, since SFs accumulate in the CT of various tissues with age, a key factor for the age-related increase in epidermal carcinoma malignancy. This suggests that these stromal cells could be important targets for therapeutic strategies aiming to block carcinoma progression, and the model indicates the proximity to the BM at the tumor location (and, likely, the number) of SFs in the CT as quantitative information that could potentially be leveraged for decision-making purposes during intervention.

<table><caption>Table 4: List of frequently used acronyms.</caption><thead><tr><td>Acronym</td><td>Meaning</td><td>Acronym</td><td>Meaning</td></tr></thead><tbody><tr><td>BM</td><td>Basement Membrane</td><td>ODE</td><td>Ordinary Differential Equation</td></tr><tr><td>CAF</td><td>Cancer-Associated Fibroblast</td><td>proMMP</td><td>MMP precursor</td></tr><tr><td>CT</td><td>Conjunctive Tissue</td><td>PDE</td><td>Partial Differential Equation</td></tr><tr><td>ECM</td><td>Extracellular Matrix</td><td>SF</td><td>Senescent Fibroblast</td></tr><tr><td>GSA</td><td>Global Sensitivity Analysis</td><td>TAM</td><td>Tumour-Associated Macrophages</td></tr><tr><td>MMP</td><td>Matrix MetalloProteinases</td><td>TIMP</td><td>Tissue-Inhibitor of MMP</td></tr><tr><td>MT1-MMP</td><td>Membrane-Type MMP</td><td>TME</td><td>Tumour Microenvironment</td></tr></tbody></table>

## Acknowledgements

This project has received funding from the European Unions Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 945298-ParisRegionFP. C.V. is a Fellow of the Paris Region Fellowship Programme, supported by the Paris Region. A.P. acknowledges support from the Labex CEMPI (ANR-11-LABX-0007-01). This work has been partially funded by the PEPS JCJC project "Conservative numerical schemes for novel structured PDE models of cancer invasion" from INSMI - CNRS.

## References

[1] A. ALBINI, Y. IWAMOTO, H. KLEINMAN, G. MARTIN, S. AARONSON, J. KOZLOWSKI, AND R. MCEWAN, *A rapid in vitro assay for quantitating the invasive potential of tumor cells*, Cancer research, 47 (1987), pp. 3239-3245.

[2] T.-M. ANNE, T.-H. TAINA, ET AL., *Levels of circulating timp-2 and mmp2-timp2 complex are decreased in squamous cervical carcinoma*, Obstetrics and Gynecology International, 2010 (2010).

[3] P. BALUK, S. MORIKAWA, A. HASKELL, M. MANCUSO, AND D. M. MCDONALD, *Abnormalities of basement membrane on blood vessels and endothelial sprouts in tumors*, The American journal of pathology, 163 (2003), pp. 1801-1815.

[4] L. BERBEN, G. FLORIS, H. WILDIERS, AND S. HATSE, *Cancer and aging: Two tightly interconnected biological processes*, Cancers, 13 (2021).
## 5.2 Modelling and interdisciplinary perspectives

**Mathematical modelling perspectives.** In this work we use mathematical modelling to represent and study the contribution of SF secretome in initial cancerous progression of carcinoma toward BM rupture and invasion. In particular, we investigate the potential impact of this partnership on the digestion of the BM, which is the first step of carcinoma invasion. While our model is able to capture numerous biological effects, there are several planned extensions that would allow the inclusion of additional biologically relevant features.

The transmission coefficient $\kappa_i (i = t, p)$ models the rate at which proteins from the CT are transmitted to and across the BM. It is computed using the diffusion coefficient $D_i$ of the protein and the width of the BM. In particular, we assumed that the diffusivity of each protein in the BM is $10^{-6}$ smaller than in the CT (see Appendix B), a coarse assumption explained by the fact that the BM is a dense sheet of collagen fibers. A more accurate description of the protein diffusivity in the BM could be achieved by integrating information on the physical properties of the BM in this parameter. Specifically, it would be relevant to use a diffusion coefficient dependent on the tortuosity of the BM, which is related to the pore sizes of the ECM network within the BM, estimated to be between 10 and 130 nm in diameter [59]. This is particularly relevant in our modelling framework because as the BM is degraded, the size of the pores increase, ultimately increasing the protein diffusivity within the BM. This effect could be implemented in our model by considering a BM density dependent transmission coefficient $\kappa_i(M) = K_i(M) \frac{D_i}{\varepsilon}$ with $0 \le K_i(M) < 1$ being the partitioning coefficient (as in [38]).

Secondly, in our model, MT1-MMP turnover is modelled as a combination of the production of monomeric MT1-MMP by cancer cells, with rate $\alpha_m$, and the internalization of monomeric, dimeric and complexes involving MT1-MMPs, via decay terms with rates $\beta_i (i = m, d, 1, 2, 3)$. While this was sufficiently detailed for the purpose of this first work, it might be an overly simplified representation of MT1-MMP turnover. As MT1-MMP turnover also plays an important role in subsequent steps of invasion in the CT, it might be relevant to adapt the model proposed in [55] to our setting.

In the present work, we only consider the digestion of the BM by MMP-2 enriched in SF secretome and activated by tumor cells. Nonetheless, degradation of the BM may also be mediated by MT1-MMPs or by other MMPs directly produced by tumor cells. This process was neglected for the moment to discriminate the effect of SFs located in the CT in the digestion of the BM. The addition of a source term of proMMP, of the form of $\alpha_p \rho_0$ ($\alpha_p$ corresponding to the rate of production of proMMP by the tumor cells, having unit nmol/cell/s) in the equation for the proMMP-2 in the BM may be sufficient to include this phenomenon.

Furthermore, while we only consider the action of MMP-2 in this work, it may be relevant to extend this framework to include the dynamics of other MMP subtypes. For instance, MMP-9 may play an important role in the digestion of the BM and the invasion process (see [57] and references therein), even though no relevant amount of MMP-9 was measured in SF secretome in [81]. In this natural extension of our work, MMP-9 dynamics may be captured with a modelling strategy analogous to that used for MMP-2.

Finally, as already demonstrated by experimental observations, SFs' secretome plays a role not only in the digestion of the BM, but also in the subsequent invasion of carcinoma cells in the CT [81]. It may, in particular, promote epithelial-to-mesenchymal transition of carcinoma cells, enhancing their motility capabilities. Our next efforts will be devoted to the development of a model of invasion of carcinoma cells inside the CT, and we will focus on capturing the cancer cell motility-boosting effect of SF secretome. Coupling a model for invasion inside the CT and the present model will lead to a complete mathematical representation of the partnership between carcinoma cells and the SFs during invasion.

**Interdisciplinary perspectives.** A limitation of our work was our ability to find precise parameter values that correspond to the specific pathophysiological context under study. Although a thorough calibration from the literature was conducted, some parameters are not found in the current literature or originate from different contexts, not specifically measured using epithelial samples. For example, reference enzyme concentrations used in this work were measured in plasma serum, due to a lack of measurements in tissues. Most probably, concentrations differ from one tissue to another. The importance of this limitation is highlighted by the results of our sensitivity analysis, which demonstrate that uncertainties on limited parameter values could tremendously impact the simulations' outputs. This will be a crucial point for our future work, where we will compare our model's numerical results to the results of biological experiments. We aim to accurately replicate the experimental setup *in silico* to achieve quantitative comparisons.

Therefore our numerical model needs to be fed by data gathered from biological models mimicking skin carcinoma progression, which will allow us to harvest critical information on specific factors involved in early cancer progression. Specifically, organotypic 3D models of skin carcinoma progression are being
[5] G. A. CABRAL-PACHECO, I. GARZA-VELOZ, C. CASTRUITA-DE LA ROSA, J. M. RAMIREZ-AHUÑA, B. A. PEREZ-ROMERO, J. F. GUERRERO-RODRIGUEZ, N. MARTINEZ-AVILA, AND M. L. MARTINEZ-FIERRO, *The roles of matrix metalloproteinases and their inhibitors in human diseases*, International Journal of Molecular Sciences, 21 (2020), p. 9739.

[6] J. M. CASTELLANO, K. I. MOSHER, R. J. ABBEY, A. A. MCBRIDE, M. L. JAMES, D. BERDNIK, J. C. SHEN, B. ZOU, X. S. XIE, M. TINGLE, ET AL., *Human umbilical cord plasma proteins revitalize hippocampal function in aged mice*, Nature, 544 (2017), pp. 488-492.

[7] M. A. CEPEDA, J. J. PELLING, C. L. EVERED, K. C. WILLIAMS, Z. FREEDMAN, I. STAN, J. A. WILLSON, H. S. LEONG, AND S. DAMJANOVSKI, *Less is more: low expression of mt1-mmp is optimal to promote migration and tumourigenesis of breast cancer cells*, Molecular cancer, 15 (2016), pp. 1-27.

[8] M. CHAPLAIN AND N. DEAKIN, *Mathematical modeling of cancer invasion: The role of membrane-bound matrix metalloproteinases*, Frontiers in Oncology, 3 (2013).

[9] M. CHAPLAIN AND G. LOLAS, *Mathematical modelling of cancer invasion of tissue: dynamic heterogeneity*, Networks and Heterogeneous Media, 1 (2006), pp. 399-439.

[10] M. A. CHAPLAIN, C. GIVERSO, T. LORENZI, AND L. PREZIOSI, *Derivation and application of effective interface conditions for continuum mechanical models of cell invasion through thin membranes*, SIAM Journal on Applied Mathematics, 79 (2019), pp. 2011-2031.

[11] M. A. CHELLAIAH, T. MA, ET AL., *Membrane localization of membrane type 1 matrix metalloproteinase by cd44 regulates the activation of pro-matrix metalloproteinase 9 in osteoclasts*, BioMed Research International, 2013 (2013).

[12] I. E. COLLIER, W. LEGANT, B. MARMER, O. LUBMAN, S. SAFFARIAN, T. WAKATSUKI, E. ELSON, AND G. I. GOLDBERG, *Diffusion of mmps on the surface of collagen fibrils: the mobile cell surface-collagen substratum interface*, PloS one, 6 (2011), p. e24029.

[13] M. DEMESTRE, G. PARKIN-SMITH, A. PETZOLD, AND A. PULLEN, *The pro and the active form of matrix metalloproteinase-9 is increased in serum of patients with amyotrophic lateral sclerosis*, Journal of neuroimmunology, 159 (2005), pp. 146-154.

[14] E. I. DERYUGINA, M. A. BOURDON, K. JUNGWIRTH, J. W. SMITH, AND A. Y. STRONGIN, *Functional activation of integrin αvβ3 in tumor cells expressing membrane-type 1 matrix metalloproteinase*, International journal of cancer, 86 (2000), pp. 15-23.

[15] E. I. DERYUGINA AND J. P. QUIGLEY, *Matrix metalloproteinases and tumor metastasis*, Cancer and metastasis reviews, 25 (2006), pp. 9-34.

[16] G. P. DIMRI, X. LEE, G. BASILE, M. ACOSTA, G. SCOTT, C. ROSKELLEY, E. E. MEDRANO, M. LINSKENS, I. RUBELJ, AND O. PEREIRA-SMITH, *A biomarker that identifies senescent human cells in culture and in aging skin in vivo.*, Proceedings of the National Academy of Sciences, 92 (1995), pp. 9363-9367.

[17] P. DOMSCHKE, D. TRUCU, A. GERISCH, AND M. A. J. CHAPLAIN, *Mathematical modelling of cancer invasion: Implications of cell adhesion variability for tumour infiltrative growth patterns*, Journal of Theoretical Biology, 361 (2014), pp. 41-60.

[18] A. DONZÉ, E. FANCHON, L. M. GATTEPAILLE, O. MALER, AND P. TRACQUI, *Robustness analysis and behavior discrimination in enzymatic reaction networks*, PloS one, 6 (2011), p. e24246.

[19] H. EMONARD, G. BELLON, L. TROEBERG, A. BERTON, A. ROBINET, P. HENRIET, E. MARBAIX, K. KIRKEGAARD, L. PATTHY, Y. EECKHOUT, ET AL., *Low density lipoprotein receptor-related protein mediates endocytic clearance of pro-mmp-2· timp-2 complex through a thrombospondin-independent mechanism*, Journal of Biological Chemistry, 279 (2004), pp. 54944-54951.

[20] H. ENDERLING AND M. AJ CHAPLAIN, *Mathematical modeling of tumor growth and treatment*, Current pharmaceutical design, 20 (2014), pp. 4934-4940.

[21] R. EYMARD, T. GALLOUËT, AND R. HERBIN, *Finite volume methods*, in Handbook of numerical analysis, Vol. VII, vol. VII of Handb. Numer. Anal., North-Holland, Amsterdam, 2000, pp. 713-1020.
developed to allow us to gather data on the impact of cellular and matricial densities, proximity of stroma cells and their secretory activity in the senescent context. Moreover, data on protein concentrations and gradients analyzed in human biopsy specimen from the clinic will also be relevant, in order to identify the histological expression patterns of critical molecules associated with different stages of cancer progression and invasion features, or in samples from patients of different ages.

Data from these experiments, and from patients, will nourish the mathematical model with quantitative information that will improve the quality and specificity of parameter values upon empirical calibration. Therefore, the concentrations of critical actors identified in the model from the sensitivity analysis will be measured in the biological model in order to replace those inferred from related literature. This will help improve the predictive capability of the model and allow the simulation of scenarios in which molecules or processes are targeted to interfere with cancer progression. In summary, the mathematical and biological model will mutually feed one another, in a synergistic numerical and experimental approach.

Carcinoma invasive growth also comprises other factors such as inflammation, collective and individual migration of tumor cells, mechanical effects such as pressure, stiffening of the tissue, and partnerships existing with other cells of the stroma (e.g. macrophages, cancer-associated fibroblasts, ...). The ultimate goal of our work will be to progressively add new cellular contributions to build a complete model of carcinoma invasion that comprises more of the previously mentioned effects. Even though such a model will be highly complex, the present work demonstrates that designing and simulating an already complex subsystem, as well as characterizing a large parameter set using the literature, is manageable.

<table><caption>Table 4: List of frequently used acronyms.</caption><thead><tr><td>Acronym</td><td>Meaning</td><td>Acronym</td><td>Meaning</td></tr></thead><tbody><tr><td>BM</td><td>Basement Membrane</td><td>ODE</td><td>Ordinary Differential Equation</td></tr><tr><td>CAF</td><td>Cancer-Associated Fibroblast</td><td>proMMP</td><td>MMP precursor</td></tr><tr><td>CT</td><td>Conjunctive Tissue</td><td>PDE</td><td>Partial Differential Equation</td></tr><tr><td>ECM</td><td>Extracellular Matrix</td><td>SF</td><td>Senescent Fibroblast</td></tr><tr><td>GSA</td><td>Global Sensitivity Analysis</td><td>TAM</td><td>Tumour-Associated Macrophages</td></tr><tr><td>MMP</td><td>Matrix MetalloProteinases</td><td>TIMP</td><td>Tissue-Inhibitor of MMP</td></tr><tr><td>MT1-MMP</td><td>Membrane-Type MMP</td><td>TME</td><td>Tumour Microenvironment</td></tr></tbody></table>

## Acknowledgements

This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 945298-ParisRegionFP. C.V. is a Fellow of the Paris Region Fellowship Programme, supported by the Paris Region. A.P. acknowledges support from the Labex CEMPI (ANR-11-LABX-0007-01). This work has been partially funded by the PEPS JCJC project "Conservative numerical schemes for novel structured PDE models of cancer invasion" from INSMI - CNRS.

## References

[1] N. J. AHERNE, A. DHAWAN, J. G. SCOTT, AND H. ENDERLING, *Mathematical oncology and its application in non melanoma skin cancer-a primer for radiation oncology professionals*, Oral oncology, 103 (2020), p. 104473.

[2] A. ALBINI, Y. IWAMOTO, H. KLEINMAN, G. MARTIN, S. AARONSON, J. KOZLOWSKI, AND R. MCEWAN, *A rapid in vitro assay for quantitating the invasive potential of tumor cells*, Cancer research, 47 (1987), pp. 3239-3245.

[3] P. M. ALTROCK, L. L. LIU, AND F. MICHOR, *The mathematics of cancer: integrating quantitative models*, Nature Reviews Cancer, 15 (2015), pp. 730-745.

[4] T.-M. ANNE, T.-H. TAINA, ET AL., *Levels of circulating timp-2 and mmp2-timp2 complex are decreased in squamous cervical carcinoma*, Obstetrics and Gynecology International, 2010 (2010).
[22] M. FANE AND A. T. WEERARATNA, *How the ageing microenvironment influences tumour progression*, Nature Reviews Cancer, 20 (2020), p. 89106.

[23] J. FERLAY, M. ERVIK, F. LAM, M. LAVERSANNE, M. COLOMBET, L. MERY, M. PIÑEROS, A. ZNAOR, I. SOERJOMATARAM, AND F. BRAY, *Global cancer observatory: Cancer today* (version 1.1), 2024. Lyon, France: International Agency for Research on Cancer. Available from: https://gco.iarc.who.int/today, accessed: (22/03/2024).

[24] W. H. FISSELL, C. L. HOFMANN, N. FERRELL, L. SCHNELL, A. DUBNISHEVA, A. L. ZYDNEY, P. D. YURCHENCO, AND S. ROY, *Solute partitioning and filtration by extracellular matrices*, American Journal of Physiology-Renal Physiology, 297 (2009), pp. F1092–F1100.

[25] O. GALLINATO, T. COLIN, O. SAUT, AND C. POIGNARD, *Tumor growth model of ductal carcinoma: from in situ phase to stroma invasion*, Journal of Theoretical Biology, 429 (2017), pp. 253–266.

[26] D. HANAHAN, *Hallmarks of cancer: new dimensions*, Cancer discovery, 12 (2022), pp. 31–46.

[27] D. HANAHAN AND R. A. WEINBERG, *The hallmarks of cancer*, cell, 100 (2000), pp. 57–70.

[28] Z. HEIDARY, J. GHAISARI, S. MOEIN, AND S. HAGHJOOY JAVANMARD, *The double-edged sword role of fibroblasts in the interaction with cancer cells; an agent-based modeling approach*, PloS one, 15 (2020), p. e0232965.

[29] P. HENRIET AND H. EMONARD, *Matrix metalloproteinase-2: Not (just) a hero of the past*, Biochimie, 166 (2019), p. 223232.

[30] D. HOSHINO, N. KOSHIKAWA, T. SUZUKI, V. QUARANTA, A. M. WeAVER, M. SEIKI, AND K. ICHIKAWA, *Establishment and validation of computational model for mt1-mmp dependent ecm degradation and intervention strategies*, PLoS computational biology, 8 (2012), p. e1002479.

[31] H. HUANG, *Matrix metalloproteinase-9 (mmp-9) as a cancer biomarker and mmp-9 biosensors: Recent advances*, Sensors, 18 (2018).

[32] B. IOOSS AND P. LEMAÎTRE, *A review on global sensitivity analysis methods*, in Uncertainty Management in Simulation-Optimization of Complex Systems: Algorithms and Applications, G. Dellino and C. Meloni, eds., Springer US, Boston, MA, 2015, pp. 101–122.

[33] R. JAYADEV AND D. R. SHERWOOD, *Basement membranes*, Current Biology, 27 (2017), pp. R207–R211.

[34] A. JIANG, K. LEHTI, X. WANG, S. J. WEISS, J. KESKI-OJA, AND D. PEI, *Regulation of membrane-type matrix metalloproteinase 1 activity by dynamin-mediated endocytosis*, Proceedings of the National Academy of Sciences, 98 (2001), pp. 13693–13698.

[35] J. JOSHI, M. D. ABNAVI, AND C. R. KOTHAPALLI, *Synthesis and secretome release by human bone marrow mesenchymal stem cell spheroids within three-dimensional collagen hydrogels: Integrating experiments and modelling*, Journal of Tissue Engineering and Regenerative Medicine, 13 (2019), pp. 1923–1937.

[36] E. D. KARAGIANNIS AND A. S. POPEL, *A theoretical model of type i collagen proteolysis by matrix metalloproteinase (mmp) 2 and membrane type 1 mmp in the presence of tissue inhibitor of metalloproteinase 2*, Journal of Biological Chemistry, 279 (2004), pp. 39105–39114.

[37] E. KIM, V. REBECCA, I. V. FEDORENKO, J. L. MESSINA, R. MATHEW, S. S. MARIA-ENGLER, D. BASANTA, K. S. SMALLEY, AND A. R. ANDERSON, *Senescent fibroblasts in melanoma initiation and progression: an integrated theoretical, experimental, and clinical approach*, Cancer research, 73 (2013), pp. 6874–6885.

[38] Y. KIM AND A. FRIEDMAN, *Interaction of tumor with its micro-environment: A mathematical model*, Bulletin of mathematical biology, 72 (2010), pp. 1029–1068.

[39] Y. KIM, J. WALLACE, F. LI, M. OSTROWSKI, AND A. FRIEDMAN, *Transformed epithelial cells and fibroblasts/myofibroblasts interaction in breast tumor: a mathematical model and experiments*, Journal of mathematical biology, 61 (2010), pp. 401–421.
[5] P. J. ASIF, C. LONGOBARDI, M. HAHNE, AND J. P. MEDEMA, *The role of cancer-associated fibroblasts in cancer invasion and metastasis*, Cancers, 13 (2021).

[6] D. J. BAKER, T. WIJSHAKE, T. TCHKONIA, N. K. LEBRASSEUR, B. G. CHILDS, B. VAN DE SLUIS, J. L. KIRKLAND, AND J. M. VAN DEURSEN, *Clearance of p16ink4a-positive senescent cells delays ageing-associated disorders*, Nature, 479 (2011), p. 232–236.

[7] P. BALUK, S. MORIKAWA, A. HASKELL, M. MANCUSO, AND D. M. MCDONALD, *Abnormalities of basement membrane on blood vessels and endothelial sprouts in tumors*, The American journal of pathology, 163 (2003), pp. 1801–1815.

[8] D. BARBOLOSI, J. CICCOLINI, B. LACARELLE, F. BARLÉSI, AND N. ANDRÉ, *Computational oncology—mathematical modelling of drug regimens for precision medicine*, Nature reviews Clinical oncology, 13 (2016), pp. 242–254.

[9] L. BARTHA AND R. EFTIMIE, *Mathematical investigation into the role of macrophage heterogeneity on the temporal and spatio-temporal dynamics of non-small cell lung cancers*, Journal of Theoretical Biology, 549 (2022), p. 111207.

[10] S. BEKISZ AND L. GERIS, *Cancer modeling: From mechanistic to data-driven approaches, and from fundamental insights to clinical applications*, Journal of Computational Science, 46 (2020), p. 101198. 20 years of computational science.

[11] I. BEN-PORATH AND R. A. WEINBERG, *When cells get stressed: an integrative view of cellular senescence*, The Journal of Clinical Investigation, 113 (2004), pp. 8–13.

[12] L. BERBEN, G. FLORIS, H. WILDIERS, AND S. HATSE, *Cancer and aging: Two tightly interconnected biological processes*, Cancers, 13 (2021).

[13] J. A. BULL AND H. M. BYRNE, *The hallmarks of mathematical oncology*, Proceedings of the IEEE, 110 (2022), pp. 523–540.

[14] G. A. CABRAL-PACHECO, I. GARZA-VELOZ, C. CASTRUITA-DE LA ROSA, J. M. RAMIREZ-ACUÑA, B. A. PEREZ-ROMERO, J. F. GUERRERO-RODRIGUEZ, N. MARTINEZ-AVILA, AND M. L. MARTINEZ-FIERRO, *The roles of matrix metalloproteinases and their inhibitors in human diseases*, International Journal of Molecular Sciences, 21 (2020), p. 9739.

[15] J. M. CASTELLANO, K. I. MOSHER, R. J. ABBEY, A. A. MCBRIDE, M. L. JAMES, D. BERDNIK, J. C. SHEN, B. ZOU, X. S. XIE, M. TINGLE, ET AL., *Human umbilical cord plasma proteins revitalize hippocampal function in aged mice*, Nature, 544 (2017), pp. 488–492.

[16] M. A. CEPEDA, J. J. PELLING, C. L. EVERED, K. C. WILLIAMS, Z. FREEDMAN, I. STAN, J. A. WILLSON, H. S. LEONG, AND S. DAMJANOVSKI, *Less is more: low expression of mt1-mmp is optimal to promote migration and tumourigenesis of breast cancer cells*, Molecular cancer, 15 (2016), pp. 1–27.

[17] M. CHAPLAIN AND N. DEAKIN, *Mathematical modeling of cancer invasion: The role of membrane-bound matrix metalloproteinases*, Frontiers in Oncology, 3 (2013).

[18] M. CHAPLAIN AND G. LOLAS, *Mathematical modelling of cancer invasion of tissue: dynamic heterogeneity*, Networks and Heterogeneous Media, 1 (2006), pp. 399–439.

[19] M. A. CHAPLAIN, C. GIVERSO, T. LORENZI, AND L. PREZIOSI, *Derivation and application of effective interface conditions for continuum mechanical models of cell invasion through thin membranes*, SIAM Journal on Applied Mathematics, 79 (2019), pp. 2011–2031.

[20] M. A. CHELLAIAH, T. MA, ET AL., *Membrane localization of membrane type 1 matrix metalloproteinase by cd44 regulates the activation of pro-matrix metalloproteinase 9 in osteoclasts*, BioMed Research International, 2013 (2013).

[21] G. CIAVOLELLA, N. DAVID, AND A. POULAIN, *Effective interface conditions for a porous medium type problem*, Interfaces Free Bound., (2024).