[40] P. LAUDAŃSKI, J. SWIATECKA, L. KOZŁOWSKI, M. LEŚNIEWSKA, M. WOJTUKIEWICZ, AND S. WOŁCZYŃSKI, *Increased serum level of membrane type 1-matrix metalloproteinase (mt1-mmp/mmp-14) in patients with breast cancer.*, Folia histochemica et cytobiologica, 48 (2010), pp. 101–103.

[41] F. LAURINO AND P. ZUNINO, *Derivation and analysis of coupled PDEs on manifolds with high dimensionality gap arising from topological model reduction*, ESAIM Math. Model. Numer. Anal., 53 (2019), pp. 2047–2080.

[42] A. LEINONEN, M. MARIYAMA, T. MOCHIZUKI, K. TRYGGVASON, AND S. T. REEDERS, *Complete primary structure of the human type iv collagen alpha 4(iv) chain. comparison with structure and expression of the other alpha (iv) chains*, Journal of Biological Chemistry, 269 (1994), pp. 26172–26177.

[43] K. V. LU, K. A. JONG, A. K. RAJASEKARAN, T. F. CLOUGHESY, AND P. S. MISCHEL, *Upregulation of tissue inhibitor of metalloproteinases (timp)-2 promotes matrix metalloproteinase (mmp)-2 activation and cell invasion in a human glioblastoma cell line*, Laboratory investigation, 84 (2004), pp. 8–20.

[44] K. M. MAK AND R. MEI, *Basement membrane type iv collagen and laminin: An overview of their biology and value as fibrosis biomarkers of liver disease*, The Anatomical Record, 300 (2017), pp. 1371–1390.

[45] N. MALAQUIN, C. VERCAMER, F. BOUALI, S. MARTIEN, E. DERUY, N. WERNERT, M. CHWASTYNIAK, F. PINET, C. ABBADIE, AND A. POURTIER, *Senescent fibroblasts enhance early skin carcinogenic events via a paracrine mmp-par-1 axis*, PLOS ONE, 8 (2013), pp. 1–17.

[46] S. MARUŠIĆ AND E. MARUŠIĆ-PALOKA, *Two-scale convergence for thin domains and its applications to some lower-dimensional models in fluid mechanics*, Asymptot. Anal., 23 (2000), pp. 23–57.

[47] S. MOMI, E. FALCINELLI, S. GIANNINI, L. RUGGERI, L. CECCETTI, T. CORAZZI, C. LIBERT, AND P. GRESELE, *Loss of matrix metalloproteinase 2 in platelets reduces arterial thrombosis in vivo*, Journal of Experimental Medicine, 206 (2009), pp. 2365–2379.

[48] M. D. MORRIS, *Factorial sampling plans for preliminary computational experiments*, Technometrics, 33 (1991), pp. 161–174.

[49] M. W. OLSON, D. C. GERVASI, S. MOBASHERY, AND R. FRIDMAN, *Kinetic analysis of the binding of human matrix metalloproteinase-2 and -9 to tissue inhibitor of metalloproteinase (timp)-1 and timp-2*, The Journal of Biological Chemistry, 272 (1997), p. 29975–29983.

[50] T. K. PALANIAPPAN, L. ŠLEKIENĖ, A.-K. JONASSON, J. GILTHORPE, AND L. GUNHAGA, *Camdelam: an in vivo approach to visualize and quantify the delamination and invasion capacity of human cancer cells*, Scientific Reports, 10 (2020), p. 10472.

[51] J. PANOVSKA AND G. PETTET, *Mathematical model for glioma cell invasion*. Unpublished, 2007.

[52] K. PAULA, T.-M. ANNE, AND T.-H. TAINA, *The sample type used affects the levels of gelatinases (mmp-2 and -9) and their inhibitors (timp-1 and -2) in circulating blood of healthy controls and breast cancer patients*, Biomarker Insights, 2 (2007), p. 117–127.

[53] F. PIANOSI, K. BEVEN, J. FREER, J. W. HALL, J. ROUGIER, D. B. STEPHENSON, AND T. WAGENER, *Sensitivity analysis of environmental models: A systematic review with practical workflow*, Environmental Modelling & Software, 79 (2016), pp. 214–232.

[54] F. PIANOSI, F. SARRAZIN, AND T. WAGENER, *A matlab toolbox for global sensitivity analysis*, Environmental Modelling & Software, 70 (2015), pp. 80–85.

[55] M. PIIPPONEN, P. RIIHILÄ, L. NISSINEN, AND V.-M. KÄHÄRI, *The role of p53 in progression of cutaneous squamous cell carcinoma*, Cancers, 13 (2021), p. 4507.

[56] B. RIBBA, O. SAUT, T. COLIN, D. BRESCH, E. GRENIER, AND J. P. BOISSEL, *A multiscale mathematical model of avascular tumor growth to investigate the therapeutic benefit of anti-invasive agents*, J. Theoret. Biol., 243 (2006), pp. 532–541.
[22] 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.

[23] J.-P. COPPE, M. BOYSEN, C. H. SUN, B. J. WONG, M. K. KANG, N.-H. PARK, P.-Y. DESPREZ, J. CAMPISI, AND A. KRTOLICA, *A Role for Fibroblasts in Mediating the Effects of Tobacco-Induced Epithelial Cell Growth and Invasion*, Molecular Cancer Research, 6 (2008), pp. 1085–1098.

[24] J.-P. COPPÉ, P.-Y. DESPREZ, A. KRTOLICA, AND J. CAMPISI, *The senescence-associated secretory phenotype: The dark side of tumor suppression*, Annual Review of Pathology: Mechanisms of Disease, 5 (2010), pp. 99–118. PMID: 20078217.

[25] 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.

[26] 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.

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

[28] 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.

[29] 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.

[30] 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.

[31] R. EFTIMIE AND L. GIBELLI, *A kinetic theory approach for modelling tumour and macrophages heterogeneity and plasticity during cancer progression*, Mathematical Models and Methods in Applied Sciences, 30 (2020), pp. 659–683.

[32] 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.

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

[34] 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.

[35] D. FABIAN AND T. FLATT, *The evolution of aging*, Nature Education Knowledge, 3 (2011), pp. 1–10.

[36] M. FANE AND A. T. WEERARATNA, *How the ageing microenvironment influences tumour progression*, Nature Reviews Cancer, 20 (2020), p. 89–106.

[37] 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).

[38] 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.
[57] E. ROIG-ROSELLO AND P. ROUSSELLE, *The human epidermal basement membrane: a shaped and cell instructive platform that aging slowly alters*, Biomolecules, 10 (2020), p. 1607.

[58] T. SAITOU, K. ITANO, D. HOSHINO, N. KOSHIKAWA, M. SEIKI, K. ICHIKAWA, AND T. SUZUKI, *Control and inhibition analysis of complex formation processes*, Theoretical Biology and Medical Modelling, 9 (2012), pp. 1–18.

[59] N. SFAKIANAKIS AND M. A. J. CHAPLAIN, *Mathematical modelling of cancer invasion: A review*, in Methods of Mathematical Oncology, T. Suzuki, C. Poignard, M. Chaplain, and V. Quaranta, eds., vol. 370, Springer Singapore, Singapore, 2021, p. 153172.

[60] L. F. SHAMPINE AND M. W. REICHELT, *The matlab ode suite*, SIAM Journal on Scientific Computing, 18 (1997), pp. 1–22.

[61] K.-I. SHIMOKAWA, M. KATAYAMA, Y. MATSUDA, H. TAKAHASHI, I. HARA, H. SATO, AND S. KANEKO, *Matrix metalloproteinase (MMP)-2 and MMP-9 activities in human seminal plasma*, Molecular Human Reproduction, 8 (2002), pp. 32–36.

[62] W. G. STETLER-STEVENSON, *The tumor microenvironment: regulation by mmp-independent effects of tissue inhibitor of metalloproteinases-2*, Cancer and Metastasis Reviews, 27 (2008), pp. 57–66.

[63] M. TAURO AND C. C. LYNCH, *Cutting to the chase: how matrix metalloproteinase-2 activity controls breast-cancer-to-bone metastasis*, Cancers, 10 (2018), p. 185.

[64] M. TOth, M. M. BERNARDO, D. C. GERVASI, P. D. SOLOWAY, Z. WANG, H. F. BIGG, C. M. OVERALL, Y. A. DeCLERCK, H. TSCHESCHE, M. L. CHER, et al., *Tissue inhibitor of metalloproteinase (timp)-2 acts synergistically with synthetic matrix metalloproteinase (mmp) inhibitors but not with timp-4 to enhance the (membrane type 1)-mmp-dependent activation of pro-mmp-2*, Journal of Biological Chemistry, 275 (2000), pp. 41415–41423.

[65] T. H. VU AND Z. WERB, *Matrix metalloproteinases: effectors of development and normal physiology*, Genes & development, 14 (2000), pp. 2123–2133.

[66] A. WATANABE, D. HOSINO, N. KOSHIKAWA, M. SEIKI, T. SUZUKI, AND K. ICHIKAWA, *Critical role of transient activity of mt1-mmp for ecm degradation in invadopodia*, PLOS Computational Biology, 9 (2013), pp. 1–13.

# A Mathematical notation and derivation of the model

## A.1 Original system

**Notation.** We consider a general function $u(x, t)$ depending on space $x \in \mathbb{R}^d$ ($d = 1, 2, 3$, we choose $d=3$ in this section) and time $t \in [0, T]$ ($T > 0$ a positive finite time horizon). We denote by $\partial_t u$ the partial derivative of $u$ with respect to time and $\nabla u = (\partial_{x_1} u, \partial_{x_2} u, \partial_{x_3} u)$ the gradient of $u$ in space. We denote by $\Delta u = \text{div}(\nabla u) = \partial^2_{x_1 x_1} u + \partial^2_{x_2 x_2} u + \partial^2_{x_3 x_3} u$ the laplacian of $u$, with $\text{div}(\vec{v}) = \partial_{x_1} v_1 + \partial_{x_2} v_2 + \partial_{x_3} v_3$ the divergence of the vector field $\vec{v} = (v_1, v_2, v_3)$.

**Full system.** In this section, we introduce the mathematical model comprising a system of PDEs describing the dynamics of the SFs' secretome in the CT, with boundary conditions modelling the transmission of these concentrations to the BM, and a system of PDEs describing the activation and degradation dynamics in the BM. Details about the modelling assumptions and the notations for the different unknowns can be found in Section 2.

We start with two domains, *i.e.* smooth bounded open subsets, $\Omega_{\text{CT}}, \Omega_{\text{BM}} \subset \mathbb{R}^d$, separated by a “membrane” denoted by $\Gamma$. We highlight that $\Gamma$ simply models the boundary separating the two domains and
[39] C. FRANTZ, K. M. STEWART, AND V. M. WEAVER, *The extracellular matrix at a glance*, Journal of Cell Science, 123 (2010), pp. 4195-4200.

[40] A. FRIEDMAN, C. HUANG, AND J. YONG, *Effective permeability of the boundary of a domain*, Communications in Partial Differential Equations, 20 (1995), pp. 59-102.

[41] 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.

[42] R. A. GATENBY AND E. T. GAWLINSKI, *A reaction-diffusion model of cancer invasion*, Cancer research, 56 (1996), pp. 5745-5753.

[43] J. GIROUD, I. BOURIEZ, H. PAULUS, A. POURTIER, F. DEBACQ-CHAINIAUX, AND O. PLUQUET, *Exploring the communication of the sasp: Dynamic, interactive, and adaptive effects on the microenvironment*, International Journal of Molecular Sciences, 24 (2023).

[44] C. GIVERSO, T. LORENZI, AND L. PREZIOSI, *Effective interface conditions for continuum mechanical models describing the invasion of multiple cell populations through thin membranes*, Applied Mathematics Letters, 125 (2022), p. 107708.

[45] L. GROSSE, N. WAGNER, A. EMELYANOV, C. MOLINA, S. LACAS-GERVAIS, K.-D. WAGNER, AND D. V. BULAVIN, *Defined p16high senescent cell types are indispensable for mouse healthspan*, Cell Metabolism, 32 (2020), pp. 87-99.e6.

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

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

[48] ——, *Hallmarks of cancer: the next generation*, cell, 144 (2011), pp. 646-674.

[49] Y. HASSONA, N. CIRILLO, K. HEESOM, E. PARKINSON, AND S. S. PRIME, *Senescent cancer-associated fibroblasts secrete active mmp-2 that promotes keratinocyte dis-cohesion and invasion*, British journal of cancer, 111 (2014), pp. 1230-1237.

[50] 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.

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

[52] T. HERNANDEZ-BOUSSARD, P. MACKLIN, E. J. GREENSPAN, A. L. GRYSHUK, E. STAHLBERG, T. SYEDA-MAHMOOD, AND I. SHMULEVICH, *Digital twins for predictive oncology will be a paradigm shift for precision cancer care*, Nature medicine, 27 (2021), pp. 2065-2066.

[53] F. HILLEN AND A. W. GRIFFIOEN, *Tumour vascularization: sprouting angiogenesis and beyond*, Cancer and Metastasis Reviews, 26 (2007), pp. 489-502.

[54] U. B. HOFMANN, J. R. WESTPHAL, G. N. VAN MUIJEN, AND D. J. RUITER, *Matrix metalloproteinases in human melanoma*, Journal of Investigative Dermatology, 115 (2000), pp. 337-344.

[55] 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.

[56] W. HU, L. FANG, R. NI, H. ZHANG, AND G. PAN, *Changing trends in the disease burden of non-melanoma skin cancer globally from 1990 to 2019 and its predicted level in 25 years*, BMC Cancer, 22 (2022).

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

[58] 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.
not a real physical membrane with a positive thickness. We have the following system of equations

$$ \partial_t \bar{c}_t = \bar{D}_t \Delta \bar{c}_t + \bar{S}_t(x) - \beta_t \bar{c}_t, \quad \text{in } \Omega_{\text{CT}} \times (0, T), \quad \text{(A.1)} $$

$$ \partial_t \bar{c}_p = \bar{D}_p \Delta \bar{c}_p + \bar{S}_p(x) - \beta_p c_p, \quad \text{in } \Omega_{\text{CT}} \times (0, T), \quad \text{(A.2)} $$

$$ \bar{c}_t = c_t, \quad \text{on } \Gamma \times (0, T), \quad \text{(A.3)} $$

$$ \bar{c}_p = c_p, \quad \text{on } \Gamma \times (0, T), \quad \text{(A.4)} $$

$$ \bar{D}_t \nabla \bar{c}_t \cdot \nu = D_t^{\text{BM}} \nabla c_t \cdot \nu, \quad \text{on } \Gamma \times (0, T), \quad \text{(A.5)} $$

$$ \bar{D}_p \nabla \bar{c}_p \cdot \nu = D_p^{\text{BM}} \nabla c_p \cdot \nu, \quad \text{on } \Gamma \times (0, T), \quad \text{(A.6)} $$

$$ \partial_t c_t = D_t^{\text{BM}} \Delta c_t + R_t(x, t) - \beta_t c_t, \quad \text{in } \Omega_{\text{BM}} \times (0, T), \quad \text{(A.7)} $$

$$ \partial_t c_p = D_p^{\text{BM}} \Delta c_p + R_p(x, t) - \beta_p c_p, \quad \text{in } \Omega_{\text{BM}} \times (0, T), \quad \text{(A.8)} $$

with $\nu$ the unit normal vector to $\Gamma$ that points from the BM to the CT. The rest of the boundary conditions on $\partial\Omega_{\text{CT}} \setminus \Gamma \times (0, T)$ and $\partial\Omega_{\text{BM}} \setminus \Gamma \times (0, T)$ are set to zero-flux boundary conditions.

Equations (A.3)-(A.6) are transmission conditions modeling the diffusion of TIMP-2 and proMMP-2 between the CT and the BM. In particular, Equations (A.3)-(A.4) model continuity of the concentration at the border between the two domains, and equations (A.5)-(A.6) model the continuity of the flux across this interface. Moreover, while the terms $\bar{S}_t(x)$ and $\bar{S}_p(x)$ in equations (A.1)-(A.2) are as defined in (2.15), the terms $R_t(x, t)$ and $R_p(x, t)$ in equations (A.7)-(A.8) include the biochemical reactions terms. Namely, we have

$$ R_t(x, 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}), \quad \text{(A.9)} $$

$$ R_p(x, t) = -k_2 c_p (c_t + c_1) + k_{-2} (c_{tp} + c_2). \quad \text{(A.10)} $$

Equations (A.1)-(A.8) are coupled, via the reaction terms $R_i(t, x)$, $i = \{t, p\}$, with the rest of the system modelling proMMP-2 activation, BM digestion and TIMP-mediated MMP inhibition: assuming the diffusion of the rest of the components is negligible, we have

$$ \partial_t M = -\gamma c_a \frac{M}{K_m + M} + r_M \left(1 - \frac{M}{M_{\text{max}}}\right)_+, \quad \text{in } \Omega_{\text{BM}} \times (0, T), $$

$$ \partial_t c_m = \alpha_m \rho^0 - k_0 c_m^2 + k_{-0} c_d - \beta_m c_m, \quad \text{in } \Omega_{\text{BM}} \times (0, T), $$

$$ \partial_t c_d = k_0 c_m^2 - k_{-0} c_d - k_1 c_d c_t + k_{-1} c_1 - \beta_d c_d, \quad \text{in } \Omega_{\text{BM}} \times (0, T), $$

$$ \partial_t c_a = k_3 c_2 - k_2 c_t c_a + k_{-2} c_{ta} - \beta_a c_a, \quad \text{in } \Omega_{\text{BM}} \times (0, T), $$

$$ \partial_t c_1 = 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 \text{in } \Omega_{\text{BM}} \times (0, T), $$

$$ \partial_t c_2 = k_2 c_1 c_p - (k_{-2} + k_3)c_2 - \beta_2 c_2, \quad \text{in } \Omega_{\text{BM}} \times (0, T), $$

$$ \partial_t c_3 = k_1 c_1 c_t - k_{-1} c_3 - \beta_3 c_2, \quad \text{in } \Omega_{\text{BM}} \times (0, T), $$

