Traction force reconstruction assessment on real three-dimensional matrices and cellular morphologies

International Journal of Engineering Science 186 (2023) 103828 Available online 20 February 2023 0020-7225/© 2023 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). Contents lists available at ScienceDirect International Journal of Engineering Science journal homepage: www.elsevier.com/locate/ijengsci Traction force reconstruction assessment on real three-dimensional matrices and cellular morphologies Alejandro Apolinar-Fernándeza, Jorge Barrasa-Fanob, Mar Cóndorb, Hans Van Oosterwyckb,c, José A. Sanz-Herreraa,∗ aEscuela Técnica Superior de Ingeniería, Universidad de Sevilla, Spain bBiomechanics section, Department of Mechanical Engineering, KU Leuven, Belgium cPrometheus division of Skeletal Tissue Engineering, KU Leuven, Belgium ARTICLE INFO Keywords: Traction force microscopy Mechanobiology Nonlinear mechanics Computational mechanics Finite element method ABSTRACT Traction force microscopy (TFM) allows to estimate tractions on the surface of cells when they mechanically interact with hydrogel substrates that mimic the extracellular matrix (ECM). The field of mechanobiology has a strong interest in using TFM in 3D in vitro models. However, there are a number of challenges that hamper the accuracy of 3D TFM and that are often bypassed. In this study, the computational efficiency and accuracy of TFM, referred to traction reconstruction from synthetically generated (control) ground truth solutions, are assessed from four different perspectives: magnitude of cellular pulling force (and hence strain level achieved in the hydrogel), effect of the complexity of the cellular morphology, accuracy and computational efficiency of forward vs inverse traction recovery methods, and the effect of incorrectly selecting a constitutive model that describes the behavior of the ECM (i.e. linear/nonlinear). The main results showed: (i) traction reconstruction is more challenging for complex cell morphologies, (ii) there is no significant impact of the magnitude of cellular pulling force on the overall reconstruction accuracy, and (iii) modeling a nonlinear hydrogel with a linear constitutive model leads to non-negligible errors (up to 80% and 30% for forward and inverse methodologies, respectively) in traction reconstruction. This study expands the characterization of the accuracy and efficiency of 3D TFM, highlighting important factors to be considered in future 3D TFM in vitro applications. 1. Introduction The mechanical interactions of cells with the extracellular matrix (ECM) are crucial for multiple physiological and pathological processes (Ingber,2003;Mammoto et al.,2013;Vogel & Sheetz,2006). Over the last decades, the field of mechanobiology has discovered that cell–ECM mechanical interactions can drive processes such as cancer invasion (Jiang et al.,2019;Kopanska et al.,2016;Kumar & Weaver,2009;Medina et al.,2019;Nguyen-Ngoc et al.,2012;Peng et al.,2019;Shields et al.,2007), angiogenesis (Elliott et al.,2015;Ingber,2002;Vaeyens et al.,2020) or stem cell differentiation (Discher et al.,2005;Engler et al., 2006;Guilak et al.,2009;Huebsch et al.,2010). In vitro models are typically used to obtain fundamental understanding of cellular mechanotransduction pathways that can become pivotal in developing tissue engineering and therapeutic techniques (Janmey et al., 2020;Jansen et al.,2015;Kim et al.,2021;Kutys & Chen,2016;Vining & Mooney,2017). ∗Correspondence to: 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] (M. Cóndor), [email protected] (H. Van Oosterwyck), [email protected] (J.A. Sanz-Herrera). https://doi.org/10.1016/j.ijengsci.2023.103828 Received 11 June 2022; Received in revised form 25 January 2023; Accepted 5 February 2023 International Journal of Engineering Science 186 (2023) 103828 2 A. Apolinar-Fernández et al. In order to extract information from these systems, researchers make use of optical microscopy combined with computational methods that include microscopy image analysis techniques and mechanical models. One such methodology is Traction Force Microscopy (TFM), which has been the main choice for calculating forces exerted by cells on the ECM for the last decades (Butler et al.,2002;Dembo & Wang,1999;Harris et al.,1980;Sabass et al.,2008). In TFM, cells are seeded on top of or embedded in a substrate or hydrogel that contains fiducial markers (often fluorescent beads) and that mimics the ECM. Microscopy images are acquired before (with mechanically active cells; stressed state) and after (relaxed state) cell removal or cell mechanical inhibition. The matrix deformations caused by cellular forces lead to changes in the position of the markers, which can be measured by means of image analysis algorithms that compare the stressed state to the relaxed state images. Finally, a constitutive model that describes the mechanical behavior of the substrate/hydrogel is used to calculate strain, stress and traction forces from the measured displacement field. Several ways of calculating traction forces from the measured displacements have been proposed in the literature. On the one hand, forward (also found in the literature as direct) methods calculate the strain field by differentiating the measured displacements, and the stress/traction field is computed by means of the constitutive relation for the ECM (Franck et al.,2011;Gjorevski et al., 2015;Toyjanova, Bar-Kochba, et al.,2014). On the other hand, inverse methods search for a traction solution that is consistent with the measured displacements while being subject to a certain constraint that limits the effect of measurement noise in the resultant traction field (Feng & Hui,2016;Legant et al.,2010). TFM is very well established for 2D in vitro cultures where cells are seeded on top of a substrate. Polyacrylamide substrates are typically chosen due to their linear elastic properties. This allows for using simple analytical formulations such as the Boussinesq solution, which heavily eases traction recovery (Huang et al.,2020;Izquierdo-Álvarez et al.,2019). However, the dimensionality of the surrounding microenvironment is crucial for cell behavior and the field considers 3D in vitro cultures as more physiologically relevant (Duval et al.,2017;Vogel & Sheetz,2006). As a result, researchers have started using other bio-mimetic materials that allow for embedding cells in a 3D environment. Synthetic materials of highly tunable properties such as Polyethylene Glycol (PEG) as well as natural materials that more closely resemble physiological conditions such as collagen or fibrin hydrogels have been primarily used (Caliari & Burdick,2016). This has lead to the development of more complex 3D TFM methodologies in the last twelve years. The approach chosen to recover cellular tractions will vary depending on two factors. First, the stress–strain relationship of the hydrogel can be linear or nonlinear. While PEG hydrogels can be tuned to display a linear constitutive response, fibrillar hydrogels are characterized by nonlinear responses and strain stiffening. Second, researchers often use the small strain theory to describe the behavior of the hydrogel. However, when deformations are large, this theory is no longer valid and the finite-strain theory should be adopted. The literature offers a vast number of approaches combining the previous two factors differently. 3D traction forces have been calculated in PEG hydrogels using linear elastic models and under a small strain elasticity theory for single cells (Legant et al., 2010) and for angiogenic sprouts (Barrasa-Fano, Shapeti, de Jong, et al.,2021;Barrasa-Fano, Shapeti, Jorge-Peñas, et al.,2021). Collagen hydrogels have been modeled as a linear elastic material (Du et al.,2018;Gjorevski et al.,2015) as well as a hyperelastic material such as Neo-Hookean (Hervas-Raluy et al.,2021;Sanz-Herrera et al.,2021;Song, Seidl, & Oberai,2020). A particular case of application in this study is the work of Steinwachs et al. (2016), who developed a specific nonlinear constitutive model that best fit the experimentally acquired (by means of e.g. shear rheology mechanical testing) strain–stress curves of a collagen hydrogel. However, the choice of a constitutive model or an elasticity theory are often taken arbitrarily either to simplify the study or to highlight its increased complexity with respect to the state of the art. There is a lack of systematic analysis on the effect of incorrectly selecting a model. Recently, Song, Dong, et al. (2020) developed a traction recovery algorithm that accounts for both material and geometric nonlinearities. Interestingly, they quantified the effect of neglecting nonlinearity by resolving cellular forces synthetically generated with a nonlinear model (ground truth) with a linear model. However, they only explored this effect for an arbitrary traction field of a given magnitude. Existing studies have not analyzed the effect of the magnitude of the cell tractions on the accuracy of TFM. In a different work, Zündel et al. (2017) performed an in silico study to analyze the effect on TFM errors of different constitutive behaviors of the hydrogel, noisy displacement fields, and focal adhesions. The analysis was established for soft hyperelastic planar substrates with finite size. The authors found an underestimation in tractions’ reconstruction for limited resolution of focal adhesions, and finite inter-bead distance. Similarly, Barrasa-Fano, Shapeti, de Jong, et al. (2021) previously reported on the severe impact of the size and distribution of focal adhesions on 3D TFM accuracy. Indeed, it is important to analyze factors that lead to traction fields that are difficult to recover. Another such factor is the complexity of the geometry of a cell. While initial studies validated their traction recovery algorithm in simplified geometries such as cylinders (Du et al.,2018), spheres (Holenstein et al.,2019) or ellipsoids (Vitale et al.,2012), later studies have used real cell geometries obtained from the segmentation of a microscopy image of a cell to make their results more representative of a real scenario (Barrasa-Fano, Shapeti, de Jong, et al.,2021;Hervas-Raluy et al., 2021;Sanz-Herrera et al.,2021;Song, Seidl, & Oberai,2020). However, these studies always perform simulations on a sole cell geometry. Cell geometries can be highly heterogeneous going from spherical to spread and protrusive. This depends on multitude of factors such as cell type, hydrogel properties (adhesiveness, degradability, stiffness, porosity, etc.), or incubation time. There are no studies that systematically report on how traction calculation accuracy is affected by the geometry of a cell. The aim of this study is to assess the accuracy and efficiency of TFM, with regard to traction reconstruction versus a (controlled) synthetic ground truth solution, from different perspectives: cellular pulling force magnitude and hence strain level achieved in the hydrogel, effect of cellular morphology complexity, accuracy and computational efficiency of forward vs inverse traction reconstruction approaches, and the effect of incorrectly selecting a constitutive model that describes the behavior of the ECM (i.e. linear/nonlinear). The paper is organized as follows: first, the theoretical bases of the methods used in this study (forward International Journal of Engineering Science 186 (2023) 103828 3 A. Apolinar-Fernández et al. and inverse) are briefly introduced. Then, the different cases of study and in silico TFM methodology and simulations are explained in Section 3. The corresponding results are detailed in Section 4, and discussed in Section 5. Finally, the most important conclusions are drawn at the end of the paper. 2. Theoretical background Traction Force Microscopy can be performed by means of two fundamentally distinct methodologies commonly referred to as the forward and the inverse methods (Schwarz & Soiné,2015). The forward method computes the cellular tractions and hydrogel stresses directly from the measured displacements, and the assumed constitutive law of the hydrogel. Alternatively, the inverse method uses this displacement field as input data in order to obtain a traction field that meets certain criteria within a minimization problem framework, based on a regularization strategy. The forward method is less computationally demanding than the inverse method, but the latter provides more precise traction field reconstructions (Barrasa-Fano, Shapeti, de Jong, et al.,2021). The methods that we use in this study are a conventional forward (direct traction computation from measured displacements via the material’s constitutive law), and an inverse formulation which searches for a new displacement field that approximates the measured one, ensuring that force equilibrium within the hydrogel is satisfied (Sanz-Herrera et al.,2021). Both the forward and inverse methods were developed and implemented in the context of both linear elasticity and finite strain hyperelasticity, under a finite element (FE) approach. 2.1. Forward method The forward method is formulated in this section within the framework of finite strain hyperelasticity. Once given the reconstructed displacements field 𝐮𝑚referring to the spatial (Eulerian) description of the medium 𝐱, the position of a point in the material description (Lagrangian) 𝐗is given by, 𝐗=𝐱−𝐮𝑚(1) Then, the deformation gradient is defined as, 𝐅=𝜕𝐱 𝜕𝐗(2) which, by using (1), yields the expression, 𝐅−1 =𝐈−𝜕𝐮𝑚 𝜕𝐱(3) being 𝐈the identity matrix. It is important to recall that the reference configuration in TFM corresponds to the deformed hydrogel and stressed cell state after which a relaxed (assumed stress free) hydrogel configuration is achieved. After a FE mesh discretization, the displacement field used in this forward method (𝐮𝑓 𝑤𝑑 ) is approximated as follows, 𝐮𝑓𝑤𝑑 ≈𝐍(𝐱)⋅𝐮𝑚 𝑖and 𝜕𝐮𝑓𝑤𝑑 𝜕𝐱≈𝜕𝐍 𝜕𝐱 ⋅𝐮𝑚 𝑖(4) where 𝐍(𝐱)denotes the matrix which contains standard shape functions associated to element interpolation, and the reconstructed displacement field 𝐮𝑚 𝑖is defined at nodal positions 𝑖. Then, using Eq. (4) into Eq. (3), the expression for the computation of the discrete inverse of the deformation gradient is given as, 𝐅−1 𝑗≈𝐈−𝜕𝐍 𝜕𝐱||||𝐱=𝐱𝑗 ⋅𝐮𝑚 𝑖(5) where derivatives are computed at the Gauss points 𝐱𝑗of the finite element mesh. Once obtained the deformation gradient, stresses are computed from this quantity by use of a certain constitutive relation. The reader is addressed to Appendix A for a detailed description of the nonlinear constitutive law employed in this study and the associated mechanical variables. Finally, tractions are computed following Cauchy’s formula as, 𝐭𝑛 𝑖=𝝈𝑖⋅𝐧𝑖(6) in which 𝐭𝑛 𝑖is the traction vector associated to the spatial configuration at nodal points 𝑖of the finite element mesh, 𝝈𝑖is the Cauchy stress tensor at nodal point 𝑖, once an averaging procedure of 𝝈𝑗from Gauss points is conducted. Finally, 𝐧𝑖denotes the outward normal vector to node 𝑖, which is numerically computed at nodal points from the surfaces of elements of the finite element mesh in the deformed (reference) configuration. 2.2. Inverse method The inverse method used in this work, PBNIM (physics-based nonlinear inverse method) (Sanz-Herrera et al.,2021), focuses on finding a new displacement field 𝐮𝑖𝑛𝑣 that approximates as closely as possible the measured one 𝐮𝑚and that fulfills the equilibrium International Journal of Engineering Science 186 (2023) 103828 4 A. Apolinar-Fernández et al. Fig. 1. Tree diagram showing the different series of cases considered in the study. condition within the medium of interest, i.e., the hydrogel in which the cell is embedded. This can be modeled as a minimization problem that includes the mechanical equilibrium constraint equation, and is mathematically described as follows, min 𝐮𝑖𝑛𝑣 (1 2‖‖𝐮𝑖𝑛𝑣 −𝐮𝑚‖‖2 2) s.t. 𝛩= 0 (7) where 𝛩generically represents the equilibrium constraint manifold where the candidate solution must lie. For the sake of compatibility with its finite element formulation, the functional 𝛩is developed in terms of the Principle of Virtual Work (PVW) equation. The reader is referred to Appendix A where this elaboration is explained in more detail in the context of its integration with ABAQUS software (used in the numerical implementation) and the chosen constitutive law. The equilibrium constraint can be integrated into the cost function (7) by means of a continuous and scalar Lagrange multiplier 𝜂as follows, min 𝐮𝑖𝑛𝑣 (1 2‖‖𝐮𝑖𝑛𝑣 −𝐮𝑚‖‖2 2+𝛩⋅𝜂)(8) Using the Gateaux derivative, (8) has its analytical minimal stationary solution at, 𝛿𝐮𝑖𝑛𝑣 =𝟎→𝐮𝑖𝑛𝑣 +𝛿𝛩 𝛿𝐮𝑖𝑛𝑣 ⋅𝜂=𝐮𝑚(9a) 𝛿𝜂 = 0 →𝛩= 0 (9b) Eqs. (9a) and (9b) are then elaborated, linearized and discretized in a FE framework. The reader is referred to Sanz-Herrera et al. (2021) for further information and details about PBNIM numerical implementation. 3. In silico TFM simulations The study presented in this article follows a holistic approach, which includes a compilation of cases up to a total of 744 in silico simulations. They are summarized in Fig. 1 and are explained next. 3.1. Ground truth cases First, a number of displacement fields were synthetically reproduced for different cell morphologies, levels of cellular traction magnitudes and different assumed mechanical behaviors of the hydrogels. These displacement fields were directly computed from FE simulations using the software ABAQUS. Cell morphologies We selected four different single cell morphologies obtained from real confocal microscopy images (voxel size: 0.57 ×0.57 ×1 μm3). Each cell displayed a different degree of morphological complexity. They were classified into 4 categories: (i) Spherical cells, with similar longitudinal and transversal lengths and few protrusions; (ii) Spread cells, whose longitudinal dimension is notably International Journal of Engineering Science 186 (2023) 103828 5 A. Apolinar-Fernández et al. Fig. 2. In silico models with selected cell morphologies: (a) star-shaped cell, (b) spread cell, (c) protrusive cell, and (d) spherical cell. The models illustrate the 3D cell embedded in the hydrogel (left part of each subfigure) and the body of the cells including location and direction of prescribed pulling forces (right part of each subfigure). longer than its transversal one; (iii) Protrusive cells, which presents an elongated morphology with protrusions; and (iv) Star-shaped cells, with no apparent dominant dimension in length, but with a large number of protrusions. Cell morphologies, as well as in silico models, domains of analysis, and loads are shown in Fig. 2. Finite element meshes were created with the Matlab toolbox TFMLAB (Barrasa-Fano, Shapeti, Jorge-Peñas, et al.,2021). Moreover, the solidity index, defined as the ratio between the cell volume and the volume of its convex hull, is used as an indicator of cell morphology complexity. Higher solidity indicates less complex (closer to a spherical shape) cell geometries. The solidity value of each cell is plotted in Fig. 3. Prescribed cellular tractions Cellular tractions are prescribed as distributed and uniform nodal forces along the tips of selected protrusions of the different analyzed cells. The direction of the forces is defined along the main axis of the protrusion, and inwardly towards the center of the cell. Fig. 2 shows the regions of application and direction of nodal forces for the different cells. Three different cellular pulling force magnitudes (small, intermediate and large) were selected in order to investigate their impact on traction reconstruction accuracy. Large cellular pulling force magnitudes were defined such that their associated maximum strains (Frobenius norm of the logarithmic strain tensor, i.e., the Euclidean norm of the logarithmic principal strains vector) ranged from 48% to 54% (see Table 1, nonlinear case) using the nonlinear model defined in Appendix A. In Table 1, the values for the large cellular pulling force magnitude prescribed for each cell lie within the range of nN found in the literature (Lekka et al.,2021). Moreover, the resulting displacements are of the order of microns, as typically seen in the laboratory in TFM (Vaeyens et al.,2020). Small force magnitudes were defined as 1/10 of the large force magnitude. Intermediate force magnitudes were defined as the mean value of the large and small force magnitudes (55% of large force). Table 1 also shows maximum strains and averaged maximum strains achieved in the different protrusions of the cells. Fig. 4 shows the ground truth strain measure indicator, i.e. the Frobenius norm of the logarithmic strain tensor on the boundary of the selected cells, for the different cases of considered cellular pulling force levels (see Table 1), assuming material nonlinearity for the hydrogel. It can be observed that the highest levels of strain are concentrated along the protrusion regions of the cells where forces are exerted (see Fig. 2). Note that the range of values shown in Fig. 4 has been chosen to provide a proper visual representation of the strain contours over the cells’ surfaces, and they do not correspond to the maximum values. Hydrogel’s constitutive behavior We selected a fibered nonlinear model to describe the nonlinear mechanical behavior of hydrogels. This model assumes a purely hyperelastic behavior and has been proven as an excellent fitting for collagen matrices in shear rheology tests (see Steinwachs International Journal of Engineering Science 186 (2023) 103828 6 A. Apolinar-Fernández et al. Fig. 3. Solidity index for selected cell morphologies. Table 1 Overall Frobenius norm of the logarithmic strain tensor (ground truth solution) for the nonlinear hydrogel case, for each protrusion of the selected cells for analysis. The values correspond to the case with the highest level of cellular pulling force magnitude (third column). Protrusion # Force [pN] Average Maximum Mean average Mean max. Maximum max. Spread 1 2.25e3 0.1870 0.5631 0.1592 0.4884 0.65312 0.1215 0.2490 3 0.1690 0.6531 Protrusive 1 0.375e3 0.1218 0.1838 0.1064 0.2412 0.6476 2 0.0795 0.1048 3 0.1103 0.1803 4 0.1125 0.2103 5 0.0814 0.1187 6 0.1331 0.6476 Spherical 10.75e30.2824 0.5607 0.2158 0.3781 0.6570 2 0.1491 0.1955 Star-shaped 1 0.67e3 0.1303 0.2065 0.1334 0.2602 0.4863 2 0.1457 0.2099 3 0.1390 0.4863 4 0.1086 0.1638 5 0.1660 0.3173 6 0.1109 0.1772 et al. (2016) and Appendix B). While previous studies have addressed the importance of viscoelasticity in cell behavior (Chaudhuri et al.,2020,2016;Vining & Mooney,2017) and implemented viscoelastic models for TFM (Toyjanova, Hannen, et al.,2014), our model neglects viscoelastic effects, as it was verified in previous studies for similar collagen hydrogels to the ones used in this paper (Steinwachs et al.,2016), that the material presents mainly an elastic response with negligible viscoelastic effects. For detailed information about the nonlinear model and its computational implementation, the reader is referred to Appendices A and B. Collagen matrices of 1.2 mg/ml concentration were prepared and tested in a shear rheology device following the protocol described in Cóndor et al. (2017). Fig. 5 shows the experimental results of shear rheology tests, as well as the best fit for the selected nonlinear model (see Appendix B). It can be observed, that the model accurately represents the real behavior of the hydrogel up to 50% of shear strain. Table 2 contains the values of the fitted nonlinear model parameters, as well as the parameters fitted for a linear isotropic model. In the latter case, Young’s modulus was computed from the combination of Poisson’s ratio and the shear International Journal of Engineering Science 186 (2023) 103828 7 A. Apolinar-Fernández et al. Fig. 4. Frobenius norm of the ground truth logarithmic strain tensor [–] on the boundary of the selected cells assuming a nonlinear matrix. (a) small cellular pulling force case; (b) intermediate cellular pulling force case; (c) large cellular pulling force case. modulus through the expression 2𝐺(1 + 𝜈). In the literature, one can find Poisson’s ratio values for collagen hydrogels ranging from 0.2 to 0.48 (see Bloom et al. (2008), Bowers et al. (2020), Du et al. (2018), Koch et al. (2012), Wang et al. (2014) and Gjorevski et al. (2015)). For this study, 0.34, a value within the range found in the literature, was chosen. The shear modulus was set as the value of the initial slope in the shear rheology experimental curve. International Journal of Engineering Science 186 (2023) 103828 8 A. Apolinar-Fernández et al. Fig. 5. Fitting of experimental curve of 1.2 mg/ml collagen hydrogel with linear elastic and selected nonlinear elastic model for shear rheology tests. The shaded blue area represents the variability of the different tested specimens. The solid blue line curve is the average of all specimens. Table 2 Fitted nonlinear and linear model parameters. Parameters 𝐸and 𝜈of the nonlinear model stand for the Young’s modulus and Poisson’s ratio, respectively (see Appendix A for the definition of the nonlinear model parameters). Model Parameter Value Nonlinear 𝜅𝑜[Pa] 156.6217 𝑑𝑜[–] 0.375 𝜆𝑠[–] 0.0111 𝑑𝑠[–] 0.0836 Linear 𝐸[Pa] 29.5076 𝜈[–] 0.340 The assumed linear and nonlinear models in the different cases of analysis consider the linear and nonlinear fittings shown in Fig. 5. On the one hand, the case of the nonlinear model uses the nonlinear versions of the forward and inverse formulations (finite strain hyperelasticity), as introduced in Section 2.1,2.2,Appendix A and Sanz-Herrera et al. (2021). On the other hand, the case of the linear model uses the linearized (infinitesimal strains) versions of the forward and inverse formulations (Barrasa-Fano, Shapeti, de Jong, et al.,2021). For the different conditions explained above, a total of 24 ground truth displacement fields (4 cells ×3 cellular force levels ×2 material behaviors) were reproduced, see Fig. 1, using the software ABAQUS. These displacement fields, and associated strain and traction variables, are referred to as the ground truth (GT) solutions. 3.2. Displacement reconstruction In a similar way to which it is done in the TFM experimental methodology, GT hydrogel displacements were sampled for each case at discrete random (bead) positions. For all the cases of analysis, we assumed an experimentally feasible bead density of approximately 0.03 beads∕μm3, which relies within the order of magnitude of experimental procedures (Barrasa-Fano, Shapeti, de Jong, et al.,2021). Table 3 shows the dimensions of the hydrogel domain, the number of beads included in the domain, and the inter-bead distance for each case of analysis. The reconstructed displacement field in the FE mesh is obtained then as a Lagrangian interpolation from bead to nodal positions, being the displacements at bead locations previously obtained by interpolation of ground truth displacements present at the nodes in the discretization, whose position does not generally coincide with that of the beads. Fig. 6 shows a schematic representation of this process. We refer to the reconstructed displacement field as the corrupted input displacement field for forward and inverse methods simulations. International Journal of Engineering Science 186 (2023) 103828 9 A. Apolinar-Fernández et al. Table 3 Size of the hydrogel domain, number of beads (𝑁beads), inter-bead distance for each cell, and size of the FE mesh. The inter-bead distance was calculated as the mean nearest neighbor Euclidean distance. 𝐿𝑋[μm] 𝐿𝑌[μm] 𝐿𝑍[μm] 𝑁beads Interbead distance [μm] FE mesh (#elements) Star-shaped 117.99 117.99 66 28280 1.799 163097 Spread 148.77 148.77 73 49728 1.796 224386 Protrusive 119.13 119.13 68 29703 1.802 147118 Spherical 92.91 92.91 66 17535 1.799 79034 Fig. 6. Schematic representation of the process to obtain corrupted displacement fields from ground truth simulations. Errors between ground truth and corrupted displacements were amplified for representation purposes. With the aim of giving statistical validity to the results, 10 realizations of bead positions for each case were performed corresponding to different instances of the reconstructed displacement fields obtained from the GT (reference) solutions. 3.3. Forward and inverse simulations The different corrupted displacement fields were used as input for the forward and inverse simulations for traction recovery. A total of 120 simulations (4 cells ×3 cellular force levels ×10 realizations) were run for the forward and inverse methods (240 simulations in total) using the linear material behavior, and the recovered displacement field of the linear material case as input. Moreover, another 240 forward/inverse simulations were carried out using the nonlinear material behavior, and the recovered displacement field of the nonlinear material case as input. And finally, a total of 240 forward/inverse simulations using the linear material behavior, and the recovered displacement field of the nonlinear material case as input, were performed. The latter serves to analyze the feasibility of assuming a linear elastic approach (or equivalently neglecting nonlinear effects) for traction reconstruction in nonlinear matrices. The simulations were performed using ABAQUS in an in-house code orchestrated by Matlab. The size of the different meshes used in the simulations is given in Table 3. The convergence of the results with the selected meshes was previously checked. 3.4. Error indicators In this paper we focus on the errors of traction recovery (either by means of the forward or inverse methods), as they are the most important variable of analysis in TFM. To this end, we define the forward and inverse traction error indicators 𝑇𝑓𝑤𝑑 𝑒,𝑇𝑖𝑛𝑣 𝑒as follows, 𝑇𝑖𝑛𝑣 𝑒= 100 ⋅(1 − corr[𝑡𝑖𝑛𝑣 pdf, 𝑡𝐺𝑇 pdf]) (10a) 𝑇𝑓𝑤𝑑 𝑒= 100 ⋅(1 − corr[𝑡𝑓𝑤𝑑 pdf , 𝑡𝐺𝑇 pdf]) (10b) where 𝑡𝑖𝑛𝑣 pdf and 𝑡𝑓𝑤𝑑 pdf are the probability density functions of the magnitudes of the tractions distribution 𝑡𝑛 𝑖(𝑖being the boundary nodes of the cell) for the inverse and forward solutions, respectively. 𝑡𝐺𝑇 pdf has the same definition but for the ground truth solution. Finally, corr[𝑓, 𝑔]is the correlation coefficients of functions 𝑓and 𝑔. International Journal of Engineering Science 186 (2023) 103828 16 A. Apolinar-Fernández et al. Fig. 12. CPU time versus mean traction error indicators (for analyzed cellular pulling force level cases) of recovered tractions solution: (a) Spread cell, (b) protrusive cell, (c) spherical cell, and (d) star-shaped cell. 5.2. Traction reconstruction accuracy referred to cellular morphology The impact of cellular morphology on the error of traction reconstruction can be seen in Fig. 11. In this figure, the traction error indicator (averaged for all the pulling force cases) is plotted versus the solidity index of the cell morphology. It is observed that there exists a strong correlation between the error and the solidity index, being more challenging to recover accurate tractions for star-shaped geometries (low solidity indices) than for spherical ones (high solidity indices), for both linear and nonlinear matrices. Indeed, complex geometries (low solidity indices) induce abrupt changes and gradients of displacements in the hydrogel domain near to the cell boundary. Barrasa-Fano, Shapeti, de Jong, et al. (2021) demonstrated that high gradients of the displacements solution (and hence referred to tractions) are difficult to capture in TFM, as seen here for increasing complexity of the cell morphology. Furthermore, Fig. 11 shows the error recovery performance as solidity index increases. Besides the superiority of the inverse method, it is shown that the error decreases with a similar trend (in a logarithmic scale) with solidity for the forward and inverse methods, for the case that the nonlinear hydrogel behavior is reconstructed assuming a nonlinear behavior of the hydrogel (Fig. 11a). On the other hand, the error decreases faster (in a logarithmic scale) with solidity for the inverse methods, for the rest of considered cases related to the assumed mechanical behavior of the hydrogel (Figs. 11b and 11c). 5.3. Traction reconstruction accuracy and computational efficiency referred to forward/inverse methodologies The outperformance of the inverse method on the accuracy of tractions reconstruction in TFM versus the forward method can be seen in Fig. 10 for all the analyzed cases. Also, these data are plotted in Fig. 11, averaged for all cellular pulling force levels, as it does not represent a significant effect, as discussed before. These figures show errors 3–5 times higher for the forward method compared to the inverse method. International Journal of Engineering Science 186 (2023) 103828 17 A. Apolinar-Fernández et al. Fig. 13. Ratio between reconstructed nonlinear and linear strain energy densities on the boundary of the selected cells assuming a nonlinear hyperelastic behavior for the matrix. (a) small cellular pulling force case; (b) intermediate cellular pulling force case; (c) large cellular pulling force case. On the contrary, the forward method is computationally less expensive than the inverse method as it includes direct and straightforward algebraic computations, provided the displacement field, either in a linear or nonlinear statement of the problem (Sanz-Herrera et al.,2021). The inverse method requires the setup of a matrix system, and an iterative process in a nonlinear International Journal of Engineering Science 186 (2023) 103828 18 A. Apolinar-Fernández et al. case (Sanz-Herrera et al.,2021). Table 4 and Fig. 12 show the CPU time of the forward and inverse methods, separated for each analyzed case and cell morphology (as it is related to the FE mesh size and hence CPU time). It is observed that higher CPU differences are always found between the forward and inverse for the nonlinear case, as it additionally requires of the mentioned iterative process. The highest CPU time differences are found for the most time consuming case (spread morphology) with a ratio of 233 s to 1849s from forward to inverse (for the case that a nonlinear behavior is reconstructed assuming a nonlinear behavior of the matrix), but with a higher traction error cost of 31.2 to 7.7% from forward to inverse, respectively (see Table 4 and Fig. 12). 5.4. Traction reconstruction accuracy and computational efficiency referred to hydrogel’s behavior Linear and nonlinear matrices were considered in the generation of ground truth synthetic solutions, and tractions were reconstructed for both cases considering their corresponding linear or nonlinear behavior. Figs. 10a and 10b show the traction error for nonlinear and linear cases, respectively. It is observed that the errors of the inverse method are similar for all the analyzed cases, whereas the forward provide slightly higher errors for the linear case, specially for the spherical cell morphology. Analogously, Zündel et al. (2017) found that although the use of a strain-stiffening material provides higher quality traction reconstructions, errors are of similar magnitude in both materials. In our case, we find an explanation for the last point in Fig. 13. One can observe that there are regions in which the linear strain energy density overcomes the nonlinear one (ratio below 1). This is due to the fact that the considered nonlinear model assumes fiber buckling and hence low stiffness, even negligible, under compression (see Appendix A), whereas the linear model present a constant stiffness under compressive/tensile stresses. Therefore it results in different values that contribute to the error metric, as Fig. 10b shows. These regions are particularly represented in the spherical cell, as illustrated in Fig. 13.Fig. 11 also shows the same trend of the error of the assumed models, but with a faster decrease of the error with solidity index for the linear model (Fig. 11b). With regards to the CPU time, the linear model runs faster than the nonlinear one for all the analyzed cases, according to Table 4 and Fig. 12. In particular, for the forward method, both the linear and nonlinear methods require only algebraic computations, although with a slightly higher CPU time cost for the nonlinear model since the evaluation of the assumed constitutive law requires a numerical integration over the unit sphere (Steinwachs et al.,2016). Nonetheless, CPU time ratios between the nonlinear and linear model approaches are lower than a 1.4 factor for all the analyzed cases ( Table 4 and Fig. 12). On the other hand, the inverse method requires an iterative procedure when the material model is assumed as nonlinear. Despite the few number of iterations until convergence for the inverse nonlinear approach (being 4.3 iterations in average for all the analyzed cellular pulling force levels for the worst case, according to Table 4), CPU time increases to ratios up to 3.3 (worst case), according to Table 4 and Fig. 12, as it increases proportionally to the number of iterations. Feasibility of linear TFM analysis in real nonlinear matrices The particular case that a nonlinear hydrogel is approached by a linear behavior for traction reconstruction, and hence neglecting nonlinear effects, was analyzed in order to provide an estimation of the feasibility of using linear modeling in TFM even though the real behavior of the matrix fits better with a nonlinear model. A linear model in TFM shows several advantages versus a nonlinear one, such as ease of numerical implementation avoiding convergence troubles, and hence no need for computational mechanics background of the user for troubleshooting; as well as better CPU time performance. Fig. 10c shows the traction error for this case, and Fig. 11c the error averaged for the considered pulling force levels. As expected, it is observed in Fig. 11c that the linear assumption for the hydrogel performs worse in traction reconstruction. Specifically, for the toughest case i.e. the starshaped morphology, the error increases from 36.8 to 74.7% for the forward, and from 9.3 to 30.5% in the inverse case; when neglecting nonlinear effects in a real nonlinear hydrogel (Figs. 11a and 11c). If we focus on the small cellular pulling force case, such that geometrical nonlinear effects may be neglected, the differences found in Figs. 10a and 10c are just a consequence of the differences of the material models used. Despite that linear and nonlinear material models could be assumed as the same for small strains according to Fig. 5, full 3D and multiaxial stress tensors (which is the most frequent stress state in the analyzed cases) show significant differences for the assumed linear and nonlinear models when computed from similar strain fields, even at small strains (specially at compression states as previously commented). On the other hand, CPU time is faster when assuming a linear model instead of a nonlinear one, as discussed before. 6. Conclusions TFM involves many multidisciplinary challenges. From the point of view of mathematical modeling, accurate material models and improved numerical methods are needed to reconstruct reliable tractions, to properly investigate the mechanobiology of in vitro models. In this sense, the impact of a number of issues of interest, such as cell morphology, cellular pulling force and assumed material model, on the traction error recovery; was investigated by means of a conventional forward approach and a version of the inverse method. The inverse method outperforms the forward for all the analyzed cases, reducing the traction error by a factor of 3–5 times. Given the order of magnitude of the traction error, the forward method is not a good approach for the analyzed cases except for spherical morphologies where acceptable errors were found for the forward (∼ 10% for the linear model and <2% for the nonlinear model). Regarding the impact of cell morphology on the performance of tractions recovery, it was found a strong correlation between the defined traction error indicator and the sphericity of the cell (solidity index). Moreover, the cellular pulling force level was proven International Journal of Engineering Science 186 (2023) 103828 19 A. Apolinar-Fernández et al. to have a residual impact on the accuracy of traction reconstruction, at least when the error indicator is defined in global terms. The accuracy of traction recovery is similar for linear and nonlinear matrices for the inverse method, but slightly more challenging for linear matrices for the forward method. Despite the advantages that a linear model offers, with better CPU time efficiency amongst other features, it is not a good assumption to represent the behavior of a nonlinear matrix, like collagen hydrogel as used in this study. Even for the inverse method, traction errors range from 15%–30% in this scenario. The knowledge acquired in the analysis performed in this study sheds light on TFM performance in a number of cases of real interest. These situations have not been sufficiently analyzed in previous works. The obtained conclusions may be useful to properly design TFM in vitro experiments, and to further investigate new methods and models in this field. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Data availability Data will be made available on request. Acknowledgments A.A.-F. and J.A.S.-H. gratefully acknowledge the Financial support from project PID2021-126051OB-C42 by the Ministerio de Ciencia e Innovación (MCI), Agencia Estatal de Investigación (AEI) and Fondo Europeo de Desarrollo Regional (FEDER); and project P20-01195 funded by the Consejería de Economía, Conocimiento, Empresas y Universidad de la Junta de Andalucía. H.V.O., J.B.- F. and M.C. gratefully acknowledge the financial support from KU Leuven Internal Funding (project C14/17/111), FWO (project G087018N, infrastructure grants G0H6316N and I009718N, postdoctoral fellowship 12ZR120N of MC) and Special Research Fund iBOF/21/083. Appendix A The aim of this appendix is to present a thorough description of the nonlinear model employed in this study (Steinwachs et al., 2016), as well as the derivation process of the corresponding fundamental mechanical relations, to be numerically implemented in the ABAQUS subroutine UMAT, according to the numerical scheme used in this paper. A brief explanation about the Finite Element Method formulation in ABAQUS is given first, in order to provide some context about the variables to be defined in the UMAT subroutine code, to implement a certain constitutive law. A.1. ABAQUS FEM theoretical background The FE analysis software ABAQUS offers the possibility of implementing any desired constitutive law into its computational process by means of the user subroutine UMAT. This subroutine is coded using the FORTRAN language. Finite element formulation The equilibrium of a deformable body can be described through the fundamental scalar expression resulting from the application of the Principle of Virtual Work in spatial configuration, 𝛿𝑊 ∶= ∫𝑣 𝝈∶𝛿𝜺𝑑𝑣 −∫𝑣 𝜌𝒃⋅𝛿𝒖𝑑𝑣 −∫𝜕𝑣 𝒕⋅𝛿𝒖𝑑𝑎 = 0 (A.1) where 𝛿𝐮denotes the virtual displacement, 𝑣the current volume, 𝜌the current density, 𝐛the volumetric forces vector, 𝛿𝝐the linear virtual strain tensor, 𝐭the traction vector resulting from the product between the Cauchy stress tensor 𝝈and the vector normal to the surface 𝜕𝑣,𝐧. A linearized form of Eq. (A.1), in the context of a Newton–Raphson procedure, can be obtained by performing the directional derivative of the virtual work in the direction of the increment 𝛥𝐮. Being 𝜙𝑘a test function, 𝛿𝑊 (𝜙𝑘, 𝛿𝒖) + 𝐷𝛿𝑊 (𝜙𝑘, 𝛿𝒖)[𝛥𝒖]=0 (A.2) ABAQUS employs the Updated Lagrangian configuration, by means of which the Principle of Virtual Work is expressed as, 𝛿𝑊 ∶= ∫𝑉 𝝉∶𝛿𝜺𝑑𝑉 −∫𝑉 𝜌𝑜𝒃⋅𝛿𝒖𝑑𝑉 −∫𝜕𝑣 𝒕⋅𝛿𝒖𝑑𝑎 = 0, 𝛿𝜺=1 2[∇𝛿𝒖+[∇𝛿𝒖]𝑇](A.3) where 𝝉denotes the Kirchhoff stress tensor. Note that the volume integrals have been now expressed in the material configuration. At the time of linearizing Eq. (A.3), ABAQUS approximates the term involving the work of internal forces by means of the Jaumann derivative or co-rotational derivative, which represents an objective quantity, 𝐷[∫𝑉 𝝉∶𝛿𝜺𝑑𝑉 ][𝛥𝒖] = ∫𝑉 𝐷[𝝉∶𝛿𝜺][𝛥𝒖]𝑑𝑉 = International Journal of Engineering Science 186 (2023) 103828 20 A. Apolinar-Fernández et al. =∫𝑉[∇𝛿𝒖∶[𝐽c𝐽−H]∶ ∇𝛿𝒖+𝝉∶[[∇𝛿𝒖]𝑇∇𝛥𝒖]]𝑑𝑉 (A.4) where c𝐽is the Jaumann elasticity tensor, which can be computed as, 𝐽c𝐽=𝐽c+H, 𝐻𝑖𝑗𝑘𝑙 =1 2[𝛿𝑗𝑘𝜏𝑖𝑙 +𝛿𝑖𝑘𝜏𝑗𝑙 +𝛿𝑗𝑙𝜏𝑖𝑘 +𝛿𝑖𝑙 𝜏𝑗𝑘](A.5) being cthe so called Fourth Cauchy Elasticity Tensor, c=1 𝐽𝜙∗[], 𝑐𝑖𝑗𝑘𝑙 =1 𝐽𝐹𝑖𝐼 𝐹𝑗𝐽 𝐹𝑘𝐾 𝐹𝑙𝐿𝐼𝐽 𝐾𝐿 (A.6) where Einstein’s summation convention has been used. In order for ABAQUS to accept a certain constitutive relation via the UMAT subroutine, 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 ABAQUS employs, 𝝈=⎡⎢⎢⎢⎢⎢⎢⎣ 𝜎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 SAEN model (Steinwachs et al.,2016) propose a model for biopolymer networks based on the combination of a micromechanical basis, i.e., focused on the individual behavior of an individual fiber; and a continuum statement. For a given sufficiently small scale corresponding to a fiber segment, deformations are considered non-affine, this is, independent of deformations of the bulk medium. However, for scales larger tan the typical interconnection distance between fibers, fiber deformations approximate macroscopic deformations 𝜆, which depends on fiber orientations and the magnitude of the considered deformations. In this way, deformations are considered to be affine for sufficiently large volumes of the material. The model considers that one fiber can exist in three different states of deformation: compression, in which the fiber becomes unstable, i.e., it buckles, and its stiffness decays exponentially towards zero; straightening, intermediate process in which the fiber recovers its native length after buckling, with constant stiffness; and stretching, a process by which the fiber deforms axially beyond its native length and it experiments an exponential stiffening. The strain energy density function is defined from the stiffness, this is, through its second derivative with respect to the macroscopic deformation measure 𝜆, 𝛹′′(𝜆) = 𝜅𝑜⎧ ⎪ ⎨ ⎪ ⎩ exp [𝜆∕𝑑𝑜] ∀𝜆 < 0 1 ∀0 ≤𝜆<𝜆𝑠 exp [[𝜆−𝜆𝑠]∕𝑑𝑠]∀𝜆≥𝜆𝑠 𝜆=‖𝑭 𝒆𝛺‖− 1 (A.8) where 𝜅𝑜is a material parameter with stress dimensions, 𝑑𝑜is the buckling dimensionless parameter, 𝑑𝑠is the strain-stiffening dimensionless parameter and 𝜆𝑠is the dimensionless parameter that defines the extension of the interval of the linear behavior, i.e., with constant stiffness. The macroscopic deformation 𝜆is defined as the variation in length of a fiber with respect its native length, 𝐞𝛺is the unit vector that represents the fiber’s orientation and 𝐅is the deformation gradient. The mechanical stress of the fiber can be defined by the expression: 𝛹′(𝜆) = ∫𝜆 0 𝛹′′(𝜆)𝑑𝜆 =𝜅𝑜⎧ ⎪ ⎨ ⎪ ⎩ 𝑑𝑜⋅[exp [𝜆∕𝑑𝑜]−1]∀𝜆 < 0 𝜆∀0 ≤𝜆<𝜆𝑠 𝜆𝑠+𝑑𝑠⋅[exp [[𝜆−𝜆𝑠]∕𝑑𝑠]− 1]∀𝜆≥𝜆𝑠 (A.9) in which integration from zero indicates that the material is not prestressed. The constitutive relation in mixed configuration can then be obtained as, 𝑷=⟨𝜕𝛹(𝜆(𝑭)) 𝜕𝑭⟩𝜴 =⟨𝜕𝛹 𝜕𝜆 𝜕𝜆 𝜕𝑭⟩𝜴 =⟨𝛹′(𝜆)[𝑭 𝒆𝛺]⊗𝒆𝛺 ‖𝑭 𝒆𝛺‖⟩𝜴 (A.10) where 𝑷denotes the First Piola–Kirchhoff stress tensor and brackets indicate averaging over all directions 𝛺in the unit sphere. From Eq. (A.10), the Cauchy stress tensor is obtained by applying a push-forward operation as follows, 𝝈=1 𝐽𝜙∗[𝑷] = 1 𝐽 𝑷 𝑭 𝑇=1 𝐽[1 4𝜋∫2𝜋 0∫𝜋 0 𝛹′(𝜆)[𝑭 𝒆𝛺]⊗𝒆𝛺 ‖𝑭 𝒆𝛺‖sin 𝜃𝑑𝜃𝑑𝜙]⋅𝑭𝑇(A.11) where the integral over the unit sphere is approximated by means of the Gauss–Legendre quadrature. International Journal of Engineering Science 186 (2023) 103828 21 A. Apolinar-Fernández et al. 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.12) association which permits the rewriting of 𝜆in terms of the Right Cauchy–Green deformation tensor, 𝜆(𝑭) = ‖𝑭 𝒆𝛺‖− 1 ⟶𝜆(𝑪)=(𝐼4)1∕2 − 1 (A.13) In this way, the constitutive relation in material configuration is obtained as, 𝑺= 2 ⟨𝜕𝛹(𝜆(𝑪)) 𝜕𝑪⟩𝜴 = 2 ⟨𝜕𝛹 𝜕𝜆 𝜕𝜆 𝜕𝑪⟩𝜴 =⟨[𝜆+ 1]−1 ⋅𝛹′(𝜆)⋅[𝒆𝛺⊗𝒆𝛺]⟩𝜴(A.14) where ⊗denotes tensor product. The Second elasticity tensor can then be retrieved from, = 2 ⟨𝜕𝑺∗(𝜆(𝑪)) 𝜕𝑪⟩𝜴 = 2 ⟨𝜕𝑺∗ 𝜕𝜆 ⊗𝜕𝜆 𝜕𝑪⟩𝜴 (A.15) in which 𝑺∗denotes the Second Piola–Kirchhoff stress tensor of a single fiber (non-homogenized) and, 𝜕𝑺∗ 𝜕𝜆 =[−[𝜆+ 1]−2 ⋅𝛹′(𝜆)+[𝜆+ 1]−1 ⋅𝛹′′(𝜆)]⋅[𝒆𝛺⊗𝒆𝛺](A.16) 𝜕𝜆 𝜕𝑪=1 2[𝜆+ 1]−1 ⋅[𝒆𝛺⊗𝒆𝛺](A.17) Appendix B We fitted the parameters that characterize the SAEN model to experimental rheological data corresponding to a simple shear experiment performed on collagen hydrogels, in order to verify that the model is able to reproduce the real material behavior. B.1. Shear rheology The hydrogel undergoes the action of a rotational rheometer of cone and plate with the aim of measuring the stress–strain relationship within a simple shear state of deformation. The deformation gradient in the case of a simple shear state of deformation is written as, 𝑭=⎡⎢⎢⎣ 1𝛾0 010 001⎤⎥⎥⎦ (B.1) where 𝛾represents the engineering shear strain. The shear stress in the test is computed from the measured torque and the setup of the test in the nominal configuration. Therefore, the shear stress resulting from the test corresponds to the shear component of the First Piola–Kirchhoff stress tensor (nominal stress). Recalling the definition given in Eq. (A.10) for this tensor, the nominal stress is then computed from the following expression, 𝑃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. 5 (main text) shows the fitting performed to the experimental curve corresponding to a hydrogel with a fiber concentration of 1.2 mg/ml. International Journal of Engineering Science 186 (2023) 103828 22 A. Apolinar-Fernández et al. References 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 Biomaterialia,126, 326–338. http: //dx.doi.org/10.1016/j.actbio.2021.03.014. Barrasa-Fano, J., Shapeti, A., Jorge-Peñas, Á., Barzegari, M., Sanz-Herrera, J. A., & Van Oosterwyck, H. (2021). TFMLAB: A MATLAB toolbox for 4D traction force microscopy. SoftwareX,15, Article 100723. http://dx.doi.org/10.1016/j.softx.2021.100723. 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. Biophysical Journal,95(8), 4077–4088. http://dx.doi.org/10.1529/biophysj.108.132738. Bowers, H. J., Fannin, E. E., Thomas, A., & Weis, J. A. (2020). Characterization of multicellular breast tumor spheroids using image data-driven biophysical mathematical modeling. Scientific Reports,10(1), 11583. http://dx.doi.org/10.1038/s41598-020-68324-4. Butler, J. P., Tolić-Nørrelykke, I. M., Fabry, B., & Fredberg, J. J. (2002). Traction fields, moments, and strain energy that cells exert on their surroundings. American Journal of Physiology - Cell Physiology,282(3), URL: http://ajpcell.physiology.org/content/282/3/C595. Caliari, S. R., & Burdick, J. A. (2016). A practical guide to hydrogels for cell culture. Nature Methods,13(5), 405–414. http://dx.doi.org/10.1038/nmeth.3839, URL: http://www.nature.com/articles/nmeth.3839. Chaudhuri, O., Cooper-White, J., Janmey, P. A., Mooney, D. J., & Shenoy, V. B. (2020). Effects of extracellular matrix viscoelasticity on cellular behaviour. Nature,584(7822), http://dx.doi.org/10.1038/nmat4489, URL: https://www.nature.com/articles/s41586-020-2612-2. Chaudhuri, O., Gu, L., Klumpers, D., Darnell, M., Bencherif, S. A., Weaver, J. C., Huebsch, N., Lee, H., Lippens, E., Duda, G. N., & Mooney, D. J. (2016). Hydrogels with tunable stress relaxation regulate stem cell fate and activity. Nature Materials,15(3), http://dx.doi.org/10.1038/nmat4489, URL: https://www.nature.com/articles/nmat4489. Cóndor, M., Steinwachs, J., Mark, C., García-Aznar, J., & Fabry, B. (2017). Traction force microscopy in 3-dimensional extracellular matrix networks. Current Protocols in Cell Biology,75(1), 10.22.1–10.22.20. http://dx.doi.org/10.1002/cpcb.24, URL: https://currentprotocols.onlinelibrary.wiley.com/doi/abs/10.1002/ cpcb.24,arXiv:https://currentprotocols.onlinelibrary.wiley.com/doi/pdf/10.1002/cpcb.24. Dembo, M., & Wang, Y. L. (1999). Stresses at the cell-to-substrate interface during locomotion of fibroblasts. Biophysical Journal,76(4), 2307–2316. http: //dx.doi.org/10.1016/S0006-3495(99)77386-8, URL: http://www.cell.com/biophysj/pdf/S0006-3495(99)77386-8.pdf. Discher, D. E., Janmey, P., & Wang, Y. L. (2005). Tissue cells feel and respond to the stiffness of their substrate. Science,310(5751), 1139–1143. http: //dx.doi.org/10.1126/science.1116995, URL: http://www.ncbi.nlm.nih.gov/pubmed/16293750. 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 Biomaterialia,67, 215–228. http://dx.doi.org/10.1016/j.actbio.2017.12.002, URL: https://www.sciencedirect.com/science/article/pii/S1742706117307584#s0065. Duval, K., Grover, H., Han, L. H., Mou, Y., Pegoraro, A. F., Fredberg, J., & Chen, Z. (2017). Modeling physiological events in 2D vs. 3D cell culture. Physiology, 32(4), 266–277. http://dx.doi.org/10.1152/physiol.00036.2016, URL: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5545611/. Elliott, H., Fischer, R. S., Myers, K. A., Desai, R. A., Gao, L., Chen, C. S., Adelstein, R. S., Waterman, C. M., & Danuser, G. (2015). Myosin II controls cellular branching morphogenesis and migration in three dimensions by minimizing cell-surface curvature. Nature Cell Biology 2014,17(2), 137–147. http://dx.doi.org/10.1038/ncb3092, URL: https://www.nature.com/articles/ncb3092. Engler, A. J., Sen, S., Sweeney, H. L., & Discher, D. E. (2006). Matrix elasticity directs stem cell lineage specification. Cell,126(4), 677–689. http://dx.doi.org/ 10.1016/J.CELL.2006.06.044, URL: https://www.sciencedirect.com/science/article/pii/S0092867406009615. Feng, X., & Hui, C. Y. (2016). Force sensing using 3D displacement measurements in linear elastic bodies. Computational Mechanics,58(1), 91–105. http: //dx.doi.org/10.1007/s00466-016-1283-1, URL: https://link.springer.com/article/10.1007/s00466-016-1283-1. Franck, C., Maskarinec, S. A., Tirrell, D. A., Ravichandran, G., & Genin, G. (2011). Three-dimensional traction force microscopy: A new tool for quantifying cell-matrix interactions. PLoS One,6(3), Article e17833. http://dx.doi.org/10.1371/journal.pone.0017833, URL: http://dx.plos.org/10.1371/journal.pone. 0017833. Gjorevski, N., Piotrowski, A. S., Varner, V. D., & Nelson, C. M. (2015). Dynamic tensile forces drive collective cell migration through three-dimensional extracellular matrices. Scientific Reports,5, 11458. http://dx.doi.org/10.1038/srep11458. Guilak, F., Cohen, D. M., Estes, B. T., Gimble, J. M., Liedtke, W., & Chen, C. S. (2009). Control of stem cell fate by physical interactions with the extracellular matrix. Cell Stem Cell,5(1), 17–26. http://dx.doi.org/10.1016/J.STEM.2009.06.016, URL: https://www.cell.com/cell-stem-cell/fulltext/S1934- 5909%2809%2900293-8. Harris, A., Wild, P., & Stopak, D. (1980). Silicone rubber substrata: a new wrinkle in the study of cell locomotion. Science,208(4440), 177–179. http: //dx.doi.org/10.1126/science.6987736, URL: http://www.sciencemag.org/cgi/doi/10.1126/science.6987736. Hervas-Raluy, S., Gomez-Benito, M. J., Borau-Zamora, C., Cóndor, M., & Garcia-Aznar, J. M. (2021). A new 3D finite element-based approach for computing cell surface tractions assuming nonlinear conditions. PLoS One,16(4 April), Article e0249018. http://dx.doi.org/10.1371/journal.pone.0249018, URL: https://dx.plos.org/10.1371/journal.pone.0249018. Holenstein, C. N., Lendi, C. R., Wili, N., & Snedeker, J. G. (2019). Simulation and evaluation of 3D traction force microscopy. Computer Methods in Biomechanics and Biomedical Engineering,22(8), 853–860. http://dx.doi.org/10.1080/10255842.2019.1599866, URL: https://www.tandfonline.com/doi/full/10.1080/10255842. 2019.1599866. Huang, Y., Gompper, G., & Sabass, B. (2020). A Bayesian traction force microscopy method with automated denoising in a user-friendly software package. Computer Physics Communications,256, Article 107313. http://dx.doi.org/10.1016/j.cpc.2020.107313,arXiv:2005.01377. Huebsch, N., Arany, P. R., Mao, A. S., Shvartsman, D., Ali, O. A., Bencherif, S. A., Rivera-Feliciano, J., & Mooney, D. J. (2010). Harnessing tractionmediated manipulation of the cell/matrix interface to control stem-cell fate. Nature Materials,9(6), 518–526. http://dx.doi.org/10.1038/nmat2732, URL: http://www.nature.com/articles/nmat2732. Ingber, D. E. (2002). Mechanical signaling and the cellular response to extracellular matrix in angiogenesis and cardiovascular physiology. Circulation Research, 91(10), 877–887. http://dx.doi.org/10.1161/01.RES.0000039537.73816.E5, URL: http://www.circresaha.org. Ingber, D. E. (2003). Mechanobiology and diseases of mechanotransduction. Annals of Medicine,35(8), 564–577. http://dx.doi.org/10.1080/07853890310016333, URL: http://www.tandfonline.com/doi/full/10.1080/07853890310016333. 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. Annals of Biomedical Engineering,47(2), 624–637. http://dx.doi.org/10.1007/s10439-018-02164-2, URL: http://link.springer.com/10.1007/s10439-018-02164-2. Janmey, P. A., Fletcher, D. A., & Reinhart-King, C. A. (2020). Stiffness sensing by cells. Physiological Reviews,100(2), 695–724. http://dx.doi.org/10.1152/physrev. 00013.2019. Jansen, K. A., Donato, D. M., Balcioglu, H. E., Schmidt, T., Danen, E. H., & Koenderink, G. H. (2015). A guide to mechanobiology: Where biology and physics meet. Biochimica et Biophysica Acta - Molecular Cell Research,1853(11), 3043–3052. http://dx.doi.org/10.1016/j.bbamcr.2015.05.007. Jiang, T., Zhao, J., Yu, S., Mao, Z., Gao, C., Zhu, Y., Mao, C., & Zheng, L. (2019). Untangling the response of bone tumor cells and bone forming cells to matrix stiffness and adhesion ligand density by means of hydrogels. Biomaterials,188, 130–143. http://dx.doi.org/10.1016/j.biomaterials.2018.10.015. International Journal of Engineering Science 186 (2023) 103828 23 A. Apolinar-Fernández et al. Kim, S., Uroz, M., Bays, J. L., & Chen, C. S. (2021). Harnessing mechanobiology for tissue engineering. Developmental Cell,56(2), 180–191. http://dx.doi.org/ 10.1016/j.devcel.2020.12.017, URL: https://linkinghub.elsevier.com/retrieve/pii/S1534580720310236. Koch, T. M., Münster, S., Bonakdar, N., Butler, J. P., & Fabry, B. (2012). 3D traction forces in cancer cell invasion. PLoS One,7(3), Article e33476. http://dx.doi.org/10.1371/journal.pone.0033476, URL: http://dx.plos.org/10.1371/journal.pone.0033476. Kopanska, K. S., Alcheikh, Y., Staneva, R., Vignjevic, D., & Betz, T. (2016). Tensile forces originating from cancer spheroids facilitate tumor invasion. PLoS One, 11(6), Article e0156442. http://dx.doi.org/10.1371/journal.pone.0156442, URL: https://dx.plos.org/10.1371/journal.pone.0156442. Kumar, S., & Weaver, V. M. (2009). Mechanics, malignancy, and metastasis: The force journey of a tumor cell. Cancer Metastasis Reviews,28(1–2), 113. http://dx.doi.org/10.1007/S10555-008-9173-4, URL: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2658728/. Kutys, M. L., & Chen, C. S. (2016). Forces and mechanotransduction in 3D vascular biology. Current Opinion in Cell Biology,42, 73–79. http://dx.doi.org/10. 1016/j.ceb.2016.04.011. 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(12), 969–971. http://dx.doi.org/10.1038/nmeth.1531, URL: http://www.nature.com/articles/nmeth.1531. Lekka, M., Gnanachandran, K., Kubiak, A., Zieliński, T., & Zemła, J. (2021). Traction force microscopy – Measuring the forces exerted by cells. Micron,150, Article 103138. http://dx.doi.org/10.1016/j.micron.2021.103138, URL: https://www.sciencedirect.com/science/article/pii/S0968432821001293. Mammoto, T., Mammoto, A., & Ingber, D. E. (2013). Mechanobiology and developmental control. Annual Review of Cell and Developmental Biology,29, 27–61. http://dx.doi.org/10.1146/annurev-cellbio-101512-122340, URL: https://www.annualreviews.org/doi/abs/10.1146/annurev-cellbio-101512-122340. Medina, S. H., Bush, B., Cam, M., Sevcik, E., DelRio, F. W., Nandy, K., & Schneider, J. P. (2019). Identification of a mechanogenetic link between substrate stiffness and chemotherapeutic response in breast cancer. Biomaterials,202, 1–11. http://dx.doi.org/10.1016/j.biomaterials.2019.02.018. Nguyen-Ngoc, K. V., Cheung, K. J., Brenot, A., Shamir, E. R., Gray, R. S., Hines, W. C., Yaswen, P., Werb, Z., & Ewald, A. J. (2012). ECM microenvironment regulates collective migration and local dissemination in normal and Malignant Mammary Epithelium. Proceedings of the National Academy of Sciences of the United States of America,109(39), E2595–E2604. http://dx.doi.org/10.1073/pnas.1212834109, URL: https://www.pnas.org/content/109/39/E2595.abstract. Peng, Y., Chen, Z., Chen, Y., Li, S., Jiang, Y., Yang, H., Wu, C., You, F., Zheng, C., Zhu, J., Tan, Y., Qin, X., & Liu, Y. (2019). ROCK isoforms differentially modulate cancer cell motility by mechanosensing the substrate stiffness. Acta Biomaterialia,88, 86–101. http://dx.doi.org/10.1016/j.actbio.2019.02.015. Sabass, B., Gardel, M. L., Waterman, C. M., & Schwarz, U. S. (2008). High resolution traction force microscopy based on experimental and computational advances. Biophysical Journal,94(1), 207–220. http://dx.doi.org/10.1529/biophysj.107.113670, URL: https://pubmed.ncbi.nlm.nih.gov/17827246/. 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(45), 10210–10222. http://dx.doi.org/10.1039/d0sm00789g, URL: https://pubs.rsc.org/en/ content/articlelanding/2021/sm/d0sm00789g. Schwarz, U. S., & Soiné, J. R. D. (2015). Traction force microscopy on soft elastic substrates: A guide to recent computational advances. Biochimica et Biophysica Acta - Molecular Cell Research,1853(11), 3095–3104. http://dx.doi.org/10.1016/j.bbamcr.2015.05.028,arXiv:1506.02394v1. Shields, J. D., Emmett, M. S., Dunn, D. B., Joory, K. D., Sage, L. M., Rigby, H., Mortimer, P. S., Orlando, A., Levick, J. R., & Bates, D. O. (2007). Chemokine-mediated migration of melanoma cells towards lymphatics - A mechanism contributing to metastasis. Oncogene,26(21), 2997–3005. http: //dx.doi.org/10.1038/sj.onc.1210114, URL: https://www.nature.com/articles/1210114. Song, D., Dong, L., Gupta, M., Li, L., Klaas, O., Loghin, A., Beall, M., Chen, C. S., & Oberai, A. A. (2020). Recovery of tractions exerted by single cells in three-dimensional nonlinear matrices. Journal of Biomechanical Engineering,142(8), http://dx.doi.org/10.1115/1.4046974. Song, D., Seidl, D. T., & Oberai, A. A. (2020). Three-dimensional traction microscopy accounting for cell-induced matrix degradation. Computer Methods in Applied Mechanics and Engineering,364, Article 112935. http://dx.doi.org/10.1016/j.cma.2020.112935. 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(2), 171–176. http://dx.doi.org/10.1038/nmeth.3685. 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. http://dx.doi.org/10.1371/journal.pone.0090976. Toyjanova, J., Hannen, E., Bar-Kochba, E., Darling, E. M., Hennan, D. L., & Franck, C. (2014). 3D viscoelastic traction force microscopy. Soft Matter,10(40), http://dx.doi.org/10.1039/c4sm01271b, URL: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4176508/. Vaeyens, M. M., 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. http://dx.doi.org/10.1007/s10456-020-09708-y, URL: http://link.springer.com/10.1007/s10456-020-09708-y. Vining, K. H., & Mooney, D. J. (2017). Mechanical forces direct stem cell behaviour in development and regeneration. Nature Reviews Molecular Cell Biology, 18(12), 728–742. http://dx.doi.org/10.1038/nrm.2017.108, URL: https://www.nature.com/articles/nrm.2017.108. Vitale, G., Preziosi, L., & Ambrosi, D. (2012). A numerical method for the inverse problem of cell traction in 3D. Inverse Problems,28(9), http://dx.doi.org/10. 1088/0266-5611/28/9/095013, URL: http://iopscience.iop.org/article/10.1088/0266-5611/28/9/095013/pdf. Vogel, V., & Sheetz, M. (2006). Local force and geometry sensing regulate cell functions. Nature Reviews Molecular Cell Biology,7(4), 265–275. http: //dx.doi.org/10.1038/nrm1890, URL: www.nature.com/reviews/molcellbio. Wang, H., Abhilash, A., Chen, C. S., Wells, R. G., & Shenoy, V. B. (2014). Long-range force transmission in fibrous matrices enabled by tension-driven alignment of fibers. Biophysical Journal,107(11), 2592–2603. http://dx.doi.org/10.1016/J.BPJ.2014.09.044, URL: https://www.sciencedirect.com/science/article/pii/ S0006349514011096?via%3Dihub. Zündel, M., Ehret, A. E., & Mazza, E. (2017). Factors influencing the determination of cell traction forces. PLoS One,12(2), http://dx.doi.org/10.1371/journal. pone.0172927, URL: https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0172927.