Inverse modeling of heterogeneous ECM mechanical properties in nonlinear 3DTFM

Contents lists available at ScienceDirect European Journal of Mechanics / A Solids journal homepage: www.elsevier.com/locate/ejmsol Full length article Inverse modeling of heterogeneous ECM mechanical properties in nonlinear 3DTFM Alejandro Apolinar-Fernándeza,b, Jorge Barrasa-Fanoc, Hans Van Oosterwyckc,d, José A. Sanz-Herreraa,b,∗ aEscuela Técnica Superior de Ingeniería, Universidad de Sevilla, Spain bInstituto de Biomedicina de Sevilla (IBIS), Spain cBiomechanics section, Department of Mechanical Engineering, KU Leuven, Belgium dPrometheus division of Skeletal Tissue Engineering, KU Leuven, Belgium A R T I C L E I N F O Keywords: Traction force microscopy Inverse methods Tikhonov regularization Mechanobiology ECM remodeling Nonlinear continuum mechanics A B S T R A C T Accurate characterization of cellular tractions is crucial for understanding cell-extracellular matrix (ECM) mechanical interactions and their implications in pathology-related situations, yet their direct measurement in experimental setups remains challenging. Traction Force Microscopy (TFM) has emerged as a key methodology to reconstruct traction fields from displacement data obtained via microscopic imaging techniques. While traditional TFM methods assume homogeneous and static ECM properties, the dynamic nature of the ECM through processes such as enzyme–induced collagen degradation or cell-mediated collagen deposition i.e. ECM remodeling, requires approaches that account for spatio-temporal evolution of ECM stiffness heterogeneity and other mechanical properties. In this context, we present a novel inverse methodology for 3DTFM, capable of reconstructing spatially heterogeneous distributions of the ECM’s stiffness. Our approach formulates the problem as a PDE-constrained inverse method which searches for both displacement and the stiffness map featured in the selected constitutive law. The elaborated numerical algorithm is integrated then into an iterative Newton–Raphson/Finite Element Method (NR/FEM) framework, bypassing the need for external iterative solvers. We validate our methodology using in silico 3DTFM cases based on real cell geometries, modeled within a nonlinear hyperelastic framework suitable for collagen hydrogels. The performance of our approach is evaluated across different noise levels, and compared versus the commonly used iterative L-BFGS algorithm. Besides the novelty of our formulation, we demonstrate the efficacy of our approach both in terms of accuracy and CPU time efficiency. 1. Introduction Mechanobiology is an active research field which investigates the relationship between mechanics and the origin and progression of multiple physiological and pathological processes. Such processes are linked to cell activity within biological tissues, and fundamental insight on the principles by which mechanical forces and conditions alter the activity of cells and their microenvironment, are of great interest to the fields of biology and medicine (Ingber, 2003; Vogel and Sheetz, 2006; Mammoto et al., 2013). Indeed, it has been found that the mechanical interaction between cells and their surrounding extracellular matrix (ECM) guides processes such as bending, stretching, and repositioning of the epithelium, which are required for tissue morphogenesis (Martin ∗Corresponding author at: Camino de los descubrimientos s/n, 41092 Seville, Spain. E-mail addresses: [email protected] (A. Apolinar-Fernández), [email protected] (J. Barrasa-Fano), [email protected] (H. Van Oosterwyck), [email protected] (J.A. Sanz-Herrera). et al., 2009; Rauzi et al., 2008). These processes also activate signaling pathways that regulate cancer invasion (Broguiere et al., 2018), angiogenesis (Ingber, 2002; Vaeyens et al., 2020), wound healing (Li et al., 2021), survival, proliferation and stem-cell differentiation (Iskratsch et al., 2014; Mui et al., 2016; Engler et al., 2006; Nelson et al., 2005; Siedlik et al., 2016; Humphrey et al., 2014), and influence other tissue or organ pathologies (Ross, 1986; Chicurel et al., 1998; Riley et al., 2002; Lammerding et al., 2004). The mechanism through which cells translate mechanical stimuli into different cellular activities is known as mechanotransduction. The intricate nature governing this phenomenon is not yet fully known, and researchers aim to unravel it by means of in vitro models that emulate the composition, organization and mechanical characteristics of human-like tissue, which may potentially boost https://doi.org/10.1016/j.euromechsol.2025.105722 Received 14 March 2025; Received in revised form 7 May 2025; Accepted 14 May 2025 European Journal of Mechanics / A Solids 114 (2025) 105722 Available online 31 May 2025 0997-7538/© 2025 The Author(s). Published by Elsevier Masson SAS. This is an open access article under the CC BY-NC-ND license ( http://creativecommons.org/licenses/by-nc-nd/4.0/ ). A. Apolinar-Fernández et al. advances in tissue engineering and therapeutic techniques (Vining and Mooney, 2017; Janmey et al., 2020; Kim et al., 2021). In this context, precise characterization of cellular tractions is essential in the study of mechanotransduction. However, the direct measurement of such tractions and forces in experimental setups is challenging, and several methodologies have been devised and elaborated in the last years (Polacheck and Chen, 2016; Wang et al., 2007). Within these, Traction Force Microscopy (TFM) is emphasized by its effectiveness and potential (Schwarz and Soiné, 2015; Style et al., 2014; Mulligan et al., 2018; Hall et al., 2013). TFM is a technique that focuses in recovering traction fields associated to specific cell-ECM interactions from measured displacement data obtained via microscopy imaging techniques in the laboratory. The general 3DTFM workflow is summarized as follows: cells are seeded and embedded in a hydrogel that mimics the ECM, which contains spatially distributed fluorescent particles (beads). Once cell activity is initiated, a set of microscopy images are captured from which cell morphology and beads motion can be retrieved. In particular, two sets of images are at least needed, one referred to the cell stressed state (reference configuration), and the second corresponding to the cell relaxed state attained after lysis. The force-induced matrix displacement field is then reconstructed by tracking beads motion between the two recorded configurations. Alternatively, previous works have used bead-free approaches that simplify the gel preparation procedure by imaging collagen fibers using second harmonic generation (Jorge-Peñas et al., 2017), confocal reflectance (Laforgue et al., 2022) or by fluorescently labeling the collagen fibers (Doyle et al., 2021). Finally, provided the measured displacement field and the constitutive model of the matrix, tractions exerted at the cell-ECM interface are retrieved using the elasticity equations. 3DTFM undoubtedly constitutes a more sophisticated methodology in comparison to its 2D counterpart, providing traction reconstructions in reliable in vitro settings that better reproduce the in vivo cell environment (Legant et al., 2010; Caliari and Burdick, 2016). 3DTFM presents, however, an increased complexity compared to 2DTFM, regarding experimental setups (Cóndor et al., 2017; Shapeti et al., 2024; Jorge-Peñas et al., 2017; Colin-York et al., 2016, 2017; Mulligan et al., 2017, 2019; Shapeti et al., 2024; BarrasaFano et al., 2021; Steinwachs et al., 2016; Broguiere et al., 2018; Vaeyens et al., 2020), imaging algorithms (Izquierdo-Álvarez et al., 2019; Huang et al., 2019; Style et al., 2014; Janke et al., 2020; Feng et al., 2014; Bar-Kochba et al., 2015; Jorge-Peñas et al., 2017, 2015), matrix mechanical characterization (Steinwachs et al., 2016; Jansen et al., 2018; Izquierdo-Álvarez et al., 2019; Javanmardi et al., 2021; Kraning-Rush et al., 2012; Toyjanova et al.; Mulligan et al., 2018), and computational procedures for traction reconstruction (BarrasaFano et al., 2021; Legant et al., 2010; Apolinar-Fernández et al., 2023, 2024a,b; Sanz-Herrera et al., 2021; Song et al., 2020a,b). Concerning the computational side of TFM, inverse methods are the most popular methodology for traction reconstruction from measured noisy displacement data. These methods state an optimization problem in which a new displacement field is searched for under the condition that it has to resemble the measured displacements, while being subjected to the fulfillment of fundamental mechanical principles (Apolinar-Fernández et al., 2024b). Inverse methods have been developed over the years and successfully implemented in multiple scenarios. In 2DTFM experiments in which linear elastic behavior is assumed for the matrix, analytical solutions can be obtained by means of Green’s functions (Schwarz and Soiné, 2015). Regularization is often applied in order to reduce the impact of noise on the results by penalizing high norm values of tractions at certain regions (Feng and Hui, 2016). Regularized inverse approaches have been also employed in 2D/3D linear cases (Makarchuk et al., 2018; Suñé-Auñón et al., 2017; Legant et al., 2010; Du et al., 2018), linear 2D/3DTFM with constraints (Ambrosi, 2006; Vitale et al., 2012a,b; Peschetola et al., 2013) (Mulligan et al., 2019), and nonlinear 2D/3D with or without constraints (Michel et al., 2013; Steinwachs et al., 2016; Dong and Oberai, 2017; Cóndor et al., 2019). In addition, constrained nonregularized methods have been proposed for 3D linear and nonlinear cases (Barrasa-Fano et al., 2021; Sanz-Herrera et al., 2021). An alternative methodology, the forward method, has also been extensively applied in TFM due to the simplicity of its implementation and its computational efficiency relative to inverse methods (Mulligan et al., 2018; Hall et al., 2013, 2012; Franck et al., 2011; Toyjanova et al., 2014). In this approach, the measured noisy displacement field is directly interpolated and transformed into a continuous field, which is differentiated to obtain the associated strain field. The latter is introduced in the constitutive equation selected for the matrix to compute the corresponding stresses and tractions. Nonetheless, the accuracy of the forward method highly relies on the quality of the input strains, specially those close to the cell surface, which is often compromised as spatial derivatives tend to amplify the noise contained in the measured displacements. Indeed, it has been shown previously that inverse methods provide substantially more accurate reconstructions than their forward counterpart (Barrasa-Fano et al., 2021; Sanz-Herrera et al., 2021; Apolinar-Fernández et al., 2023, 2024a). An important limitation of the referred works is that these TFM studies assumed the mechanical properties of the ECM as spatially homogeneous and static without considering remodeling. Actually, the ECM is composed of a complex interconnected network of biopolymers, providing structural support for cell activity, that serves as a platform for the diffusion of bio-chemicals within tissues. However, the ECM possesses a highly dynamic nature and it undergoes a continuous remodeling process that is crucial for the regulation of diverse cellular behaviors (Daley et al., 2008). In particular, cells are able to express enzymes, such as the proteolytic matrix metalloproteinases (MMPs), that can locally degrade the collagen microstructure of the ECM (Kleiner and Stetler-Stevenson, 1999; Keating et al., 2017; Trappmann et al., 2017; Jerrell and Parekh, 2014; Khetan et al., 2013; Vincent and Engler, 2013). MMP-induced degradation facilitates cancer growth and spread by cancer cell migration and invasion of surrounding tissues (Werb, 1997; López-Otín and Overall, 2002). In addition, cancer cells deposit significant amounts of collagen in their microenvironment, resulting in localized stiffened regions of the matrix (Zhao et al., 2024) that act as walls that shield them from drugs that are used during treatments (Calvo et al., 2013; Liu et al., 2019). As a result, the assumption of homogeneity of the ECM properties is not justified in multiple scenarios, therefore compromising the accuracy of traction reconstruction when heterogeneity is disregarded. In fact, the relevance of considering ECM heterogeneity in measuring stresses has been emphasized in several experimental studies (Gjorevski and Nelson, 2012; Bloom et al., 2008). On the other hand, reconstruction of heterogeneous ECM patterns may have also a great interest in a number of different biophysical problems (Zhang et al., 2023). Nonetheless, direct measure of ECM mechanical properties constitutes a difficult challenge, mainly due to cell-induced alterations concentrated in the proximity of cells’ surface (Keating et al., 2017). Motivated by this, some authors have proposed innovative inverse formulations to recover material stiffness distributions in 3DTFM experiments. Chen et al. (2019) implemented a mixed forward/inverse methodology for the reconstruction of ECM stiffness profiles from a measured displacement field. Although interesting, its applicability is limited as it was formulated within the context of linear elasticity. Also, the authors showed stiffness reconstructions directly from measured displacements and, like the majority of forward methods, the associated errors may potentially be non negligible. Song et al. (2020a) proposed a nonlinear inverse formulation for reconstructing ECM stiffness profiles. However, the region in which heterogeneous stiffness values were assumed was small and located in the vicinity of the cell surface. Moreover, the solution approach relies on the use of an external iterative solver, particularly the L-BFGS algorithm, which is potentially inefficient as will be demonstrated in the present study. In addition, they considered a Neo-Hookean model for the ECM, which European Journal of Mechanics / A Solids 114 (2025) 105722 2 A. Apolinar-Fernández et al. takes into account geometric nonlinearities due to large deformations, but it does not properly characterize the specific nonlinear behavior of the fibered materials usually employed in 3DTFM. Khang et al. (2023) implemented a methodology conceptually similar to that of Song et al. (2020a), but in this study the region in which heterogeneity can be reconstructed is not limited. The authors modeled the ECM as a NeoHookean material as well, and employed the L-BFGS algorithm to solve their formulated inverse problem. Besides, the authors did not verify the accuracy of reconstruction of the proposed methodology with in silico cases in which noise in the displacement data is added, as they only assumed noiseless ground truth cases. On the other hand, the study carried out by Mei et al. (2023), methodologically similar to Song et al. (2020a), Khang et al. (2023), outlined a computational methodology for the 2D inverse reconstruction of the shear modulus of cell nuclei (nuclear elastography) by means of a specific version of the adjoint method. They presented an in silico TFM study in which they assess on the accuracy of reconstruction of the methodology for different levels of noise added to the ground truth displacements. As in the referred works, the authors employed the iterative solver L-BFGS. In a previous study, we presented a novel 3DTFM formulation which accounts for MMP-induced ECM degradation (Apolinar-Fernández et al., 2024a). This work establishes a multiscale approach in which 3D degradation profiles are obtained from the solution of a system of partial differential equations. The obtained degradation values are then used to calibrate the parameters of a fiber microstructure generator (SanzHerrera et al., 2024), that fits a continuum (macroscale) nonlinear hyperelastic model in order to consider heterogeneous properties in the 3DTFM workflow. The applicability of this approach is however limited considering that empirical data of ECM degradation are needed to calibrate model parameters of the degradation model. In this framework, we present in this paper a novel methodology to reconstruct spatially heterogeneous distributions of the mechanical parameters of the ECM in a 3DTFM context. The problem is formulated as a PDE-constrained inverse method in which new displacements and the stiffness heterogeneity, are retrieved from the input of the standard measured displacements employed in TFM. The solution to the global constrained problem is approached by iteratively solving two coupled subproblems in which the equilibrium constraint is substituted by a Tikhonov regularization term acting on the hydrogel forces (‘‘weak’’ enforcement of the equilibrium). Differently to other similar inverse approaches, in which iterative methods are used to minimize the cost function, our methodology is developed in the context of a coupled iterative Newton–Raphson/Finite Element Method framework (in which nonlinear elasticity is considered), that does not rely on external solvers. In order to test the methodology, we generated in silico 3DTFM cases from a ground truth reference solution, using the geometry of a real cell, and the mechanical behavior of the ECM modeled via a nonlinear hyperelastic law devised for collagen hydrogels. The performance of the proposed methodology is evaluated versus the iterative solver L-BFGS method. Also, we show 3DTFM reconstructions and results corresponding to three different cases referred to the level of noise added to the ground truth displacements solution. The paper is structured in the following sections: Section 2 presents a detailed description of the mathematical formulation and computational implementation of the proposed inverse methodology. Section 3 illustrates the process of generating the ground truth reference case, the definition of the noise added to the ground truth displacements, and the error metrics used in assessing the accuracy of reconstruction. Section 4 shows and discusses the main results obtained from the study. Finally, the paper ends with some concluding remarks featured in Section 5. 2. Mathematical formulation In this section, the mathematical and computational aspects of the method are presented in detail. The aim of the formulation is to reconstruct the heterogeneous stiffness profile 𝜅𝑜(𝑥, 𝑦, 𝑧) by means of a PDE-constrained inverse method, in the context of a FE framework in which the equilibrium condition is imposed. The parameter 𝜅𝑜 is the one associated to the stiffness of the material in the nonlinear hyperelastic law selected to reproduce the ECM’s mechanical behavior (Section 3.2). The main inverse problem is divided into two coupled subproblems that are solved iteratively similarly to a fixed point procedure. Both subproblems contain a Tikhonov regularization term that serves as a weak form of imposing the equilibrium condition. The two subproblems are integrated and solved within a Newton– Raphson/Finite Element Method (NR/FEM) formulation. We first start by providing details about the nonlinear hyperelastic model selected to represent the mechanical behavior of the hydrogel. 2.1. Mechanical characterization of the hydrogel The mechanical characterization of the ECM-mimicking collagen hydrogel is carried out by fitting the model parameters to a specific nonlinear hyperelastic model in shear rheology experiments. The selected nonlinear model was used in previous publications by the authors (Apolinar-Fernández et al., 2023, 2024a,b), and it was originally presented by Steinwachs et al. in Steinwachs et al. (2016), where it was proven to be an excellent model to reproduce the mechanical behavior of collagen hydrogels. In its application in Steinwachs et al. (2016), Böhringer et al. (2024) to similar hydrogels to the one we have considered in this study, the authors verified that no plastic nor viscoelastic behavior was present. Indeed, regarding the latter, in TFM it is common to assume that the ECM relaxation process after cell lysis is quasi-static (Mulligan et al., 2018; Hall et al., 2013), having been reported to typically take a time within the range of 30–120 min (Legant et al., 2010; Steinwachs et al., 2016; Song et al., 2020b; Böhringer et al., 2024). On the other hand, the relaxation time characteristic of hydrogels similar to the one considered here is of around 3 s (Steinwachs et al., 2016). This discrepancy in orders of magnitude between the time of viscoelastic response of the material and that of the drug-induced ECM relaxation reinforces the quasi-static assumption. Moreover, we assume that degradation only takes place during cell activity, and that matrix properties do not change during the ECM unloading process. As a consequence, under these conditions, a hyperelastic model is considered to be valid for fully characterizing the ECM’s behavior in the context of this study. Specifically, the model considers the mechanical behavior of a single fiber, that can be subjected to three different states of deformation: buckling (𝜆 < 0), in which the fiber is unstable and its stiffness decays exponentially; straightening (0≤𝜆 < 𝜆𝑠), in which the fiber recovers its natural length (𝜆𝑠) as it straightens, being the stiffness kept as constant; and stretching (𝜆≥𝜆𝑠), in which the fiber stretches further than its natural length and its stiffness increases exponentially. 𝜆 denotes the macroscopic stretch of a fiber. Accordingly, the strain energy density function associated to a single fiber is defined from the stiffness through its second derivative with respect to 𝜆, 𝛹′′(𝜆) = 𝜅𝑜⎧ ⎪ ⎨ ⎪ ⎩ exp [𝜆∕𝑑𝑜] ∀𝜆 < 0 1 ∀0 ≤𝜆 < 𝜆𝑠 exp [[𝜆−𝜆𝑠]∕𝑑𝑠]∀𝜆≥𝜆𝑠 𝜆=‖𝐅⋅𝐞𝛺‖2− 1 (1) where 𝜅𝑜 is the material stiffness, 𝑑𝑜 is the non-dimensional buckling parameter, 𝑑𝑠 the non-dimensional strain-stiffening (dimensionless) parameter, and 𝜆𝑠 the non-dimensional parameter that defines the extension of the interval within which the fiber behaves linearly, i.e. with constant stiffness. Moreover, 𝐞𝛺 is the unit vector that represents the fiber’s orientation and 𝐅 is the deformation gradient tensor. The macroscopic behavior of the hydrogel is obtained by averaging fiber contributions to the strain energy in all directions. Thus, the constitutive law in mixed configuration is obtained as, 𝐏=⟨𝜕𝛹 (𝜆(𝐅)) 𝜕𝐅⟩𝜴 =⟨𝜕𝛹 𝜕𝜆 𝜕𝜆 𝜕𝐅⟩𝜴 =⟨𝛹′(𝜆)[𝐅𝐞𝛺]⊗𝐞𝛺 ‖𝐅𝐞𝛺‖2⟩𝜴 (2) European Journal of Mechanics / A Solids 114 (2025) 105722 3 A. Apolinar-Fernández et al. Fig. 1. Hydrogel characterization: Synthetic curve (solid line) reproducing the shear rheology behavior (markers) of a collagen hydrogel with 1.2 mg/ml density after fitting the associated parameters to the experimental data. Table 1 Virgin (non-degraded, homogeneous ECM) values of the nonlinear model parameters (Eq. (1)). Parameter Value 𝜅𝑜[Pa] 217 𝑑𝑜[–] 0.375 𝜆𝑠[–] 0.0426 𝑑𝑠[–] 0.0387 where 𝐏 denotes the First Piola–Kirchhoff stress tensor and the brackets indicate the average over all directions 𝛺 in the unit sphere. The model is integrated into the NR/FEM computational framework described in Sections 2.2–2.6 by means of an ABAQUS UMAT subroutine. The reader is referred to Appendix A of this paper, where a thorough description of the model formulation and its computational implementation is presented. For this study, the parameters in Eq. (1) were fitted to reproduce the shear rheology curves obtained from real collagen hydrogels. Fig. 1 shows the shear rheology curve of 1.2 mg/ml collagen hydrogel, together with the curve reproduced with the nonlinear model. The values of the parameters resulting from the fitting process are given in Table 1. These values correspond to the virgin state (non-degraded, homogeneous) of the ECM. These rheological data are identical to those utilized in a previous paper by the authors (Apolinar-Fernández et al., 2024b). The reader is referred to Appendix B of this paper for more details about the curve fitting process. 2.2. Global problem formulation The global problem is written to find both the new displacement field 𝐮 and the heterogeneous stiffness distribution 𝜅𝑜 so that the former resemble the measured displacement field 𝐮∗ and the latter behaves smoothly, provided that the solution to the problem satisfies the equilibrium of forces (nil internal body forces within the hydrogel domain for a self-balanced cell), min 𝐮,𝜅𝑜(1 2‖‖𝐮−𝐮∗‖‖2 2+𝛼 2∫𝛺||∇𝜅𝑜||2𝑑𝑣) s.t. 𝜣=𝟎 (3) where 𝛼 is the regularization parameter, 𝛺 represents the hydrogel domain and 𝛩 the equilibrium constraint manifold in which the candidate solution must lie. Notations ‖⋅‖2 2 and ∫|⋅|2𝑑 indicate the squared L2-norm. This problem has a unique solution, given that the material parameter to be retrieved from its solution, 𝜅𝑜 (Eq. (1)), is involved as a linear scaling factor acting on a nominal virgin value of the material stiffness (Table 1) within the chosen constitutive formulation. This means that, for a given displacement field 𝐮, an explicit expression for 𝜅𝑜 can be extracted from the equilibrium equations. Of all the possible regularized displacement fields 𝐮 that satisfy the imposed boundary conditions and constraints, for a given value of the regularization parameter 𝛼, only one of them is the closest to the measured displacement field 𝐮∗ that minimizes the objective function. Thus, the solution to the problem in Eq. (3) corresponds to a unique pair (𝐮, 𝜅𝑜). This is not guaranteed if more than one material parameter are retrieved simultaneously, or if the dependency of the constitutive law with the parameter in question is more complex (anisotropy, inelastic effects). Similar discussions can be found in Chen et al. (2019), Song et al. (2020a), Khang et al. (2023). Normally, in order to solve this problem, the equilibrium constraint would be introduced into the cost function through a Lagrange multiplier field 𝜂 in the following way, min 𝐮,𝜅𝑜(1 2‖‖𝐮−𝐮∗‖‖2 2+𝜼⋅𝜣+𝛼 2∫𝛺||∇𝜅𝑜||2𝑑𝑣)(4) The methodology presented in this paper is based on the substitution of the ‘‘strong’’ constraint, represented by 𝜼⋅𝜣, by a term that regularizes and minimizes the nodal forces within the hydrogel domain. We refer to this as a ‘‘weak’’ imposition of the equilibrium condition, and it is motivated to keep the numerical workflow of the algorithm, as will be seen in the next sections. The global problem (4) is subdivided into two additional coupled problems that are solved simultaneously: the forward and inverse subproblems, associated to the terms 𝛼 2∫𝛺||∇𝜅𝑜||2𝑑𝑣 and 1 2‖𝐮−𝐮∗‖2 2, respectively; both approximating the condition represented by 𝜼⋅𝜣 with the aforementioned regularization term. The details about the implementation of both subproblems are described thoroughly in the next sections. 2.3. Forward subproblem The forward subproblem is formulated for a given displacement field 𝐮 as input data. As in a conventional forward method, this displacement field is fixed during the solution of this subproblem. At the first global iteration, which involves the forward and inverse subproblems, its value is fixed to the measured displacements field 𝐮∗. The problem (4) is then stated as, min 𝜅𝑜(1 2∫𝛺||∇𝜅𝑜||2𝑑𝑣 + 𝛼𝜅𝑜 2∫𝛺|𝐛|2𝑑𝑣 +𝑔(𝜅𝑜))(5) in which 𝛼𝜅𝑜 is a regularization parameter associated to the forward subproblem, and 𝐛 represents the body forces acting in the ECM’s interior domain. Overall, this regularization term establishes the equilibrium condition through 𝛼𝜅𝑜 as 𝐛 must vanish within the hydrogel domain. 𝑔(𝜅𝑜) is a sigmoid (activation) function included to avoid negative values of 𝜅𝑜 during the iterative process. The problem stated in Eq. (5) has a minimum stationary solution at, 𝜇(𝜅𝑜) + 𝛼𝜅𝑜𝜔(𝜅𝑜) + 𝜐(𝜅𝑜) = 0 (6) where, 𝜇(𝜅𝑜) = 1 2 𝜕 𝜕𝜅𝑜∫𝛺||∇𝜅𝑜||2𝑑𝑣 𝜔(𝜅𝑜) = 1 2 𝜕 𝜕𝜅𝑜[∫𝛺|𝐛|2𝑑𝑣] 𝜐(𝜅𝑜) = 𝜕𝑔(𝜅𝑜) 𝜕𝜅𝑜 (7) European Journal of Mechanics / A Solids 114 (2025) 105722 4 A. Apolinar-Fernández et al. Eqs. (7) are developed into their discretized form within a Finite Element (FE) framework. In this sense, field variables are discretized as 𝜅𝑜≈𝐍𝑒⋅𝜿𝑖 𝒐, 𝛥𝜅𝑜≈𝐍𝑒⋅𝛥𝜿𝑖 𝒐, 𝐛≈𝐍𝑒⋅𝐛𝑖 and 𝑔≈𝐍𝑒⋅𝒈𝑖, where 𝐍𝑒 is the matrix that contains the shape functions of the finite element (𝑒). Further, 𝜿𝑖 𝒐, 𝛥𝜿𝑖 𝒐, 𝐛𝑖 and 𝒈𝑖 are discrete vectors which contain the discrete values of 𝜅𝑜, 𝛥𝜅𝑜, 𝐛 and 𝑔 at node positions 𝑖, respectively, in the FE element (𝑒). Introducing this discretization into Eqs. (7) yields, 1 2∫𝛺||∇𝜅𝑜||2𝑑𝑣 ⟶1 2‖‖𝐁⋅𝜿𝒐‖‖2 2 𝛼𝜅𝑜 2∫𝛺|𝐛|2𝑑𝑣 ⟶ 𝛼𝜅𝑜 2‖‖𝐅𝑏‖‖2 2 (8) 𝐁 is the global gradient matrix after assembly, and 𝜿𝒐 the global vector with discrete values of 𝜅𝑜 at each node of the FE mesh. 𝐅𝑏 is the global vector, after assembly, with discrete values corresponding to the body forces in the internal ECM domain. Functions 𝜇, 𝜔 and 𝜐 in Eq. (7) now take the discretized expressions, 𝝁(𝜿𝒐) = 1 2 𝜕 𝜕𝜿𝒐‖‖𝐁⋅𝜿𝒐‖‖2 2=𝐁⊺⋅𝐁⋅𝜿𝒐(9) 𝝎(𝜿𝒐) = 1 2 𝜕 𝜕𝜿𝒐‖‖𝐅𝑏‖‖2 2=[𝜕𝐅𝑏 𝜕𝜿𝒐]⊺ ⋅𝐅𝑏(10) 𝝊(𝜿𝒐) = 𝜕𝑔 𝜕𝜿𝒐 =𝒈′(𝜿𝒐)(11) with 𝒈′ a nodal global vector which contains the discrete values of the derivative of function 𝑔 evaluated at 𝜿𝒐. These equations are linearized as follows: 𝝁(𝜿𝒐+𝛥𝜿𝒐) ≅ [𝝁]𝜿𝒐+[𝜕𝝁 𝜕𝜿𝒐]𝜿𝒐 ⋅𝛥𝜿𝒐=𝐁⊺⋅𝐁⋅𝜿𝒐+𝐁⊺⋅𝐁⋅𝛥𝜿𝒐(12) 𝝎(𝜿𝒐+𝛥𝜿𝒐) ≅ [𝝎]𝜿𝒐+[𝜕𝝎 𝜕𝜿𝒐]𝜿𝒐 ⋅𝛥𝜿𝒐≅[[𝜕𝐅𝑏 𝜕𝜿𝒐]⊺ ⋅𝐅𝑏]𝜿𝒐 + +[[𝜕𝐅𝑏 𝜕𝜿𝒐]⊺ ⋅ 𝜕𝐅𝑏 𝜕𝜿𝒐]𝜿𝒐 ⋅𝛥𝜿𝒐(13) 𝝊(𝜿𝒐+𝛥𝜿𝒐) ≅ [𝝊]𝜿𝒐+[𝜕𝝊 𝜕𝜿𝒐]𝜿𝒐 ⋅𝛥𝜿𝒐=[𝒈′]𝜿𝒐+[𝒈′′]𝜿𝒐⋅𝛥𝜿𝒐(14) with 𝒈′′ a nodal global vector which contains the discrete values of the second derivative of function 𝑔 evaluated at 𝜿𝒐. Note that the second derivative of 𝐅𝑏 with respect to 𝜿𝒐 was neglected in Eq. (13), which does not impact the accuracy of converged solutions, only affecting the convergence rate of the iterative scheme. On the other hand, using the spatial version of the discretized principle of virtual work (PVW), the following relation can be written 𝐅𝑏=𝐃(𝐮)⋅𝜿𝒐(15) Global matrix 𝐃 is obtained after assembly of the PVW which depends on the displacements, which are known and fixed in this forward subproblem. Note that 𝐅𝑏 is linearly related with 𝜿𝒐 since stresses are linear with 𝜅𝑜 as is seen in the selected hyperelastic law (see Section 2.1). Using Eq. (15) in Eqs. (12), (13) and (14); Eq. (6) yields, [𝐁⊺⋅𝐁+𝛼𝜅𝑜 𝐃⊺⋅𝐃+𝒈′′ 𝑘]⋅𝛥𝜿𝒐= = − [𝐁⊺⋅𝐁⋅𝜿𝑜,𝑘 +𝛼𝜅𝑜 𝐃⊺⋅𝐃⋅𝜿𝑜,𝑘 +𝒈′ 𝑘](16) Remarkably, nonlinearity is introduced in Eq. (16) through the activation function 𝑔. Therefore, this equation is iteratively solved until convergence is achieved. Note that vector 𝜿𝑜 is defined in the nodes of the mesh according to the developed FE formulation. As a result, it is mathematically impossible to satisfy the strong form of the equilibrium condition as it would involve the fulfillment of 3 equations per node with just one unknown variable (the stiffness) per node. This overdetermined system of PDEs may not have a solution. Nonetheless, we know there is at least one distribution of 𝜅𝑜 that fulfills the equations, which is the one utilized for the generation of the ground truth scenarios, and that is thus a solution to the problem. For a similar discussion, the reader is referred to Song et al. (2020a). Consequently, this condition is imposed in our formulation through parameter 𝛼𝜅𝑜 in Eq. (16). The solution obtained from this problem is corrected at each iteration by the updated displacements obtained in the inverse subproblem and viceversa, until both solutions converge into the one corresponding to the global problem. 2.4. Inverse subproblem The inverse subproblem is formulated under the assumption that 𝜅𝑜 is known and fixed, as it results from the forward subproblem. Hence, the problem is stated as a normal unconstrained inverse problem with regularization in which displacements are recomputed and the material properties are known (although inhomogeneous). The problem can be written as follows, min 𝐮(1 2‖‖𝐮−𝐮∗‖‖2 2+𝛼𝑢 2∫𝛺|𝐛|2𝑑𝑣)(17) in which 𝛼𝑢 is the regularization term parameter associated to the inverse subproblem, and 𝐛 is defined analogously to the forward subproblem. A very similar formulation to the inverse subproblem was previously described by the authors in Appendix A of Apolinar-Fernández et al. (2024b). Similarly to the previous section, Eq. (17) has a minimum stationary solution using the Gateaux derivative with respect to 𝐮, 𝐮+𝛼𝑢𝝓(𝐮) = 𝐮∗(18) where 𝝓(𝐮) is defined as a nonlinear function that depends implicitly on 𝐮 through 𝐛, 𝝓(𝐮) = 1 2 𝜕 𝜕𝐮[∫𝛺|𝐛|2𝑑𝑣](19) 𝝓 is linearized in the vicinity of 𝐮 in the following fashion, 𝝓(𝐮+𝛥𝐮) ≅ 𝝓(𝐮) + 𝐷𝝓(𝐮)[𝛥𝐮] = [𝝓]𝐮+[𝜕𝝓 𝜕𝐮]𝐮 ⋅𝛥𝐮(20) being 𝐷𝝓(𝐮)[𝛥𝐮] the directional derivative of 𝝓(𝐮) in the direction of an arbitrary small increment 𝛥𝐮. The linearized version of Eq. (18) is then written as, 𝐮+𝛼𝑢[𝝓(𝐮) + 𝐷𝝓(𝐮)[𝛥𝐮]]=𝐮∗(21) which will be discretized following a Finite Element approach and solved iteratively using a Newton–Raphson procedure. As in the forward subproblem, the field variables are discretized such that 𝐮≈𝐍𝑒⋅𝐮𝑖, 𝛥𝐮≈𝐍𝑒⋅𝛥𝐮𝑖 and 𝐮∗≈𝐍𝑒⋅𝐮∗,𝑖. Besides, 𝐮𝑖, 𝛥𝐮𝑖 and 𝐮∗,𝑖 are vectors whose components correspond to the discrete values of 𝐮, 𝛥𝐮 and 𝐮∗ at node positions 𝑖, respectively, in the FE element (𝑒). Similarly to what was described in Eq. (8), 𝛼𝑢 2∫𝛺|𝐛|2𝑑𝑣 ⟶ 𝛼𝑢 2‖‖𝐅𝑏‖‖2 2(22) and 𝝓(𝐮) is rewritten as, 𝝓(𝐮) = 1 2 𝜕 𝜕𝐮‖‖𝐅𝑏‖‖2 2=[𝜕𝐅𝑏 𝜕𝐮]⊺ ⋅𝐅𝑏(23) which is linearized in the proximity of 𝐮 in the following fashion, 𝝓(𝐮+𝛥𝐮) ≅ 𝝓(𝐮) + 𝐷𝝓(𝐮)[𝛥𝐮] = [[𝜕𝐅𝑏 𝜕𝐮]⊺ ⋅𝐅𝑏]𝐮 + +[[𝜕𝐅𝑏 𝜕𝐮]⊺ ⋅ 𝜕𝐅𝑏 𝜕𝐮]𝐮 ⋅𝛥𝐮(24) European Journal of Mechanics / A Solids 114 (2025) 105722 5 A. Apolinar-Fernández et al. Fig. 2. Workflow of the iterative process: Convergence of the forward subproblem, provided a fixed displacement vector 𝐮, yields the discrete stiffness vector 𝜿𝒐, which is the input quantity of the inverse subproblem. Convergence of the inverse subproblem, provided a fixed 𝜿𝒐, yields the discrete displacement vector 𝐮, which is the input quantity of the forward subproblem at the next global iteration. The global iterative procedure is stopped when the global convergence is achieved. in which the second derivative of 𝐅𝑏 with respect to 𝐮 was neglected as it was done in Eq. (13). Finally, after assembly, Eq. (24) is introduced in Eq. (21) and the resulting linearized global system yields the expression, [𝟏+𝛼𝑢⋅𝐊⊺ 𝑙⋅𝐊𝑙]⋅𝐮𝑙+1 =𝐮∗+𝛼𝑢⋅𝐊⊺ 𝑙⋅[𝐊𝑙⋅𝐮𝑙−𝐅𝑏,𝑙](25) where 𝟏 is the unit matrix, 𝐮𝑙+1 and 𝐮𝑙 are the vectors corresponding to the nodal values of the reconstructed displacements at the current (𝑙+1), and previous (𝑙) iterations, respectively. 𝐅𝑏,𝑙 and 𝐊𝑙 are the global vector of internal forces and the tangent stiffness matrix (𝜕𝐅𝑏∕𝜕𝐮), respectively, both evaluated at the 𝑙-th iteration. The coupled algorithm of solving both the forward and inverse subproblems within the global iterative procedure is explained in the next section. 2.5. Global iteration description Both forward and inverse subproblems are solved sequentially within a global iterative loop that is stopped when the global convergence criterion ‖𝜿𝒐,m+1−𝜿𝒐,m‖2 ‖𝜿𝒐,m‖2≤TOL is achieved (𝑚 being the counter of the global loop). On the one hand, the forward subproblem iterates with respect to a particular local convergence criterion defined for ‖𝜿𝒐,k+1−𝜿𝒐,k‖2 ‖𝜿𝒐,k‖2≤TOL (𝑘 being the local counter of iterations of the forward subproblem). On the other hand, the inverse subproblem proceeds up to convergence of ‖𝐮𝑖 l+1−𝐮𝑖 l‖2 ‖𝐮𝑖 l‖2≤TOL (𝑙 being the local counter of iterations of the inverse subproblem). TOL is the relative convergence tolerance, which is set as 1e-5 in this work. A schematic representation of this process is illustrated in Fig. 2. 2.6. Calibration of the regularization parameters The ‘‘weak’’ enforcement of the equilibrium condition through regularization in the forward and inverse subproblems yields regularization parameters 𝛼𝜅𝑜 and 𝛼𝑢. The calibration procedure for both parameters is based on a modified version of the frequently used L-curve method (Vitale et al., 2012a; Peschetola et al., 2013; Michel et al., 2013; Peschetola, 2011; Hansen, 2001). First, in the forward subproblem, 𝛼𝜅𝑜 is taken as large as possible to enforce equilibrium while at the same time keeping the reconstructed stiffness values in the non-negative range. Once 𝛼𝜅𝑜 is calibrated under this criterion, a sweep of simulations for different values of 𝛼𝑢 is performed to reproduce the L-curve. In order to do so, the regularization term ‖∇𝜅𝑜‖2 2 is represented against the data mismatch term ‖𝐮−𝐮∗‖2 2 for the simulated values of 𝛼𝑢. Then, the optimal value is determined as the one that is closest to the point of maximum curvature of the resulting curve. The reader is addressed to Apolinar-Fernández et al. (2024b) for a more detailed description of the L-curve method applied to the inverse subproblem. Note that 𝛼𝑢 is the only one parameter calibrated with the L-curve method, and it is not required to be repeatedly fitted for multiple values of 𝛼𝜅𝑜. 3. Setup of the synthetic 3D TFM cases This section describes the process for generating the ground truth case that will be taken as a reference for the assessment of the validity of the results. The displacements required as an input in the methodology are then taken from this ground truth case. First, information about the real 3DTFM experiment and the employed cell geometry are given. Afterwards, the configuration of the ground truth scenario, and the process of generating the input displacement fields for the in silico TFM reconstructions are explained. Finally, the error metrics used for the assessment of the performance of the methodology are introduced. 3.1. Acquisition of microscopy images, segmentation and generation of the FE mesh The in vitro setting from which the geometry of the cell considered for this study was extracted is outlined in previous papers published European Journal of Mechanics / A Solids 114 (2025) 105722 6 A. Apolinar-Fernández et al. Fig. 3. Setup of the ground truth reference case: (a) Problem domain containing both the hydrogel and the cell. (b) 3D tetrahedral Finite Element mesh. (c) Ground truth forces (arrows) prescribed on the tips of the cell’s protrusions. by our lab (Apolinar-Fernández et al., 2024b; Blázquez-Carmona et al., 2024). In this section, a brief description of the generation process of the finite element mesh corresponding to the cell geometry concerned is given. The selected cell geometry corresponds to a human breast cancer cell (line MDA-MB-231). For the associated experiment, the cell was embedded in a 1.2 mg/ml collagen-based hydrogel prepared with a mix of rat tail and bovine skin collagen, following the protocol described in Cóndor et al. (2017). 24 h after polymerization, confocal images of the cell were captured, featuring a voxel size of 0.227 𝑥 0.227 𝑥 1.48 μm3. The obtained images were then enhanced by a contrast stretching operation and subsequently segmented, resulting in a 3D binarized image of the cell geometry. Finally, the MatLab toolbox Iso2Mesh (Qianqian and Boas, 2009) was employed for the meshing of the voxelized cell boundary geometry resulting from the segmentation process. The finite element mesh is composed of 4-noded tetrahedra with linear interpolation, resulting in a total of 34666 nodes and 199691 elements, which includes the hydrogel domain and the cell in its center (Fig. 3a). 3.2. Generation of the ground-truth cases The aim of the methodology presented in this paper is to reconstruct stiffness profiles from standard input displacement data in the context of TFM. Therefore, different heterogeneous stiffness profiles 𝜅𝑜 are prescribed in the hydrogel domain of the ground truth reference case presented in Fig. 3a. The FE mesh of this case of study is shown in Fig. 3b. These profiles simulate degradation effects of the substances secreted by the cell and remodeling on the collagen hydrogel microstructure. In particular, two different stiffness heterogeneities were simulated, hereby referred to as tip-localized matrix degradation pattern (first pattern) and diffuse matrix remodeling pattern (second pattern). The tip-localized degradation pattern is generated as an exponential decay, following a sigmoid-type function in the hydrogel domain from the emission spot that is located at the tip of one of the cell’s protrusions. Fig. 4a shows a view of the mid cross-section of the hydrogel domain in which the referred ground truth degradation pattern is outlined. Fig. 4b illustrates a projection of the cross-section on the XY-plane, and Fig. 4c plots the ground-truth stiffness profile with respect to the distance from the emission point (𝛾) along the green line represented in Fig. 4b. Secondly, the diffuse matrix remodeling pattern is shown in Fig. 5. Differently to the tip-localized matrix degradation pattern, this one outlines the process of matrix remodeling which induces stiffening of the ECM around the vicinity of the cell’s surface. This can be attributed to the process of cell-mediated collagen deposition, followed by matrix degradation. Then, nominal (virgin) hydrogel stiffness is assumed far from the cell location. This kind of stiffness distribution is typically seen in in vitro setups (Calvo et al., 2013; Liu et al., 2019; Zhao et al., 2024; Iordan et al., 2010). Afterwards, for the two considered stiffness profiles, nodal forces are prescribed on the tips of the cell’s protrusions and null normal displacement boundary conditions are applied to the hydrogel boundaries. The values of the forces were calibrated such that the magnitude of the associated deformations is aligned with our observations in the laboratory. Fig. 3c illustrates the regions where these cellular forces are prescribed. The field variables of interest (stresses, strains, tractions, displacements) are retrieved from the solution to the corresponding nonlinear elasticity problem simulated by means of the FE software Simulia ABAQUS, and using the described mechanical behavior of the hydrogel through an external UMAT subroutine. Together with the material stiffness distribution, the recovery of these variables using the presented inverse methodology, are compared versus the ground truth solution for the different levels of noise considered in the study. The recovery of the ground truth solutions corresponding to the 2 stiffness patterns (Figs. 4and 5) are compared using our proposed methodology versus the frequently used L-BFGS iterative algorithm (Section 4.1). 3.3. Input displacement fields and error metric In order to generate the input data of the inverse algorithm described in Section 2, the ground truth displacements (obtained as explained in Section 3.3) are ‘‘corrupted’’ by adding a Gaussian noise to the value of each degree of freedom of the FEM-discretized ground truth displacement vector, 𝐮∗∼(𝐮GT, 𝜁)(26) with 𝐮∗ being the measured (simulated) displacement vector and 𝐮GT the ground truth displacement vector. (𝐮GT, 𝜁) is a normal distribution with mean the values of vector 𝐮GT, and 𝜁 its standard deviation with respect to the ground truth displacements. Three different levels of noise have been considered for assessment: 2% (𝜁= 0.000292), 5% (𝜁= 0.000741) and 10% (𝜁= 0.001475) error in displacements. On the other hand, in order to quantify the accuracy of the method in regards to material stiffness reconstruction at the different noise levels employed in this study, a specific error metric is defined as follows: 𝑆NL e= 100 ⋅‖𝜿NL 𝒐−𝜿GT 𝒐‖2 ‖𝜿𝐺𝑇 𝒐‖2 (27) European Journal of Mechanics / A Solids 114 (2025) 105722 7 A. Apolinar-Fernández et al. Fig. 4. Material stiffness 𝜅𝑜 contours for the ground truth case (tip-localized matrix degradation pattern): (a) 3D view of mid cross-section of the hydrogel domain outlining spatial material stiffness distribution. (b) XY-plane projection of the mid cross section. (c) Stiffness plotted with respect to the distance from the emission point (cell’s tip) highlighted as a green line drawn in subfigure (b). where 𝜿NL 𝒐 represents the vector of the (interpolated) stiffness at hydrogel points along the reference line (𝛾) shown in Fig. 4b for a certain noise level (NL), and 𝜿𝐺𝑇 𝒐 being the corresponding ground truth case. 4. Results and discussion In this section, the main results of the study are presented. First, the proposed methodology is compared with an iterative algorithm usually employed in the context of the application on hand. After this, the results corresponding to the recovery of variables performed for synthetic 3DTFM cases are presented, for three different levels of noise added to the ground truth solution. The resulting data were postprocessed in different ways in order to assess the accuracy that the methodology provides. In particular, its performance is evaluated based on the accuracy of reconstruction for stiffness, tractions and displacements. In terms of computer performance, the number of iterations taken in each case, until local and global convergence are achieved, are discussed. Both local convergence and the global convergence criteria were established as a relative change (tolerance) between iterations of 1e-5 (as described in Section 2.5). 4.1. Noise-free comparison between NR/FEM and L-BFGS algorithms In order to establish a comparison between the proposed methodology and the frequently used Limited-memory Broyden–Fletcher– Goldfarb–Shanno (L-BFGS) iterative algorithm (Song et al., 2020a; Khang et al., 2023; Mei et al., 2023), noise-free recovery of ground truth solutions are performed. The L-BFGS algorithm can be quite straightforwardly implemented as it just requires the cost function (Eq. (5)) and its gradient (Eq. (6)), the latter being utilized to identify the direction of steepest descent, and also to form an estimate of the Hessian matrix (second derivative), which is not explicitly defined, as opposed to the methodology proposed in this paper (NR/FEM). Two degradation patterns are considered for reconstruction given in Figs. 4and 5. As displacements are considered as noise-free, just the forward subproblem is required for the comparison between solvers. Figs. 6and 7 show the reconstructed ground truth stiffness profiles for the tip-localized matrix degradation and diffuse matrix remodeling patterns, respectively, using the L-BFGS iterative solver (Figs. 6a and 7a), and the NR/FEM proposed methodology (Figs. 6b and 7b). These results are qualitatively compared to the ground truth solution in Figs. 6c and 7c; and quantitatively in (Figs. 6d and 7d). On the other hand Table 2 provides, for each methodology, the values of the defined performance metrics for reconstruction of each degradation pattern. Qualitatively, it can be observed in Figs. 6and 7 that both methodologies are able to successfully recover the ground truth solution, being the NR/FEM reconstruction overall better than the L-BFGS solver. Nevertheless, Table 2 reveals that, although L-BFGS takes a tenth of the time per iteration that NR/FEM, as no inversion of any matrix is required in this algorithm, the number of total iterations is significantly bigger for the L-BFGS. In fact, the predefined number of iterations (2000) was achieved before the predefined tolerance for convergence (the same for both methodologies) was reached for the tip-localized pattern. As a result, the total time of execution of L-BFGS relative European Journal of Mechanics / A Solids 114 (2025) 105722 8 A. Apolinar-Fernández et al. Fig. 5. Material stiffness 𝜅𝑜 contours for the ground truth case (diffuse matrix remodeling pattern): (a) 3D view of mid cross-section of the hydrogel domain outlining spatial material stiffness distribution. (b) XY-plane projection of the mid cross section. (c) Stiffness plotted with respect to the distance from the emission point (cell’s tip) highlighted as a green line drawn in subfigure (b). Table 2 Performance metrics for both the L-BFGS algorithm and the proposed methodology introduced in this paper (NR/FEM) (in downward direction): Number of iterations, CPU time per iteration and total CPU time taken for the reconstruction. The calibrated regularization parameters of the forward subproblem 𝛼𝜅𝑜 are given for each degradation pattern. Tip-localized Diffuse 𝛼𝜅𝑜=5e4 𝛼𝜅𝑜=1e3 L-BFGS NR/FEM L-BFGS NR/FEM 𝑵iter [-] 2000 49 1967 90 𝑻iter [s] 1.13 8.16 1.17 8.41 𝑻total [s] 2256 400 2311 757 to NR/FEM is nearly sixfold in the case of the tip-localized matrix degradation pattern, and around threefold in the case of the diffuse matrix remodeling pattern. Also, it has to be taken into account that these times of execution correspond to the forward subproblem. The full iteration procedure (both forward and inverse subproblems, see Fig. 2) is required when noisy displacement fields are considered. Therefore, the tangent stiffness matrix has to be updated and assembled with the new displacement field during iterations. As a consequence, given the discrepancy between the number of iterations performed by L-BFGS and NR/FEM, the total CPU time will heavily increase, highlighting that the proposed NR/FEM scheme outperforms the L-BFGS iterative solver. 4.2. Results of reconstruction with noisy displacements Results of reconstructed variables, following the proposed methodology introduced in Section 2, are given in this section. Different levels of noise in the input displacements are considered, using Eq. (26), for the tip-localized degradation pattern case. Figs. 8, 9and 10 show a visual representation of the global reconstructed stiffness profiles for noise levels 2%, 5% and 10%, respectively. The ground truth solution is also included in these figures for comparison purposes. It can be observed in Figs. 8a, 9a and 10a that, in general terms, the recovered profiles are similar to that of the ground truth case. Moreover, it is shown in Figs. 8c, 9c and 10c that the stiffness distribution progressively struggles to capture the values located in the vicinity of the emission point on the cell protrusion as noise level increases. Moreover, Table 3 quantitatively provides the error metric 𝑆𝑒 (Eq. (27)) of the stiffness reconstruction, showing that this error steadily increases from 5.31% to 8.56%, when the added noise in the displacement fields increases from 2% to 10%. In addition, Table 3 provides the number of global and local iterations required for achieving convergence of the method. It can be seen that, as noise level increases, the convergence of the iterative processes becomes more difficult, requiring generally more iterations. Particularly, the amount of both local iterations and that of the global iterations steadily increases with noise level, this tendency being substantially emphasized for the 10% noise level case. In the latter, the number of iterations required for the convergence European Journal of Mechanics / A Solids 114 (2025) 105722 9 A. Apolinar-Fernández et al. Fig. 12. Traction magnitude contours for each of the analyzed cases: (a) 2% noise. (b) 5% noise. (c) 10% noise. (d) Ground Truth solution. European Journal of Mechanics / A Solids 114 (2025) 105722 16 A. Apolinar-Fernández et al. In the UMAT subrountine, it is necessary to explicitly define both the Cauchy stress tensor and the tangent stiffness matrix DDSDDE as a function of the deformation gradient i.e. 𝐽c𝐽. Using the inherent symmetries of both quantities, these need to be expressed in the particular version of the Voigt notation that is utilized by ABAQUS, 𝝈=⎡⎢⎢⎢⎢⎢⎢⎣ 𝜎11 𝜎22 𝜎33 𝜎12 𝜎13 𝜎23 ⎤⎥⎥⎥⎥⎥⎥⎦ 𝐃𝐃𝐒𝐃𝐃𝐄 =⎡⎢⎢⎢⎢⎢⎢⎣ 𝐷1111 𝐷1122 𝐷1133 𝐷1112 𝐷1113 𝐷1123 𝐷2211 𝐷2222 𝐷2233 𝐷2212 𝐷2213 𝐷2223 𝐷3311 𝐷3322 𝐷3333 𝐷3312 𝐷3313 𝐷3323 𝐷1211 𝐷1222 𝐷1233 𝐷1212 𝐷1213 𝐷1223 𝐷1311 𝐷1322 𝐷1333 𝐷1312 𝐷1313 𝐷1323 𝐷2311 𝐷2322 𝐷2333 𝐷2312 𝐷2313 𝐷2323 ⎤⎥⎥⎥⎥⎥⎥⎦(A.7) A.2. Computational implementation of the nonlinear hyperelastic model A general description of the nonlinear hyperelastic model considered in this paper is given in Section 2.1. Here we gather the expressions and mathematical elaborations corresponding to the rest of the variables required for its implementation within the NR/FEM computational framework that were not desribed in Section 2.1. From the expression for the mechanical stiffness of a single fiber (Eq. (1)), its mechanical stress can be defined as, 𝛹′(𝜆) = ∫𝜆 0 𝛹′′(𝜆)𝑑𝜆 =𝜅𝑜⎧ ⎪ ⎨ ⎪ ⎩ 𝑑𝑜⋅[exp [𝜆∕𝑑𝑜]−1]∀𝜆 < 0 𝜆∀0 ≤𝜆 < 𝜆𝑠 𝜆𝑠+𝑑𝑠⋅[exp [[𝜆−𝜆𝑠]∕𝑑𝑠]− 1]∀𝜆≥𝜆𝑠 (A.8) in which integration from zero indicates that the material is not prestressed. From the constitutive relation in mixed configuration (Eq. (2)), the Cauchy stress tensor is obtained by applying a pushforward operation as follows, 𝝈=1 𝐽𝜙∗[𝑷] = 1 𝐽𝑷 𝑭 𝑇 =1 𝐽[1 4𝜋∫2𝜋 0∫𝜋 0 𝛹′(𝜆)[𝑭 𝒆𝛺]⊗𝒆𝛺 ‖𝑭 𝒆𝛺‖sin 𝜃𝑑𝜃𝑑𝜙]⋅𝑭𝑇(A.9) where the integral over the unit sphere is approximated by means of the Gauss–Legendre quadrature. The definition of the tangent stiffness matrix DDSDDE depends ultimately on the definition of the Second elasticity tensor , which can be transformed using Eqs. (A.5) and (A.6) to accommodate to the ABAQUS format. The fourth invariant of the Right Cauchy–Green deformation tensor is described as, 𝐼4=𝒆𝑇 𝛺⋅𝑪⋅𝒆𝛺=𝒆𝑇 𝛺⋅𝑭𝑇𝑭⋅𝒆𝛺=‖𝑭 𝒆𝛺‖2= (𝜆+ 1)2(A.10) which enables the rewriting of 𝜆 in terms of the Right Cauchy–Green deformation tensor, 𝜆(𝑭) = ‖𝑭 𝒆𝛺‖− 1 ⟶𝜆(𝑪) = (𝐼4)1∕2 − 1 (A.11) As a consequence, the constitutive relation in the reference or material configuration is obtained as, 𝑺= 2 ⟨𝜕𝛹 (𝜆(𝑪)) 𝜕𝑪⟩𝜴 = 2 ⟨𝜕𝛹 𝜕𝜆 𝜕𝜆 𝜕𝑪⟩𝜴 =⟨[𝜆+ 1]−1 ⋅𝛹′(𝜆)⋅[𝒆𝛺⊗𝒆𝛺]⟩𝜴 (A.12) where ⊗ denotes tensor product. The Second elasticity tensor can then be retrieved from, = 2 ⟨𝜕𝑺∗(𝜆(𝑪)) 𝜕𝑪⟩𝜴 = 2 ⟨𝜕𝑺∗ 𝜕𝜆 ⊗𝜕𝜆 𝜕𝑪⟩𝜴 (A.13) in which 𝑺∗ denotes the Second Piola–Kirchhoff stress tensor of a single fiber (non-homogenized) and, 𝜕𝑺∗ 𝜕𝜆 =[−[𝜆+ 1]−2 ⋅𝛹′(𝜆)+[𝜆+ 1]−1 ⋅𝛹′′(𝜆)]⋅[𝒆𝛺⊗𝒆𝛺](A.14) 𝜕𝜆 𝜕𝑪=1 2[𝜆+ 1]−1 ⋅[𝒆𝛺⊗𝒆𝛺](A.15) Appendix B The parameters governing the hyperelastic model considered in this paper were fitted to rheological data corresponding to a simple shear experiment performed on a 1.2 mg/ml collagen hydrogel. B.1. Shear rheology The hydrogel undergoes the action of a rotational rheometer of cone and plate, providing measurements of the stress–strain relationship in the context of a simple shear state of deformation. The deformation gradient in the case of simple shear receives the expression, 𝑭=⎡⎢⎢⎣ 1𝛾0 010 001⎤⎥⎥⎦ (B.1) where 𝛾 represents the engineering shear strain. In the experiment, the shear stress is retrieved from the measured torque and the setup of the test in the nominal configuration. Thus, the measured shear stress resulting from the experiment corresponds to the shear component of the First Piola–Kirchhoff stress tensor (nominal stress). Recalling the definition given in Eq. (2), the nominal stress is then computed as follows, 𝑃12 =𝑑𝑊 (𝑭(𝛾)) 𝑑𝐹12 =𝑑𝑊 (𝑭(𝛾)) 𝑑𝛾 =⟨𝜕𝑤(𝜆) 𝜕𝜆 𝜕𝜆 𝜕𝑭∶𝜕𝑭 𝜕𝛾 ⟩𝜴 = =⟨𝑤′(𝜆)⋅trace[𝜕𝜆 𝜕𝑭⋅[𝜕𝑭 𝜕𝛾 ]𝑇]⟩𝜴 (B.2) being, 𝜕𝑭 𝜕𝛾 =⎡⎢⎢⎣ 010 000 000⎤⎥⎥⎦ (B.3) B.2. Fitting The fitting of the model parameters have been carried out by means of the MATLAB function lsqnonlin, which employs the nonlinear version of least squares minimization. Fig. 1 (main text) shows the fitting performed to the experimental curve corresponding to a hydrogel with a fiber concentration of 1.2 mg/ml. Data availability Data will be made available on request. References Ambrosi, D., 2006. Cellular traction as an inverse problem. SIAM J. Appl. Math. 66, 2049–2060. European Journal of Mechanics / A Solids 114 (2025) 105722 17 A. Apolinar-Fernández et al. Apolinar-Fernández, A., Barrasa-Fano, J., Cóndor, M., Van Oosterwyck, H., SanzHerrera, J.A., 2023. Traction force reconstruction assessment on real threedimensional matrices and cellular morphologies. Internat. J. Engrg. Sci. 186, 10, 3828. Apolinar-Fernández, A., Barrasa-Fano, J., Oosterwyck, H.V., Sanz-Herrera, J.A., 2024a. Multiphysics modeling of 3d traction force microscopy with application to cancer cell-induced degradation of the extracellular matrix. Eng. Comput.. Apolinar-Fernández, A., Blázquez-Carmona, P., Ruiz-Mateos, R., Barrasa-Fano, J., Van Oosterwyck, H., Reina-Romo, E., Sanz-Herrera, J.A., 2024b. Regularization techniques and inverse approaches in 3D Traction Force Microscopy. Int. J. Mech. Sci. 283, 10, 9592. Bar-Kochba, E., Toyjanova, J., Andrews, E.J., Kim, K.-S., Franck, C., 2015. A fast iterative digital volume correlation algorithm for large deformations. Exp. Mech. 55, 261–274. Barrasa-Fano, J., Shapeti, A., de Jong, J., Ranga, A., Sanz-Herrera, J.A., 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. Blázquez-Carmona, P., Ruiz-Mateos, R., Barrasa-Fano, J., Shapeti, A., MartínAlfonso, J.E., Domínguez, J., Van Oosterwyck, H., Reina-Romo, E., SanzHerrera, J.A., 2024. Quantitative atlas of collagen hydrogels reveals mesenchymal cancer cell traction adaptation to the matrix nanoarchitecture. Acta Biomater. 185, 281–295. Bloom, R.J., George, J.P., Celedon, A., Sun, S.X., Wirtz, D., 2008. Mapping local matrix remodeling induced by a migrating tumor cell using three-dimensional multiple-particle tracking. Biophys. J. 95, 4077–4088. Böhringer, D., Cóndor, M., Bischof, L., Czerwinski, T., Gampl, N., Ngo, P.A., Bauer, A., Voskens, C., López-Posadas, R., Franze, K., Budday, S., Mark, C., Fabry, B., Gerum, R., 2024. Dynamic traction force measurements of migrating immune cells in 3d biopolymer matrices. Nat. Phys. 20, 1816–1823. Broguiere, N., Isenmann, L., Hirt, C., Ringel, T., Placzek, S., Cavalli, E., Ringnalda, F., Villiger, L., Züllig, R., Lehmann, R., Rogler, G., Heim, M.H., Schüler, J., ZenobiWong, M., Schwank, G., 2018. Growth of epithelial organoids in a defined hydrogel. Adv. Mater. 30, 180, (1621). Caliari, S.R., Burdick, J.A., 2016. A practical guide to hydrogels for cell culture. Nature Methods 13, 405–414. Calvo, F., Ege, N., Grande-Garcia, A., Hooper, S., Jenkins, R.P., Chaudhry, S.I., Harrington, K., Williamson, P., Moeendarbary, E., Charras, G., Sahai, E., 2013. Mechanotransduction and YAP-dependent matrix remodelling is required for the generation and maintenance of cancer-associated fibroblasts. Nature Cell Biol. 15, 637–646. Chen, S., Xu, W., Kim, J., Nan, H., Zheng, Y., Sun, B., Jiao, Y., 2019. Novel inverse finite-element formulation for reconstruction of relative local stiffness in heterogeneous extra-cellular matrix and traction forces on active cells. Phys. Biology 16, 03, (6002). Chicurel, M.E., Chen, C.S., Ingber, D.E., 1998. Cellular control lies in the balance of forces. Curr. Opin. Cell Biol. 10, 232–239. Colin-York, H., Eggeling, C., Fritzsche, M., 2017. Dissection of mechanical force in living cells by super-resolved traction force microscopy. Nat. Protoc. 12, 783–796. Colin-York, H., Shrestha, D., Felce, J.H., Waithe, D., Moeendarbary, E., Davis, S.J., Eggeling, C., Fritzsche, M., 2016. Super-resolved traction force microscopy (stfm). Nano Lett. 16, 2633–2638, PMID: 2692 (3775). Cóndor, M., Mark, C., Gerum, R.C., Grummel, N.C., Bauer, A., García-Aznar, J.M., Fabry, B., 2019. Breast cancer cells adapt contractile forces to overcome steric hindrance. Biophys. J. 116, 1305–1312. Cóndor, M., Steinwachs, J., Mark, C., García-Aznar, J., Fabry, B., 2017. Traction force microscopy in 3-dimensional extracellular matrix networks. Curr. Protoc. Cell Biology 75, 10. Daley, W.P., Peters, S.B., Larsen, M., 2008. Extracellular matrix dynamics in development and regenerative medicine. J. Cell Sci. 121, 255–264. Dong, L., Oberai, A.A., 2017. Recovery of cellular traction in three-dimensional nonlinear hyperelastic matrices. Comput. Methods Appl. Mech. Engrg. 314, 296–313, Special Issue on Biological Systems Dedicated to William S. Klug.. Doyle, A.D., Sykora, D.J., Pacheco, G.G., Kutys, M.L., Yamada, K.M., 2021. 3D mesenchymal cell migration is driven by anterior cellular contraction that generates an extracellular matrix prestrain. Dev. Cell 56, 826–841. Du, Y., Herath, S.C., Wang, Q.G., Asada, H., Chen, P.C., 2018. Determination of Green’s function for three-dimensional traction force reconstruction based on geometry and boundary conditions of cell culture matrices. Acta Biomater. 67, 215–228. Engler, A.J., Sen, S., Sweeney, H.L., Discher, D.E., 2006. Matrix elasticity directs stem cell lineage specification. Cell 126, 677–689. Feng, X., Hall, M.S., Wu, M., Hui, C.-Y., 2014. An adaptive algorithm for tracking 3d bead displacements: application in biological experiments. Meas. Sci. Technol. 25, 05. Feng, X., Hui, C.Y., 2016. Force sensing using 3D displacement measurements in linear elastic bodies. Comput. Mech. 58, 91–105. Franck, C., Maskarinec, S.A., Tirrell, D.A., Ravichandran, G., Genin, G., 2011. Threedimensional traction force microscopy: A new tool for quantifying cell-matrix interactions. PLoS One 6, e17833. Gjorevski, N., Nelson, C.M., 2012. Mapping of mechanical strains and stresses around quiescent engineered three-dimensional epithelial tissues. Biophys. J. 103, 152–162. Hall, M.S., Long, R., Feng, X., Huang, Y., Hui, C.-Y., Wu, M., 2013. Toward single cell traction microscopy within 3D collagen matrices. Exp. Cell Res. 319, 2396–2408. Hall, M.S., Long, R., Hui, C.-Y., Wu, M., 2012. Mapping Three-Dimensional Stress and Strain Fields within a Soft Hydrogel Using a Fluorescence Microscope. Biophys. J. 102, 2241–2250. Hansen, P., 2001. The l-curve and its use in the numerical treatment of inverse problems. the Advances in computational bioengineering series, In: Computational Inverse Problems in Electrocardiography 2001, vol. 5. Huang, Y., Schell, C., Huber, T.B., Şimşek, A.N., Hersch, N., Merkel, R., Gompper, G., Sabass, B., 2019. Traction force microscopy with optimized regularization and automated bayesian parameter selection for comparing cells. Sci. Rep. 9, 539. Humphrey, J.D., Dufresne, E.R., Schwartz, M.A., 2014. Mechanotransduction and extracellular matrix homeostasis. Nature Rev. Mol. Cell Biol. 15, 802–812. Ingber, D.E., 2002. Mechanical signaling and the cellular response to extracellular matrix in angiogenesis and cardiovascular physiology. Ingber, D.E., 2003. Mechanobiology and diseases of mechanotransduction. Iordan, A., Duperray, A., Gérard, A., Grichine, A., Verdier, C., 2010. Breakdown of cell-collagen networks through collagen remodeling. Biorheology 47, 277–295. Iskratsch, T., Wolfenson, H., Sheetz, M.P., 2014. Appreciating force and shape – the rise of mechanotransduction in cell biology. Nature Rev. Mol. Cell Biol. 15, 825–833. Izquierdo-Álvarez, A., Vargas, D.A., Jorge-Peñas, Á., Subramani, R., Vaeyens, M.M., Van Oosterwyck, H., 2019. Spatiotemporal analyses of cellular tractions describe subcellular effect of substrate stiffness and coating. Ann. Biomed. Eng. 47, 624–637. Janke, T., Schwarze, R., Bauer, K., 2020. Part2track: A matlab package for double frame and time resolved particle tracking velocimetry. SoftwareX 11, 10. Janmey, P.A., Fletcher, D.A., Reinhart-King, C.A., 2020. Stiffness sensing by cells. Physiol. Rev. 100, 695–724. Jansen, K., Licup, A., Sharma, A., Rens, R., MacKintosh, F., Koenderink, G., 2018. The role of network architecture in collagen mechanics. Biophys. J. 114, 2665–2678. Javanmardi, Y., Colin-York, H., Szita, N., Fritzsche, M., Moeendarbary, E., 2021. Quantifying cell-generated forces: Poisson’s ratio matters. Commun. Phys. 4, 237. Jerrell, R.J., Parekh, A., 2014. Cellular traction stresses mediate extracellular matrix degradation by invadopodia. Acta Biomater. 10, 1886–1896. Jorge-Peñas, A., Bové, H., Sanen, K., Vaeyens, M.-M., Steuwe, C., Roeffaers, M., Ameloot, M., Van Oosterwyck, H., 2017. 3D full-field quantification of cellinduced large deformations in fibrillar biomaterials by combining non-rigid image registration with label-free second harmonic generation. Biomaterials 136, 86–97. Jorge-Peñas, A., Izquierdo-Alvarez, A., Aguilar-Cuenca, R., Vicente-Manzanares, M., Garcia-Aznar, J.M., Van Oosterwyck, H., de Juan-Pardo, E.M., Ortiz-de Solorzano, C., Muñoz-Barrutia, A., 2015. Free form deformation–based image registration improves accuracy of traction force microscopy. PLoS One 10, 1–22. Keating, M., Kurup, A., Alvarez-Elizondo, M., Levine, A., Botvinick, E., 2017. Spatial distributions of pericellular stiffness in natural extracellular matrices are dependent on cell-mediated proteolysis and contractility. Acta Biomater. 57, 304–312. Khang, A., Steinman, J., Tuscher, R., Feng, X., Sacks, M.S., 2023. Estimation of aortic valve interstitial cell-induced 3D remodeling of poly(ethylene glycol) hydrogel environments using an inverse finite element approach. Acta Biomater. 160, 123–133. Khetan, S., Guvendiren, M., Legant, W.R., Cohen, D.M., Chen, C.S., Burdick, J.A., 2013. Degradation-mediated cellular traction directs stem cell fate in covalently crosslinked three-dimensional hydrogels. Nat. Mater. 12, 458–465. Kim, S., Uroz, M., Bays, J.L., Chen, C.S., 2021. Harnessing mechanobiology for tissue engineering. Dev. Cell 56, 180–191. Kleiner, D.E., Stetler-Stevenson, W.G., 1999. Matrix metalloproteinases and metastasis. Cancer Chemother. Pharmacol. 43, S42–S51. Kraning-Rush, C.M., Califano, J.P., Reinhart-King, C.A., 2012. Cellular traction stresses increase with increasing metastatic potential. PLoS One 7, 1–10. Laforgue, L., Fertin, A., Usson, Y., Verdier, C., Laurent, V.M., 2022. Efficient deformation mechanisms enable invasive cancer cells to migrate faster in 3d collagen networks. Sci. Rep. 12, 7867. Lammerding, J., Kamm, R.D., Lee, R.T., 2004. Mechanotransduction in cardiac myocytes. Ann. New York Acad. Sci. 1015, 53–70. Legant, W.R., Miller, J.S., Blakely, B.L., Cohen, D.M., Genin, G.M., Chen, C.S., 2010. Measurement of mechanical tractions exerted by cells in three-dimensional matrices. Nature Methods 7, 969–971. Li, D., Colin-York, H., Barbieri, L., Javanmardi, Y., Guo, Y., Korobchevskaya, K., Moeendarbary, E., Li, D., Fritzsche, M., 2021. Astigmatic traction force microscopy (atfm). Nat. Commun. 12, 2168. Liu, T., Zhou, L., Li, D., Andl, T., Zhang, Y., 2019. Cancer-associated fibroblasts build and secure the tumor microenvironment. Front. Cell Dev. Biology 7. López-Otín, C., Overall, C.M., 2002. Protease degradomics: A new challenge for proteomics. Nature Rev. Mol. Cell Biol. 3, 509–519. Makarchuk, S., Beyer, N., Gaiddon, C., Grange, W., Hébraud, P., 2018. Holographic traction force microscopy. Sci. Rep. 8, 3038. Mammoto, T., Mammoto, A., Ingber, D.E., 2013. Mechanobiology and developmental control. Annu. Rev. Cell. Dev. Biol. 29, 27–61. European Journal of Mechanics / A Solids 114 (2025) 105722 18 A. Apolinar-Fernández et al. Martin, A.C., Kaschube, M., Wieschaus, E.F., 2009. Pulsed contractions of an actin–myosin network drive apical constriction. Nature 457, 495–499. Mei, Y., Feng, X., Jin, Y., Kang, R., Wang, X., Zhao, D., Ghosh, S., Neu, C.P., Avril, S., 2023. Cell nucleus elastography with the adjoint-based inverse solver. Comput. Methods Programs Biomed. 242, 10, (7827). Michel, R., Peschetola, V., Vitale, G., tienne, J., Duperray, A., Ambrosi, D., Preziosi, L., Verdier, C., 2013. Mathematical framework for traction force microscopy. ESAIM: Proc. 42, 61–83. Mui, K.L., Chen, C.S., Assoian, R.K., 2016. The mechanical regulation of integrin– cadherin crosstalk organizes cells, signaling and forces. J. Cell Sci. 129, 1093–1100. Mulligan, J.A., Bordeleau, F., Reinhart-King, C.A., Adie, S.G., 2017. Measurement of dynamic cell-induced 3d displacement fields in vitro for traction force optical coherence microscopy. Biomed. Opt. Express 8, 1152–1171. Mulligan, J.A., Bordeleau, F., Reinhart-King, C.A., Adie, S.G., 2018. Traction force microscopy for noninvasive imaging of cell forces. Adv. Exp. Med. Biol. 1092, 319–349. Mulligan, J., Feng, X., A. S.G, 2019. Quantitative reconstruction of time-varying 3d cell forces with traction force optical coherence microscopy. Sci. Rep. 9, 4086. Nelson, C.M., Jean, R.P., Tan, J.L., Liu, W.F., Sniadecki, N.J., Spector, A.A., Chen, C.S., 2005. Emergent patterns of growth controlled by multicellular form and mechanics. Proc. Natl. Acad. Sci. 102, 1. Peschetola, V., 2011. Détermination des forces de traction au cours de la migration de cellules cancéreuses sur des gels. (Ph.D. thesis). Université de Grenoble. Peschetola, V., Laurent, V.M., Duperray, A., Michel, R., Ambrosi, D., Preziosi, L., Verdier, C., 2013. Time-dependent traction force microscopy for cancer cells as a measure of invasiveness. Cytoskeleton 70, 201–214. Polacheck, W.J., Chen, C.S., 2016. Measuring cell-generated forces: A guide to the available tools. Nature Methods 13, 415–423. Qianqian, F., Boas, D., 2009. Tetrahedral mesh generation from volumetric binary and grayscale images. IEEE Int. Symp. Biomed. Imaging: From Nano To Macro 1142–1145. Rauzi, M., Verant, P., Lecuit, T., Lenne, P.-F., 2008. Nature and anisotropy of cortical forces orienting Drosophila tissue morphogenesis. Nature Cell Biol. 10, 1401–1410. Riley, G.P., Curry, V., DeGroot, J., van El, B., Verzijl, N., Hazleman, B.L., Bank, R.A., 2002. Matrix metalloproteinase activities and their relationship with collagen remodelling in tendon pathology. Matrix Biol. 21, 185–195. Ross, R., 1986. The pathogenesis of atherosclerosis – an update. N. Engl. J. Med. 314, 488–500. Sanz-Herrera, J., 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. Sanz-Herrera, J.A., 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, 1, (0210)–1 (0222). Schwarz, U.S., Soiné, J.R.D., 2015. Traction force microscopy on soft elastic substrates: A guide to recent computational advances. Biochim. Biophys. Acta - Mol. Cell Res. 1853, 3095–3104. Shapeti, A., Barrasa-Fano, J., Fattah, A.R.A., de Jong, J., Sanz-Herrera, J.A., Pezet, M., Assou, S., de Vet, E., Elahi, S.A., Ranga, A., Faurobert, E., Van Oosterwyck, H., 2024. Force-mediated recruitment and reprogramming of healthy endothelial cells drive vascular lesion growth. Nat. Commun. 15, 8660. Siedlik, M.J., Varner, V.D., Nelson, C.M., 2016. Pushing, pulling, and squeezing our way to understanding mechanotransduction. Methods 94, 4–12. Song, D., Dong, L., Gupta, M., Li, L., Klaas, O., Loghin, A., Beall, M., Chen, C.S., Oberai, A.A., 2020b. Recovery of tractions exerted by single cells in three-dimensional nonlinear matrices. J. Biomech. Eng. 142. Song, D., Seidl, D.T., Oberai, A.A., 2020a. Three-dimensional traction microscopy accounting for cell-induced matrix degradation. Comput. Methods Appl. Mech. Engrg. 364, 11, 2935. Steinwachs, J., Metzner, C., Skodzek, K., Lang, N., Thievessen, I., Mark, C., Münster, S., Aifantis, K.E., Fabry, B., 2016. Three-dimensional force microscopy of cells in biopolymer networks. Nature Methods 13, 171–176. Style, R.W., Boltyanskiy, R., German, G.K., Hyland, C., MacMinn, C.W., Mertz, A.F., Wilen, L.A., Xu, Y., Dufresne, E.R., 2014. Traction force microscopy in physics and biology. Soft Matter 10, 4047–4055. Suñé-Auñón, A., Jorge-Peñas, A., Aguilar-Cuenca, R., Vicente-Manzanares, M., Van Oosterwyck, H., Muñoz-Barrutia, A., 2017. Full l1-regularized traction force microscopy over whole cells. BMC Bioinformatics 18, 365. 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, 1–12. Toyjanova, J., Hannen, E., Bar-Kochba, E., Darling, E.M., Hennan, D.L., Franck, C., 3D viscoelastic traction force microscopy. Soft Matter 10 (40), 2014. Trappmann, B., Baker, B.M., Polacheck, W.J., Choi, C.K., Burdick, J.A., Chen, C.S., 2017. Matrix degradability controls multicellularity of 3D cell migration. Nat. Commun. 8, 371. Vaeyens, M.M., Jorge-Peñs, 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, 315–324. Vincent, L.G., Engler, A.J., 2013. Post-degradation forces kick in. Nat. Mater. 12, 384–386. Vining, K.H., Mooney, D.J., 2017. Mechanical forces direct stem cell behaviour in development and regeneration. Nature Rev. Mol. Cell Biol. 18, 728–742. Vitale, G., Preziosi, L., Ambrosi, D., 2012a. A numerical method for the inverse problem of cell traction in 3D. Inverse Problems 28. Vitale, G., Preziosi, L., Ambrosi, D., 2012b. Force traction microscopy: An inverse problem with pointwise observations. J. Math. Anal. Appl. 395, 788–801. Vogel, V., Sheetz, M., 2006. Local force and geometry sensing regulate cell functions. Wang, J., Zhang, S., Jin, Y., Qin, G., Yu, L., Zhang, J., 2007. Elevated levels of platelet– monocyte aggregates and related circulating biomarkers in patients with acute coronary syndrome. Int. J. Cardiol. 115, 361–365. Werb, Z., 1997. Ecm and cell surface proteolysis: Regulating cellular ecology. Cell 91, 439–442. Zhang, Q., An, Z.-Y., Jiang, W., Jin, W.-L., He, X.-Y., 2023. Collagen code in tumor microenvironment: Functions. Mol. Mech. Ther. Implic. Biomed. Pharmacother. 166, 115390. Zhao, F., Zhang, M., Nizamoglu, M., Kaper, H.J., Brouwer, L.A., Borghuis, T., Burgess, J.K., Harmsen, M.C., Sharma, P.K., 2024. Fibroblast alignment and matrix remodeling induced by a stiffness gradient in a skin-derived extracellular matrix hydrogel. Acta Biomater. 182, 67–80. European Journal of Mechanics / A Solids 114 (2025) 105722 19