$$ \partial_t c_{tp} = k_2 c_t c_p - k_{-2} c_{tp} - \beta_{tp} c_{tp}, \quad \text{in } \Omega_{\text{BM}} \times (0, T), $$

$$ \partial_t c_{ta} = k_2 c_t c_a - k_{-2} c_{ta} - \beta_{ta} c_{ta}, \quad \text{in } \Omega_{\text{BM}} \times (0, T). $$

## A.2 Dimension reduction of the BM

In order to improve computational tractability of our model, as the BM is a very thin structure compared to the rest of the tissue [3], we reduce the dimension of $\Omega_{\text{BM}}$, which then becomes a subset of $\mathbb{R}^{d-1}$ (whereas $\Omega_{\text{CT}} \subset \mathbb{R}^d$). We do this following a simple and formal dimension reduction procedure, inspired by ideas from [41] to derive 1D-3D coupled problems.

**Smallness assumptions and decomposition.** We introduce a small parameter $0 < \varepsilon \ll 1$ describing the thickness of the BM, i.e. $\Omega_{\text{BM}} = \Gamma \times [0, \varepsilon]$. Assuming $\Gamma$ is a flat surface, we consider that $\forall x \in \Omega_{\text{BM}}$, $x$ is defined by $x := (x_1, x_2, x_3)$ with $(x_1, x_3) \in \Gamma$ and $x_2 = \varepsilon \bar{x}_2$, $(\bar{x}_2 \in (0, 1))$, coherently with the coordinate system adopted in Section 2 and illustrated in Figure 2. We decompose the concentrations $c_i$ as

$$ c_i(x_1, x_2, x_3, t) = \hat{c}_i(x_1, x_3, t) + \tilde{c}_i(x_1, x_2, x_3, t) \quad \text{for } i = \{t, p, d, tp, ta, 1, 2, 3\}, \quad \text{(A.11)} $$

where $\hat{c}_i(x_1, x_3, t)$ corresponds to the average of the concentration in the $x_2$ direction, i.e.

$$ \hat{c}_i = \frac{1}{\varepsilon} \int_0^\varepsilon c_i \, dx_2 = \int_0^1 c_i(\varepsilon \bar{x}_2) \, d\bar{x}_2, \quad \text{(A.12)} $$
[59] R. JAYADEV AND D. R. SHERWOOD, *Basement membranes*, Current Biology, 27 (2017), pp. R207–R211.

[60] 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.

[61] 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.

[62] 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.

[63] O. KEDEM AND A. KATCHALSKY, *Thermodynamic analysis of the permeability of biological membranes to non-electrolytes*, Biochimica et biophysica Acta, 27 (1958), pp. 229–246.

[64] 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.

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

[66] 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.

[67] F. KLEINHANS, *Membrane permeability modeling: Kedem–katchalsky vs a two-parameter formalism*, Cryobiology, 37 (1998), pp. 271–289.

[68] H. KNÚTSDÓTTIR, E. PÁLSSON, AND L. EDELSTEIN-KESHET, *Mathematical model of macrophage-facilitated breast cancer cells invasion*, Journal of theoretical biology, 357 (2014), pp. 184–199.

[69] A. KRTOLICA, S. PARRINELLO, S. Lockett, P.-Y. DESPREZ, AND J. CAMPISI, *Senescent fibroblasts promote epithelial cell growth and tumorigenesis: A link between cancer and aging*, Proceedings of the National Academy of Sciences, 98 (2001), pp. 12072–12077.

[70] J. K. KULAR, S. BASU, AND R. I. SHARMA, *The extracellular matrix: Structure, composition, age-related differences, tools for analysis and applications for tissue engineering*, Journal of Tissue Engineering, 5 (2014), p. 204173141455711.

[71] E. LACONI, F. MARONGIU, AND J. DEGREGORI, *Cancer as a disease of old age: changing mutational and microenvironmental landscapes*, British journal of cancer, 122 (2020), pp. 943–952.

[72] H. LARONHA AND J. CALDEIRA, *Structure and function of human matrix metalloproteinases*, Cells, 9 (2020), p. 1076.

[73] P. LAUDAŃSKI, J. SWIATECKA, L. KOZŁOWSKI, M. LEŚNIEWSKA, M. WOJTUKIEWICZ, AND S. WOŁCZYŃSKI, *Increased serum level of membrane type 1-matrix metalloproteinase (mt1-mmp/mmp-14) in patients with breast cancer*, Folia histochemica et cytobiologica, 48 (2010), pp. 101–103.

[74] F. LAURINO AND P. ZUNINO, *Derivation and analysis of coupled PDEs on manifolds with high dimensionality gap arising from topological model reduction*, ESAIM Math. Model. Numer. Anal., 53 (2019), pp. 2047–2080.

[75] A. LEINONEN, M. MARIYAMA, T. MOCHIZUKI, K. TRYGGVASON, AND S. T. REEDERS, *Complete primary structure of the human type iv collagen alpha 4(iv) chain. comparison with structure and expression of the other alpha (iv) chains*, Journal of Biological Chemistry, 269 (1994), pp. 26172–26177.
and $\tilde{c}_i(x_1, x_2, x_3, t)$ is a fluctuation around this average. Integrating both sides of the decomposition (A.11) along the width of the BM (i.e. the $x_2$ direction) and using the definition (A.12) of the average $\hat{c}_i$, we obtain that the fluctuation must satisfy

$$ \int_0^\epsilon \tilde{c}_i \, dx_2 = 0. $$

We assume that the fluctuations are small and that the concentration fields evolve smoothly within the BM, i.e. spatial derivatives of the fluctuations in the BM are small.

Under assumption (A.11) and definitions (A.9)-(A.10), we have that both sources $R_i$ ($i = \{t, p\}$) in (A.7) and (A.8) also admit a decomposition in the form (A.11).

**Simplified dynamics in the BM.** Since we assumed that the fluctuation of the concentrations within the BM – i.e. $\tilde{c}_i$ in (A.11) – are small, we simplify the dynamics in the BM by focusing on the evolution of the averages of the concentration fields in the BM along the $x_2$ direction (i.e. along the width of the BM) – i.e. $\hat{c}_i$ defined in (A.12). This will allow us to drop the dependency on $x_2$ of the dynamics, providing an approximation that is well suited for a model in which the dimensionality of the BM is reduced to $d-1$. Given that the thickness of the BM is very small, we expect that the diffusion of the enzymes in the $x_2$ direction will have a strong impact on the average concentration $\hat{c}_i$, while we assume diffusion in the $(x_1, x_3)$ plane (i.e. tangential to $\Gamma$) is negligible.

Integrating Equations (A.7)–(A.8) with respect to $x_2$, dividing by $\epsilon$, and neglecting the diffusion in the tangential plane $(x_1, x_3)$, we obtain for any $(x_1, x_3) \in \Gamma$ and $t \in (0, T)$

$$ \frac{1}{\epsilon} \int_{0}^{\epsilon} \partial_t c_i \, dx_2 - \frac{D_i^{\text{BM}}}{\epsilon} \left[ \nabla c_i \cdot \nu \bigg|_{x_2=\epsilon} - \nabla c_i \cdot \nu \bigg|_{x_2=0} \right] = \frac{1}{\epsilon} \int_{0}^{\epsilon} R_i - \beta_i c_i \, dx_2, $$

where we used the divergence theorem, with $\nu$ being the outward normal to $\partial\Omega_{\text{BM}}$. Since we assumed zero-flux boundary conditions on $\partial\Omega_{\text{BM}} \setminus \Gamma \times (0, T)$ and using the definition (A.12), we have

$$ \partial_t \hat{c}_i + \frac{D_i^{\text{BM}}}{\epsilon} \partial_{x_2} c_i(x_1, 0, x_3, t) \cdot \nu = \hat{R}_i - \beta_i \hat{c}_i, \quad (A.13) $$

for $i = \{t, p\}$. We note that the directional derivative $\partial_{x_2} c_i(x_1, 0, x_3, t) \cdot \nu$ also appears on the right-hand-side of (A.5)-(A.6).

**Approximation of concentration gradients at the interface.** We now construct an approximation of $\partial_{x_2} c_i(x_1, 0, x_3, t)$ for $i = \{t, p\}$. Since $\epsilon$ is small, assuming $c_i$ is sufficiently smooth to allow a Taylor series approximation, we consider the following expansion around $x_2 = 0$

$$ c_i(x_1, \epsilon, x_3, t) = c_i(x_1, 0, x_3, t) + \epsilon \partial_{x_2} c_i(x_1, 0, x_3, t) + \mathcal{O}(\epsilon^2). $$

Using this, we know that the following holds

$$ \partial_{x_2} c_i(x_1, 0, x_3, t) = \frac{c_i(x_1, \epsilon, x_3, t) - c_i(x_1, 0, x_3, t)}{\epsilon} + \mathcal{O}(\epsilon). $$

Using the continuity of concentrations (A.3)-(A.4) at the interface $\Gamma$, and the decomposition (A.11) at $x_2 = \epsilon$, we rewrite this as

$$ \partial_{x_2} c_i(x_1, 0, x_3, t) = \frac{\hat{c}_i(x_1, x_3, t) - \bar{c}_i(x_1, 0, x_3, t)}{\epsilon} + \frac{\tilde{c}_i(x_1, \epsilon, x_3, t)}{\epsilon} + \mathcal{O}(\epsilon). $$

We approximate $\tilde{c}_i(x_1, \epsilon, x_3, t)$ as

$$ \tilde{c}_i(x_1, \epsilon, x_3, t) = (C_{\epsilon}^i - 1)(\hat{c}_i(x_1, x_3, t) - \bar{c}_i(x_1, 0, x_3, t)), $$

with constant $C_{\epsilon}^i > 0$ for $\epsilon > 0$. Thus, for $\epsilon$ small, we have

$$ \partial_{x_2} c_i(x_1, 0, x_3) \approx C_{\epsilon}^i \frac{\hat{c}_i(x_1, x_3) - \bar{c}_i(x_1, 0, x_3)}{\epsilon} \quad \text{for } i = \{t, p\}. \quad (A.14) $$

This approximation is justified by the fact that the $c_i$ ($i = \{t, p\}$) satisfy a reaction-diffusion equation in a thin domain, with Dirichlet boundary conditions on $\Gamma$ and zero-flux boundary conditions at the other end of the domain. At equilibrium, this leads to a concentration profile along the $x_2$ direction analogous to the one shown in Figure 9, where we simulated this scenario in 1D (with $\Gamma$ the left boundary) and tested the approximation (A.14) for different values of $C_{\epsilon}^i$. We clearly have that $C_{\epsilon}^i > 0$ for $\epsilon > 0$ (analogous profiles were obtained for decreasing values of the domain size).
[76] Z. LIAO, H. L. YEO, S. W. WONG, AND Y. ZHAO, *Cellular senescence: Mechanisms and therapeutic potential*, Biomedicines, 9 (2021).

[77] K. V. LU, K. A. JONG, A. K. RAJASEKARAN, T. F. CLOUGHESY, AND P. S. MISCHEL, *Upregulation of tissue inhibitor of metalloproteinases (timp)-2 promotes matrix metalloproteinase (mmp)-2 activation and cell invasion in a human glioblastoma cell line*, Laboratory investigation, 84 (2004), pp. 8–20.

[78] I. LÜHR, A. FRIEDL, T. OVERATH, A. THOLEY, T. KUNZE, F. HILPERT, S. SEBENS, N. ARNOLD, F. RÖSEL, H.-H. OBERG, N. MAASS, C. MUNDHENKE, W. JONAT, AND M. BAUER, *Mammary fibroblasts regulate morphogenesis of normal and tumorigenic breast epithelial cells by mechanical and paracrine signals*, Cancer Letters, 325 (2012), pp. 175–188.

[79] G. MAHLBACHER, L. T. CURTIS, J. LOWENGRUB, AND H. B. FRIEBOES, *Mathematical modeling of tumor-associated macrophage interactions with the cancer microenvironment*, Journal for immunotherapy of cancer, 6 (2018), pp. 1–17.

[80] K. M. MAK AND R. MEI, *Basement membrane type iv collagen and laminin: An overview of their biology and value as fibrosis biomarkers of liver disease*, The Anatomical Record, 300 (2017), pp. 1371–1390.

[81] N. MALAQUIN, C. VERCAMER, F. BOUALI, S. MARTIEN, E. DERUY, N. WERNERT, M. CHWASTYNIAK, F. PINET, C. ABBADIE, AND A. POURTIER, *Senescent fibroblasts enhance early skin carcinogenic events via a paracrine mmp-par-1 axis*, PLOS ONE, 8 (2013), pp. 1–17.

[82] S. MARUŠIĆ AND E. MARUŠIĆ-PALOKA, *Two-scale convergence for thin domains and its applications to some lower-dimensional models in fluid mechanics*, Asymptot. Anal., 23 (2000), pp. 23–57.

[83] S. MOMI, E. FALCINELLI, S. GIANNINI, L. RUGGERI, L. CECCHETTI, T. CORAZZI, C. LIBERT, AND P. GRESELE, *Loss of matrix metalloproteinase 2 in platelets reduces arterial thrombosis in vivo*, Journal of Experimental Medicine, 206 (2009), pp. 2365–2379.

[84] M. D. MORRIS, *Factorial sampling plans for preliminary computational experiments*, Technometrics, 33 (1991), pp. 161–174.

[85] G. MURPHY, *Tissue inhibitors of metalloproteinases*, Genome biology, 12 (2011), pp. 1–7.

[86] M. W. OLSON, D. C. GERVASI, S. MOBASHERY, AND R. FRIDMAN, *Kinetic analysis of the binding of human matrix metalloproteinase-2 and -9 to tissue inhibitor of metalloproteinase (timp)-1 and timp-2*, The Journal of Biological Chemistry, 272 (1997), p. 29975–29983.

[87] M. R. OWEN AND J. A. SHERRATT, *Mathematical modelling of macrophage dynamics in tumours*, Mathematical Models and Methods in Applied Sciences, 9 (1999), pp. 513–539.

[88] T. K. PALANIAPPAN, L. ŠLEKIENĖ, A.-K. JONASSON, J. GILTHORPE, AND L. GUNHAGA, *Camdelam: an in vivo approach to visualize and quantify the delamination and invasion capacity of human cancer cells*, Scientific Reports, 10 (2020), p. 10472.

[89] J. PANOVSKA AND G. PETTET, *Mathematical model for glioma cell invasion*. Unpublished, 2007.

[90] K. PAULA, T.-M. ANNE, AND T.-H. TAINA, *The sample type used affects the levels of gelatinases (mmp-2 and -9) and their inhibitors (timp-1 and -2) in circulating blood of healthy controls and breast cancer patients*, Biomarker Insights, 2 (2007), p. 117–127.

[91] F. PIANOSI, K. BEVEN, J. FREER, J. W. HALL, J. ROUGIER, D. B. STEPHENSON, AND T. WAGENER, *Sensitivity analysis of environmental models: A systematic review with practical workflow*, Environmental Modelling & Software, 79 (2016), pp. 214–232.

[92] F. PIANOSI, F. SARRAZIN, AND T. WAGENER, *A matlab toolbox for global sensitivity analysis*, Environmental Modelling & Software, 70 (2015), pp. 80–85.

[93] M. PHPPONEN, P. RIIHILÄ, L. NISSINEN, AND V.-M. KÄHÄRI, *The role of p53 in progression of cutaneous squamous cell carcinoma*, Cancers, 13 (2021), p. 4507.
<figure><img src="image_11.png" /><figcaption>Figure 9: Illustration of the best $C_{\varepsilon}$ to approximate the derivative. We compute the analytical solution of the equation $\partial_t c = D\Delta c - \beta c$ in the 1D domain $[0, \varepsilon]$ at steady-state, with $\beta = 10^{-4}$ and $D = \varepsilon^2$, using $\varepsilon = 10^{-6}$ (a, left) and $\varepsilon = 10^{-2}$ (b, right). The equation is supplemented with Dirichlet boundary condition on the left boundary $c(t, 0) = \bar{c} = 0.1$ and zero-flux boundary condition on the right boundary. The steady state of the problem, i.e. $c^{\infty}(x_2)$ satisfying $D\Delta c^{\infty} = \beta c^{\infty}$, is plotted in a solid blue line. The dashed black line marks $\hat{c}^{\infty}$, i.e. the average of $c^{\infty}$ in the domain. The red dashed-dotted line corresponds to the derivative of $c^{\infty}$ at $x_2 = 0$. The colored dashed lines display different approximations of this derivative using the formula $\partial_{x_2}c^{\infty}(0) \approx C_{\varepsilon}(\hat{c}^{\infty} - \bar{c})/\varepsilon$, i.e. the equivalent of (A.14), for different values of $C_{\varepsilon}$: $C_{\varepsilon} = 1$ (yellow dotted line), $C_{\varepsilon} = 2$ (green dotted line), $C_{\varepsilon} = 3$ (blue dotted line), and $C_{\varepsilon} = 4$ (cyan dotted line). Analogous results were obtained varying the grid size, the parameter $\bar{c}$, and $\varepsilon$ maintaining the relation $D = \varepsilon^2$. Similar results were obtained varying $\beta \in [0, 10^{-3}]$, although larger values deteriorate the quality of the derivative approximation with $C_{\varepsilon} = 3$. This deterioration is also observed varying $\varepsilon$ or $D$ independently of each other, hence breaking the proportionality relation $D \propto \varepsilon^2$.</figcaption></figure>

**Dimension reduction.** Substituting the approximation (A.14) in (A.13), we obtain the equations for the evolution of the average quantities

$$ \partial_t \hat{c}_i + \frac{D_i^{\text{BM}}}{\varepsilon^2} C_{\varepsilon}^i (\hat{c}_i - \bar{c}_i(x_1, 0, x_3)) = \hat{R}_i - \beta_i \hat{c}_i, \quad (A.15) $$

with $i = \{t, p\}$. Analogously, approximation (A.14) can be introduced in (A.5) and (A.6).

As this quantity is constant along the $x_2$ axis and its computation only involves the knowledge of the concentrations on $\Gamma$, we can reduce the dimension of $\Omega_{\text{BM}}$ and use only $\Gamma$ as the "BM" domain.

