Multiphysics modeling of 3D traction force microscopy with application to cancer cell-induced degradation of the extracellular matrix

Vol.:(0123456789) Engineering with Computers https://doi.org/10.1007/s00366-024-02017-8 ORIGINAL ARTICLE Multiphysics modeling of3D traction force microscopy withapplication tocancer cell‑induced degradation oftheextracellular matrix AlejandroApolinar‑Fernández1· JorgeBarrasa‑Fano2· HansVanOosterwyck2,3· JoséA.Sanz‑Herrera1 Received: 29 September 2023 / Accepted: 21 June 2024 © The Author(s) 2024 Abstract 3D Traction Force Microscopy (3DTFM) constitutes a powerful methodology that enables the computation of realistic forces exerted by cells on the surrounding extracellular matrix (ECM). The ECM is characterized by its highly dynamic structure, which is constantly remodeled in order to regulate most basic cellular functions and processes. Certain pathological processes, such as cancer and metastasis, alter the way the ECM is remodeled. In particular, cancer cells are able to invade its surrounding tissue by the secretion of metalloproteinases that degrade the extracellular matrix to move and migrate towards different tissues, inducing ECM heterogeneity. Typically, 3DTFM studies neglect such heterogeneity and assume homogeneous ECM properties, which can lead to inaccuracies in traction reconstruction. Some studies have implemented ECM degradation models into 3DTFM, but the associated degradation maps are defined in an ad hoc manner. In this paper, we present a novel multiphysics approach to 3DTFM with evolving mechanical properties of the ECM. Our modeling considers a system of partial differential equations based on the mechanisms of activation of diffusive metalloproteinase MMP2 by membrane-bound metalloproteinase MT1-MMP. The obtained ECM density maps in an ECM-mimicking hydrogel are then used to compute the heterogeneous mechanical properties of the hydrogel through a multiscale approach. We perform forward and inverse TFM simulations both accounting for and omitting degradation, and results are compared to ground truth reference solutions in which degradation is considered. The main conclusions resulting from the study are: (i) the inverse methodology yields results that are significantly more accurate than those provided by the forward methodology; (ii) ignoring ECM degradation results in a considerable overestimation of tractions and non negligible errors in all analyzed cases. Keywords Traction force microscopy· Metalloproteinase-induced ECM degradation· Mechanobiology· Nonlinear mechanics· Computational mechanics· Finite element method 1 Introduction The field of mechanobiology studies the influence of mechanical phenomena on the origin and development of processes associated with biological cells and tissues. In the last years, researchers have put efforts on understanding how physiological forces and mechanical conditions of the cell can alter its cellular function, and can change the structure and composition of the extracellular matrix (ECM). In this context, it has been repeatedly found that the mechanical interactions of cells with the extracellular matrix are pivotal in the development of multiple physiological and pathological processes [1]. Such processes include cancer invasion [2], angiogenesis [3], stem cell differentiation [4] as well as other tissue or organ pathologies [5]. However, the mechanisms by which cells sense forces and * José A. Sanz-Herrera [email protected] Alejandro Apolinar-Fernández [email protected] Jorge Barrasa-Fano jorge.barrasaf[email protected] Hans Van Oosterwyck [email protected] 1 Escuela Técnica Superior de Ingeniería, Universidad de Sevilla, Camino de los descubrimientos s/n, Seville41092, Spain 2 Biomechanics section, Department ofMechanical Engineering, KU Leuven, Leuven, Belgium 3 Prometheus division ofSkeletal Tissue Engineering, KU Leuven, Leuven, Belgium Engineering with Computers convert such mechanical signals into biological responses, i.e. mechanotransduction, remain unclear. In order to obtain fundamental understanding of cellular mechanotransduction pathways that could lead to the development of innovative tissue engineering and therapeutic techniques, invitro models with human-like tissue composition, organization, and mechanical behavior are devised [6]. The extracellular matrix is a highly dynamic structure that exists under the influence of a constant remodeling process, where ECM components are deposited, degraded, or otherwise modified. It constitutes a major component of a cell microenvironment, taking part in most basic cell behaviors, from cell proliferation, adhesion and migration, to cell differentiation and death [7]. In the context of continuously changing cellular activity, remodeling of the extracellular matrix represents a crucial mechanism whereby diverse cellular behaviors are regulated [8]. Indeed, normal development, physiology, and robustness of organ systems require ECM dynamics to be tightly regulated. Alteration or dysfunction of such regulation leads to congenital defects and diseases, including cancer. The study of the process of cancer mestastasis, in particular invasive cellular migration towards surrounding tissue by means of degradation of the collagen microstructure that supports the ECM, is of great interest. Indeed, cells are able to induce ECM degradation through the secretion of matrix metalloproteinases (MMPs), a type of proteolytic enzymes which have the capacity to degrade (potentially) all components of the surrounding tissue [9]. MMP-induced degradation hence facilitates cancer growth and spread by virtue of the new available space generated, as well as by other proteins that are released by the degraded tissue [10, 11]. Cellular mechanical interactions with the ECM are typically quantified by combining optical microscopy with computational methods, including microscopy image analysis techniques and numerical models. Over the last decades, Traction Force Microscopy (TFM) has been one of the preferred methodologies for the computation of cell forces exerted on the ECM [12–15]. In TFM, cultured cells in artificial substrates with known mechanical properties, such as hydrogels, are seeded with a number of fiducial markers (often fluorescent beads). These beads track the motion generated within the substrate by a mechanically active cell from a reference configuration. Microscopy imaging is used to capture the stressed state of the hydrogel (during cell activity), as well as the relaxed state of hydrogel (after mechanical inhibition of the cell). Then, image processing techniques are used for acquiring the deformations from the tracked motion field. Finally, strain, stress and traction forces can be calculated by employing a constitutive model that describes the mechanical behavior of the substrate/hydrogel. TFM has been first formulated in the context of 2D invitro cultures, where cells are seeded on the top of a substrate, and in which simple analytical formulations that facilitate traction reconstruction are available [16, 17]. Nonetheless, cell microenvironment dimensionality heavily influences cell behavior and 3D invitro cultures are therefore established as more physiologically relevant [1, 18, 19]. The field of 3D TFM has lead to the development of complex and advanced methodologies during the last decade. There are two main approaches to traction force calculation in TFM: the forward and inverse methods. On the one hand, forward methods, that prescribe all or a subset of the measured displacements, to obtain the associated strain field from which to retrieve stresses and tractions in a direct manner, through the selected constitutive relation for the ECM [20–26]. On the other hand, inverse methods search for a traction solution that is coherent with the measured displacements while being subjected to a certain constraint/regularization that limits the effect of displacements measurement noise in the resultant traction field [27–32]. The majority of studies involving 3D TFM have been carried out under the assumption that the ECM has spatially homogeneous mechanical properties. However, such assumption becomes questionable, especially for longer term experiments (multiple hours/days) in which cells have time to locally degrade the surrounding matrix via MMPs, and where the resulting mechanical heterogeneity may impact on the accuracy of traction recovery. Several experimental studies have highlighted the importance of ECM heterogeneity when measuring cellular stresses [23, 33]. This motivates the development of new 3D TFM approaches taking into account the heterogeneity of ECM mechanical properties as well as their evolution in time. As already mentioned, there are few TFM studies in which heterogeneity of the ECM is considered. In [34] the authors implement a methodology for solving 3DTFM cases while considering ECM heterogeneity, although not due to ECM degradation. Nevertheless, the applicability of the methodology is limited, since they consider small strains and assume the material to behave as linear elastic. Also, the authors utilize a forward method to perform traction reconstruction, which, as it was concluded in other studies [29, 30], exhibits a significantly reduced efficacy when compared to inverse methods. In contrast, our methodology accounts for finite strains and nonlinear constitutive behavior of the ECM; and tractions are recovered by means of both forward and inverse methods. Song etal. [35] proposed an in silico study in which the authors perform 3D TFM while accounting for heterogeneity of matrix mechanical properties, considering nonlinearity of the ECM mechanical behavior and conducting traction reconstruction through an inverse method. However, the law that correlates ECM degradation with the material parameters is imposed by the authors, which may not represent a realistic case. On the contrary, our approach considers the simulation of realistic Engineering with Computers time-dependent ECM density maps by implementing a thorough mathematical model based on enzymatic reactions in 3D TFM within the framework of MMP-induced ECM degradation. Moreover, our approach explicitly takes into consideration the ECM microstructure to establish the relationship between ECM degradation and corresponding material parameters. There have been several attempts of developing mathematical models that involve MMP-induced ECM degradation. More specifically, models regarding the activation of MMP-2 by MT1-MMP can be seen in [36–38]. Deakin and Chaplain [39] devised a model based on a simplified version of the MMP2 activation process in which the spatial diffusion of species was incorporated. However, the model is referred to cell populations in 2D planar substrates and is not used in the context of complex 3D individual cell geometries. In this study, we present a complete 3D TFM framework which considers ECM degradation retrieved from the solution of a set of equations belonging to a previously described MMP-induced degradation mathematical model [39]. The model was adapted for its application to individual complex cell geometries in 3D. Degradation is then correlated with a nonlinear hyperelastic model, that was calibrated with synthetic shear rheology tests corresponding to different levels of ECM degradation. The shear rheology curves were obtained from computational (virtual) simple shear tests applied to synthetically generated fibered matrices. Shear rheology tests corresponding to real collagen hydrogels tested in the laboratory, were used to fit and calibrate such a virtual tool. As a result, TFM accuracy is assessed in these heterogeneous matrices referred to traction reconstruction versus a (controlled) synthetic ground truth solution in which degradation of the matrix is included. The analysis of cellular forces in several scenarios is then conducted: under the presence and absence of ECM degradation, and forward / inverse reconstructions for each case of study. The next sections include: (i) a brief description of the study, (ii) the mathematical model for MMP-induced degradation of the ECM and its computational implementation, (iii) the virtual shear rheology tests and their corresponding synthetic data using our fibered matrix generator tool, and (iv) the multiscale assembly of these methodologies in our forward/ inverse TFM simulations. Finally, the obtained results are discussed and the main conclusions of the study are highlighted at the end of the paper. 2 Description ofthestudy The present study follows a weakly coupled multiscale approach. Two separated scales (macro- and microscales) are defined. At the macroscale, ECM density maps are extracted from the resolution of a mathematical model implemented to simulate metalloproteinase-induced collagen degradation. The equations are solved in a 3D (macro) domain, composed of a real cell geometry embedded within an ECM. The macroscopic mechanical behavior of the degraded material macro (Gauss) point is given through a synthetic shear rheology curve, corresponding to the level of ECM density of this macro point. This curve is obtained by means of the representation of the underlying fibered microstructure, using a matrix generator and a homogenization scheme developed by the authors. The parameters involved in these numerical algorithms are calibrated with real shear rheology data, resulting from collagen hydrogels tested in the laboratory. At the macroscale, the space-dependent mechanical behavior of the matrix is treated by using a macroscopic nonlinear hyperelastic model, in which the parameters are fitted heterogeneously with the referred synthetically obtained shear rheology data at each macroscopic point. The use of a previously calibrated (offline) macroscopic nonlinear hyperelastic model to represent the macroscopic mechanical behavior of the ECM through the underlying microstructure turns the approach into a weakly coupled multiscale approach, and it avoids expensive FE 2 simulations (i.e. multiscale finite element coupled analyses at 2 different scales of observation). The heterogeneous macroscopic mechanical model of the ECM, obtained as referred above, is used to perform TFM simulations. First, a ground truth (baseline case) solution is built for an individual cell embedded in an ECM that is subjected to metalloproteinase-induced collagen degradation. Then, an artificial noise is introduced in the recovered displacement field of the ground truth solution. This perturbed displacement field is used to recover tractions in two different scenarios: (i) assuming a homogeneous ECM, and (ii) assuming the degraded and heterogeneous ECM at different timepoints. The results of these two cases are then compared and discussed versus the ground truth solution. Figure1 shows a schematic outline of the study. 3 Modeling metalloproteinase‑driven ECM degradation In this section we present the mathematical model for simulating collagen degradation of the extracellular matrix due to its interaction with diffusive metalloproteinase MMP-2 and membrane-bound metalloproteinase MT1-MMP. The model has been adapted from the one proposed by [39], applied in our case for a single cell within a 3D geometry. More details about the original implementation of the model can be found in [40]. Engineering with Computers 3.1 Phenomenological description The model considers two different degradation pathways associated with the two species of metalloproteinases involved in the remodeling process: the highly localized tissue degradation produced by MT1-MMP, as well as the diffusive degradation corresponding to the activity of MMP-2. MT1-MMP is produced by the cell and is bound to its membrane. MMP-2 is activated by MT1- MMP through its interaction with the tissue inhibitor TIMP2, which consequently acts, dually, as an inhibitor of both species, and as an activator of MMP-2. The mechanism that mediates this process is illustrated in Fig.2a, as well as outlined in more detail in Fig.2b. Fundamentally, TIMP2 intercepts MT1-MMP, inhibiting its activity. The resulting complex MT1-MMP:TIMP2 interacts with MMP-2 proenzyme (pMMP-2), an inactive precursor of MMP-2. As a consequence, the complex MT1- MMP:TIMP2:pMMP-2 is formed, which, upon interacting with MT1-MMP, triggers the production of the enzyme MMP-2. Species MT1-MMP, TIMP2, and pMMP-2 are produced, while species MMP-2 and the intermediate complexes described above are formed as a result of the interactions between the former. In addition, the intermediate complexes can dissociate into their components. 3.2 Mathematical model The model is based on a simplified version of the MMP-2 activation mechanism in which the intermediate phases regarding the formation of complexes MT1-MMP:TIMP2, MT1-MMP:TIMP2:pMMP-2 and MT1-MMP:TIMP2:pMMP- 2:MT1-MMP, are summarized into one single phase containing a generic intermediate complex formed from the combination of MT1-MMP and TIMP2. These two complexes, upon interacting with MT1-MMP, trigger the production of MMP-2 and MT1-MMP, as illustrated by the diagram in Fig.2c. This generic complex is able to dissociate into MT1-MMP and TIMP2. The model is composed of five nonlinear and coupled PDEs of the reaction-diffusion type, (1a) 𝜕v 1 𝜕t =−𝛿1c[v2+𝛿2v3]+𝜇v[1−c−v1 ] (1b) 𝜕v 2 𝜕t =∇⋅[Dv2∇v2]−k21v4v2+k22v3v5−𝛽v2 v 2 (1c) 𝜕v 3 𝜕t =−k31v4v3+k32v5−𝛽v3 v3+𝛼v3 c[1+v1 ] Fig. 1 Outline of the study and its different stages Engineering with Computers Fig. 2 a Illustration of the MMP-2 activation process. Created with BioRender. b Schematic representation of the MMP-2 activation process. c Schematic representation of the simplified MMP-2 activation process Engineering with Computers where c represents the cancer density, v1 the ECM density, and v2 , v3 , v4 and v5 the concentration of MMP-2, MT1- MMP, TIMP2 and the generic complex, respectively. The positive terms at the right side of the equations correspond to the interactions between species/complexes that lead to the production of the respective variables, whilst the negative terms represent the consumption due to inhibition/dissociation of the involved variables. Apart from the referred mechanisms of production/consumption, v3 and v4 are also spontaneously produced at a rate depending on the level of cancer density. Both c and v1 range between values 0 and 1, 0 meaning total absence/deterioration of cancer/ECM, and 1 meaning complete presence/integrity of cancer/ECM, respectively. v1 is reduced (ECM degradation) by v2 and v3 , and is increased (ECM regeneration) in regions where c is greater than 0 and v1 is less than 1. Only v2 and v4 diffuse within the domain, since MMP-2 and TIMP2 are not bound to the cell surface (see Fig.2a). The variable of interest in the context of this study is v1 (ECM density), that will be related in Sect.4 to the heterogenous mechanical properties of the hydrogel. Variable c (cancer density) is conceived as a fixed stationary function that indicates (value of 1) the regions in which MMPs are emitted, and therefore degradation is initiated, as well as the regions in which degradation is not prescribed (value of 0). This is further explained in Sect.3.3.2. The equations presented in (1a–1f) are dimensionless quantities. 3.3 Computational implementation The model was implemented in COMSOL Multiphysics, which offers a powerful platform for implementing systems of PDEs in problems with complex 3D geometries. 3.3.1 Geometry andmesh We extracted the geometry of a real cell from confocal microscopy images obtained in the laboratory. This geometry was then represented in COMSOL, and embedded within a 3D cartesian domain. The cell volume was subtracted from the domain and the resulting geometry (a rectangular prism that represents the ECM with a cell-shaped void in the middle) was meshed. Figure3 shows a 3D representation of the domain of analysis and the resulting mesh. (1d) 𝜕v 4 𝜕t =∇⋅[Dv4∇v4]−k41v4v2−k42v4v3+k43v5+𝛼v4 c (1e) 𝜕v 5 𝜕t =k51v4v3−k52v5v3−k53v 5 (1f) c=c(x,y,z) 3.3.2 Definition ofc, boundary andinitial conditions We consider the cancer density (c) quantity as a binary variable taking the value of 1 where MMPs are emitted in the cell, i.e. certain regions on the cell boundary of the domain, and 0 elsewhere. Thus, the ECM degradation process starts in the prescribed regions (named MMP-emitting regions) along the cell surface, which in the case of this study correspond to the tips of the cell protrusions (from which MMPs are emitted). As a consequence, this variable is implemented as a prescribed field variable constant in time, that is, a time independent function c(x,y,z) which takes the value 1 in the MMP-emitting regions and 0 at any other point. Since the problem is continuous and in order to avoid convergence problems, the transition between 0 and 1 was made smooth with continuous smoothed step functions ( C2 functions) in the direction normal to the cell surface, and in the axial direction of the cell protrusion, respectively. The variable v1 was set as free in the cell surface, i.e., its flux was set to zero (homogeneous Neumann boundary condition). Regarding the rest of variables, Dirichlet boundary conditions were defined by assigning to the variables a value of 1 at MMP- emitting regions on the cell surface and 0 elsewhere (the transition between the two values was made smooth as in the case of c). All variables fluxes were set to zero at the external boundaries of the domain as homogeneous Neumann boundary conditions, assuming that the domain is isolated. The initial condition for v1 was set as 1 (intact ECM), and as c for the rest of variables. The cancer density variable (c) definition, as well as the regions in which it is applied (and equivalently, the previously described boundary conditions) are represented in Fig.3. 3.3.3 Results A degradation process throughout 28 dimensionless time steps (from t=0 to t=27) was simulated. The values of the parameters that govern the proposed degradation model (Eqs. (1a–1f)) were fitted in order to obtain a maximum ECM degradation level of 50% (ECM density level of 50%) at the end of the degradation process at t=27 . Since the scope of this study is computational, the calibrated parameters do not necessarily coincide with values found in the literature, although some of them remain unchanged with respect to the ones used in [39, 40], which are either extracted from the literature or estimated by the authors. Table1 contains the fitted nondimensional values for the parameters. Parameters 𝛿1 , 𝛿2 , 𝜇v and 𝛼v3 were calibrated by means of a trial and error procedure. Dv2 and Dv4 were modified only with respect to their non-dimensional values according to the dimensions of our particular problem, that are different to those presented in [39, 40]. The resultant ECM density maps throughout the simulated Engineering with Computers degradation process are shown in Fig.4. Only four representative time snapshots are shown: t=0, 9, 18 and 27 (dimensionless time). We calibrated these degradation levels with the experimental measurements presented in [41]. In this work, the authors provide data corresponding to the volume of DQ collagen (which becomes fluorescent while being degraded) per cell within a period of 24h. By computing the volume of each tetrahedron of the finite element mesh of our problem, we obtain the total degraded volume at t=9 (equivalent to approximately 1 day, according to the dimensionless definition of this variable [39]), performing the sum of the degraded volume fractions of the tetrahedra for the corresponding ECM density map. The computed quantity was 8456 𝜇m3 , which lies within the range given in the referred study ( 1000 −10000 𝜇 m3∕cell ). The referred experimental study was established at similar circumstances than the ones presented in our study, although not strictly the same. Therefore, this calibration should be understood as a qualitative estimation to show the applicability of the degradation model, without compromising its validity. Figure5 plots variables v2 , v3 , v4 and v5 of the degradation model at t=27 . Fig. 3 Different steps of the computational implementation of the model in COMSOL Multiphysics: geometry, mesh, MMP-emitting regions (labelled in blue), and definition of prescribed field variable cancer density (c(x,y,z)), which triggers the emission of MMPs in the selected MMP-emitting regions Table 1 Calibrated nondimensional values for the parameters used in the simulations Parameter Dimensionless value 𝛿1 0.02 𝛿2 0.1 𝜇v 0.002 Dv2 129 k21 5 k22 0.195 𝛽v2 0.1 𝛼v3 1.55 k31 27.4 k32 2 𝛽 v 3 0.1 Dv4 129 𝛼 v 4 4 k41 5 k42 27.4 k43 2 k51 27.4 k52 0.195 k53 2 Engineering with Computers Fig. 4 ECM (dimensionless) density maps (XY cross midsection) for 4 different timepoints: t=0 , 9, 18 and 27 Fig. 5 Variables v2 , v3 , v4 and v5 (XY cross midsection) at t = 27 . These variables correspond to the concentrations of MMP-2, MT1-MMP, TIMP2 and the generic intermediate complex, respectively Engineering with Computers 4 Mechanical characterization ofthehydrogel includingdegradation In this section we describe the multiscale approach to obtain the (virtual) macroscopic behavior of fibered hydrogels, as well as the tool developed for producing randomly generated interconnected fiber network microstructures. This methodology is used to generate synthetic shear rheology tests corresponding to different levels of ECM density (variable v1 in the degradation model previously described, corresponding v1=1 and v1=0 to 100% and 0% of the initial ECM density, respectively). The reader is referred to [42] for details of this approach. We only consider heterogeneity of the mechanical properties of the medium due to matrix degradation, and therefore all the synthetically generated fiber microstructures are considered to behave homogeneously, at the initial time of analysis. 4.1 Fibered matrix generator andshear rheology tests The generator creates randomly distributed and interconnected fiber networks and is aimed at replicating the microstructure of fibered materials, such as collagen hydrogels. Fibers are constructed from sequences of nonlinear beam finite elements. The main input parameters of the generator are the fiber volumetric density, length of a single fiber and the curvature of a single fiber [42]. These parameters were calibrated by inspecting Second-Harmonic Generation (SHG) images of real collagen hydrogels obtained in the laboratory. Values representative of the average were retrieved from the segmentation of a set of fibers, ensuring that they remained within the ranges found in the literature. Figure6 shows a comparison between a real collagen hydrogel microstructure and the synthetically generated one. Figure6 (left) shows a SHG image of a real collagen hydrogel, which correspond to a cross section of a 3D sample of the hydrogel. Equivalently, Fig.6 (right) shows a cross section of the reproduced 3D microstructure using our fiber matrix generator. Table2 includes the values employed in the matrix generation along with those that are typically present in the literature. The generated matrices are then subjected to computational (virtual) simple shear tests, by using the referred multiscale approach [42]. 1.2 mg/ml collagen type I hydrogels were prepared following the protocol described in [47]. Briefly, we mixed 400 μ L rat tail collagen (collagen R, 1.2 mg/mL; Matrix Bioscience), 400 μ L bovine skin collagen (collagen G, 4mg/mL; Matrix Bioscience), 191 μ L NaHCO3 (23mg/mL; Supplier), 191 μ L 10x DMEM (Biochrom), 802 μ L distilled water, and 14.37 μ L NaOH (1M) to adjust the pH to 10. The resultant 2mL final collagen solution is pipetted into a 35-mm petri dish and polymerized in an incubator at 37 ◦ C, 95% relative humidity, and 5% CO2 Fig. 6 Real and generated (synthetic) microstructures of 1.2 mg/ml collagen hydrogel Table 2 Calibrated parameters for matrix generation and virtual shear rheology tests Symbol Definition Value Range found in the literature Efiber Elastic modulus of a single fiber 4.3 [MPa] Tens-hundreds of MPa [43, 44] Dfiber Diameter of a single fiber 600 [nm] Hundreds of nm [43, 45] Lfiber Length of a single fiber 10 [ 𝜇 m] Tens of 𝜇 m [46] Engineering with Computers and inverse methods. Note that in Fig.14a, errors remain within the same order of magnitude through the different timepoints, since all these cases reconstruct tractions assuming the real heterogeneous properties of the ECM (although for the degradation patterns associated to the different considered timepoints). On the other hand, Fig.14b shows that errors increase from t=0 (which corresponds to traction reconstruction using homogeneous and virgin ECM properties) when ECM degradation evolves in time and increases. This trend means that traction reconstruction errors increase as ECM properties deviates from the assumed homogeneous and virgin ECM properties for traction reconstruction (conventional TFM). 6 Conclusions In this study, we implemented and solved a complex transient system of partial differential equations, to model metalloproteinase-mediated ECM degradation within a 3D domain involving the geometry of a real cell. This type of degradation is the one that promotes cancer cells to migrate and progress within the extracellular space. As far as degradation induces heterogeneity and softening of the mechanical properties of the ECM, accounting for this phenomenon is critical to accurately compute cellular forces, as for instance, following the approach exposed in this work. As part of our future work, we plan to calibrate and validate the proposed model of degradation with invitro models of cancer cells, which allows to obtain heterogeneous ECM Fig. 13 Ground truth and both forward and inverse conventional (homogeneous) TFM reconstruction results (magnitude of displacements, Frobenius norm of the logarithmic strain tensor, and magnitude of tractions) neglecting ECM degradation ( t=27 ). Insets in the figures represent a detail of the variables localized in the tips of the cell Engineering with Computers Fig. 14 Traction error metric magnitude corresponding to the considered timepoints of the study, for both forward (red boxes) and inverse (blue boxes) formulations t = 0t = 9t = 18 t = 27 10 15 20 25 30 35 40 Traction error [%] DEG-DEG Inverse method Forward method t = 0t = 9t = 18 t = 27 10 20 30 40 50 60 70 Traction error [%] DEG-NODEG Inverse method Forward method Engineering with Computers properties that are needed to reconstruct tractions following TFM techniques. The concept lies in the tuning of the parameters of the model with different real degradation patterns obtained in the laboratory. With a sufficient amount of data and validation, it is assumed that experimental ECM degradation measurements will be greatly reduced, complemented and accurately reproduced by the model. Indeed, experimentally, degradation due to proteinase activity can be visualized and characterized in both space and time by means of live-imaging with DQ-collagen [52]. DQ-colla- gen is initially non-fluorescent, but as it starts to become degraded, it releases fluorescent particles that can be imaged by means of confocal microscopy. This can be used to characterize the rate (speed) of degradation and, furthermore, to tune the degradation model parameters. In addition, only MT1-MMP and MMP-2 -induced degradation was considered in this work. This constitutes a limitation since, in reality, cancer cells are able to secrete a great variety of additional substances that influence the ECM remodeling process (other MMPs, as well as proteases pertaining to other families, such as ADAMs/ADAMTSs, among others). The relation between the degradation model outcomes (i.e. ECM density) and corresponding ECM mechanical properties was established through a multiscale methodology that provides virtual rheology tests for different ECM densities. In this context, we proved the efficacy of our fibered matrix generator tool at the time of providing synthetic rheology data corresponding to generated collagen microstructures that mimic real collagen hydrogels. The shear rheology curves obtained were accurately fitted by the chosen hyperelastic nonlinear model (used to reconstruct macroscopic tractions), which allowed for the retrieval of the relationship between ECM density and material model parameter 𝜅o . Nonetheless, a homogeneous behavior of the hydrogel was considered at the initial time of analysis, as commonly considered in TFM. This limitation can be addressed in future works as an application of our proposed computational framework. In this case, the heterogeneous mechanical behavior of each material macroscopic point of the domain can be obtained through the initial matrix density distribution of fibers, similarly to ECM degradation. For that purpose, an accurate reconstruction of the fibered microstructure of the hydrogel is required for the whole domain of analysis (including the embedded cell), which can be obtained, for instance, from SHG microscopy or fiber staining images. Moreover, at the time of reproducing the collagen hydrogel microstructure, constant geometric parameters (length and curvature of fibers) representative of the average values found in the examined images were considered for all fibers in the synthetic microstructures. In the future application of the methodology to a real invitro experiment, a reconstruction of the exact material microstructure resulting from the segmentation of the 3D images and subsequent meshing will be employed, for a more realistic representation of the assessed cell microenvironment. Regarding the TFM results, there are several interesting conclusions that can be drawn from this study. First, from a simulation point of view, displacements and strains are underestimated at the same time that tractions become higher when neglecting degradation. This is the result of collagen degradation in the vicinity of the boundary of the cell, which induces a local reduction of the stiffness (softening) of the ECM in the location in which the cell exerts forces. Second, with respect to performance of both forward and inverse methodologies at traction reconstruction, the clear outperformance of the inverse method is concluded, as it was previously found in different applications [30, 50]. Finally, it is demonstrated that neglecting ECM degradation (as in conventional TFM) results in a significant overestimation of tractions and non negligible errors in all analyzed cases, but specially for advanced degradation stages. The study of the effect of degradation in traction reconstruction may therefore be more relevant for intermediate and long time scales, and not much for short time scales. However, we believe it is also relevant to study the mechanical behavior of cells under the effect of matrix degradation, to analyze the associated transient dynamics during the whole process and its effect on the mechanical behavior of cells. This suggests that the presented approach may be an available methodology to accurately reconstruct tractions in TFM in the presence of metalloproteinase-mediated ECM degradation. Acknowledgements A.A.-F. and J.A.S.-H. gratefully acknowledge the Financial support from MCIN/AEI/10.13039/501100011033 [PID2021-126051OB-C42]; and project P20-01195 funded by the Consejería de Economía, Conocimiento, Empresas y Universidad de la Junta de Andalucía. J.B-F. gratefully acknowledges the support from Research Foundation--Flanders (FWO, Grant 1259223N). Funding Funding for open access publishing: Universidad de Sevilla/ CBUA. Data availability The data employed in the study, as well as additional information and details about any of the parts that compose the described methodology are available from the corresponding author upon reasonable request. Declarations Conflict of interest All authors declare that they have no Conflict of interest. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will Engineering with Computers need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. References 1. Vogel V, Sheetz M (2006) Local force and geometry sensing regulate cell functions. https:// doi. org/ 10. 1038/ nrm18 90. www. nature. com/ revie ws/ molce llbio 2. Kumar S, Weaver VM (2009) Mechanics, malignancy, and metastasis: the force journey of a tumor cell. Cancer Metastasis Rev 28(1–2):113. https:// doi. org/ 10. 1007/ S10555- 008- 9173-4 3. Ingber DE (2002) Mechanical signaling and the cellular response to extracellular matrix in angiogenesis and cardiovascular physiology . https:// doi. org/ 10. 1161/ 01. RES. 00000 39537. 73816. E5 . http://www.circresaha.org 4. Engler AJ, Sen S, Sweeney HL, Discher DE (2006) Matrix elasticity directs stem cell lineage specification. Cell 126(4):677–689. https:// doi. org/ 10. 1016/J. CELL. 2006. 06. 044 5. Ross R (1986) The pathogenesis of atherosclerosis–an update. New Engl J Med 314(8):488–500. https:// doi. org/ 10. 1056/ NEJM1 98602 20314 0806 6. Vining KH, Mooney DJ (2017) Mechanical forces direct stem cell behaviour in development and regeneration. Nature Publishing Group . https:// doi. org/ 10. 1038/ nrm. 2017. 108 . https://www. nature.com/articles/nrm.2017.108 7. Hynes RO (2009) The extracellular matrix: Not just pretty fibrils. Science 326(5957), 1216–1219 https:// doi. org/ 10. 1126/ scien ce. 11760 09 https://www.science.org/doi/pdf/10.1126/ science.1176009 8. Daley WP, Peters SB, Larsen M (2008) Extracellular matrix dynamics in development and regenerative medicine. J Cell Sci 121(3):255–264. https:// doi. org/ 10. 1242/ jcs. 006064 (https:// journ als. biolo gis ts. com/ jcs/ artic lepdf/ 121/3/ 255/ 15024 98/ 255. pdf) 9. Kleiner DE, Stetler-Stevenson WG (1999) Matrix metalloproteinases and metastasis. Cancer Chemother Pharmacol 43(1):42–51 10. Werb Z (1997) Ecm and cell surface proteolysis: regulating cellular ecology. Cell 91:439–442 11. López-Otín C, Overall CM (2002) Protease degradomics: a new challenge for proteomics. Nat Rev Mol Cell Biol 3:509–519 12. Harris A, Wild P, Stopak D (1980) Silicone rubber substrata: a new wrinkle in the study of cell locomotion . https:// doi. org/ 10. 1126/ scien ce. 69877 36 . http://www.sciencemag.org/cgi/ doi/10.1126/science.6987736 13. Dembo M, Wang YL (1999) Stresses at the cell-to-substrate interface during locomotion of fibroblasts. Biophys J 76(4):2307–2316. https:// doi. org/ 10. 1016/ S0006- 3495(99) 77386-8 14. Butler JP, Tolić-Nørrelykke IM, Fabry B, Fredberg JJ (2002) Traction fields, moments, and strain energy that cells exert on their surroundings. Am J Physiol Cell Physiol 282(3):2 15. Sabass B, Gardel ML, Waterman CM, Schwarz US (2008) High resolution traction force microscopy based on experimental and computational advances. Biophys J 94(1):207–20. https:// doi. org/ 10. 1529/ bioph ysj. 107. 113670 16. Izquierdo-Álvarez A, Vargas DA, Jorge-Peñas Á, Subramani R, Vaeyens MM, Van Oosterwyck H (2019) Spatiotemporal analyses of cellular tractions describe subcellular effect of substrate stiffness and coating. Ann Biomed Eng 47(2):624–637. https:// doi. org/ 10. 1007/ s10439- 018- 02164-2 17. Huang Y, Gompper G, Sabass B (2020) A Bayesian traction force microscopy method with automated denoising in a user-friendly software package. Comput Phys Commun 256:107313 https:// doi. org/ 10. 1016/j. cpc. 2020. 107313arXiv: 2005. 01377 18. Duval K, Grover H, Han LH, Mou Y, Pegoraro AF, Fredberg J, Chen Z (2017) Modeling physiological events in 2D vs. 3D cell culture. Am Physiol Soc. https:// doi. org/ 10. 1152/ physi ol. 00036. 2016 . /pmc/articles/PMC5545611/?report=abstracthttps://www. ncbi.nlm.nih.gov/pmc/articles/PMC5545611/ 19. Caliari SR, Burdick JA (2016) A practical guide to hydrogels for cell culture. Nat Methods 13(5):405–414. https:// doi. org/ 10. 1038/ nmeth. 3839 20. Franck C, Maskarinec SA, Tirrell DA, Ravichandran G, Genin G (2011) Three-dimensional traction force microscopy: a new tool for quantifying cell-matrix interactions. PLoS ONE 6(3):17833. https:// doi. org/ 10. 1371/ journ al. pone. 00178 33 21. Toyjanova J, Bar-Kochba E, López-Fagundo C, Reichner J, Hoffman-Kim D, Franck C (2014) High resolution, large deformation 3D traction force microscopy. PLoS ONE 9(4):1–12. https:// doi. org/ 10. 1371/ journ al. pone. 00909 76 22. Gjorevski N, Piotrowski AS, Varner VD, Nelson CM (2015) Dynamic tensile forces drive collective cell migration through three-dimensional extracellular matrices. Sci Rep 5:11458. https:// doi. org/ 10. 1038/ srep1 1458 23. Gjorevski N, Nelson CM (2012) Mapping of mechanical strains and stresses around quiescent engineered three-dimensional epithelial tissues. Biophys J 103(1):152–162. https:// doi. org/ 10. 1016/j. bpj. 2012. 05. 048 24. Maskarinec SA, Franck C, Tirrell DA, Ravichandran G (2009) Quantifying cellular traction forces in three dimensions. Proc Natl Acad Sci 106(52):22108–22113 https:// doi. org/ 10. 1073/ pnas. 09045 65106 https://www.pnas.org/doi/pdf/10.1073/ pnas.0904565106 25. Yang Z, Lin J-S, Chen J, Wang JH-C (2006) Determining substrate displacement and cell traction fields-a new approach. J Theor Biol 242(3):607–616. https:// doi. org/ 10. 1016/j. jtbi. 2006. 05. 005 26. Tang X, Tofangchi A, Anand SV, Saif TA (2014) A novel cell traction force microscopy to study multi-cellular system. PLoS Comput Biol 10(6):1003631 27. Legant WR, Miller JS, Blakely BL, Cohen DM, Genin GM, Chen CS (2010) Measurement of mechanical tractions exerted by cells in three-dimensional matrices. Nat Methods 7(12):969–971. https:// doi. org/ 10. 1038/ nmeth. 1531 28. Feng X, Hui CY (2016) Force sensing using 3D displacement measurements in linear elastic bodies. Comput Mech 58(1):91– 105. https:// doi. org/ 10. 1007/ s00466- 016- 1283-1 29. Sanz-Herrera JA, Barrasa-Fano J, Cóndor M, Van Oosterwyck H (2021) Inverse method based on 3D nonlinear physically constrained minimisation in the framework of traction force microscopy. Soft Matter 17(45):10210–10222. https:// doi. org/ 10. 1039/ d0sm0 0789g 30. Apolinar-Fernández A, Barrasa-Fano J, Cóndor M, Van Oosterwyck H, Sanz-Herrera JA (2023) Traction force reconstruction assessment on real three-dimensional matrices and cellular morphologies. Int J Eng Sci 186:103828. https:// doi. org/ 10. 1016/j. ijeng sci. 2023. 103828 31. Barrasa-Fano J, Shapeti A, Jorge-Peñas Barzegari M, Sanz-Her- rera JA, Van Oosterwyck H (2021) Tfmlab: A matlab toolbox for 4d traction force microscopy. SoftwareX 15, 100723 https:// doi. org/ 10. 1016/j. softx. 2021. 100723 32. Song D, Hugenberg N, Oberai AA (2019) Three-dimensional traction microscopy with a fiber-based constitutive model. Comput Methods Appl Mech Eng 357:112579. https:// doi. org/ 10. 1016/j. cma. 2019. 112579 33. Bloom RJ, George JP, Celedon A, Sun SX, Wirtz D (2008) Mapping local matrix remodeling induced by a migrating tumor cell using three-dimensional multiple-particle tracking. Biophys J 95(8):4077–4088. https:// doi. org/ 10. 1529/ bioph ysj. 108. 132738 34. Chen S, Xu W, Kim J, Nan H, Zheng Y, Sun B, Jiao Y (2019) Novel inverse finite-element formulation for reconstruction of Engineering with Computers relative local stiffness in heterogeneous extra-cellular matrix and traction forces on active cells. Phys Biol 16(3):036002. https:// doi. org/ 10. 1088/ 1478- 3975/ ab0463 35. Song D, Seidl DT, Oberai AA (2020) Three-dimensional traction microscopy accounting for cell-induced matrix degradation. Comput Methods Appl Mech Eng 364:112935. https:// doi. org/ 10. 1016/j. cma. 2020. 112935 36. Karagiannis ED, Popel AS (2004) 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. J Biol Chem 279:39105–39114 37. Donzé A, Fanchon E, Gattepaille LM, Maler O, Tracqui P (2011) Robustness analysis and behavior discrimination in enzymatic reaction networks. PLoS ONE 6(9):1–16. https:// doi. org/ 10. 1371/ journ al. pone. 00242 46 38. Hoshino D, Koshikawa N, Suzuki T, Quaranta V, Weaver AM, Seiki M, Ichikawa K (2012) Establishment and validation of computational model for mt1-mmp dependent ecm degradation and intervention strategies. PLoS Comput Biol 8(4):1–10. https:// doi. org/ 10. 1371/ journ al. pcbi. 10024 79 39. Deakin N, Chaplain M (2013) Mathematical modeling of cancer invasion: the role of membrane-bound matrix metalloproteinases. Front Oncol. https:// doi. org/ 10. 3389/ fonc. 2013. 00070 40. Deakin NE (2015) Ph.D. Thesis: Mathematical modelling of cancer growth and spread: The role of matrix metalloproteinases. University of Dundee 41. Anguiano M, Morales X, Castilla C, Pena AR, Ederra C, Martínez M, Ariz M, Esparza M, Amaveda H, Mora M, Movilla N, Aznar JMG, Cortés-Domínguez I, Ortiz-de-Solorzano C (2020) The use of mixed collagen-matrigel matrices of increasing complexity recapitulates the biphasic role of cell adhesion in cancer cell migration: Ecm sensing, remodeling and forces at the leading edge of cancer invasion. PLoS ONE 15(1):1–29. https:// doi. org/ 10. 1371/ journ al. pone. 02200 19 42. Sanz-Herrera JA, Apolinar-Fernandez A, Jimenez-Aires A, Perez- Alcantara P, Dominguez J, Reina-Romo E (2024) Multiscale characterization of the mechanics of curved fibered structures with application to biological materials. bioRxiv https:// doi. org/ 10. 1101/ 2024. 01. 09. 574800 https://www.biorxiv.org/content/early/ 2024/01/09/2024.01.09.574800.full.pdf 43. Malandrino A, Trepat X, Kamm RD, Mak M (2019) Dynamic filopodial forces induce accumulation, damage, and plastic remodeling of 3d extracellular matrices. PLoS Comput Biol 15(4):1–26. https:// doi. org/ 10. 1371/ journ al. pcbi. 10066 84 44. Varma S, Orgel JPRO, Schieber JD (2016) Nanomechanics of type i collagen. Biophys J 111(1):50–56. https:// doi. org/ 10. 1016/j. bpj. 2016. 05. 038 45. Yang L, van der Werf KO, Fitié CFC, Bennink ML, Dijkstra PJ, Feijen J (2008) Mechanical properties of native and cross-linked type i collagen fibrils. Biophys J 94(6):2204–2211. https:// doi. org/ 10. 1529/ bioph ysj. 107. 111013 46. Harjanto D, Maffei JS, Zaman MH (2011) Quantitative analysis of the effect of cancer invasiveness and collagen concentration on 3d matrix remodeling. PLoS ONE 6(9):1–8. https:// doi. org/ 10. 1371/ journ al. pone. 00248 91 47. Cóndor M, Steinwachs J, Mark C, García-Aznar JM, Fabry B (2017) Traction force microscopy in 3-dimensional extracellular matrix networks. Current Protocols in Cell Biology 75(1):10– 221102220 https:// doi. org/ 10. 1002/ cpcb. 24 https://currentprotocols.onlinelibrary.wiley.com/doi/pdf/10.1002/cpcb.24 48. Steinwachs J, Metzner C, Skodzek K, Lang N, Thievessen I, Mark C, Münster S, Aifantis KE, Fabry B (2016) Three-dimensional force microscopy of cells in biopolymer networks. Nat Methods 13(2):171–176. https:// doi. org/ 10. 1038/ nmeth. 3685 49. Schwarz US, Soiné JRD (2015) Traction force microscopy on soft elastic substrates: A guide to recent computational advances. Biochim Biophys Acta Mol Cell Res 1853(11):3095–3104. https:// doi. org/ 10. 1016/j. bbamcr. 2015. 05. 028arXiv: 1506. 02394 v1 50. Barrasa-Fano J, Shapeti A, Jong J, Ranga A, Sanz-Herrera JA, Van Oosterwyck H (2021) Advanced in silico validation framework for three-dimensional traction force microscopy and application to an invitro model of sprouting angiogenesis. Acta Biomater 126:326–338. https:// doi. org/ 10. 1016/j. actbio. 2021. 03. 014 51. Vaeyens MM, Jorge-Peñas A, Barrasa-Fano J, Steuwe C, Heck T, Carmeliet P, Roeffaers M, Van Oosterwyck H (2020) Matrix deformations around angiogenic sprouts correlate to sprout dynamics and suggest pulling activity. Angiogenesis 23(3):315– 324. https:// doi. org/ 10. 1007/ s10456- 020- 09708-y 52. Sameni M, Olive M, Moin K, Sloane B (2008) Visualizing protease activity in living cells: from two dimensions to four dimensions. Curr Protoc Cell Biol 4:4–20. https:// doi. org/ 10. 1002/ 04711 43030. cb042 0s39 Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.