**Physical considerations and transmission coefficients.** While in a rigorous dimension reduction procedure one needs to consider the limit $\varepsilon \to 0$, we keep the $\varepsilon$-dependency in (A.15) in order to incorporate information on the thickness of the BM into the rates at which different enzymes cross this membrane.

We consider a BM width $\varepsilon = 2 \times 10^{-6}$ dm, the average value measured by new experimental techniques as reported in [57]. The diffusion coefficients of TIMP-2 and proMMP-2 within the ECM of the CT, i.e. $\bar{D}_i$, estimated in [12] are of about $1.29 \times 10^{-6}$ dm$^2$/s. In the absence of empirical measurements or estimation of the diffusivity of these molecules within the BM, we consider that the much higher density of ECM in the BM compared to that in the CT drastically reduces molecular diffusion rates of these molecules in the BM, and we assume that $D_i^{\text{BM}} \approx \mathcal{O}(\varepsilon^2/s)$.

We note that while the value of $C_{\varepsilon}^i$ may vary as $\varepsilon$ or $D_i^{\text{BM}}$ vary independently, varying $\varepsilon$ while maintaining the proportionality relation $D_i^{\text{BM}} \propto \varepsilon^2$ does not change the value of $C_{\varepsilon}^i$ providing the best approximation using (A.14) (cf. Figure 9, for which up to 16 orders of magnitude were tested for $\varepsilon$).

We introduce the following parameters: the diffusive permeability of the molecule denoted by the index $i$ across the BM (crossing the BM in the $x_2$ direction) $\kappa_i$, and the corresponding transmission rate of the molecule to the BM $\hat{\kappa}_i$, defined as

$$ \kappa_i = \frac{C_{\varepsilon}^i D_i^{\text{BM}}}{\varepsilon} \quad \text{and} \quad \hat{\kappa}_i = \frac{\kappa_i}{\varepsilon} = \frac{C_{\varepsilon}^i D_i^{\text{BM}}}{\varepsilon^2}, \quad i = \{t, p\}. \quad (A.16) $$
[94] B. RIBBA, O. SAUT, T. COLIN, D. BRESCH, E. GRENIER, AND J. P. BOISSEL, *A multiscale mathematical model of avascular tumor growth to investigate the therapeutic benefit of anti-invasive agents*, J. Theoret. Biol., 243 (2006), pp. 532-541.

[95] R. C. ROCKNE, A. HAWKINS-DAARUD, K. R. SWANSON, J. P. SLUKA, J. A. GLAZIER, P. MACKLIN, D. A. HORMUTH, A. M. JARRETT, E. A. LIMA, J. T. ODEN, ET AL., *The 2019 mathematical oncology roadmap*, Physical biology, 16 (2019), p. 041005.

[96] E. ROIG-ROSELLO AND P. ROUSSELLE, *The human epidermal basement membrane: a shaped and cell instructive platform that aging slowly alters*, Biomolecules, 10 (2020), p. 1607.

[97] M. ROSE, *Evolutionary Biology of Aging*, Oxford University Press, 1994.

[98] T. SAITOU, K. ITANO, D. HOSHINO, N. KOSHIKAWA, M. SEIKI, K. ICHIKAWA, AND T. SUZUKI, *Control and inhibition analysis of complex formation processes*, Theoretical Biology and Medical Modelling, 9 (2012), pp. 1-18.

[99] N. SFAKIANAKIS AND M. A. J. CHAPLAIN, *Mathematical modelling of cancer invasion: A review*, in *Methods of Mathematical Oncology*, T. Suzuki, C. Poignard, M. Chaplain, and V. Quaranta, eds., vol. 370, Springer Singapore, Singapore, 2021, p. 153-172.

[100] L. F. SHAMPINE AND M. W. REICHELT, *The matlab ode suite*, SIAM Journal on Scientific Computing, 18 (1997), pp. 1-22.

[101] K.-I. SHIMOKAWA, M. KATAYAMA, Y. MATSUDA, H. TAKAHASHI, I. HARA, H. SATO, AND S. KANEKO, *Matrix metalloproteinase (MMP)-2 and MMP-9 activities in human seminal plasma*, Molecular Human Reproduction, 8 (2002), pp. 32-36.

[102] W. G. STETLER-STEVSON, *The tumor microenvironment: regulation by mmp-independent effects of tissue inhibitor of metalloproteinases-2*, Cancer and Metastasis Reviews, 27 (2008), pp. 57-66.

[103] S. SUVEGES, R. EFTIMIE, AND D. TRUCU, *Directionality of macrophages movement in tumour invasion: A multiscale moving-boundary approach*, Bulletin of Mathematical Biology, 82 (2020), pp. 1-50.

[104] M. TAURO AND C. C. LYNCH, *Cutting to the chase: how matrix metalloproteinase-2 activity controls breast-cancer-to-bone metastasis*, Cancers, 10 (2018), p. 185.

[105] M. TOOTH, M. M. BERNARDO, D. C. GERVASI, P. D. SOLOWAY, Z. WANG, H. F. BIGG, C. M. OVERALL, Y. A. DECLERCK, H. TSCHESCHE, M. L. CHER, ET AL., *Tissue inhibitor of metalloproteinase (timp)-2 acts synergistically with synthetic matrix metalloproteinase (mmp) inhibitors but not with timp-4 to enhance the (membrane type 1)-mmp-dependent activation of pro-mmp-2*, Journal of Biological Chemistry, 275 (2000), pp. 41415-41423.

[106] T. H. VU AND Z. WEBB, *Matrix metalloproteinases: effectors of development and normal physiology*, Genes & development, 14 (2000), pp. 2123-2133.

[107] T. WANG, F. NOTTA, R. NAVAB, J. JOSEPH, E. IBRAHIMOV, J. XU, C.-Q. ZHU, A. BORGIDA, S. GALLINGER, AND M.-S. TSAO, *Senescent Carcinoma-Associated Fibroblasts Upregulate IL8 to Enhance Prometastatic Phenotypes*, Molecular Cancer Research, 15 (2017), pp. 3-14.

[108] H. E. V. WART AND H. BIRKEDAL-HANSEN, *The cysteine switch: a principle of regulation of metalloproteinase activity with potential applicability to the entire matrix metalloproteinase gene family*, Proceedings of the National Academy of Sciences, 87 (1990), pp. 5578-5582.

[109] A. WATANABE, D. HOSINO, N. KOSHIKAWA, M. SEIKI, T. SUZUKI, AND K. ICHIKAWA, *Critical role of transient activity of mt1-mmp for ecm degradation in invadopodia*, PLOS Computational Biology, 9 (2013), pp. 1-13.

[110] M. C. WHITE, D. M. HOLMAN, J. E. BOEHM, L. A. PEIPINS, M. GROSSMAN, AND S. J. HENLEY, *Age and cancer risk: a potentially modifiable relationship*, American journal of preventive medicine, 46 (2014), pp. S7-S15.
**Remark 3.** We emphasise that our method is formal and we do not send $\varepsilon \to 0$, but only retain the dynamics of the average concentration along the direction selected for dimensionality reduction, assuming that $\varepsilon$, i.e. the BM width, is very small.

**Remark 4.** While we adopt a formal method for the dimension reduction of the BM, rigorous methods yielding analogous results have been developed – we refer the interested reader to the method of two-scale convergence for dimension reduction problems, e.g. see [46]. However, we also mention that in other models modelling the effect of thin membranes, transmission is modelled as Kedem-Katchalsky conditions, in which the flux is computed from the permeability and the concentration difference between the sides of the membrane. Here, using this formal dimension reduction, we retain an information about the average of these quantities in the BM in a simple manner.

**Summary of the full reduced model.** Altogether, using the bar notation to denote concentrations in the CT $\Omega_{CT}$ and the hat notation to denote the concentrations in the BM $\Omega_{BM} \equiv \Gamma$ (which in now an open subset of $\mathbb{R}^{d-1}$), we have the model

$$
\begin{align*}
\partial_t \bar{c}_t &= \bar{D}_t \Delta \bar{c}_t + \bar{S}_t(x) - \beta_t \bar{c}_t, && \text{in } \Omega_{CT} \times (0, T), \\
\partial_t \bar{c}_p &= \bar{D}_p \Delta \bar{c}_p + \bar{S}_p(x) - \beta_p \bar{c}_p, && \text{in } \Omega_{CT} \times (0, T), \\
-\bar{D}_t \nabla \bar{c}_t \cdot \nu &= \kappa_t (\bar{c}_t - \hat{c}_t), && \text{on } \Gamma \times (0, T), \\
-\bar{D}_p \nabla \bar{c}_p \cdot \nu &= \kappa_p (\bar{c}_p - \hat{c}_p), && \text{on } \Gamma \times (0, T), \\
\partial_t \hat{c}_t &= \hat{\kappa}_t (\bar{c}_t - \hat{c}_t) + \hat{R}_t - \beta_t \hat{c}_t, && \text{in } \Gamma \times (0, T), \\
\partial_t \hat{c}_p &= \hat{\kappa}_p (\bar{c}_p - \hat{c}_p) + \hat{R}_p - \beta_p \hat{c}_p, && \text{in } \Gamma \times (0, T), \\
\partial_t \hat{M} &= -\gamma \hat{c}_a \frac{\hat{M}}{K_m + \hat{M}} + r_M \left(1 - \frac{\hat{M}}{M_{\text{max}}}\right)_+, && \text{in } \Gamma \times (0, T), \\
\partial_t \hat{c}_m &= \alpha_m \hat{\rho}^0 - k_0 \hat{c}_m^2 + k_{-0} \hat{c}_d - \beta_m \hat{c}_m, && \text{in } \Gamma \times (0, T), \\
\partial_t \hat{c}_d &= k_0 \hat{c}_m^2 - k_{-0} \hat{c}_d - k_1 \hat{c}_d \hat{c}_t + k_{-1} \hat{c}_1 - \beta_d \hat{c}_d, && \text{in } \Gamma \times (0, T), \\
\partial_t \hat{c}_a &= k_3 \hat{c}_2 - k_2 \hat{c}_t \hat{c}_a + k_{-2} \hat{c}_{ta} - \beta_a \hat{c}_a, && \text{in } \Gamma \times (0, T), \\
\partial_t \hat{c}_1 &= k_1 \hat{c}_d \hat{c}_t - \hat{c}_1 [k_1 \hat{c}_t + k_2 \hat{c}_p] + k_{-1} (\hat{c}_3 - \hat{c}_1) + (k_{-2} + k_3) \hat{c}_2 - \beta_1 \hat{c}_1, && \text{in } \Gamma \times (0, T), \\
\partial_t \hat{c}_2 &= k_2 \hat{c}_1 \hat{c}_p - (k_{-2} + k_3) \hat{c}_2 - \beta_2 \hat{c}_2, && \text{in } \Gamma \times (0, T), \\
\partial_t \hat{c}_3 &= k_1 \hat{c}_1 \hat{c}_t - k_{-1} \hat{c}_3 - \beta_3 \hat{c}_3, && \text{in } \Gamma \times (0, T), \\
\partial_t \hat{c}_{tp} &= k_2 \hat{c}_t \hat{c}_p - k_{-2} \hat{c}_{tp} - \beta_{tp} \hat{c}_{tp}, && \text{in } \Gamma \times (0, T), \\
\partial_t \hat{c}_{ta} &= k_2 \hat{c}_t \hat{c}_a - k_{-2} \hat{c}_{ta} - \beta_{ta} \hat{c}_{ta}, && \text{in } \Gamma \times (0, T).
\end{align*}
$$

We note that in Section 2 the hat is dropped for convenience, while $\hat{\kappa}_i(\bar{c}_i - \hat{c}_i) = S_i(x, t)$ and the remaining biochemical reaction terms introduced in equations (2.4) and (2.5) are denoted by $\hat{R}_i$.

## B Computation details of biologically relevant parameter values

Most baseline parameter values reported in Table 1 are found directly in the literature, short of a change of unit to fit our reference system when required. These come from experimental measurements, were computed using a computational model (different from ours), or were estimated to fit a relevant range. As we use multiple sources to find a priori values for the parameters, we refine our parameter set using optimization techniques to ensure that the outputs of our model are biologically consistent. Below we give details of the references used and computations performed to estimate parameter values which were not found in the literature (marked with a * in Table 1).

**BM carrying capacity.** The BM carrying capacity $M_{\text{max}}$ represents the maximum ECM density in the BM. We consider the BM matrigel density used in *in vitro* invasion assays reported in [1], whose authors estimated a protein concentration of about 10 mg/ml. Considering a molecular weight in the typical weight range of collagen IV, the main component of the ECM found in the BM, i.e. about 160kDa [42], we compute

$$M_{\text{max}} = \frac{10 \text{ mg/ml}}{160 \times 10^3 \text{ Da}} \approx 6.2510^{-5} \text{ mol/dm}^3 = 62.5 \times 10^3 \text{ nM}.$$
[111] E. C. WOENNE, W. LEDERLE, S. ZWICK, M. PALMOWSKI, H. KRELL, W. SEMMLER, M. M. MUELLER, AND F. KIESSLING, *Mmp inhibition blocks fibroblast-dependent skin cancer invasion, reduces vascularization and alters vegf-a and pdgf-bb expression*, Anticancer Research, 30 (2010), pp. 703-711.

[112] C. WU, G. LORENZO, D. A. HORMUTH, E. A. B. F. LIMA, K. P. SLAVKOVA, J. C. DI-CARLO, J. VIROSTKO, C. M. PHILLIPS, D. PATT, C. CHUNG, AND T. E. YANKEELOV, *Integrating mechanism-based modeling with biomedical imaging to build practical digital twins for clinical oncol-ogy*, Biophysics Reviews, 3 (2022), p. 021304.

# A Mathematical notation and derivation of the model

## A.1 Original system

**Notation.** We consider a general function $u(x, t)$ depending on space $x \in \mathbb{R}^d$ ($d = 1, 2, 3$, we choose $d = 3$ in this section) and time $t \in [0, T]$ ($T > 0$ a positive finite time horizon). We denote by $\partial_t u$ the partial derivative of $u$ with respect to time and $\nabla u = (\partial_{x_1} u, \partial_{x_2} u, \partial_{x_3} u)$ the gradient of $u$ in space. We denote by $\Delta u = \text{div}(\nabla u) = \partial_{x_1}^2 u + \partial_{x_2}^2 u + \partial_{x_3}^2 u$ the laplacian of $u$, with $\text{div}(\vec{v}) = \partial_{x_1} v_1 + \partial_{x_2} v_2 + \partial_{x_3} v_3$ the divergence of the vector field $\vec{v} = (v_1, v_2, v_3)$.

**Full system.** In this section, we introduce the mathematical model comprising a system of PDEs describing the dynamics of the SFs' secretome in the CT, with boundary conditions modelling the transmission of these concentrations to the BM, and a system of PDEs describing the activation and degradation dynamics in the BM. Details about the modelling assumptions and the notations for the different unknowns can be found in Section 2.

We start with two domains, i.e. smooth bounded open subsets, $\Omega_{\text{CT}}, \Omega_{\text{BM}} \subset \mathbb{R}^d$, separated by a "membrane" denoted by $\Gamma$. We highlight that $\Gamma$ simply models the boundary separating the two domains and not a real physical membrane with a positive thickness. We have the following system of equations

$$ \partial_t \bar{c}_t = \bar{D}_t \Delta \bar{c}_t + \bar{S}_t(x) - \beta_t \bar{c}_t, \quad \text{in } \Omega_{\text{CT}} \times (0, T), \quad (A.1) $$

$$ \partial_t \bar{c}_p = \bar{D}_p \Delta \bar{c}_p + \bar{S}_p(x) - \beta_p c_p, \quad \text{in } \Omega_{\text{CT}} \times (0, T), \quad (A.2) $$

$$ \bar{c}_t = c_t, \quad \text{on } \Gamma \times (0, T), \quad (A.3) $$

$$ \bar{c}_p = c_p, \quad \text{on } \Gamma \times (0, T), \quad (A.4) $$

$$ \bar{D}_t \nabla \bar{c}_t \cdot \nu = D_t^{\text{BM}} \nabla c_t \cdot \nu, \quad \text{on } \Gamma \times (0, T), \quad (A.5) $$

$$ \bar{D}_p \nabla \bar{c}_p \cdot \nu = D_p^{\text{BM}} \nabla c_p \cdot \nu, \quad \text{on } \Gamma \times (0, T), \quad (A.6) $$

$$ \partial_t c_t = D_t^{\text{BM}} \Delta c_t + R_t(x, t) - \beta_t c_t, \quad \text{in } \Omega_{\text{BM}} \times (0, T), \quad (A.7) $$

$$ \partial_t c_p = D_p^{\text{BM}} \Delta c_p + R_p(x, t) - \beta_p c_p, \quad \text{in } \Omega_{\text{BM}} \times (0, T), \quad (A.8) $$

with $\nu$ the unit normal vector to $\Gamma$ that points from the BM to the CT. The rest of the boundary conditions on $\partial \Omega_{\text{CT}} \setminus \Gamma \times (0, T)$ and $\partial \Omega_{\text{BM}} \setminus \Gamma \times (0, T)$ are set to zero-flux boundary conditions.

Equations (A.3)-(A.6) are transmission conditions modeling the diffusion of TIMP-2 and proMMP-2 between the CT and the BM. In particular, Equations (A.3)-(A.4) model continuity of the concentration at the border between the two domains, and equations (A.5)-(A.6) model the continuity of the flux across this interface. Moreover, while the terms $\bar{S}_t(x)$ and $\bar{S}_p(x)$ in equations (A.1)-(A.2) are as defined in (2.15), the terms $R_t(x, t)$ and $R_p(x, t)$ in equations (A.7)-(A.8) include the biochemical reactions terms. Namely, we have

$$ R_t(x, 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}), \quad (A.9) $$

$$ R_p(x, t) = -k_2 c_p (c_t + c_1) + k_{-2} (c_{tp} + c_2). \quad (A.10) $$

Equations (A.1)-(A.8) are coupled, via the reaction terms $R_i(t, x)$, $i = \{t, p\}$, with the rest of the system modelling proMMP-2 activation, BM digestion and TIMP-mediated MMP inhibition: assuming
Physiological levels of TIMP-2, proMMP-2, active MMP-2 and TIMP-2/MMP-2 complex.

We use as reference, the serum levels reported in [52] for TIMP-2, proMMP-2, active MMP-2 and TIMP-2/MMP-2 complex, which are in good agreement with the values reported by recent studies. Knowing the molecular masses of TIMP-2 (21 kDa, see e.g. [62] and references therein), proMMP-2 (72 kDa [61]), active MMP-2 (67 kDa [61]), and the complex TIMP-2/MMP-2 (88kDa, sum of the two free proteins), we compute the physiological molar concentrations

$$c_t^{ph} = \frac{201.8 \text{ ng/cm}^3}{21 \times 10^3 \text{ Da}} \approx 9.6 \text{ nmol/dm}^3 = 9.6 \text{nM}$$

$$c_p^{ph} = \frac{1259.9 \text{ ng/cm}^3}{72 \times 10^3 \text{ Da}} \approx 18 \text{nM},$$

$$c_a^{ph} = \frac{29.4 \text{ ng/cm}^3}{67 \times 10^3 \text{ Da}} \approx 0.44 \text{nM},$$

$$c_{ta}^{ph} = \frac{1456.4 \text{ ng/cm}^3}{88 \times 10^3 \text{ Da}} \approx 17 \text{nM}.$$

For the level of TIMP-2/proMMP-2 $c_{tp}^*$, we use the value reported in [2] and obtain

$$c_{tp}^{ph} = \frac{625 \times 10^3 \text{ ng/L}}{93 \times 10^3 \text{ Da}} \approx 6.72 \text{nM}.$$

For the reduced model, detailed in Appendix A, we assume equilibrium level of MT1-MMP dimers are given. With analogous calculations to those reported above, we take this to be $\bar{c}_d = 8.55/(63 \times 10^3) \approx 0.14$ nM, from normal serum concentration of MT1-MMP given in [40] and molecular mass in [11].

**Remark 5.** In the review article [29], the authors survey the different values of the level of MMP-2, proMMP-2, TIMP-2 and the complex TIMP-2/MMP-2. We observe that, from one study to another, discrepancies in serum or plasma levels of the proteins are found, which could be due to different techniques used to measure these quantities but also due to differences in the cohort. The chosen reference publication for the baseline parameter values used here ensures higher consistency in origin and measurement technique among these values. Moreover, as already mentioned, these values are in good agreement with those reported by recent studies.

**Remark 6.** 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 measurements in tissue samples, we use these values for the moment.

**Physiological production rate of proMMP-2 and TIMP-2 in the CT.** We assume the physiological concentrations of TIMP-2 and pro-MMP2 for a healthy patient are equivalent to the equilibria $c_i^*$ ($i = \{t, p\}$) of their dynamics in the CT (2.13)–(2.12) in the absence of SFs. This yields the relation

$$c_i^{ph} = \bar{c}_i^* = \frac{s_i^{ph}}{\beta_i}, \quad i = \{t, p\}.$$

Using the values for decay rates and physiological levels reported in Table 1, we find

$$s_p^{ph} = 1.09 \times 10^{-2} \text{nM/s},$$

$$s_t^{ph} = 4.28 \times 10^{-4} \text{nM/s}.$$

**Molecular decay rates.** Decay rates $\beta_i$ in System (2.1)–(2.12) represent molecular natural decay, but may also capture the effect of molecule degradation and internalization by other actors, e.g. internalization of membrane bound MT1-MMP [34], or internalization of TIMP-2/proMMP-2 [19]. We assume the decay rates of MT1-MMP dimers and complexes match the one of MT1-MMP monomers, for which we use the rate of MT1-MMP internalization estimated in [30].

We assume physiological levels of TIMP-2/proMMP-2 and TIMP-2/MMP-2 complexes correspond to their concentrations at the equilibrium of system (2.1)–(2.13) in the absence of SFs. We estimate their decay rates from the equilibrium of equations (2.10)–(2.11), using the reaction rates indicated in Table 1, obtaining

$$\beta_{tp} = \frac{k_2 c_t^* c_p^*}{c_{tp}^*} - k_{-2} = \frac{k_2 c_t^{ph} c_p^{ph}}{c_{tp}^{ph}} - k_{-2} = 3.60 \text{s}^{-1},$$

$$\beta_{ta} = \frac{k_2 c_t^* c_a^*}{c_{ta}^*} - k_{-2} = \frac{k_2 c_t^{ph} c_a^{ph}}{c_{ta}^{ph}} - k_{-2} = 0.035 \text{s}^{-1}.$$
the diffusion of the rest of the components is negligible, we have

$$ \partial_t M = -\gamma c_a \frac{M}{K_m + M} + r_M \left(1 - \frac{M}{M_{\max}}\right)_+ , \quad \text{in } \Omega_{BM} \times (0, T), $$

$$ \partial_t c_m = \alpha_m \rho^0 - k_0 c_m^2 + k_{-0} c_d - \beta_m c_m, \quad \text{in } \Omega_{BM} \times (0, T), $$

$$ \partial_t c_d = k_0 c_m^2 - k_{-0} c_d - k_1 c_d c_t + k_{-1} c_1 - \beta_d c_d, \quad \text{in } \Omega_{BM} \times (0, T), $$

$$ \partial_t c_a = k_3 c_2 - k_2 c_t c_a + k_{-2} c_{ta} - \beta_a c_a, \quad \text{in } \Omega_{BM} \times (0, T), $$

$$ \partial_t c_1 = 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 \text{in } \Omega_{BM} \times (0, T), $$

$$ \partial_t c_2 = k_2 c_1 c_p - (k_{-2} + k_3)c_2 - \beta_2 c_2, \quad \text{in } \Omega_{BM} \times (0, T), $$

$$ \partial_t c_3 = k_1 c_1 c_t - k_{-1} c_3 - \beta_3 c_2, \quad \text{in } \Omega_{BM} \times (0, T), $$

$$ \partial_t c_{tp} = k_2 c_t c_p - k_{-2} c_{tp} - \beta_{tp} c_{tp}, \quad \text{in } \Omega_{BM} \times (0, T), $$

$$ \partial_t c_{ta} = k_2 c_t c_a - k_{-2} c_{ta} - \beta_{ta} c_{ta}, \quad \text{in } \Omega_{BM} \times (0, T). $$

## A.2 Dimension reduction of the BM

In order to improve computational tractability of our model, as the BM is a very thin structure compared to the rest of the tissue [7], we reduce the dimension of $\Omega_{BM}$, which then becomes a subset of $\mathbb{R}^{d-1}$ (whereas $\Omega_{CT} \subset \mathbb{R}^d$). We do this following a simple and formal dimension reduction procedure, inspired by ideas from [74] to derive 1D-3D coupled problems.

**Smallness assumptions and decomposition.** We introduce a small parameter $0 < \varepsilon \ll 1$ describing the thickness of the BM, i.e. $\bar{\Omega}_{BM} = \Gamma \times [0, \varepsilon]$. Assuming $\Gamma$ is a flat surface, we consider that $\forall x \in \bar{\Omega}_{BM}$, $x$ is defined by $x := (x_1, x_2, x_3)$ with $(x_1, x_3) \in \Gamma$ and $x_2 = \varepsilon \bar{x}_2$, $(\bar{x}_2 \in (0, 1))$, coherently with the coordinate system adopted in Section 2 and illustrated in Figure 2. We decompose the concentrations $c_i$ as

$$ c_i(x_1, x_2, x_3, t) = \hat{c}_i(x_1, x_3, t) + \tilde{c}_i(x_1, x_2, x_3, t) \quad \text{for } i = \{t, p, d, tp, ta, 1, 2, 3\}, \quad (A.11) $$

where $\hat{c}_i(x_1, x_3, t)$ corresponds to the average of the concentration in the $x_2$ direction, i.e.

$$ \hat{c}_i = \frac{1}{\varepsilon} \int_0^\varepsilon c_i \, dx_2 = \int_0^1 c_i(\varepsilon \bar{x}_2) \, d\bar{x}_2, \quad (A.12) $$

and $\tilde{c}_i(x_1, x_2, x_3, t)$ is a fluctuation around this average. Integrating both sides of the decomposition (A.11) along the width of the BM (i.e. the $x_2$ direction) and using the definition (A.12) of the average $\hat{c}_i$, we obtain that the fluctuation must satisfy

$$ \int_0^\varepsilon \tilde{c}_i \, dx_2 = 0. $$

We assume that the fluctuations are small and that the concentration fields evolve smoothly within the BM, i.e. spatial derivatives of the fluctuations in the BM are small.

Under assumption (A.11) and definitions (A.9)-(A.10), we have that both sources $R_i$ ($i = \{t, p\}$) in (A.7) and (A.8) also admit a decomposition in the form (A.11).

**Simplified dynamics in the BM.** Since we assumed that the fluctuation of the concentrations within the BM – i.e. $\tilde{c}_i$ in (A.11) – are small, we simplify the dynamics in the BM by focusing on the evolution of the averages of the concentration fields in the BM along the $x_2$ direction (i.e. along the width of the BM) – i.e. $\hat{c}_i$ defined in (A.12). This will allow us to drop the dependency on $x_2$ of the dynamics, providing an approximation that is well suited for a model in which the dimensionality of the BM is reduced to $d-1$. Given that the thickness of the BM is very small, we expect that the diffusion of the enzymes in the $x_2$ direction will have a strong impact on the average concentration $\hat{c}_i$, while we assume diffusion in the $(x_1, x_3)$ plane (i.e. tangential to $\Gamma$) is negligible.

Integrating Equations (A.7)–(A.8) with respect to $x_2$, dividing by $\varepsilon$, and neglecting the diffusion in the tangential plane $(x_1, x_3)$, we obtain for any $(x_1, x_3) \in \Gamma$ and $t \in (0, T)$

$$ \frac{1}{\varepsilon} \int_0^\varepsilon \partial_t c_i \, dx_2 - \frac{D_i^{BM}}{\varepsilon} \left[ \nabla c_i \cdot \nu \bigg|_{x_2=\varepsilon} - \nabla c_i \cdot \nu \bigg|_{x_2=0} \right] = \frac{1}{\varepsilon} \int_0^\varepsilon R_i - \beta_i c_i \, dx_2, $$
The other decay rates, assumed to model natural decay, are computed from the half-life $t_{i,1/2}$ of the molecule via the formula $\beta_i = \ln(2)/t_{i,1/2}$. We know from [6] that the half-life of TIMP-2 in blood is $t_{t,1/2} = 4.33$ hours, and compute

$$ \beta_t = \frac{\ln(2)}{4.33 \times 3600} \approx 4.45 \times 10^{-5} \text{s}^{-1}. $$

In [47], the authors reported a half-life in the plasma of $t_{p,1/2} = 19$ min for human recombinant proMMP-2 injected in mice deficient in MMP-2. As reported in [13], active MMP-9 has a shorter half-life compared to the pro and inactive enzyme. Hence, we assume that $t_{a,1/2} = 0.5t_{p,1/2} = 9.5$ min. This gives the coefficients

$$ \beta_p = \frac{\ln(2)}{19 \times 60} \approx 6.08 \times 10^{-4} \text{s}^{-1}, $$

and

$$ \beta_a = \frac{\ln(2)}{9.5 \times 60} \approx 1.22 \times 10^{-3} \text{s}^{-1}. $$

**Diffusive permeabilities and transmission rates at the BM.** As detailed in Appendix A, we consider a BM of width $\varepsilon = 2 \times 10^{-6}$ dm, estimated for the human epidermal BM via new measuring techniques, as reported in [57]. We consider diffusion coefficients of TIMP-2 and proMMP-2 within the ECM of the CT estimated in [12], i.e. $\bar{D}_i = 1.29 \times 10^{-6}$ dm$^2$/s. Given the much higher density of ECM in the BM compared to that in the CT, we take the diffusivity of TIMP-2 and proMMP-2 in the BM to be $D_i^{BM} = 1.29 \times 10^{-12}$ dm$^2$/s. This is consistent with the assumption that $D_i^{BM} \approx O(\varepsilon^2/s)$ (while $\bar{D}_i$ is of the order of $\varepsilon$). Finally, to estimate the diffusive permeability of TIMP-2 and proMMP-2 in the BM, $\kappa_i$, and their transmission rate across the BM, $\hat{\kappa}_i$, using definitions (A.16), we take $C_i^e = 3$. The best value to use for this parameter, introduced in the approximation of the concentration gradient in (A.14), was investigated numerically (see Figure 9) Altogether this gives

$$ \kappa_i = \frac{3 \times 1.29 \times 10^{-12} \text{ dm}^2}{2 \times 10^{-6} \text{ dm}} \approx 1.935 \times 10^{-6} \text{ dm/s}, $$

and

$$ \hat{\kappa}_i = \frac{3 \times 1.29 \times 10^{-12} \text{ dm}^2}{4 \times 10^{-12} \text{ dm}^2} \approx 0.9675 \text{ s}^{-1}. $$

## C Reduced system and analytical results

We reduce System (2.1)-(2.11), comprising 11 coupled nonlinear equations, to an analytically tractable model that retains the key biological processes under consideration. We then study the stability of this reduced system and derive necessary conditions for BM rupture, i.e. tumour invasion.

**Simplifying assumptions.** We make the following simplifying assumptions:

SA1. the inhibition of proMMP-2 and active MMP-2 by TIMP-2 can not reverse (i.e. we take $k_{-2}c_{tp} = k_{-2}c_{ta} \approx 0$ in Equations (2.4)-(2.6) and (2.10)-(2.11));

SA2. the activation and inhibition of proMMP-2 are very fast compared to the rest of the dynamics (hence, we use the quasi steady-state assumption in Equations (2.7)-(2.11));

SA3. stable MT1-MMP expression by cancer cells is already in place and MT1-MMP dimerisation occurs rapidly, so that MT1-MMP dimers quickly stabilize to some equilibrium value $\bar{c}_d$, which is kept constant henceforth;

SA4. the natural decay of proMMP-2, TIMP-2 and active MMP-2 inside the BM is negligible compared to the faster activation/degradation dynamics they are involved in (i.e. we take $\beta_p = \beta_t = \beta_a \approx 0$ in Equations (2.4)-(2.6)).
where we used the divergence theorem, with $\nu$ being the outward normal to $\partial\Omega_{BM}$. Since we assumed zero-flux boundary conditions on $\partial\Omega_{BM} \setminus \Gamma \times (0, T)$ and using the definition (A.12), we have

$$ \partial_t \hat{c}_i + \frac{D_i^{BM}}{\varepsilon} \partial_{x_2} c_i(x_1, 0, x_3, t) \cdot \nu = \hat{R}_i - \beta_i \hat{c}_i, \quad (A.13) $$

for $i = \{t, p\}$. We note that the directional derivative $\partial_{x_2} c_i(x_1, 0, x_3, t) \cdot \nu$ also appears on the right-hand side of (A.5)-(A.6).

**Approximation of concentration gradients at the interface.** We now construct an approximation of $\partial_{x_2} c_i(x_1, 0, x_3, t)$ for $i = \{t, p\}$. Since $\varepsilon$ is small, assuming $c_i$ is sufficiently smooth to allow a Taylor series approximation, we consider the following expansion around $x_2 = 0$

$$ c_i(x_1, \varepsilon, x_3, t) = c_i(x_1, 0, x_3, t) + \varepsilon \partial_{x_2} c_i(x_1, 0, x_3, t) + \mathcal{O}(\varepsilon^2). $$

Using this, we know that the following holds

$$ \partial_{x_2} c_i(x_1, 0, x_3, t) = \frac{c_i(x_1, \varepsilon, x_3, t) - c_i(x_1, 0, x_3, t)}{\varepsilon} + \mathcal{O}(\varepsilon). $$

Using the continuity of concentrations (A.3)-(A.4) at the interface $\Gamma$, and the decomposition (A.11) at $x_2 = \varepsilon$, we rewrite this as

$$ \partial_{x_2} c_i(x_1, 0, x_3, t) = \frac{\hat{c}_i(x_1, x_3, t) - \bar{c}_i(x_1, 0, x_3, t)}{\varepsilon} + \frac{\tilde{c}_i(x_1, \varepsilon, x_3, t)}{\varepsilon} + \mathcal{O}(\varepsilon). $$

We approximate $\tilde{c}_i(x_1, \varepsilon, x_3, t)$ as

$$ \tilde{c}_i(x_1, \varepsilon, x_3, t) = (C_{\varepsilon}^i - 1)(\hat{c}_i(x_1, x_3, t) - \bar{c}_i(x_1, 0, x_3, t)), $$

with constant $C_{\varepsilon}^i > 0$ for $\varepsilon > 0$. Thus, for $\varepsilon$ small, we have

$$ \partial_{x_2} c_i(x_1, 0, x_3) \approx C_{\varepsilon}^i \frac{\hat{c}_i(x_1, x_3) - \bar{c}_i(x_1, 0, x_3)}{\varepsilon} \quad \text{for } i = \{t, p\}. \quad (A.14) $$

This approximation is justified by the fact that the $c_i$ ($i = \{t, p\}$) satisfy a reaction-diffusion equation in a thin domain, with Dirichlet boundary conditions on $\Gamma$ and zero-flux boundary conditions at the other end of the domain. At equilibrium, this leads to a concentration profile along the $x_2$ direction analogous to the one shown in Figure 9, where we simulated this scenario in 1D (with $\Gamma$ the left boundary) and tested the approximation (A.14) for different values of $C_{\varepsilon}^i$. We clearly have that $C_{\varepsilon}^i > 0$ for $\varepsilon > 0$ (analogous profiles were obtained for decreasing values of the domain size).

**Dimension reduction.** Substituting the approximation (A.14) in (A.13), we obtain the equations for the evolution of the average quantities

$$ \partial_t \hat{c}_i + \frac{D_i^{BM}}{\varepsilon^2} C_{\varepsilon}^i (\hat{c}_i - \bar{c}_i(x_1, 0, x_3)) = \hat{R}_i - \beta_i \hat{c}_i, \quad (A.15) $$

with $i = \{t, p\}$. Analogously, approximation (A.14) can be introduced in (A.5) and (A.6).

As this quantity is constant along the $x_2$ axis and its computation only involves the knowledge of the concentrations on $\Gamma$, we can reduce the dimension of $\Omega_{BM}$ and use only $\Gamma$ as the "BM" domain.

**Physical considerations and transmission coefficients.** While in a rigorous dimension reduction procedure one needs to consider the limit $\varepsilon \to 0$, we keep the $\varepsilon$-dependency in (A.15) in order to incorporate information on the thickness of the BM into the rates at which different enzymes cross this membrane.

We consider a BM width $\varepsilon = 2 \times 10^{-6}$ dm, the average value measured by new experimental techniques as reported in [96]. The diffusion coefficients of TIMP-2 and proMMP-2 within the ECM of the CT, i.e. $\bar{D}_i$, estimated in [22] are of about $1.29 \times 10^{-6}$ dm$^2$/s. In the absence of empirical measurements or estimation of the diffusivity of these molecules within the BM, we consider that the much higher density of ECM in the BM compared to that in the CT drastically reduces molecular diffusion rates of these molecules in the BM, and we assume that $D_i^{BM} \approx \mathcal{O}(\varepsilon^2/s)$.
**Reduced model.** Under these assumptions, System (2.1)-(2.11) reduces to the System (2.25)-(2.27), i.e.

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

$$ \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, $$

$$ \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, $$

where we have introduced the notation $\bar{c}_i \equiv \bar{c}_i(x, t)|_\Gamma$ ($i = \{t, p\}$) for the TIMP-2 and proMMP-2 concentration at the BM coming from the CT, and equation (2.1) for the digestion of the BM, complemented by the quantities obtained from the quasi-state assumptions SA2., i.e.

$$ c_1 = \frac{k_1 \bar{c}_d}{k_{-1}} c_t, \quad c_2 = \frac{k_2 k_1 \bar{c}_d}{k_{-1}(k_{-2} + k_3)} c_t c_p, \quad c_3 = \left( \frac{k_1}{k_{-1}} \right)^2 \bar{c}_d c_t^2, \quad c_{tp} = \frac{k_2 c_t c_p}{\beta_{tp}}, \quad c_{ta} = \frac{k_2 c_t c_a}{\beta_{ta}}. $$

In the absence of spatial movement in the BM, the reduced system holds pointwise in $x \in \Gamma$.

## C.1 Analytical results of the reduced system

We first analyse the steady state of the ODE system (2.25)-(2.27) and its linear stability. We then exploit these results to derive sufficient conditions for BM rupture while analysing Equation (2.24).

### C.1.1 Steady state

For simplicity, we assume the TIMP-2 and proMMP-2 dynamics in the CT have reached an equilibrium and denote their equilibrium value at the BM by

$$ \bar{c}_i = \bar{c}_i^* \equiv \bar{c}_i^*(x), \quad i = \{t, p\}, \quad x \in \Gamma. \qquad (C.1) $$

**Proposition 7.** If $\bar{c}_t^* > 0$, denoting

$$ K_R = k_2 + \frac{k_1 k_2 k_3 \bar{c}_d}{k_{-1}(k_{-2} + k_3)}, \qquad (C.2) $$

the steady state solution of the ODE system (2.25)-(2.27) under assumption (C.1) is

$$ \left\{ \begin{aligned} c_t^* &= \frac{\hat{\kappa}_t \bar{c}_t^*}{\hat{\kappa}_t + K_R c_p^*}, \\ c_p^* &= \sqrt{\left( \frac{\hat{\kappa}_t \bar{c}_t^* - \hat{\kappa}_p \bar{c}_p^*}{2\hat{\kappa}_p} + \frac{\hat{\kappa}_t}{2K_R} \right)^2 + \frac{\hat{\kappa}_t \bar{c}_t^*}{K_R}} - \left( \frac{\hat{\kappa}_t \bar{c}_t^* - \hat{\kappa}_p \bar{c}_p^*}{2\hat{\kappa}_p} + \frac{\hat{\kappa}_t}{2K_R} \right), \\ c_a^* &= \frac{k_1 k_3 \bar{c}_d}{k_{-1}(k_{-2} + k_3)} c_p^*. \end{aligned} \right. \qquad (C.3) $$

If $\bar{c}_t^* = 0$, the system (2.25)-(2.27) under assumption (C.1) admits infinitely many steady states in the form

$$ c_t^* = 0, \quad c_p^* = \bar{c}_p^*, \quad c_a^* = C_a, \qquad (C.4) $$

with $C_a \ge 0$ for biological relevance. From the non-negativity of parameters we have that $c_t^*, c_p^*, c_a^* \ge 0$.

**Proof.** The steady state $Y^* = (c_t^*, c_p^*, c_a^*)$ of system (2.25)-(2.27), under assumption (C.1), satisfies

$$ \left\{ \begin{aligned} \hat{\kappa}_t \bar{c}_t^* - k_2 c_t^* (c_p^* + c_a^*) - \hat{\kappa}_t c_t^* &= 0, \\ \hat{\kappa}_p \bar{c}_p^* - K_R c_p^* c_t^* - \hat{\kappa}_p c_p^* &= 0, \\ (K_R - k_2) c_t^* c_p^* - k_2 c_t^* c_a^* &= 0, \end{aligned} \right. \qquad (C.5) $$

with $K_R$ defined in (C.2). Rearranging (C.5)$_1$ we obtain

$$ c_t^* = \frac{\hat{\kappa}_t \bar{c}_t^*}{\hat{\kappa}_t + k_2(c_p^* + c_a^*)}. \qquad (C.6) $$
<figure><img src="image_10.png" /><figcaption>Figure 9: Illustration of the best $C_{\varepsilon}$ to approximate the derivative. We compute the analytical solution of the equation $\partial_t c = D\Delta c - \beta c$ in the 1D domain $[0, \varepsilon]$ at steady-state, with $\beta = 10^{-4}$ and $D = \varepsilon^2$, using $\varepsilon = 10^{-6}$ (a, left) and $\varepsilon = 10^{-2}$ (b, right). The equation is supplemented with Dirichlet boundary condition on the left boundary $c(t, 0) = \bar{c} = 0.1$ and zero-flux boundary condition on the right boundary. The steady state of the problem, i.e. $c^{\infty}(x_2)$ satisfying $D\Delta c^{\infty} = \beta c^{\infty}$, is plotted in a solid blue line. The dashed black line marks $\hat{c}^{\infty}$, i.e. the average of $c^{\infty}$ in the domain. The red dashed-dotted line corresponds to the derivative of $c^{\infty}$ at $x_2 = 0$. The colored dashed lines display different approximations of this derivative using the formula $\partial_{x_2}c^{\infty}(0) \approx C_{\varepsilon}(\hat{c}^{\infty} - \bar{c})/\varepsilon$, i.e. the equivalent of (A.14), for different values of $C_{\varepsilon}$: $C_{\varepsilon} = 1$ (yellow dotted line), $C_{\varepsilon} = 2$ (green dotted line), $C_{\varepsilon} = 3$ (blue dotted line), and $C_{\varepsilon} = 4$ (cyan dotted line). Analogous results were obtained varying the grid size, the parameter $\bar{c}$, and $\varepsilon$ maintaining the relation $D = \varepsilon^2$. Similar results were obtained varying $\beta \in [0, 10^{-3}]$, although larger values deteriorate the quality of the derivative approximation with $C_{\varepsilon} = 3$. This deterioration is also observed varying $\varepsilon$ or $D$ independently of each other, hence breaking the proportionality relation $D \propto \varepsilon^2$.</figcaption></figure>

We note that while the value of $C_{\varepsilon}^i$ may vary as $\varepsilon$ or $D_i^{BM}$ vary independently, varying $\varepsilon$ while maintaining the proportionality relation $D_i^{BM} \propto \varepsilon^2$ does not change the value of $C_{\varepsilon}^i$ providing the best approximation using (A.14) (cf. Figure 9, for which up to 16 orders of magnitude were tested for $\varepsilon$).

We introduce the following parameters: the diffusive permeability of the molecule denoted by the index $i$ across the BM (crossing the BM in the $x_2$ direction) $\kappa_i$, and the corresponding transmission rate of the molecule to the BM $\hat{\kappa}_i$, defined as

$$ \kappa_i = \frac{C_{\varepsilon}^i D_i^{BM}}{\varepsilon} \quad \text{and} \quad \hat{\kappa}_i = \frac{\kappa_i}{\varepsilon} = \frac{C_{\varepsilon}^i D_i^{BM}}{\varepsilon^2}, \quad i = \{t, p\}. \qquad (A.16) $$

**Remark 3.** We emphasise that our method is formal and we do not send $\varepsilon \to 0$, but only retain the dynamics of the average concentration along the direction selected for dimensionality reduction, assuming that $\varepsilon$, i.e. the BM width, is very small.

**Remark 4.** While we adopt a formal method for the dimension reduction of the BM, rigorous methods yielding analogous results have been developed – we refer the interested reader to the method of two-scale convergence for dimension reduction problems, e.g. see [82]. However, we also mention that in other models modelling the effect of thin membranes, transmission is modelled as Kedem-Katchalsky conditions, in which the flux is computed from the permeability and the concentration difference between the sides of the membrane. Here, using this formal dimension reduction, we retain an information about the average of these quantities in the BM in a simple manner.

**Summary of the full reduced model.** Altogether, using the bar notation to denote concentrations in the CT $\Omega_{CT}$ and the hat notation to denote the concentrations in the BM $\Omega_{BM} \equiv \Gamma$ (which in now an
Thus, if $\overline{c_t}^* > 0$ we have that $c_t^* \neq 0$ and from (C.5)$_3$ we retrieve

$$c_a^* = \frac{(K_R - k_2)}{k_2} c_p^* = \frac{k_1 k_3 \bar{c}_d}{k_{-1}(k_{-2} + k_3)} c_p^*, \qquad (C.7)$$

i.e. (C.3)$_3$. Substituting (C.7) in (C.6) we obtain

$$c_t^* = \frac{\hat{\kappa}_t \overline{c_t}^*}{\hat{\kappa}_t + K_R c_p^*}, \qquad (C.8)$$

i.e. (C.3)$_1$. Substituting (C.8) in (C.5)$_2$, multiplying by $(\hat{\kappa}_t + K_R c_p^*)$ and rearranging, we obtain

$$[\hat{\kappa}_p K_R] (c_p^*)^2 + [\hat{\kappa}_p \hat{\kappa}_t + K_R (\hat{\kappa}_t \overline{c_t}^* - \hat{\kappa}_p \overline{c_p}^*)] c_p^* - \hat{\kappa}_p \hat{\kappa}_t \overline{c_p}^* = 0. \qquad (C.9)$$

Thus we have

$$c_p^* = \sqrt{\left(\frac{\hat{\kappa}_t \overline{c_t}^* - \hat{\kappa}_p \overline{c_p}^*}{2\hat{\kappa}_p} + \frac{\hat{\kappa}_t}{2K_R}\right)^2 + \frac{\hat{\kappa}_t \overline{c_t}^*}{K_R} - \left(\frac{\hat{\kappa}_t \overline{c_t}^* - \hat{\kappa}_p \overline{c_p}^*}{2\hat{\kappa}_p} + \frac{\hat{\kappa}_t}{2K_R}\right)},$$

i.e. (C.3)$_2$, where we have chosen the largest root of (C.9) to ensure non-negativity of the proMMP-2 concentration $c_p^*$, for biological relevance, which indeed we have from the non-negativity of the parameter values. Having $c_p^* \ge 0$, we also have that $c_a^* \ge 0$ and $c_t^* \ge 0$ from (C.7) and (C.8), respectively.

On the other hand, if $\overline{c}_t^* = 0$, we automatically have from (C.6) that $c_t^* = 0$. Then, from (C.5)$_2$ we have that $c_p^* = \overline{c}_p^*$, while equation (C.5)$_3$ becomes an identity and any value of $c_a^* = C_a \ge 0$ qualifies as a biologically relevant steady state. Altogether, this gives (C.4). $\square$

### C.1.2 Condition for BM rupture

We consider rupture of the BM to have occurred if the BM density $M$ goes below a critical value $M_{crit}$, low enough to allow cancer cells to pass through. While this may temporarily happen in short finite time, we focus on the case in which BM rupture is sustained in the long-term, i.e.

$$M^*(x) < M_{\text{crit}} < M_{\text{max}} \quad \text{for some } x \in \Gamma, \qquad (C.10)$$

where $M^*$ is the BM density at equilibrium.

**Proposition 8.** *Assuming the system evolves according to the dynamics described by (2.24)–(2.27), under (C.1) and assuming $\overline{c}_t^*(x) > 0$ for some $x \in \Gamma$, long-term sustained BM rupture (C.10) will occur if*

$$C_{RUP}(x) > 0 \quad \text{for some } x \in \Gamma, \qquad (C.11)$$

where $C_{RUP}$ is defined as

$$C_{RUP}(x) = \frac{k_1 k_3 \bar{c}_d}{k_{-1}(k_{-2} + k_3)} c_p^*(x) - \frac{r_M}{\gamma} \left(1 + \frac{K_M}{M_{\text{crit}}}\right) \left(1 - \frac{M_{\text{crit}}}{M_{\text{max}}}\right), \qquad (C.12)$$

with $c_p^*$ given in (C.3)$_2$.

*Proof.* We know from Proposition (7) that at the $x \in \Gamma$ where $\overline{c}_t^*(x) > 0$, the steady state concentrations of TIMP-2, proMMP-2 and active MMP-2 are given by (C.3). For BM rupture (C.10) to occur the following must hold

$$\left. \frac{dM}{dt} \right|_{c_a^*, M_{\text{crit}}} < 0,$$

which from (2.24) gives

$$-\gamma c_a^* \frac{M_{\text{crit}}}{K_M + M_{\text{crit}}} + r_M \left(1 - \frac{M_{\text{crit}}}{M_{\text{max}}}\right) < 0.$$

Rearranging this and plugging in $c_a^*$ from (C.3)$_3$, we obtain (C.11) under definition (C.12). $\square$
open subset of $\mathbb{R}^{d-1}$), we have the model

$$ \partial_t \bar{c}_t = \bar{D}_t \Delta \bar{c}_t + \bar{S}_t(x) - \beta_t \bar{c}_t, \quad \text{in } \Omega_{\text{CT}} \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 \text{in } \Omega_{\text{CT}} \times (0, T), $$

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

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

$$ \partial_t \hat{c}_t = \hat{\kappa}_t (\bar{c}_t - \hat{c}_t) + \hat{R}_t - \beta_t \hat{c}_t, \quad \text{in } \Gamma \times (0, T), $$

$$ \partial_t \hat{c}_p = \hat{\kappa}_p (\bar{c}_p - \hat{c}_p) + \hat{R}_p - \beta_p \hat{c}_p, \quad \text{in } \Gamma \times (0, T), $$

$$ \partial_t \hat{M} = -\gamma \hat{c}_a \frac{\hat{M}}{K_m + \hat{M}} + r_M \left( 1 - \frac{\hat{M}}{M_{\text{max}}} \right)_+ , \quad \text{in } \Gamma \times (0, T), $$

$$ \partial_t \hat{c}_m = \alpha_m \hat{\rho}^0 - k_0 \hat{c}_m^2 + k_{-0} \hat{c}_d - \beta_m \hat{c}_m, \quad \text{in } \Gamma \times (0, T), $$

$$ \partial_t \hat{c}_d = k_0 \hat{c}_m^2 - k_{-0} \hat{c}_d - k_1 \hat{c}_d \hat{c}_t + k_{-1} \hat{c}_1 - \beta_d \hat{c}_d, \quad \text{in } \Gamma \times (0, T), $$

$$ \partial_t \hat{c}_a = k_3 \hat{c}_2 - k_2 \hat{c}_t \hat{c}_a + k_{-2} \hat{c}_{ta} - \beta_a \hat{c}_a, \quad \text{in } \Gamma \times (0, T), $$

$$ \partial_t \hat{c}_1 = k_1 \hat{c}_d \hat{c}_t - \hat{c}_1 [k_1 \hat{c}_t + k_2 \hat{c}_p] + k_{-1} (\hat{c}_3 - \hat{c}_1) + (k_{-2} + k_3) \hat{c}_2 - \beta_1 \hat{c}_1, \quad \text{in } \Gamma \times (0, T), $$

$$ \partial_t \hat{c}_2 = k_2 \hat{c}_1 \hat{c}_p - (k_{-2} + k_3) \hat{c}_2 - \beta_2 \hat{c}_2, \quad \text{in } \Gamma \times (0, T), $$

$$ \partial_t \hat{c}_3 = k_1 \hat{c}_1 \hat{c}_t - k_{-1} \hat{c}_3 - \beta_3 \hat{c}_3, \quad \text{in } \Gamma \times (0, T), $$

$$ \partial_t \hat{c}_{tp} = k_2 \hat{c}_t \hat{c}_p - k_{-2} \hat{c}_{tp} - \beta_{tp} \hat{c}_{tp}, \quad \text{in } \Gamma \times (0, T), $$

$$ \partial_t \hat{c}_{ta} = k_2 \hat{c}_t \hat{c}_a - k_{-2} \hat{c}_{ta} - \beta_{ta} \hat{c}_{ta}, \quad \text{in } \Gamma \times (0, T). $$

We note that in Section 2 the hat is dropped for convenience, while $\hat{\kappa}_i(\bar{c}_i - \hat{c}_i) = S_i(x, t)$ and the remaining biochemical reaction terms introduced in equations (2.4) and (2.5) are denoted by $\hat{R}_i$.

## B Computation details of biologically relevant parameter values

Most baseline parameter values reported in Table 1 are found directly in the literature, short of a change of unit to fit our reference system when required. These come from experimental measurements, were computed using a computational model (different from ours), or were estimated to fit a relevant range. As we use multiple sources to find a priori values for the parameters, we refine our parameter set using optimization techniques to ensure that the outputs of our model are biologically consistent. Below we give details of the references used and computations performed to estimate parameter values which were not found in the literature (marked with a * in Table 1).

**BM carrying capacity.** The BM carrying capacity $M_{\text{max}}$ represents the maximum ECM density in the BM. We consider the BM matrigel density used in *in vitro* invasion assays reported in [2], whose authors estimated a protein concentration of about 10 mg/ml. Considering a molecular weight in the typical weight range of collagen IV, the main component of the ECM found in the BM, i.e. about 160kDa [75], we compute

$$ M_{\text{max}} = \frac{10 \text{ mg/ml}}{160 \times 10^3 \text{ Da}} \approx 6.2510^{-5} \text{ mol/dm}^3 = 62.5 \times 10^3 \text{ nM.} $$

**Physiological levels of TIMP-2, proMMP-2, active MMP-2 and TIMP-2/MMP-2 complex.**

We use as reference, the serum levels reported in [90] for TIMP-2, proMMP-2, active MMP-2 and TIMP-2/MMP-2 complex, which are in good agreement with the values reported by recent studies. Knowing the molecular masses of TIMP-2 (21 kDa, see e.g. [102] and references therein), proMMP-2 (72 kDa [101]), active MMP-2 (67 kDa [101]), and the complex TIMP-2/MMP-2 (88kDa, sum of the two free proteins),
# D Details of the numerical scheme for System (2.1)–(2.13)

In this section, we describe the discretization of our system that combines a cell-centered finite volume spatial discretization and an explicit Runge-Kutta time integration method. We describe separately the case $dim(\Omega) = 2$, $dim(\Gamma) = 1$ and the case $dim(\Omega) = 1$, $dim(\Gamma) = 0$.

We describe the discretization of the prototypal system for TIMP-2

$$ \partial_t \bar{c}_t - D_t \Delta \bar{c}_t = S_t(x, t) - \beta_t \bar{c}_t, \qquad \Omega \times (0, T), \qquad (D.1) $$

$$ -D_t \nabla \bar{c}_t \cdot \nu = \kappa_t (\bar{c}_t - c_t), \qquad \Gamma \times (0, T), \qquad (D.2) $$

$$ -D_t \nabla \bar{c}_t \cdot \nu = 0, \qquad \partial \Omega \setminus \Gamma \times (0, T), \qquad (D.3) $$

$$ \partial_t c_t = \hat{\kappa}_t (\bar{c}_t - c_t) + R_t - \beta_t c_t, \qquad \Gamma \times (0, T). \qquad (D.4) $$

## D.1 Spatial discretization

**Case dim($\Omega$) = 2, dim($\Omega_{BM}$) = 1.** We assume that the domain $\Omega = [0, L] \times [0, L]$ is a square of side length $L$. Let $\mathcal{T} = (K_{ij})_{i=1,...,N_{x_1};j=1,...,N_{x_2}}$ (with $N_{x_1}, N_{x_2}$ the number of nodes in the $x_1$ and $x_2$ direction respectively) be an admissible mesh for $\Omega$ in the sense of [21, Proposition 9.2]. For simplicity, we use equal square control volumes (we also use the denomination "cells" for the control volumes), i.e. $h_i = h = k_j$ for $i = 1, ..., N_{x_1}$, $j = 1, ..., N_{x_2}$ ($h_i, k_j$ being the side lengths of the control volume $K_{i,j}$ in the $x_1$ and $x_2$ directions respectively) and $N_{x_1} = N_{x_2} = N_x$.

We approximate the concentration $\bar{c}_t$ with its cell-averages stored at the cell centers i.e.

$$ \bar{c}_t(x, t) \approx \bar{c}_{t,i,j} = \frac{1}{|K_{i,j}|} \int_{K_{i,j}} \bar{c}_t \, dx. $$

We formulate the finite volume discretization of equation (D.1) in any control-volume $K_{i,j}$: we have the semi-discrete problem

$$ \frac{d\bar{c}_{t,i,j}}{dt} - \frac{1}{|K_{i,j}|} \int_{\partial K_{i,j}} D_t \nabla \bar{c}_{t,i,j} \cdot \nu \, ds = \frac{1}{|K_{i,j}|} \int_{K_{i,j}} S_t \, dx - \beta_t \bar{c}_{t,i,j}. $$

We discretize the flux in an interior cell (i.e. for any $i, j = 2, ..., N_x - 1$) in the standard way using two-point flux approximation

$$ \int_{\partial K_{ij}} D_t \nabla \bar{c}_{t,i,j} \cdot \nu \, ds = h_i (F_{i,j+\frac{1}{2}} - F_{i,j-\frac{1}{2}}) + k_j (F_{i+\frac{1}{2},j} - F_{i-\frac{1}{2},j}), $$

with e.g.

$$ F_{i,j+\frac{1}{2}} = D_t \frac{\bar{c}_{t,i,j+1} - \bar{c}_{t,i,j}}{\frac{1}{2}(k_{j+1} + k_j)}. $$

Furthermore, we approximate the source term as its average on each cell as the concentration, i.e.

$$ \frac{1}{|K_{i,j}|} \int_{K_{i,j}} S_t \, dx = S_{t,i,j}. $$

Altogether, for each interior cell $K_{i,j}$, we have

$$ \frac{d\bar{c}_{t,i,j}}{dt} - \frac{1}{|K_{i,j}|} \left( h_i \left( F_{i,j+\frac{1}{2}} - F_{i,j-\frac{1}{2}} \right) + k_j \left( F_{i+\frac{1}{2},j} - F_{i-\frac{1}{2},j} \right) \right) = S_{t,i,j} - \beta_t \bar{c}_{t,i,j}. $$

As we assumed square cells, the equation simplifies to

$$ \frac{d\bar{c}_{t,i,j}}{dt} - \frac{D_t}{h^2} (\bar{c}_{t,i,j+1} - 4\bar{c}_{t,i,j} + \bar{c}_{t,i,j-1} + \bar{c}_{t,i+1,j} + \bar{c}_{t,i-1,j}) = S_{t,i,j} - \beta_t \bar{c}_{t,i,j}. $$

For a cell touching the boundary $\partial \Omega \setminus \Gamma$, the corresponding flux e.g. $F_{i-1,j}$ for a cell touching the left boundary of $\Omega$, we have $F_{i-1,j} = 0$.

For a cell touching the boundary $\Gamma$, the discretization is different to take into account the Robin boundary condition. We first describe how Equation (D.4) is discretized. The one dimensional domain $\Gamma$ is partitioned into $N_x$ one dimensional cells $C_i$. The TIMP-2 concentration in each of these cells is approximated by its average such that we have

$$ \frac{dc_{t,i}}{dt} = \hat{\kappa}_t (\bar{c}_{t,i,N_y} - c_{t,i}) + R_{t,i} - \beta_t c_{t,i}, $$
we compute the physiological molar concentrations

$$c_t^{\text{ph}} = \frac{201.8 \text{ ng/cm}^3}{21 \times 10^3 \text{ Da}} \approx 9.6 \text{ nmol/dm}^3 = 9.6 \text{nM}$$

$$c_p^{\text{ph}} = \frac{1259.9 \text{ ng/cm}^3}{72 \times 10^3 \text{ Da}} \approx 18 \text{nM},$$

$$c_a^{\text{ph}} = \frac{29.4 \text{ ng/cm}^3}{67 \times 10^3 \text{ Da}} \approx 0.44 \text{nM},$$

$$c_{ta}^{\text{ph}} = \frac{1456.4 \text{ ng/cm}^3}{88 \times 10^3 \text{ Da}} \approx 17 \text{nM}.$$

For the level of TIMP-2/proMMP-2 $c_{tp}^*$, we use the value reported in [4] and obtain

$$c_{tp}^{\text{ph}} = \frac{625 \times 10^3 \text{ ng/L}}{93 \times 10^3 \text{ Da}} \approx 6.72 \text{nM}.$$

For the reduced model, detailed in Appendix A, we assume equilibrium level of MT1-MMP dimers are given. With analogous calculations to those reported above, we take this to be $\bar{c}_d = 8.55/(63 \times 10^3) \approx 0.14$ nM, from normal serum concentration of MT1-MMP given in [73] and molecular mass in [20].

**Remark 5.** In the review article [51], the authors survey the different values of the level of MMP-2, proMMP-2, TIMP-2 and the complex TIMP-2/MMP-2. We observe that, from one study to another, discrepancies in serum or plasma levels of the proteins are found, which could be due to different techniques used to measure these quantities but also due to differences in the cohort. The chosen reference publication for the baseline parameter values used here ensures higher consistency in origin and measurement technique among these values. Moreover, as already mentioned, these values are in good agreement with those reported by recent studies.

**Remark 6.** 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 measurements in tissue samples, we use these values for the moment.

**Physiological production rate of proMMP-2 and TIMP-2 in the CT.** We assume the physiological concentrations of TIMP-2 and pro-MMP2 for a healthy patient are equivalent to the equilibria $c_i^*$ ($i = \{t, p\}$) of their dynamics in the CT (2.13)–(2.12) in the absence of SFs. This yields the relation

$$c_i^{ph} = \bar{c}_i^* = \frac{s_i^{ph}}{\beta_i}, \quad i = \{t, p\}.$$

Using the values for decay rates and physiological levels reported in Table 1, we find

$$s_p^{ph} = 1.09 \times 10^{-2} \text{nM/s},$$

$$s_t^{ph} = 4.28 \times 10^{-4} \text{nM/s}.$$

**Molecular decay rates.** Decay rates $\beta_i$ in System (2.1)–(2.12) represent molecular natural decay, but may also capture the effect of molecule degradation and internalization by other actors, e.g. internalization of membrane bound MT1-MMP [60], or internalization of TIMP-2/proMMP-2 [32]. We assume the decay rates of MT1-MMP dimers and complexes match the one of MT1-MMP monomers, for which we use the rate of MT1-MMP internalization estimated in [55].

We assume physiological levels of TIMP-2/proMMP-2 and TIMP-2/MMP-2 complexes correspond to their concentrations at the equilibrium of system (2.1)–(2.13) in the absence of SFs. We estimate their decay rates from the equilibrium of equations (2.10)–(2.11), using the reaction rates indicated in Table 1, obtaining

$$\beta_{tp} = \frac{k_2 c_t^* c_p^*}{c_{tp}^*} - k_{-2} = \frac{k_2 c_t^{ph} c_p^{ph}}{c_{tp}^{ph}} - k_{-2} = 3.60 \text{s}^{-1},$$

$$\beta_{ta} = \frac{k_2 c_t^* c_a^*}{c_{ta}^*} - k_{-2} = \frac{k_2 c_t^{ph} c_a^{ph}}{c_{ta}^{ph}} - k_{-2} = 0.035 \text{s}^{-1}.$$
for all $i = 1, \dots, N_x$, in which we approximated the integral of the reaction term $R_t$ over the cell $C_i$, is approximated by its average denoted $R_{t,i}$.

Assuming that the BM, $\Gamma$, corresponds to the top boundary of $\Omega$. For any $K_{i,j}$ with $i = 1, \dots, N_x$ and $j = N_y$, we have

$$F_{i,j+\frac{1}{2}} = D_t \nabla c_{t,i,j+\frac{1}{2}} = -\kappa_t (\bar{c}_{t,i,j} - c_{t,i}).$$

**Case dim($\Omega$) = 1, dim($\Omega_{BM}$) = 0.** We assume that the BM, denoted by the boundary $\Gamma$, is located at the right-hand side of the domain, i.e. $x = L$. As in the previous subsection, we use the finite volume method and partition the domain $\Omega = [0, L]$ ($L$ being the size of the domain) in equally sized control volumes $\mathcal{T} = (K_i)_{i=1,\dots,N_x}$. The size of each cell is denoted $h$. We approximate the concentration $\bar{c}_t$ in each cell with its average $c_{t,i}$. The semi-discrete form of Equation (D.1) is

$$\frac{d\bar{c}_{t,i}}{dt} - \frac{D_t}{h} (F_{i+\frac{1}{2}} - F_{i-\frac{1}{2}}) = S_{t,i} - \beta_t \bar{c}_{t,i},$$

for $i = 2, \dots, N_x - 1$. The flux is defined by

$$F_{i+\frac{1}{2}} = \frac{\bar{c}_{t,i+1} - \bar{c}_{t,i}}{h}, \quad F_{i-\frac{1}{2}} = \frac{\bar{c}_{t,i} - \bar{c}_{t,i-1}}{h}.$$

for the leftmost cell, i.e. $i = 1$, the flux $F_{i-\frac{1}{2}} = 0$ For the rightmost cell, the flux $F_{i+\frac{1}{2}} = -\kappa_t (\bar{c}_{t,i} - c_t)$, with $c_t$ being the concentration in the BM given by the ODE (D.4).

## D.2 Time integration scheme

Our semi-discrete system is solved using the ode15s function of the MATLAB ODE suite [60]. This function is a variable step, variable order integration method used to integrate stiff problems (see [60] for details about the algorithm).
The other decay rates, assumed to model natural decay, are computed from the half-life $t_{i,1/2}$ of the molecule via the formula $\beta_i = \ln(2)/t_{i,1/2}$. We know from [15] that the half-life of TIMP-2 in blood is $t_{t,1/2} = 4.33$ hours, and compute

$$ \beta_t = \frac{\ln(2)}{4.33 \times 3600} \approx 4.45 \times 10^{-5} \text{s}^{-1}. $$

In [83], the authors reported a half-life in the plasma of $t_{p,1/2} = 19$ min for human recombinant proMMP-2 injected in mice deficient in MMP-2. As reported in [25], active MMP-9 has a shorter half-life compared to the pro and inactive enzyme. Hence, we assume that $t_{a,1/2} = 0.5t_{p,1/2} = 9.5$ min. This gives the coefficients

$$ \beta_p = \frac{\ln(2)}{19 \times 60} \approx 6.08 \times 10^{-4} \text{s}^{-1}, $$

and

$$ \beta_a = \frac{\ln(2)}{9.5 \times 60} \approx 1.22 \times 10^{-3} \text{s}^{-1}. $$

**Diffusive permeabilities and transmission rates at the BM.** As detailed in Appendix A, we consider a BM of width $\varepsilon = 2 \times 10^{-6}$ dm, estimated for the human epidermal BM via new measuring techniques, as reported in [96]. We consider diffusion coefficients of TIMP-2 and proMMP-2 within the ECM of the CT estimated in [22], i.e. $\bar{D}_i = 1.29 \times 10^{-6} \text{dm}^2/\text{s}$. Given the much higher density of ECM in the BM compared to that in the CT, we take the diffusivity of TIMP-2 and proMMP-2 in the BM to be $D_i^{BM} = 1.29 \times 10^{-12} \text{dm}^2/\text{s}$. This is consistent with the assumption that $D_i^{BM} \approx \mathcal{O}(\varepsilon^2/\text{s})$ (while $\bar{D}_i$ is of the order of $\varepsilon$). Finally, to estimate the diffusive permeability of TIMP-2 and proMMP-2 in the BM, $\kappa_i$, and their transmission rate across the BM, $\hat{\kappa}_i$, using definitions (A.16), we take $C_i^i = 3$. The best value to use for this parameter, introduced in the approximation of the concentration gradient in (A.14), was investigated numerically (see Figure 9) Altogether this gives

$$ \kappa_i = \frac{3 \times 1.29 \times 10^{-12} \text{dm}^2}{2 \times 10^{-6} \text{dm}} \approx 1.935 \times 10^{-6} \text{dm/s}, $$

and

$$ \hat{\kappa}_i = \frac{3 \times 1.29 \times 10^{-12} \text{dm}^2}{4 \times 10^{-12} \text{dm}^2} \approx 0.9675 \text{s}^{-1}. $$

## C Reduced system and analytical results

We reduce System (2.1)-(2.11), comprising 11 coupled nonlinear equations, to an analytically tractable model that retains the key biological processes under consideration. We then study the stability of this reduced system and derive necessary conditions for BM rupture, i.e. tumour invasion.

**Simplifying assumptions.** We make the following simplifying assumptions:

SA1. the inhibition of proMMP-2 and active MMP-2 by TIMP-2 can not reverse (i.e. we take $k_{-2}c_{tp} = k_{-2}c_{ta} \approx 0$ in Equations (2.4)-(2.6) and (2.10)-(2.11));

SA2. the activation and inhibition of proMMP-2 are very fast compared to the rest of the dynamics (hence, we use the quasi steady-state assumption in Equations (2.7)-(2.11));

SA3. stable MT1-MMP expression by cancer cells is already in place and MT1-MMP dimerisation occurs rapidly, so that MT1-MMP dimers quickly stabilize to some equilibrium value $\bar{c}_d$, which is kept constant henceforth;

SA4. the natural decay of proMMP-2, TIMP-2 and active MMP-2 inside the BM is negligible compared to the faster activation/degradation dynamics they are involved in (i.e. we take $\beta_p = \beta_t = \beta_a \approx 0$ in Equations (2.4)-(2.6)).
E Supplementary figures

<figure><img src="image_12.png" /><figcaption>Figure S1: Concentration of active MMP-2 at steady state of the reduced ordinary differential equation (ODE) system (2.25)–(2.27), for varying MT1-MMP dimer concentrations $\bar{c}_d$ corresponding to 50% (left), 75% (center) and 200% (right) of the physiological value $c_d^{ph}$ in Table 1. The remaining details are analogous to those of Figure 3.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure S2: BM density for the healthy (green), tumour (blue) and tumour+SF (red) test cases predicted by the full model over the course of 21 days (3 weeks). The remaining details are analogous to those of Figure 4 (right plot).</figcaption></figure>
**Reduced model.** Under these assumptions, System (2.1)-(2.11) reduces to the System (3.2)-(3.4), i.e.

$$
\begin{align*}
\frac{dc_t}{dt} &= \hat{\kappa}_t \bar{c}_t - k_2 c_t (c_p + c_a) - \hat{\kappa}_t c_t, \\
\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, \\
\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,
\end{align*}
$$

where we have introduced the notation $\bar{c}_i \equiv \bar{c}_i(x, t)|_\Gamma$ ($i = \{t, p\}$) for the TIMP-2 and proMMP-2 concentration at the BM coming from the CT, and equation (2.1) for the digestion of the BM, complemented by the quantities obtained from the quasi-state assumptions SA2., i.e.

$$
c_1 = \frac{k_1 \bar{c}_d}{k_{-1}} c_t, \quad c_2 = \frac{k_2 k_1 \bar{c}_d}{k_{-1}(k_{-2} + k_3)} c_t c_p, \quad c_3 = \left(\frac{k_1}{k_{-1}}\right)^2 \bar{c}_d c_t^2, \quad c_{tp} = \frac{k_2 c_t c_p}{\beta_{tp}}, \quad c_{ta} = \frac{k_2 c_t c_a}{\beta_{ta}}.
$$

In the absence of spatial movement in the BM, the reduced system holds pointwise in $x \in \Gamma$.

C.1 Analytical results of the reduced system

We first analyse the steady state of the ODE system (3.2)-(3.4) and its linear stability. We then exploit
these results to derive sufficient conditions for BM rupture while analysing Equation (3.1).

C.1.1 Steady state

For simplicity, we assume the TIMP-2 and proMMP-2 dynamics in the CT have reached an equilibrium
and denote their equilibrium value at the BM by

$$
\bar{c}_i = \bar{c}_i^* \equiv \bar{c}_i^*(x), \quad i = \{t, p\}, \quad x \in \Gamma. \tag{C.1}
$$

**Proposition 7.** If $\bar{c}_t^* > 0$, *denoting*

$$
K_R = k_2 + \frac{k_1 k_2 k_3 \bar{c}_d}{k_{-1}(k_{-2} + k_3)}, \qquad (C.2)
$$

the steady state solution of the ODE system (3.2)-(3.4) under assumption (C.1) is

$$
\left\{
\begin{aligned}
c_t^* &= \frac{\hat{\kappa}_t \bar{c}_t^*}{\hat{\kappa}_t + K_R c_p^*}, \\
c_p^* &= \sqrt{\left( \frac{\hat{\kappa}_t \bar{c}_t^* - \hat{\kappa}_p \bar{c}_p^*}{2\hat{\kappa}_p} + \frac{\hat{\kappa}_t}{2K_R} \right)^2 + \frac{\hat{\kappa}_t \bar{c}_t^*}{K_R}} - \left( \frac{\hat{\kappa}_t \bar{c}_t^* - \hat{\kappa}_p \bar{c}_p^*}{2\hat{\kappa}_p} + \frac{\hat{\kappa}_t}{2K_R} \right), \\
c_a^* &= \frac{k_1 k_3 \bar{c}_d}{k_{-1}(k_{-2} + k_3)} c_p^*.
\end{aligned}
\right.
\quad (C.3)
$$

If $\bar{c}_t^* = 0$, the system (3.2)-(3.4) under assumption (C.1) admits infinitely many steady states in the form

$$
c_t^* = 0, \quad c_p^* = \bar{c}_p^*, \quad c_a^* = C_a, \tag{C.4}
$$

with $C_a \ge 0$ for biological relevance. From the non-negativity of parameters we have that $c_t^*, c_p^*, c_a^* \ge 0$.

Proof. The steady state $Y^* = (c_t^*, c_p^*, c_a^*)$ of system (3.2)-(3.4), under assumption (C.1), satisfies

$$
\left\{
\begin{array}{l}
\hat{\kappa}_t \bar{c}_t^* - k_2 c_t^* (c_p^* + c_a^*) - \hat{\kappa}_t c_t^* = 0, \\
\\
\hat{\kappa}_p \bar{c}_p^* - K_R c_p^* c_t^* - \hat{\kappa}_p c_p^* = 0, \\
\\
(K_R - k_2) c_t^* c_p^* - k_2 c_t^* c_a^* = 0,
\end{array}
\right.
\qquad (C.5)
$$

with $K_R$ defined in (C.2). Rearranging (C.5)$_1$ we obtain

$$
c_t^* = \frac{\hat{\kappa}_t \bar{c}_t^*}{\hat{\kappa}_t + k_2(c_p^* + c_a^*)}. \qquad (C.6)
$$
<figure><img src="image_14.png" /><figcaption>Figure S3: Molecular concentrations in the BM (represented here as the $x_1$ axis) at $t = 0$ days (dotted line), $t \approx 6.79$ hours (dashed line) and $t = 7$ days (solid line) predicted by the full model for the healthy test case. The remaining details are analogous to those of Figure 4 (healthy test case). We observe a decay of both TIMP-2 and proMMP-2 concentrations and a stabilization to steady values. The active MMP-2 concentration decays to 0 exponentially fast from natural decay and inhibition by TIMP-2. As no tumor cells are present in this test case, MT1-MMP concentration decays to 0 and activation of proMMP-2 is null.</figcaption></figure>
Thus, if $\overline{c}_t^* > 0$ we have that $c_t^* \neq 0$ and from (C.5)$_3$ we retrieve

$$c_a^* = \frac{(K_R - k_2)}{k_2} c_p^* = \frac{k_1 k_3 \bar{c}_d}{k_{-1}(k_{-2} + k_3)} c_p^*, \qquad (C.7)$$

i.e. (C.3)$_3$. Substituting (C.7) in (C.6) we obtain

$$c_t^* = \frac{\hat{\kappa}_t \overline{c}_t^*}{\hat{\kappa}_t + K_R c_p^*}, \qquad (C.8)$$

i.e. (C.3)$_1$. Substituting (C.8) in (C.5)$_2$, multiplying by $(\hat{\kappa}_t + K_R c_p^*)$ and rearranging, we obtain

$$[\hat{\kappa}_p K_R] (c_p^*)^2 + [\hat{\kappa}_p \hat{\kappa}_t + K_R (\hat{\kappa}_t \overline{c}_t^* - \hat{\kappa}_p \overline{c}_p^*)] c_p^* - \hat{\kappa}_p \hat{\kappa}_t \overline{c}_p^* = 0. \qquad (C.9)$$

Thus we have

$$c_p^* = \sqrt{\left(\frac{\hat{\kappa}_t \overline{c}_t^* - \hat{\kappa}_p \overline{c}_p^*}{2\hat{\kappa}_p} + \frac{\hat{\kappa}_t}{2K_R}\right)^2 + \frac{\hat{\kappa}_t \overline{c}_t^*}{K_R} - \left(\frac{\hat{\kappa}_t \overline{c}_t^* - \hat{\kappa}_p \overline{c}_p^*}{2\hat{\kappa}_p} + \frac{\hat{\kappa}_t}{2K_R}\right)},$$

i.e. (C.3)$_2$, where we have chosen the largest root of (C.9) to ensure non-negativity of the proMMP-2 concentration $c_p^*$, for biological relevance, which indeed we have from the non-negativity of the parameter values. Having $c_p^* \ge 0$, we also have that $c_a^* \ge 0$ and $c_t^* \ge 0$ from (C.7) and (C.8), respectively.

On the other hand, if $\overline{c}_t^* = 0$, we automatically have from (C.6) that $c_t^* = 0$. Then, from (C.5)$_2$ we have that $c_p^* = \overline{c}_p^*$, while equation (C.5)$_3$ becomes an identity and any value of $c_a^* = C_a \ge 0$ qualifies as a biologically relevant steady state. Altogether, this gives (C.4). □

### C.1.2 Condition for BM rupture

We consider rupture of the BM to have occurred if the BM density $M$ goes below a critical value $M_{crit}$, low enough to allow cancer cells to pass through. While this may temporarily happen in short finite time, we focus on the case in which BM rupture is sustained in the long-term, i.e.

$$M^*(x) < M_{\text{crit}} < M_{\text{max}} \quad \text{for some } x \in \Gamma, \qquad (C.10)$$

where $M^*$ is the BM density at equilibrium.

**Proposition 8.** *Assuming the system evolves according to the dynamics described by (3.1)-(3.4), under (C.1) and assuming $\overline{c}_t^*(x) > 0$ for some $x \in \Gamma$, long-term sustained BM rupture (C.10) will occur if*

$$C_{RUP}(x) > 0 \quad \text{for some } x \in \Gamma, \qquad (C.11)$$

where $C_{RUP}$ is defined as

$$C_{RUP}(x) = \frac{k_1 k_3 \bar{c}_d}{k_{-1}(k_{-2} + k_3)} c_p^*(x) - \frac{r_M}{\gamma} \left(1 + \frac{K_M}{M_{\text{crit}}}\right) \left(1 - \frac{M_{\text{crit}}}{M_{\text{max}}}\right), \qquad (C.12)$$

with $c_p^*$ given in (C.3)$_2$.

*Proof.* We know from Proposition (7) that at the $x \in \Gamma$ where $\overline{c}_t^*(x) > 0$, the steady state concentrations of TIMP-2, proMMP-2 and active MMP-2 are given by (C.3). For BM rupture (C.10) to occur the following must hold

$$\left. \frac{dM}{dt} \right|_{c_a^*, M_{\text{crit}}} < 0,$$

which from (3.1) gives

$$-\gamma c_a^* \frac{M_{\text{crit}}}{K_M + M_{\text{crit}}} + r_M \left(1 - \frac{M_{\text{crit}}}{M_{\text{max}}}\right) < 0.$$

Rearranging this and plugging in $c_a^*$ from (C.3)$_3$, we obtain (C.11) under definition (C.12). □
<figure><img src="image_15.png" /><figcaption>Figure S4: Molecular concentrations at the BM (represented here as the $x_1$ axis) at $t = 0$ days (dotted line), $t \approx 6.79$ hours (dashed line) and $t = 7$ days (solid line) predicted by the full model for the tumour test case. The remaining details are analogous to those of Figure 4 (tumour test case).</figcaption></figure>
# D Details of the numerical scheme for System (2.1)–(2.13)

In this section, we describe the discretization of our system that combines a cell-centered finite volume spatial discretization and an explicit Runge-Kutta time integration method. We describe separately the case $dim(\Omega) = 2$, $dim(\Gamma) = 1$ and the case $dim(\Omega) = 1$, $dim(\Gamma) = 0$.

We describe the discretization of the prototypal system for TIMP-2

$$ \partial_t \bar{c}_t - D_t \Delta \bar{c}_t = S_t(x, t) - \beta_t \bar{c}_t, \qquad \Omega \times (0, T), \qquad (D.1) $$

$$ -D_t \nabla \bar{c}_t \cdot \nu = \kappa_t (\bar{c}_t - c_t), \qquad \Gamma \times (0, T), \qquad (D.2) $$

$$ -D_t \nabla \bar{c}_t \cdot \nu = 0, \qquad \partial \Omega \setminus \Gamma \times (0, T), \qquad (D.3) $$

$$ \partial_t c_t = \hat{\kappa}_t (\bar{c}_t - c_t) + R_t - \beta_t c_t, \qquad \Gamma \times (0, T). \qquad (D.4) $$

## D.1 Spatial discretization

**Case dim($\Omega$) = 2, dim($\Omega_{BM}$) = 1.** We assume that the domain $\Omega = [0, L] \times [0, L]$ is a square of side length $L$. Let $\mathcal{T} = (K_{ij})_{i=1,...,N_{x_1};j=1,...,N_{x_2}}$ (with $N_{x_1}, N_{x_2}$ the number of nodes in the $x_1$ and $x_2$ direction respectively) be an admissible mesh for $\Omega$ in the sense of [34, Proposition 9.2]. For simplicity, we use equal square control volumes (we also use the denomination "cells" for the control volumes), i.e. $h_i = h = k_j$ for $i = 1, ..., N_{x_1}$, $j = 1, ..., N_{x_2}$ ($h_i, k_j$ being the side lengths of the control volume $K_{i,j}$ in the $x_1$ and $x_2$ directions respectively) and $N_{x_1} = N_{x_2} = N_x$.

We approximate the concentration $\bar{c}_t$ with its cell-averages stored at the cell centers i.e.

$$ \bar{c}_t(x, t) \approx \bar{c}_{t, i, j} = \frac{1}{|K_{i,j}|} \int_{K_{i,j}} \bar{c}_t \, dx. $$

We formulate the finite volume discretization of equation (D.1) in any control-volume $K_{i,j}$: we have the semi-discrete problem

$$ \frac{d\bar{c}_{t,i,j}}{dt} - \frac{1}{|K_{i,j}|} \int_{\partial K_{i,j}} D_t \nabla \bar{c}_{t,i,j} \cdot \nu \, ds = \frac{1}{|K_{i,j}|} \int_{K_{i,j}} S_t \, dx - \beta_t \bar{c}_{t,i,j}. $$

We discretize the flux in an interior cell (i.e. for any $i, j = 2, ..., N_x - 1$) in the standard way using two-point flux approximation

$$ \int_{\partial K_{ij}} D_t \nabla \bar{c}_{t,i,j} \cdot \nu \, ds = h_i (F_{i,j+\frac{1}{2}} - F_{i,j-\frac{1}{2}}) + k_j (F_{i+\frac{1}{2},j} - F_{i-\frac{1}{2},j}), $$

with e.g.

$$ F_{i,j+\frac{1}{2}} = D_t \frac{\bar{c}_{t,i,j+1} - \bar{c}_{t,i,j}}{\frac{1}{2}(k_{j+1} + k_j)}. $$

Furthermore, we approximate the source term as its average on each cell as the concentration, i.e.

$$ \frac{1}{|K_{i,j}|} \int_{K_{i,j}} S_t \, dx = S_{t,i,j}. $$

Altogether, for each interior cell $K_{i,j}$, we have

$$ \frac{d\bar{c}_{t,i,j}}{dt} - \frac{1}{|K_{i,j}|} \left( h_i \left( F_{i,j+\frac{1}{2}} - F_{i,j-\frac{1}{2}} \right) + k_j \left( F_{i+\frac{1}{2},j} - F_{i-\frac{1}{2},j} \right) \right) = S_{t,i,j} - \beta_t \bar{c}_{t,i,j}. $$

As we assumed square cells, the equation simplifies to

$$ \frac{d\bar{c}_{t,i,j}}{dt} - \frac{D_t}{h^2} (\bar{c}_{t,i,j+1} - 4\bar{c}_{t,i,j} + \bar{c}_{t,i,j-1} + \bar{c}_{t,i+1,j} + \bar{c}_{t,i-1,j}) = S_{t,i,j} - \beta_t \bar{c}_{t,i,j}. $$

For a cell touching the boundary $\partial \Omega \setminus \Gamma$, the corresponding flux e.g. $F_{i-1,j}$ for a cell touching the left boundary of $\Omega$, we have $F_{i-1,j} = 0$.

For a cell touching the boundary $\Gamma$, the discretization is different to take into account the Robin boundary condition. We first describe how Equation (D.4) is discretized. The one dimensional domain $\Gamma$ is partitioned into $N_x$ one dimensional cells $C_i$. The TIMP-2 concentration in each of these cells is approximated by its average such that we have

$$ \frac{dc_{t,i}}{dt} = \hat{\kappa}_t (\bar{c}_{t,i,N_y} - c_{t,i}) + R_{t,i} - \beta_t c_{t,i}, $$
<figure><img src="image_16.png" /><figcaption>Figure S5: Molecular concentrations at the BM (represented here as the $x_1$ axis) at $t = 0$ days (dotted line), $t \approx 6.79$ hours (dashed line) and $t = 7$ days (solid line) predicted by the full model for the tumour+SF test case. The remaining details are analogous to those of Figure 4 (tumour+SF test case).</figcaption></figure>
for all $i = 1, \dots, N_x$, in which we approximated the integral of the reaction term $R_t$ over the cell $C_i$, is approximated by its average denoted $R_{t,i}$.

Assuming that the BM, $\Gamma$, corresponds to the top boundary of $\Omega$. For any $K_{i,j}$ with $i = 1, \dots, N_x$ and $j = N_y$, we have

$$F_{i,j+\frac{1}{2}} = D_t \nabla c_{t,i,j+\frac{1}{2}} = -\kappa_t (\bar{c}_{t,i,j} - c_{t,i}).$$

**Case dim($\Omega$) = 1, dim($\Omega_{BM}$) = 0.** We assume that the BM, denoted by the boundary $\Gamma$, is located at the right-hand side of the domain, i.e. $x = L$. As in the previous subsection, we use the finite volume method and partition the domain $\Omega = [0, L]$ ($L$ being the size of the domain) in equally sized control volumes $T = (K_i)_{i=1,\dots,N_x}$. The size of each cell is denoted $h$. We approximate the concentration $\bar{c}_t$ in each cell with its average $c_{t,i}$. The semi-discrete form of Equation (D.1) is

$$\frac{d\bar{c}_{t,i}}{dt} - \frac{D_t}{h} (F_{i+\frac{1}{2}} - F_{i-\frac{1}{2}}) = S_{t,i} - \beta_t \bar{c}_{t,i},$$

for $i = 2, \dots, N_x - 1$. The flux is defined by

$$F_{i+\frac{1}{2}} = \frac{\bar{c}_{t,i+1} - \bar{c}_{t,i}}{h}, \quad F_{i-\frac{1}{2}} = \frac{\bar{c}_{t,i} - \bar{c}_{t,i-1}}{h}.$$

for the leftmost cell, i.e. $i = 1$, the flux $F_{i-\frac{1}{2}} = 0$ For the rightmost cell, the flux $F_{i+\frac{1}{2}} = -\kappa_t (\bar{c}_{t,i} - c_t)$, with $c_t$ being the concentration in the BM given by the ODE (D.4).

## D.2 Time integration scheme

Our semi-discrete system is solved using the ode15s function of the MATLAB ODE suite [100]. This function is a variable step, variable order integration method used to integrate stiff problems (see [100] for details about the algorithm).
# Mathematical modelling of the contribution of senescent fibroblasts to basement membrane digestion during carcinoma invasion

## Supplementary material – Literature Review

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

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

# 1 Biological context

## 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 [34, 35].

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 [25].

The risk of being diagnosed for a carcinoma is related to ageing [9, 77]. 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 [42]. 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 [39], from lack of nutrients and/or space. Regulation by the environment may 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 [65, 63]. 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 [53, 77]. 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.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 [26, 52].
E Supplementary figures

<figure><img src="image_11.png" /><figcaption>Figure S1: Concentration of active MMP-2 at steady state of the reduced ordinary differential equation (ODE) system (3.2)–(3.4), for varying MT1-MMP dimer concentrations $\bar{c}_d$ corresponding to 50% (left), 75% (center) and 200% (right) of the physiological value $c_d^{ph}$ in Table 1. The remaining details are analogous to those of Figure 3.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure S2: BM density for the healthy (green), tumour (blue) and tumour+SF (red) test cases predicted by the full model over the course of 21 days (3 weeks). The remaining details are analogous to those of Figure 4 (right plot).</figcaption></figure>
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 [54]. 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 [75]. 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 [73]. In cancer, however, the physiological ratio between TIMPs and proMMPs may be disrupted, favoring tumor invasion. TIMPs are commonly found in the TME [18], and produced by stromal cells such as fibroblasts (see e.g. [11]).

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 [72]. 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 [40, 54, 61]. 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 [76]. 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 [60]. 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 [26]. Among the stromal cells in the TME, senescent fibroblasts (SFs) have been observed to secrete pro- and active MMP1 and MMP-2 [60].

## 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 [8, 55]. Senescent cells accumulate with age in various tissues as in the skin [19]. Among many others, cancer is an age-related disease [9, 77]. In association with the cell-cycle arrest, inducing senescence of cancer cells have 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 [57], 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 [33], 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 [24, 16, 17, 36, 46, 51, 60]

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) [3] that similarly secrete proMMP-2 on top of many growth factors, cytokines such as the interleukin IL-8 [74], ECM components and other enzymes that affect the TME similarly to the senescent fibroblasts SASP. In line with the antagonistic pleiotropy hypothesis [68, 23], 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 [4, 32]. Thus, SFs may promote carcinogenesis and carcinoma progression in aged organisms. The SASP have numerous possible impacts leading to different therapeutic concepts [30, 55]. 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
<figure><img src="image_13.png" /><figcaption>Figure S3: Molecular concentrations in the BM (represented here as the $x_1$ axis) at $t = 0$ days (dotted line), $t \approx 6.79$ hours (dashed line) and $t = 7$ days (solid line) predicted by the full model for the healthy test case. The remaining details are analogous to those of Figure 4 (healthy test case). We observe a decay of both TIMP-2 and proMMP-2 concentrations and a stabilization to steady values. The active MMP-2 concentration decay to 0 exponentially fast from natural decay and inhibition by TIMP-2. As no tumor cells are present in this test case, MT1-MMP concentration decays to 0 and activation of proMMP-2 is null.</figcaption></figure>
matrix digestion or via the remodelling the latter (stiff tracks) [3, 78]. 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.

# 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, 2, 5, 7, 10, 67] 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 [38, 79]. 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 [22, 70] 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 [29]. 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.

## 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 [6, 21, 50, 58, 62, 71]. 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 coworkers [6, 71] 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 [48, 47] or agent-based [37] models, the dynamics stemming from the molecular crosstalk between cancer cells and fibroblasts, leading to fibroblast differentiation into myofibroblasts [48, 47] or CAFs [37], which then enhanced cancer cell proliferation and motility. Kim et al. [48], 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 [47]. 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. [28], 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. [46], 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.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 [56].
<figure><img src="image_14.png" /><figcaption>Figure S4: Molecular concentrations at the BM (represented here as the $x_1$ axis) at $t = 0$ days (dotted line), $t \approx 6.79$ hours (dashed line) and $t = 7$ days (solid line) predicted by the full model for the tumour test case. The remaining details are analogous to those of Figure 4 (tumour test case).</figcaption></figure>
## 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 [44]. 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 [69], 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 [44] was extended by Donze et al. [20] 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. [41] 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 [12] 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 [13] which considered ECM digestion by cancer cells to be regulated by the plasminogen activation system. The model in [12], 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 [64] 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. [43], 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 [28] 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.

## 2.3 Models of cancer invasion mediated by the digestion of the basal membrane

Ribba et al. [66] 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 [59], 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 2.1 [47, 48], 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 [48] 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. [27] 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. [28]. In their model, the authors apply nonlinear Kedem-Katchalsky transmission conditions [45, 49] 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 [28] is assumed to increase over time under the proteolytic action of MMPs, as previously described in Section 2.1. Chaplain et al. [14] presented a formal derivation of biophysically-consistent transmission conditions which, compared to those in [28], allow the cell volume fraction to also be discontinuous across the membrane, while ensuring mass conservation. These conditions, which were
<figure><img src="image_15.png" /><figcaption>Figure S5: Molecular concentrations at the BM (represented here as the $x_1$ axis) at $t = 0$ days (dotted line), $t \approx 6.79$ hours (dashed line) and $t = 7$ days (solid line) predicted by the full model for the tumour+SF test case. The remaining details are analogous to those of Figure 4 (tumour+SF test case).</figcaption></figure>
later extended to the case of multiple populations in [31] and rigorously derived in [15], can be regarded as a nonlinear generalisation of the classical Kedem-Katchalsky interface conditions [45]. In [14], 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.

## 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] P. M. ALTROCK, L. L. LIU, AND F. MICHOR, *The mathematics of cancer: integrating quantitative models*, Nature Reviews Cancer, 15 (2015), pp. 730-745.

[3] 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).

[4] 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. 232236.

[5] 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.

[6] 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.

[7] 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.

[8] 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.

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

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

[11] 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.

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

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

[14] 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.

[15] G. CIAVOLELLA, N. DAVID, AND A. POULAIN, *Effective interface conditions for a porous medium type problem*, Interfaces Free Bound., (2024).
[16] J.-P. COPPE, M. BOYSEN, C. H. SUN, B. J. WONG, M. K. KANG, N.-H. PARK, P.-Y. DESPREZ, J. CAMPISI, AND A. KRTOLICA, *A Role for Fibroblasts in Mediating the Effects of Tobacco-Induced Epithelial Cell Growth and Invasion*, Molecular Cancer Research, 6 (2008), pp. 1085-1098.

[17] J.-P. COPPÉ, P.-Y. DESPREZ, A. KRTOLICA, AND J. CAMPISI, *The senescence-associated secretory phenotype: The dark side of tumor suppression*, Annual Review of Pathology: Mechanisms of Disease, 5 (2010), pp. 99-118. PMID: 20078217.

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

[19] 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.

[20] 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.

[21] R. EFTIMIE AND L. GIBELLI, *A kinetic theory approach for modelling tumour and macrophages heterogeneity and plasticity during cancer progression*, Mathematical Models and Methods in Applied Sciences, 30 (2020), pp. 659-683.

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

[23] D. FABIAN AND T. FLATT, *The evolution of aging*, Nature Education Knowledge, 3 (2011), pp. 1-10.

[24] M. FANE AND A. T. WEERARATNA, *How the ageing microenvironment influences tumour progression*, Nature Reviews Cancer, 20 (2020), p. 89106.

[25] 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).

[26] C. FRANTZ, K. M. STEWART, AND V. M. WEAVER, *The extracellular matrix at a glance*, Journal of Cell Science, 123 (2010), pp. 4195-4200.

[27] A. FRIEDMAN, C. HUANG, AND J. YONG, *Effective permeability of the boundary of a domain*, Communications in Partial Differential Equations, 20 (1995), pp. 59-102.

[28] 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.

[29] R. A. GATENBY AND E. T. GAWLINSKI, *A reaction-diffusion model of cancer invasion*, Cancer research, 56 (1996), pp. 5745-5753.

[30] J. GIROUD, I. BOURIEZ, H. PAULUS, A. POURTIER, F. DEBACQ-CHAINIAUX, AND O. PLUQUET, *Exploring the communication of the sasp: Dynamic, interactive, and adaptive effects on the microenvironment*, International Journal of Molecular Sciences, 24 (2023).

[31] C. GIVERSO, T. LORENZI, AND L. PREZIOSI, *Effective interface conditions for continuum mechanical models describing the invasion of multiple cell populations through thin membranes*, Applied Mathematics Letters, 125 (2022), p. 107708.

[32] L. GROSSE, N. WAGNER, A. EMELYANOV, C. MOLINA, S. LACAS-GERVAIS, K.-D. WAGNER, AND D. V. BULAVIN, *Defined p16high senescent cell types are indispensable for mouse healthspan*, Cell Metabolism, 32 (2020), pp. 87-99.e6.

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

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

[35] ——, *Hallmarks of cancer: the next generation*, cell, 144 (2011), pp. 646-674.
[36] Y. HASSONA, N. CIRILLO, K. HEESOM, E. PARKINSON, AND S. S. PRIME, *Senescent cancer-associated fibroblasts secrete active mmp-2 that promotes keratinocyte dis-cohesion and invasion*, British journal of cancer, 111 (2014), pp. 1230-1237.

[37] 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.

[38] T. HERNANDEZ-BOUSSARD, P. MACKLIN, E. J. GREENSPAN, A. L. GRYSHUK, E. STAHLBERG, T. SYEDA-MAHMOOD, AND I. SHMULEVICH, *Digital twins for predictive oncology will be a paradigm shift for precision cancer care*, Nature medicine, 27 (2021), pp. 2065-2066.

[39] F. HILLEN AND A. W. GRIFFIOEN, *Tumour vascularization: sprouting angiogenesis and beyond*, Cancer and Metastasis Reviews, 26 (2007), pp. 489-502.

[40] U. B. HOFMANN, J. R. WESTPHAL, G. N. VAN MUIJEN, AND D. J. RUITER, *Matrix metalloproteinases in human melanoma*, Journal of Investigative Dermatology, 115 (2000), pp. 337-344.

[41] 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.

[42] W. HU, L. FANG, R. NI, H. ZHANG, AND G. PAN, *Changing trends in the disease burden of non-melanoma skin cancer globally from 1990 to 2019 and its predicted level in 25 years*, BMC Cancer, 22 (2022).

[43] 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.

[44] 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.

[45] O. KEDEM AND A. KATCHALSKY, *Thermodynamic analysis of the permeability of biological membranes to non-electrolytes*, Biochimica et biophysica Acta, 27 (1958), pp. 229-246.

[46] 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.

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

[48] 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.

[49] F. KLEINHANS, *Membrane permeability modeling: Kedem-katchalsky vs a two-parameter formalism*, Cryobiology, 37 (1998), pp. 271-289.

[50] H. KNÚTSDÓTTIR, E. PÁLSSON, AND L. EDELSTEIN-KESHET, *Mathematical model of macrophage-facilitated breast cancer cells invasion*, Journal of theoretical biology, 357 (2014), pp. 184-199.

[51] A. KRTOLICA, S. PARRINELLO, S. LOCKETT, P.-Y. DESPREZ, AND J. CAMPISI, *Senescent fibroblasts promote epithelial cell growth and tumorigenesis: A link between cancer and aging*, Proceedings of the National Academy of Sciences, 98 (2001), pp. 12072-12077.

[52] J. K. KULAR, S. BASU, AND R. I. SHARMA, *The extracellular matrix: Structure, composition, age-related differences, tools for analysis and applications for tissue engineering*, Journal of Tissue Engineering, 5 (2014), p. 204173141455711.
[53] E. LACONI, F. MARONGIU, AND J. DEGREGORI, *Cancer as a disease of old age: changing mutational and microenvironmental landscapes*, British journal of cancer, 122 (2020), pp. 943-952.

[54] H. LARONHA AND J. CALDEIRA, *Structure and function of human matrix metalloproteinases*, Cells, 9 (2020), p. 1076.

[55] Z. LIAO, H. L. YEO, S. W. WONG, AND Y. ZHAO, *Cellular senescence: Mechanisms and therapeutic potential*, Biomedicines, 9 (2021).

[56] K. V. LU, K. A. JONG, A. K. RAJASEKARAN, T. F. CLOUGHESY, AND P. S. MISCHEL, *Upregulation of tissue inhibitor of metalloproteinases (timp)-2 promotes matrix metalloproteinase (mmp)-2 activation and cell invasion in a human glioblastoma cell line*, Laboratory investigation, 84 (2004), pp. 8-20.

[57] I. LÜHR, A. FRIEDL, T. OVERATH, A. THOLEY, T. KUNZE, F. HILPERT, S. SEBENS, N. ARNOLD, F. RÖSEL, H.-H. OBERG, N. MAASS, C. MUNDHENKE, W. JONAT, AND M. BAUER, *Mammary fibroblasts regulate morphogenesis of normal and tumorigenic breast epithelial cells by mechanical and paracrine signals*, Cancer Letters, 325 (2012), pp. 175-188.

[58] G. MAHLBACHER, L. T. CURTIS, J. LOWENGRUB, AND H. B. FRIEBOES, *Mathematical modeling of tumor-associated macrophage interactions with the cancer microenvironment*, Journal for immunotherapy of cancer, 6 (2018), pp. 1-17.

[59] K. M. MAK AND R. MEI, *Basement membrane type iv collagen and laminin: An overview of their biology and value as fibrosis biomarkers of liver disease*, The Anatomical Record, 300 (2017), pp. 1371-1390.

[60] N. MALAQUIN, C. VERCAMER, F. BOUALI, S. MARTIEN, E. DERUY, N. WERNERT, M. CHWASTYNIAK, F. PINET, C. ABBADIE, AND A. POURTIER, *Senescent fibroblasts enhance early skin carcinogenic events via a paracrine mmp-par-1 axis*, PLOS ONE, 8 (2013), pp. 1-17.

[61] G. MURPHY, *Tissue inhibitors of metalloproteinases*, Genome biology, 12 (2011), pp. 1-7.

[62] M. R. OWEN AND J. A. SHERRATT, *Mathematical modelling of macrophage dynamics in tumours*, Mathematical Models and Methods in Applied Sciences, 9 (1999), pp. 513-539.

[63] T. K. PALANIAPPAN, L. ŠLEKIENĖ, A.-K. JONASSON, J. GILTHORPE, AND L. GUNHAGA, *Camdelam: an in vivo approach to visualize and quantify the delamination and invasion capacity of human cancer cells*, Scientific Reports, 10 (2020), p. 10472.

[64] J. PANOVSKA AND G. PETTET, *Mathematical model for glioma cell invasion*. Unpublished, 2007.

[65] M. PIIPPONEN, P. RIIHILÄ, L. NISSINEN, AND V.-M. KÄHÄRI, *The role of p53 in progression of cutaneous squamous cell carcinoma*, Cancers, 13 (2021), p. 4507.

[66] B. RIBBA, O. SAUT, T. COLIN, D. BRESCH, E. GRENIER, AND J. P. BOISSEL, *A multiscale mathematical model of avascular tumor growth to investigate the therapeutic benefit of anti-invasive agents*, J. Theoret. Biol., 243 (2006), pp. 532-541.

[67] R. C. ROCKNE, A. HAWKINS-DAARUD, K. R. SWANSON, J. P. SLUKA, J. A. GLAZIER, P. MACKLIN, D. A. HORMUTH, A. M. JARRETT, E. A. LIMA, J. T. ODEN, ET AL., *The 2019 mathematical oncology roadmap*, Physical biology, 16 (2019), p. 041005.

[68] M. ROSE, *Evolutionary Biology of Aging*, Oxford University Press, 1994.

[69] T. SAITOU, K. ITANO, D. HOSHINO, N. KOSHIKAWA, M. SEIKI, K. ICHIKAWA, AND T. SUZUKI, *Control and inhibition analysis of complex formation processes*, Theoretical Biology and Medical Modelling, 9 (2012), pp. 1-18.

[70] N. SFAKIANAKIS AND M. A. J. CHAPLAIN, *Mathematical modelling of cancer invasion: A review*, in Methods of Mathematical Oncology, T. Suzuki, C. Poignard, M. Chaplain, and V. Quaranta, eds., vol. 370, Springer Singapore, Singapore, 2021, p. 153172.

[71] S. SUVEGES, R. EFTIMIE, AND D. TRUCU, *Directionality of macrophages movement in tumour invasion: A multiscale moving-boundary approach*, Bulletin of Mathematical Biology, 82 (2020), pp. 1-50.
[72] M. TAURO AND C. C. LYNCH, *Cutting to the chase: how matrix metalloproteinase-2 activity controls breast-cancer-to-bone metastasis*, Cancers, 10 (2018), p. 185.

[73] T. H. VU AND Z. WERB, *Matrix metalloproteinases: effectors of development and normal physiology*, Genes & development, 14 (2000), pp. 2123-2133.

[74] T. WANG, F. NOTTA, R. NAVAB, J. JOSEPH, E. IBRAHIMOV, J. XU, C.-Q. ZHU, A. BORGIDA, S. GALLINGER, AND M.-S. TSAO, *Senescent Carcinoma-Associated Fibroblasts Upregulate IL8 to Enhance Prometastatic Phenotypes*, Molecular Cancer Research, 15 (2017), pp. 3-14.

[75] H. E. V. WART AND H. BIRKEDAL-HANSEN, *The cysteine switch: a principle of regulation of metalloproteinase activity with potential applicability to the entire matrix metalloproteinase gene family.*, Proceedings of the National Academy of Sciences, 87 (1990), pp. 5578-5582.

[76] A. WATANABE, D. HOSINO, N. KOSHIKAWA, M. SEIKI, T. SUZUKI, AND K. ICHIKAWA, *Critical role of transient activity of mt1-mmp for ecm degradation in invadopodia*, PLOS Computational Biology, 9 (2013), pp. 1-13.

[77] M. C. WHITE, D. M. HOLMAN, J. E. BOEHM, L. A. PEIPINS, M. GROSSMAN, AND S. J. HENLEY, *Age and cancer risk: a potentially modifiable relationship*, American journal of preventive medicine, 46 (2014), pp. S7-S15.

[78] E. C. WOENNE, W. LEDERLE, S. ZWICK, M. PALMOWSKI, H. KRELL, W. SEMMLER, M. M. MUELLER, AND F. KIESSLING, *Mmp inhibition blocks fibroblast-dependent skin cancer invasion, reduces vascularization and alters vegf-a and pdgf-bb expression*, Anticancer Research, 30 (2010), pp. 703-711.

[79] C. WU, G. LORENZO, D. A. HORMUTH, E. A. B. F. LIMA, K. P. SLAVKOVA, J. C. DICARLO, J. VIROSTKO, C. M. PHILLIPS, D. PATT, C. CHUNG, AND T. E. YANKEELOV, *Integrating mechanism-based modeling with biomedical imaging to build practical digital twins for clinical oncology*, Biophysics Reviews, 3 (2022), p. 021304.