Advanced nonlinear inverse formulations in 3D Traction Force Microscopy: Novel computational Implementations for traction Reconstruction in in vitro Problems

Author: Alejandro Apolinar Fernández Supervisor: Prof. José Antonio Sanz Herrera March 2025 Advanced nonlinear inverse formulations in 3D Traction Force Microscopy : novel computational implementations for traction reconstruction in in vitro problems Ph.D. Thesis “main” — 2025/2/25 — 15:02 — page i— #1 Alejandro Apolinar Fernández Advanced nonlinear inverse formulations in 3D Traction Force Microscopy novel computational implementations for traction reconstruction in in vitro problems Ph.D. Thesis Universidad de Sevilla March 2025 “main” — 2025/2/25 — 15:02 — page ii — #2 “main” — 2025/2/25 — 15:02 — page i— #3 Universidad de Sevilla Escuela Técnica Superior de Ingeniería Tesis Doctoral: Advanced nonlinear inverse formulations in 3D Traction Force Microscopy: novel computational implementations for traction reconstruction in in vitro problems Autor: Alejandro Apolinar Fernández Director: Prof. José Antonio Sanz Herrera El tribunal nombrado para juzgar la Tesis arriba indicada, compuesto por los siguientes doctores: Presidente: Vocales: Secretario: acuerdan otorgarle la calificación de: El Secretario del Tribunal Fecha: Marzo de 2025 “main” — 2025/2/25 — 15:02 — page ii — #4 “main” — 2025/2/25 — 15:02 — page iii — #5 A mis abuelos, Alejandro y Gracia. “main” — 2025/2/25 — 15:02 — page iv — #6 “main” — 2025/2/25 — 15:02 — page v— #7 Resumen La presente tesis doctoral desarrolla el trabajo realizado por Alejandro Apolinar Fernández, alumno del Programa de Doctorado en Ingeniería Mecánica y de Organización Industrial de la Universidad de Sevilla, durante el periodo comprendido entre los años 2020 y 2024. El núcleo del trabajo expuesto se encuentra dentro del marco de la Mecanobiología, área de conocimiento que persigue investigar la relación entre el origen y desarrollo de multitud de procesos biológicos, así como desentrañar la intrincada naturaleza mecánica de los eventos que los integran. En este sentido, se pretende caracterizar el desarrollo de patologías y otros procesos fisiológicos a través del estudio de las interacciones mecánicas entre células y su microentorno. En concreto, el trabajo se enmarca en el estudio de la técnica 3D Traction Force Microscopy (3D TFM). Se trata de una técnica híbrida experimental y computacional con la que es posible reconstruir el campo de tracciones existente en la superficie de las células estudiadas a través de la medición de las deformaciones que las mismas provocan en la matriz extracelular (ECM). Como parte del experimento, se realizan mediciones de los desplazamientos resultantes de la actividad de una célula cultivada en un material biomimético, que en el caso de este trabajo se trata de un hidrogel de colágeno que imita las propiedades mecánicas de la matriz, empleándose para ello técnicas específicas de microscopía confocal. Estas mediciones se obtienen del seguimiento de las posiciones de un conjunto de partículas fluorescentes, llamadas beads, distribuidas en el interior del volumen del hidrogel, entre los dos estados de interés: el estado de referencia en el que la célula se encuentra activa y ejerciendo fuerzas sobre el hidrogel; y el estado final en el que la actividad de la célula es inhibida y el gel retorna a su configuración no deformada. Una vez obtenidos los desplazamientos, estos pueden ser utilizados para calcular tracciones celulares a través de la ley constitutiva seleccionada para representar el comportamiento mecánico del hidrogel. El objetivo último de la metodología es el de facilitar el estudio de la conexión entre los campos de tracciones reconstruidos con procesos biológicos específicos, como pueden ser la adhesión celular, migración celular o la mecanotransducción, entre otros. El entendimiento profundo de estos mecanismos, en especial el de la mecanotransducción, resulta crucial en el desarrollo de terapias para el tratamiento de diferentes cuadros patológicos, como el cáncer, la aterosclerosis, la fibrosis, el asma, la osteoporosis y el fallo cardíaco, entre otros. En TFM, resulta de especial importancia el realizar una caracterización precisa de la cambiante naturaleza de la estructura que compone “main” — 2025/2/25 — 15:02 — page vi — #8 la ECM, la cual experimenta un proceso de remodelación constante que contribuye a regular gran parte de las funciones y procesos celulares elementales. Algunas de estas alteraciones en las propiedades mecánicas de la ECM son originadas por la actividad de las células embebidas y han sido relacionadas con el desarrollo del cáncer y la metástasis. En concreto, las células tumorales son capaces de invadir el tejido cercano a través de la secreción de metaloproteinasas, enzimas proteolíticas que degradan las fibras de colágeno de la ECM que las contiene, pudiendo así facilitar el proceso de migración celular. Este mecanismo modifica la manera en la que la remodelación de la ECM se lleva a cabo e induce heterogeneidades en sus propiedades mecánicas. El presente documento se estructura en ocho capítulos. Los tres primeros son de carácter introductorio. Particularmente, el primero presenta los contextos teóricos general y específico, esto es, una descripción general de la mecanobiología y de la (3D) TFM, respectivamente. El segundo capítulo da una breve introducción a los fundamentos de la teoría de la mecánica del continuo no lineal y la hiperelasticidad, así como presenta en detalle el modelo hiperelástico no lineal empleado para caracterizar los hidrogeles de colágeno considerados en los estudios realizados. Seguidamente, el tercer capítulo expande la teoría de los métodos inversos para la reconstrucción de tracciones, enfatizando la formulación generalizada de la cual se extraen los distintos esquemas de reconstrucción estudiados en los trabajos ejecutados. La formulación se elabora en el contexto de un esquema combinado de Newton-Raphson/Método de los Elementos Finitos (NR/FEM) que proporciona soluciones convergidas en pocas iteraciones. De entre estos métodos se destaca el Physics-based Nonlinear Inverse Method (PBNIM), desarrollado previamente por el autor, el cual se basa en un problema de minimización con restricciones donde se busca un nuevo campo de desplazamientos, que sea lo más cercano posible al medido en el laboratorio, pero que satisfaga la restricción del equilibrio mecánico en la totalidad del dominio analizado. De hecho, la formulación generalizada resulta de añadir un término de regularización de Tikhonov de orden cero al método PBNIM y distinguiendo entre si se considera o no la restricción del equilibrio, pudiéndose de ella obtener métodos restringidos, regularizados, e híbridos entre ambos enfoques. La regularización penaliza los valores elevados de la norma del vector de fuerzas, y su influencia en el proceso de regularización se controla a través de un parámetro que ha de ser calibrado previamente. El cuarto capítulo representa el primer bloque de trabajo llevado a cabo como parte de la tesis doctoral. Este capítulo describe un estudio integral en el que se analiza la precisión con la que “main” — 2025/2/25 — 15:02 — page xiii — #15 regularization). The solutions obtained for both methodologies are studied for three different deformation cases, and for four cell morphologies of varying complexity (higher or lower level of sphericity) corresponding to real cells used in the laboratory. It is concluded that the inverse method is superior to the forward method in all cases, that the reconstruction of tractions becomes more challenging for more complex (non-spherical) cell morphologies, and that the magnitude of the deformations studied does not significantly impact the results. Furthermore, it is confirmed that a linear approximation of the collagen hydrogel’s behavior is not adequate in the cases of analysis considered, and that the errors are not significantly affected by geometrical non-linearities. In the fifth chapter, the second block of work of the thesis is developed, where a synthetic 3D TFM problem is solved with a multiscale approach, in which the temporal evolution of the spatially heterogeneous profiles of the matrix properties due to the degradation phenomenon is taken into consideration. This process is mediated by the secretion of metalloproteinase enzymes by the cell. For this purpose, the degradation process is first modelled by means of a system of five partial differential equations of the reaction-diffusion type, which is solved over time and in three dimensions. Once the degradation profiles are obtained, they are related to the volumetric density of collagen fibres of the hydrogel selected as an artificial substitute for the ECM, and the dependence of the parameters governing the hyperelastic constitutive law considered (Chapter 2) on the spatial distribution of densities is obtained. For this purpose a multiscale fibered matrix generator model, previously developed by the author for another study, is used, which provides the macroscopic behavior of a fiber matrix as a result of the homogenization of the individual contributions to the total strain energy corresponding to the fibers composing the microstructure concerned. The model allows the generation of data derived from synthetic rheological experiments which are used to adjust the parameters of the selected constitutive law according to the different levels of degradation present in the corresponding profiles. Finally, the in silico 3D TFM cases necessary to complete the study are proposed and solved using the same reconstruction methods (forward and inverse) as in the previous chapter. It is concluded that it is not possible to ignore the degradation process without incurring large errors in the reconstruction of tractions, especially for medium and large time scales, and that cellular tractions increase in the case where degradation is neglected. Again, the PBNIM method yields, in each case, substantially more accurate results than the forward method to which it is compared. The sixth “main” — 2025/2/25 — 15:02 — page xiv — #16 chapter makes full use of the generalized formulation presented in Chapter 3, and presents a study that compares the performance of five inverse methodologies in traction reconstruction for an in silico case of 3D TFM. Specifically, the effect of using reconstruction methods with constraints, regularization, or a combination of both concepts in different ways; in terms of accuracy, efficiency (CPU time) and the relative implementation complexity inherent to each method is studied. The formulations employed in the study are the following: mixed PBNIM method (constrained, regularised); Unconstrained 1 method (unconstrained, regularization applied to the cell surface); Unconstrained 2 method (unconstrained, regularization applied to the whole study domain); PBNIM method (constrained, nonregularised); and Forward Modified method (partially constrained, nonregularised). Simulations are carried out for three different levels of noise present in the measured displacements. The results show that, by applying regularization and constraints (based on the fulfilment of fundamental principles; Mixed PBNIM), the best reconstruction of the traction profile is obtained, simultaneously ensuring an optimal estimation of the maximum traction, at the cost of high CPU time and low efficiency. Unconstrained methods (Unconstrained 1 , Unconstrained 2) exhibit high computational efficiency and results of varying quality, with regularization of all forces of the domain giving the best and most physically consistent results in the absence of constraints. The latter suggests the possibility of using such regularization as a weak approximation of the equilibrium constraint. At the same time, regularization-based methods introduce the difficulty of calibrating the associated parameter, which is usually done under subjective criteria. On the other hand, non-regularised but constrained methods (PBNIM, Forward Modified) could provide a good compromise between accuracy and efficiency, while avoiding pre-processing and calibration of the regularization parameter. The importance of considering not only the quality of the traction reconstruction, but also the efficiency and complexity of the implementation (intrusiveness) when selecting an appropriate inverse method for TFM analysis is emphasised. The seventh chapter is framed within a context analogous to that of Chapter 5, i.e. the heterogeneous character of the mechanical properties of the ECM due to the remodeling process induced by the cell under analysis. A new inverse methodology for 3D TFM is presented, capable of reconstructing spatially heterogeneous distributions of ECM stiffness without, unlike the study in Chapter 5, the need for the remodeling phenomenon to be explicitly modelled. This approach formulates the problem as a PDE-constrained inverse method that “main” — 2025/2/25 — 15:02 — page xv — #17 searches for both the displacements and the map of stiffnesses characterising the selected constitutive law. The numerical algorithm developed is then integrated into an iterative NR/FEM framework, avoiding the need for external iterative solvers (that do not require the jacobian matrix of the system). The methodology is validated by using in silico 3D TFM cases based on real cell geometries, modelled in a non-linear hyperelastic framework suitable for collagen hydrogels. The performance of the technique is evaluated through different noise levels and compared to the commonly used iterative L-BFGS methodology. In addition to the originality of the formulation, the effectiveness of the method is demonstrated both in terms of accuracy and efficiency in terms of CPU time. Finally, the eighth and last chapter recapitulates the main conclusions drawn from the work as a whole. In particular, the conclusions described in each chapter are reorganised in terms of the original contributions of the thesis, and then interrelated with each other. Likewise, following these conclusions, the set of future adjustments and improvements of the developed techniques that are planned to be carried out as part of future work related to the exposed thematic, as well as part of other applications in different fields, are described. “main” — 2025/2/25 — 15:02 — page xvi — #18 “main” — 2025/2/25 — 15:02 — page xvii — #19 Agradecimientos En primer lugar, deseo dar gracias al Prof. José Antonio Sanz Herrera, director de esta tesis doctoral, por la siempre acertada guía de mis esfuerzos, que ha ejercido en todo momento durante su consecución. Ha sabido, de manera eficaz, transmitirme sus conocimientos y me ha facilitado los recursos necesarios para que haya podido realizar este trabajo, tan deseado por mí. También he de agradecer la extensa y grata colaboración del Dr. Jorge Barrasa-Fano y el Prof. Hans Van Oosterwyck, de la KU Leuven. La mayor parte del trabajo aquí presentado se ha llevado a cabo haciendo uso de los datos experimentales que ellos han obtenido en su laboratorio, concretamente todas las curvas de reología de los hidrogeles de colágeno considerados en los estudios y las geometrías de las células empleadas en los Capítulos 4 y 5. Además, han contribuido sustancialmente a la consecución de los trabajos en cuestión, ofreciendo perspectivas valiosísimas para su perfeccionamiento. Especialmente, agradezco a Jorge el haber dedicado parte de su tiempo en mi formación al respecto del uso de los paquetes de software de procesamiento de medios típicos de Traction Force Microscopy, así como a él mismo haber hecho cálculos que fueron necesarios en el post-procesado de algunos de los resultados de los estudios. Seguidamente, agradezco al Dr. Pablo Blázquez Carmona y Raquel RuizMateos Brea, del grupo de investigación MecBioLab de la Universidad de Sevilla, a parte de su amistad, su preciada colaboración en la realización de parte del trabajo de la tesis. Al igual que con Jorge y Hans, los últimos dos estudios utilizan geometrías de células obtenidas por ellos en experimentos de 3DTFM realizados en su laboratorio. Pablo, además, contribuyó a la escritura del estudio expuesto en el Capítulo 6. Asimismo, quiero dar las gracias al resto del equipo que conforma MecBioLab: Ana Carrasco Mantis, Juan José Toscano Angulo, Elías Núñez Ortega, María Esther Reina Romo y Jaime Domínguez Abascal; por su apoyo, afable trato y amistad. Agradezco muchísimo a mis tíos, Luis y Gracia, el haberme acogido y tratado no distinto a como se trataría a un hijo. Gracias a ellos, mi estancia en Sevilla ha sido mucho más cálida y esperanzadora, permitiéndome ver a mi abuela, también Gracia, casi todas las semanas. Me he sentido como en casa. Por último, doy las gracias, más que a nadie, a mis padres, Gustavo y Verónica. Porque es a su imperecedero apoyo, y a su amor incondicional, a los que debo el haber conseguido superar con éxito todos los retos a los que me he enfrentado a lo largo de mi vida. “main” — 2025/2/25 — 15:02 — page xviii — #20 “main” — 2025/2/25 — 15:02 — page xix — #21 Acknowledgements First of all, I would like to thank Prof. José Antonio Sanz Herrera, the supervisor of this PhD thesis, for guiding my efforts with great care at all times during its completion. He has efficiently transmitted his knowledge to me and has provided me with the necessary resources to enable me to carry out this work, which I had so much hoped for. I would also like to thank Dr. Jorge Barrasa-Fano and Prof. Hans Van Oosterwyck of the KU Leuven for their extensive and valuable collaboration. Most of the work presented here has been carried out using the experimental data they have obtained in their laboratory, namely all the rheology curves corresponding to the collagen hydrogels considered in the studies as well as the cell geometries employed in Chapters 4 and 5. I am especially grateful to Jorge for dedicating part of his time to my training in the use of the media processing software packages common in Traction Force Microscopy, as well as to he himself for having made calculations that were necessary in the post-processing of some of the results of the studies. Next, I would like to thank Dr. Pablo Blázquez Carmona and Raquel Ruiz-Mateos Brea, from the MecBioLab research group at the University of Seville, for their friendship and helpful collaboration in carrying out part of the work of the thesis. As with Jorge and Hans, the last two studies use cell geometries obtained by them in 3DTFM experiments carried out in their laboratory. Pablo also contributed to the writing of the study presented in Chapter 6. I would also like to thank the rest of the MecBioLab team: Ana Carrasco Mantis, Juan José Toscano Angulo, Elías Núñez Ortega, María Esther Reina Romo and Jaime Domínguez Abascal for their support, kindness and friendship. I am very grateful to my uncle and aunt, Luis and Gracia, for taking me in and treating me no differently than one would treat a son. Thanks to them, my stay in Seville has been much warmer and more hopeful, allowing me to see my grandmother, also Gracia, almost every week. I have felt at home. Finally, I am grateful, more than anyone else, to my parents, Gustavo and Verónica. It is to their undying support and unconditional love that I owe my success in overcoming all the challenges I have faced throughout my life. “main” — 2025/2/25 — 15:02 — page xx — #22 xx “main” — 2025/2/25 — 15:02 — page xxi — #23 Contents 0. Motivation and structure of the thesis 1 0.1. Motivation and objectives . . . . . . . . . . . . . . . . . . . 1 0.2. Structure of the document . . . . . . . . . . . . . . . . . . . 2 0.3. Scientific production . . . . . . . . . . . . . . . . . . . . . . 3 0.3.1. Peer-reviewed publications . . . . . . . . . . . . . . . 3 0.3.2. Unpublished papers . . . . . . . . . . . . . . . . . . 4 0.4. Conferences........................... 4 0.5. Courses ............................. 7 0.6. Funding............................. 7 1. Introduction 9 1.1. Mechanobiology......................... 9 1.2. Traction Force Microscopy . . . . . . . . . . . . . . . . . . . 15 2. Nonlinear continuum mechanics and hyperelasticity 23 2.1. Theoretical framework . . . . . . . . . . . . . . . . . . . . . 23 2.1.1. Fundamental measures of stress and deformation in nonlinear continuum mechanics . . . . . . . . . . . . 24 2.1.2. Hyperelasticity . . . . . . . . . . . . . . . . . . . . . 26 2.2. Theoretical description and computational implementation oftheSAENmodel....................... 28 2.2.1. ABAQUS FEM theoretical background . . . . . . . 28 2.2.2. Computational implementation of SAEN model . . . 30 2.2.3. Alternative formulation of the elasticity tensor . . . 32 2.3. Mechanical characterization process with the SAEN model . 34 2.3.1. Shear rheology . . . . . . . . . . . . . . . . . . . . . 34 2.3.2. Fitting.......................... 35 xxi “main” — 2025/2/25 — 15:02 — page xxii — #24 3. Inverse methods in Traction Force Microscopy 39 3.1. Conceptual definition of the Physics-Based Nonlinear InverseMethod.......................... 40 3.2. Generalized inverse formulation . . . . . . . . . . . . . . . . 40 3.2.1. Unconstrained methods . . . . . . . . . . . . . . . . 42 3.2.2. Constrained methods . . . . . . . . . . . . . . . . . . 43 3.2.3. Finite element discretization . . . . . . . . . . . . . 46 4. Holistic 3D TFM study with real cell morphologies 51 4.1. Introduction........................... 51 4.2. Theoretical background . . . . . . . . . . . . . . . . . . . . 52 4.2.1. Forward method . . . . . . . . . . . . . . . . . . . . 53 4.2.2. Inverse method . . . . . . . . . . . . . . . . . . . . . 54 4.3. In silico TFMsimulations................... 55 4.3.1. Ground truth cases . . . . . . . . . . . . . . . . . . . 56 4.3.2. Displacement reconstruction . . . . . . . . . . . . . . 63 4.3.3. Forward and inverse simulations . . . . . . . . . . . 63 4.3.4. Error indicators . . . . . . . . . . . . . . . . . . . . . 64 4.4. Results.............................. 65 4.5. Discussion............................ 66 4.5.1. Traction reconstruction accuracy referred to cellular pulling force magnitude . . . . . . . . . . . . . . . . 71 4.5.2. Traction reconstruction accuracy referred to cellular morphology....................... 72 4.5.3. Traction reconstruction accuracy and computational efficiency referred to forward/inverse methodologies 73 4.5.4. Traction reconstruction accuracy and computational efficiency referred to hydrogel’s behavior . . . . . . . 74 4.6. Conclusions........................... 75 5. Multiscale framework with heterogeneous matrices 81 5.1. Introduction........................... 81 5.2. Description of the study . . . . . . . . . . . . . . . . . . . . 83 5.3. Modeling metalloproteinase-driven ECM degradation . . . . 84 5.3.1. Phenomenological description . . . . . . . . . . . . . 84 5.3.2. Mathematical model . . . . . . . . . . . . . . . . . . 85 5.3.3. Computational implementation . . . . . . . . . . . . 88 5.4. Mechanical characterization of the hydrogel including degradation.............................. 90 5.4.1. Fibered matrix generator and shear rheology tests . 93 xxii “main” — 2025/2/25 — 15:02 — page xxix — #31 5.9. The progressive effect of degradation on the ECM’s microstructure illustrated: different degraded ECM microstructures, synthetically obtained using the developed fiber generator (see the Appendix for more details). . . . . . 97 5.10. Macroscopic characterization of the degraded material properties of the matrix: fitting of the macroscopic hyperelastic nonlinear model (black dashed lines) to the different virtual shear rheology curves synthetically generated for different levels of ECM density (percentages, coloredsolidlines). ........................ 99 5.11. Preparation of the ground truth case (forces): (left) Cell geometry within the hydrogel domain. (Right) Prescribed contractile (self-balanced) cellular forces on specific regions of the cell boundary. . . . . . . . . . . . . . . . . . . 101 5.12. Preparation of the ground truth case (stiffness): distribution maps of material model parameter κoin a XY cross midsection of the hydrogel (for the considered timepoints, see Figure 5.5) used in the simulations. . . . . . . . . . . . . 103 5.13. In silico 3D TFM results (DEG-DEG): ground truth and both forward and inverse TFM reconstruction results (magnitude of displacements, Frobenius norm of the logarithmic strain tensor, and magnitude of tractions) considering ECM degradation (t= 27). Insets in the figures represent a detail of the variables localized in the tips of the cell. ...............................105 5.14. In silico 3D TFM results (DEG-NODEG): ground truth and both forward and inverse conventional (homogeneous) TFM reconstruction results (magnitude of displacements, Frobenius norm of the logarithmic strain tensor, and magnitude of tractions) neglecting ECM degradation (t= 27). Insets in the figures represent a detail of the variables localized in the tips of the cell. . . . . . . . . . . . . . 106 5.15. In silico 3D TFM results: traction error metric magnitude corresponding to the considered timepoints of the study, for both forward (red boxes) and inverse (blue boxes) formulations. ..........................107 xxix “main” — 2025/2/25 — 15:02 — page xxx — #32 A.1. Schematics of the fiber generator: (a) Definition of the discrete points xi(for n= 3). (b) Different individual randomly-curved fiber in 3D representation, for input parameters n= 3, 5 and 10; and Lf/L = 1.11, 1.21 and 1.4 for blue, green and magenta fibers, respectively. . . . . . . . 111 A.2. Validation of the single fiber computational implementation: Exact solution of a nonlinear cantilever rod subjected to a point dimensionless bending moment ˆ M= ML EI versus its finite element implementation. The analytical solution is an arc with dimensionless radius ˆ R= 1/ˆ M[1]. The rod was discretized into 150 elements, and the moment step was ∆ˆ M=π/500......................113 A.3. Schematics of the proposed multiscale formulation of fibered matrices: relationship between the microand macroscales through homogenization of the RVE-associated variables. ............................115 6.1. Process of data acquisition regarding cell geometry and ECM mechanical behavior from an in vitro experiment: (a) In vitro model of breast cancer cells embedded in a collagen hydrogel. (b) Generation of the 3D cell geometry through the segmentation of confocal images. Scale bar is 2 µm; (c) Finite Element model of the cell embedded in the surrounding ECM; (d) Fitting of the parameters governing a nonlinear hyperelastic model to a curve resulting from a real shear rheology experiment performed on a 1.2 mg/ml collagen hydrogel (ECM model). . . . . . . . . . . . 124 6.2. Preparation of the ground truth scenario from which the input measured displacements are generated: (a) Domain of analysis including the cell embedded in the ECMmimicking hydrogel. (b) FE model showing the tetrahedral mesh of the problem. (c) Prescribed forces (represented by arrows) acting on the tips of the cell protrusions. . . . . . . 126 6.3. Values of the error metrics defined for the analyzed reconstructed variables: (a) Traction error metric. (b) Displacements error metric. . . . . . . . . . . . . . . . . . . 131 xxx “main” — 2025/2/25 — 15:02 — page xxxi — #33 6.4. Traction field magnitude contours represented on the cell surface for the 5 analyzed inverse methods, together with the ground truth solution, corresponding to the high noise case. Magnitude in Pascals: (a) Mixed method (constrained and regularized). (b) Unconstrained 1 method (unconstrained and regularized, cell surface nodes only). (c) Unconstrained 2 method (unconstrained and regularized, cell surface and hydrogel nodes). (d) PBNIM (constrained and non-regularized). (e) Forward Modified method (constrained and non-regularized). (f) Ground truth solution. . . . . . . . . . . . . . . . . . . . 132 6.5. Displacement field magnitude contours represented on the cell surface for the 5 analyzed inverse methods, together with the ground truth solution, corresponding to the high noise case. Magnitude in microns: (a) Mixed method (constrained and regularized). (b) Unconstrained 1 method (unconstrained and regularized, cell surface nodes only). (c) Unconstrained 2 method (unconstrained and regularized, cell surface and hydrogel nodes). (d) PBNIM (constrained and non-regularized). (e) Forward Modified method (constrained and non-regularized). (f) Ground truth solution. . . . . . . . . . . . . . . . . . . . 133 6.6. Nodal forces magnitude in the interior nodes of the hydrogel domain, with respect to their distance to the cell surface, for the 5 analyzed inverse methods corresponding to the high noise case: (a) Mixed method (constrained and regularized). (b) Unconstrained 1 method (unconstrained and regularized, cell surface nodes only). (c) Unconstrained 2 method (unconstrained and regularized, cell surface and hydrogel nodes). (d) PBNIM (constrained and non-regularized). (e) Forward Modified method (constrained and non-regularized). (f) Ground truth solution. Force values are presented in a dimensionless fashion dividing by the norm of the full nodal forces vector (interior and boundary nodes) of the ground truth case. . . . . . . . 134 xxxi “main” — 2025/2/25 — 15:02 — page xxxii — #34 6.7. Optimal value for the regularization parameter (αopt) in the case of the Mixed method, high noise case: (a) L-curve procedure in which the criterion of maximum curvature for the location of the optimal point is employed. (b) Minimum traction error indicator procedure in which the defined error indicator (Eq. (6.2)) is represented as a function of the regularization parameter α. .........136 7.1. Flow diagram of the iterative process: Convergence of the forward subproblem, provided a fixed displacement vector u, yields the discrete stiffness vector κo, which is the input quantity of the inverse subproblem. Convergence of the inverse subproblem, provided a fixed κo, yields the discrete displacement vector u, 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. . . . . . . . . . . . . . . . . . . . . . . 151 7.2. 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. . . 154 7.3. Material stiffness κocontours 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). . . . . . 155 7.4. Material stiffness κocontours 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). . . . . . . . . . . . . 159 xxxii “main” — 2025/2/25 — 15:02 — page xxxiii — #35 7.5. Material stiffness κocontours for the reconstructed ground truth case (tip-localized matrix degradation pattern): (a) 3D view of mid cross-section of the hydrogel domain outlining reconstructed spatial material stiffness distribution (L-BFGS). (b) 3D view of mid cross-section of the hydrogel domain outlining reconstructed spatial material stiffness distribution (NR/FEM). (c) Ground Truth spatial material stiffness distribution. (d) Comparison of reconstructed stiffness profiles with the ground truth profile along γ(seeFigure7.3). .......................160 7.6. Material stiffness κocontours for the reconstructed ground truth case (diffuse matrix remodeling pattern): (a) 3D view of mid cross-section of the hydrogel domain outlining reconstructed spatial material stiffness distribution (L-BFGS). (b) 3D view of mid cross-section of the hydrogel domain outlining reconstructed spatial material stiffness distribution (NR/FEM). (c) Ground Truth spatial material stiffness distribution. (d) Comparison of reconstructed stiffness profiles with the ground truth profile along γ(seeFigure7.4). .......................161 7.7. Material stiffness κocontours for the 10% noise case: (a) Perspective cross section view of the profile obtained from an unique execution of the forward subproblem. (b) Perspective cross section view of the profile obtained from a complete inverse reconstructoin. . . . . . . . . . . . . . . 162 7.8. Material stiffness κocontours for the 2% noise case: (a) 3D view of mid cross-section of the hydrogel domain outlining reconstructed spatial material stiffness distribution. (b) 3D view of mid cross-section of the hydrogel domain outlining the ground truth spatial material stiffness distribution. (c) Reconstructed stiffness profile compared with the ground truth profile along γ(see Figure 7.3). . . . . . . 165 7.9. Material stiffness κocontours for the 5% noise case: (a) 3D view of mid cross-section of the hydrogel domain outlining reconstructed spatial material stiffness distribution. (b) 3D view of mid cross-section of the hydrogel domain outlining the ground truth spatial material stiffness distribution. (c) Reconstructed stiffness profile compared with the ground truth profile along γ(see Figure 7.3). . . . . . . 166 xxxiii “main” — 2025/2/25 — 15:02 — page xxxiv — #36 7.10. Material stiffness κocontours for the 10% noise case: (a) 3D view of mid cross-section of the hydrogel domain outlining reconstructed spatial material stiffness distribution. (b) 3D view of mid cross-section of the hydrogel domain outlining the ground truth spatial material stiffness distribution. (c) Reconstructed stiffness profile compared with the ground truth profile along γ(see Figure 7.3). . . . . . . 167 7.11. Displacements magnitude contours for each of the analyzed cases: (a) 2% noise. (b) 5% noise. (c) 10% noise. (d) Ground Truth solution. . . . . . . . . . . . . . . . 168 7.12. Traction magnitude contours for each of the analyzed cases: (a) 2% noise. (b) 5% noise. (c) 10% noise. (d) Ground Truth solution. . . . . . . . . . . . . . . . . . . 169 xxxiv “main” — 2025/2/25 — 15:02 — page xxxv — #37 List of Tables 4.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). . . . . . . . . . . . . . . . 59 4.2. Fitted nonlinear and linear model parameters. Parameters Eand νof the nonlinear model stand for the Young’s modulus and Poisson’s ratio, respectively (see Section 2 in Chapter 2 for the definition of the nonlinear model parameters). 61 4.3. Size of the hydrogel domain, number of beads (Nbeads), 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. . . . . . . . . . . . . . . . . . 63 4.4. Mean CPU time and mean traction error indicators (averaged over the three analyzed cellular pulling force level cases) of recovered traction solution assuming either linear or nonlinear behaviors of the matrix. Traction errors were calculated with respect to the ground truth solution for the case of the considered nonlinear matrix. . . . . . . . . . . . 66 5.1. Calibrated nondimensional values for the parameters used in the simulations. . . . . . . . . . . . . . . . . . . . . . . . 93 5.2. Calibrated parameters for matrix generation and virtual shear rheology tests. . . . . . . . . . . . . . . . . . . . . . . 95 5.3. Fitted values of the model parameters for the macroscopic hyperelastic nonlinear model in the absence of degradation (non-degraded ECM). . . . . . . . . . . . . . . . . . . . . . 98 5.4. Parameter κo(single fiber stiffness) of the macroscopic hyperelastic nonlinear model, for different levels of ECM density (percentages). . . . . . . . . . . . . . . . . . . . . . . . 98 xxxv “main” — 2025/2/25 — 15:02 — page xxxvi — #38 xxxvi 6.1. Fitted nonlinear model parameters. . . . . . . . . . . . . . . 125 6.2. Values of the defined error metrics (tractions and displacements, percentages) for each method and for each noise level case. The values of αcorrespond to the optimal regularization parameters obtained via the L-curve method. . . . . . 131 6.3. Characteristic CPU times of each analyzed inverse method, for the high noise case. From left to right columns: considered inverse method, CPU time per Newton-Raphson iteration, total number of iterations until convergence, total CPU time. All the computations were performed in a standard laptop (AMD Ryzen 7 4800H 2.90 GHz, 16GB RAM) 135 7.1. Virgin (non-degraded, homogeneous ECM) values of the nonlinear model parameters (Eq. 2.26). . . . . . . . . . . . 153 7.2. Performance metrics for both the L-BFGS algorithm and the methodology introduced in this work (NR/FEM) (in downward direction): Number of iterations, CPU time per iteration, total CPU time taken for the reconstruction. The calibrated values of the regularization parameter of the forward subproblem ακoare given for each degradation pattern.162 7.3. Values of the defined error metrics (stiffness, percentages) for each noise level case. The last three columns contain, respectively, the values of the number of global iterations, local iterations of the forward subproblem and local iterations of the inverse subproblem. The values of ακoand αucorrespond to the optimal regularization parameters obtained via the modified L-curve method described in Section 2.4.................................170 “main” — 2025/2/25 — 15:02 — page 1 — #39 Chapter 0 Motivation and structure of the thesis 0.1. Motivation and objectives This doctoral thesis offers a profound investigation into the computational intricacies of the 3D Traction Force Microscopy (3D TFM) technique, focusing on specific areas of study that had otherwise not been previously explored by other researchers in the past. The work herein presented delves into important topics in TFM, which are condensed into the fulfilment of the following objectives: (i) to study the general aspects of importance in in vitro 3D TFM experiments with respect to their influence on reconstruction results; (ii) to develop novel inverse methodologies for traction reconstruction in nonlinear 3D TFM, in regards to both theoretical and computational aspects; and (iii) to design computational frameworks for 3D TFM in which the mechanical behavior of the extracellular matrix (ECM) is adequately characterized, that is, its heterogeneous and evolving nature is properly modeled and/or retrieved. With regard to the first, associated matters are dealt with holistically and in greater detail than in the previous research conducted by other authors. These include: cell morphology; material and geometrical nonlinearities; noise in the input measured displacements; magnitude of cell forces; selection criteria for adequate traction reconstruction schemes; and computational cost of the algorithms concerned. In the fulfilment of the second objective, the study of constraints and regularization in the nonlinear inverse approaches employed for traction reconstruction, as well as the design of novel methods based on combinations of those two aspects, elaborated within the “main” — 2025/2/25 — 15:02 — page 2 — #40 Chapter 0. Motivation and structure of the thesis 2 context of an original combined Newton-Raphson/Finite Element Method (NR/FEM) computational implementation; are carried out. The NR/FEM numerical scheme distinguishes the methods herein devised from other approaches found in the literature, in which iterative methods that do not require the definition of the jacobian matrix of the system are employed. With respect to the third objective, it is widely known that cancer cells are able to secret sundry proteolytic compounds that alter the arrangement and structural integrity of the collagen fiber-constructed structure of the ECM. However, most TFM studies assume the matrix’s mechanical behavior to be homogeneous, and non-evolving. More than half of the work featured in this document is dedicated to the development of sophisticated computational frameworks aimed at facilitating the realization of traction reconstructions in situations in which heterogeneity and evolution of the mechanical properties of the matrix play a relevant role. The contributions resulting from the different studies featured in this work are of great interest in the field of 3D TFM, in that they are novel and bring significant advancements to the computational paradigms most commonly employed in these applications. Moreover, they provide new insight into the topics previously researched in other publications. 0.2. Structure of the document The document is composed of a total of eight chapters. It is structured as follows: Chapter 1 offers an introduction into the theoretical background of the work, as well as features a review of the state-of-the-art of the thesis’s principal subjects, mechanobiology and the Traction Force Microscopy technique. It provides a thorough description of the concept behind mechanobiological phenomena, as well as the main agents involved. After this, it elaborates on the general aspects of TFM, as well as delves deeper into its three-dimensional application, regarding both experimental and computational considerations. Chapters 2 and 3 present the fundamental mathematical concepts associated to the constitutive and numerical frameworks of the developed methodologies. The former gives a concise introduction to nonlinear continuum mechanics and its most important magnitudes, as well as to the topic of hyperelasticity, and explains the mathematical details of the nonlinear hyperelastic model selected for its use all throughout the thesis. The latter describes thoroughly the philosophical and mathematical principles upon which the inverse methodologies employed in this work are based, in particular concerning the general- “main” — 2025/2/25 — 15:02 — page 9 — #47 Chapter 1 Introduction 1.1. Mechanobiology Mechanobiology is a scientific discipline that focuses on the research of the influence mechanical stimuli possess on the development of the processes characterizing the activity of biological cells and the properties of the associated tissues and organs. Indeed, fundamental insight on the principles through which the latter are altered by mechanical forces and conditions associated to cells and their microenvironment is of great interest to the fields of biology and medicine. Unlike medicine, which typically searches for the genetic or biochemical origin of disease, mechanobiology suggests that changes in the mechanical behavior of cells and the structure of the extracellular matrix (ECM) contributes in regulating a vast collection of physiological and pathological processes [2–4]. It has been found that the mechanical interaction between cells and the surrounding extracellular matrix (ECM) guide functions such as bending, stretching, and repositioning of the epithelium, which are required for the regulation of tissue morphogenesis in embryonic development [5–8]; as well as it activates signalling pathways that regulate cancer invasion [9–13], angiogenesis [14, 15], wound healing [16], survival, proliferation, tissue homeostasis and stem-cell differentiation [17–24]; and influence other tissue or organ pathologies [25– 28]. Ultimately, these processes entail important implications in human health, driving the onset and progression of a multitude of disease scenarios which provide substantial decreases in quality of life and are potentially life-threatening, such as atherosclerosis, fibrosis, asthma, osteoporosis, and heart failure, among others [2, 29–32]. During the last last decades, the acknowledgement of the importance “main” — 2025/2/25 — 15:02 — page 10 — #48 Chapter 1. Introduction 10 of mechanical phenomena in cellular response has been conducive to the establishment of this emergent field of research, which constitutes a powerful platform for the integration of the disciplines of biology, physics, chemistry, mathematics and engineering. Traditionally, the most popular term "biomechanics" has been employed to describe the mechanical component of the study of biological systems. As fundamental understanding of cell biology increased, mechanobiology naturally evolved from biomechanics by integrating crucial elements of molecular and cellular biology. Consequently, the modern term "mechanobiology" constitutes, in current times, the most accurate descriptor of an interdisciplinary field of study that deals with the comprehensive cascade of biological events initiated or mediated by mechanical forces. Mechanobiology includes the study of physiology at multiple scales, namely the cellular, molecular and tissue scales. The study of mechanobiology is of great importance in the search for fundamental understanding of cell function and, in essence, its implications on human health. Regarding its current state of development, characterized by high activity, mechanobiology is considered as an emerging field. Nonetheless, conceptually similar forms to the latter can be traced back for more than a century. Biologist D’Arcy Thompson suggested in its 1917 book "On Growth and Form" that physico-mathematical principles could be employed for the modeling of structure and growth patterns in animals and plants, as well as that the morphology of living organisms may be influenced by mechanical states [33]. In the context of the evolutionary paradigm that was, at that time, universally accepted for the description of the dynamics of physiological changes in living beings, and in which the consideration of the role played by physical laws was largely ignored, Thompson’s work provided valuable insight into the reality of biological phenomena, planting the seed towards the sprouting of novel lines of interdisciplinary investigation. Technology ultimately limited the feasibility of addressing Thompson’s hypotheses, which paused the development of the newly conceived field of research for decades. It was not until the 1980s that essential technological advancements, which built the platform for the establishment of mechanobiology as a robust field of research, appeared, including tools for the measurement of mechanical-related variables, fluorescent imaging techniques, and nanomicro manufacturing. Most notably, the invention of the Atomic Force Microscope (AFM) and the optical tweezers played a paramount part in mechanobiological investigation within in vitro contexts. The AFM en- “main” — 2025/2/25 — 15:02 — page 11 — #49 Chapter 1. Introduction 11 abled the characterization of the morphological and mechanical aspects of complex biological systems by means of the application and measurement of forces with values contained within the pico- (pN) and micro-Newton (muN) ranges in precise spatially delimited regions [34]. Optical tweezers, originally called single-beam gradient force trap, provided the possibility of meticulous physical interaction with microscopic and sub-microscopic objects such as atoms, nanoparticles and droplets; them being manipulated in a way similar to how normal tweezers are used. This is achieved by using a highly focused laser beam to exert repulsive and attractive forces in order to hold and move the specified object. Optical tweezers have been effectively utilized in the study of mechanotransduction, cell micro-rheology, membrane mechanics, as well as cell-cell interactions [35]. The subsequent decade would see important improvements in the technologies associated to the fabrication of matrix microenvironments and microscopy imaging that led to the conception of the Traction Force Microscopy (TFM) technique. TFM allowed researchers to more accurately determine traction forces and investigate their impact on key biological processes such as embryogenesis, inflammation, cancer metastasis and angiogenesis [36, 37]. Moreover, the recent development of super-resolution microscopy and single-molecule imaging techniques has provided the means for visualization of mechanobiological dynamics at a molecular scale, within which the respective interactions of nucleic acids and proteins can be analyzed [38, 39]. Mechanosensing is the process by which cells detect mechanical signals from their microenvironment, effectively sensing and responding to changes in the mechanical stiffness of the ECM. Cells probe the rigidity of their surrounding microenvironment by actively interacting with the ECM via the exertion of traction forces. They perform this interaction with the embedding substrate via transmembrane proteins called integrins [40]. The intricate nature of the mechanosensing process remains rather obscure to this day, and research is focused on unveiling how cells are able to sense matrix stiffness by traction force exertion, as well as the translation of the gathered mechanical information into specific cellular activities related to growth, motility and differentiation. The latter process is termed mechanotransduction. A prominent motive for the complexity of this process is the elevated multiplicity of agents that participate in it. More particularly, experimental reserach has identified a large group of mechanosensors and mechanotransducers that operate together within the cell’s mechanical information processing chain [41]. These include: paxillin [42], vinculin [43, 44], talin [45], p130CAS [46, 47], integrins [48, 49], the actin cytoskele- “main” — 2025/2/25 — 15:02 — page 12 — #50 Chapter 1. Introduction 12 ton (CSK) [50–52] and mechanosensitive ion channels [53]. Additionally, mechanotransduction is the mechanism by which cells are able to translate mechanical stimuli into a differentiated response. Resulting from specific cell mechanotransduction pathways, a vast set of phenomena associated with physiological changes emerges. Some of these are apparent, such as the increase in muscle mass due to physical exercise (or decrease when neglecting it) and the effective change in mass in musculoskeletal tissues, such as in bones, ligaments and tendons, as stated in Wolff’s and David’s laws, respectively, when subject to varying degrees of mechanical loading. The translation of pressure waves reaching the hair cells of the inner ear into auditory signals that can be interpreted in the brain constitutes a good example of mechanotransduction [54]. The heart’s regular activity, which is carried out during such extended periods of time, is sustained by a soft tissue structure that evolves continually due to hypertrophy of the composing tissue, which in pathological circumstances will alternatively entail the potentially catastrophic event of heart failure [55]. Mechanosensing and mechanotransduction processes involve both intra and extra-cellular components. The key participants in the process are three: integrins, ECM proteins and the cell cytoskeleton (CSK). These elements are responsible for regulating gene expression through the begetting of multiple downstream signaling cascades, which ultimately determines the fate and behavior of cells. Integrins are transmembrane glycoproteins that establishes the physical connection between cells’ CSKs in cell-cell adhesion, as well as between CSK to the ECM in cell-matrix adhesion. They gather in localized regions of the cell surface, forming clusters known as "focal adhesions" (FAs). Numerous types of integrins exist, and one cell generally exhibits multiple different types on its surface. FA formation is crucial for enabling mechanosensing. The CSK is the protein structure that ensures cell overall morphology is preserved, providing the required mechanical strength. Its function allows cells withstand external forces while undergoing controlled deformation processes in a variety of possible configurations [56]. The CSK is constructed by an interwoven network of protein filaments, which are of three types in the case of the cells of mammals: actin, microtubules (MTs) and intermediate filaments (IFs). Generally, actin and IFs are the main contributors to cell overall stiffness regarding tensile and shear loads, whereas MFs more specifically provide resistance to compression [57]. These net- “main” — 2025/2/25 — 15:02 — page 13 — #51 Chapter 1. Introduction 13 works become stiffer when subjected to increasing loads, that is, they exhibit strain-stiffening nonlinear mechanical behavior, which effectively limits the amount of deformation cells and epithelial tissues can be subjected to [58–60]. The extracellular matrix (ECM) is composed of an intricate protein microstructure that enables cell adhesion, providing mechanical support and constraint necessary for all cell activity. It constitutes the medium within which cells operate, and serves as a platform for the diffusion of bio-chemicals in tissues. Within the ECM, a rich variety of growth factors, cytokines and proteolytic enzymes that intervene in cell-ECM interactions are contained. There are two main types of ECM. The first one, the basement membranes, correspond to thin structures that serve as a twodimensional substrate to which polarized cells, e.g. epithelial and endothelial cells, adhere. Basement membranes are composed of laminin, type IV collagen, nidogen and heparan sulfate proteoglycnans [61]. The second one correspond to the connective tissue. It provides a three-dimensional structure composed mainly of fibers correponding to varying mixes of different types (I, II, III and V) of collagen; as well as of proteoglychans and glycosaminoglycans [62]. These collagen microstructures present different configurations in regards to the geometrical characteristics of the fibers, depending on the native tissue’s biomechanical function. A notable example of this is the compositional differences existing between stiff tissues such as tendons, which exhibit thick and aligned fibers that contribute to tensile strength, and flexible tissue like the cornea, which presents a set of thin fibers organized in interweaving orthogonal layers, creating mesh-like disposition that provides optical transparency. The ECM is a highly dynamic structure that exists under the influence of a constant remodeling process, where ECM components are deposited, degraded, or otherwise modified. It constitutes a major component of a cell microenvironment, taking part in most basic cell behaviors, from cell proliferation, adhesion and migration, to cell differentiation and death [63]. In the context of continuously changing cellular activity, remodeling of the extracellular matrix represents a crucial mechanism whereby diverse cellular behaviors are regulated [64]. Indeed, normal development, physiology, and robustness of organ systems require ECM dynamics to be tightly regulated. Alteration or dysfunction of such regulation leads to congenital defects and diseases, including cancer. As a matter of fact, disturbance of the normal mechanical conditions of the ECM is implicated in the onset and progression of cancer and fibrosis [65]. Survivavility of tumor cells and their capacity for proliferation are “main” — 2025/2/25 — 15:02 — page 14 — #52 Chapter 1. Introduction 14 positively correlated to drastic increases in stiffness, cancer tissues being found to become up to 10-fold stiffer than healthy tissues [65–67]. Regarding the study of the process of cancer mestastasis, in particular invasive cellular migration towards surrounding tissue by means of degradation of the collagen microstructure that supports the ECM, is of great interest. Cells are able to induce ECM degradation through the expression of matrix metalloproteinases (MMPs), a type of proteolytic enzymes which have the capacity to degrade (potentially) all components of the surrounding tissue [68–73]. MMP-induced degradation hence facilitates cancer growth and spread by virtue of the new available space generated, as well as by other proteins that are released by the degraded tissue [74, 75]. Furthermore, cancer cells deposit siginificant amounts of collagen in their microenvironment, resulting in localized stiffened regions of the matrix that serve different purposes, e.g. aiding in the cicatrization process during wound healing [76]. These regions have also been found to act as walls that shield cancer cells from drugs administered during treatments [77, 78]. Researchers aim at unraveling these mechanisms by employing in vitro models that emulate the composition, organization and mechanical characteristics of human-like tissue, which have provided important advancements in tissue engineering and therapeutic techniques [79–81]. Mechanobiology indeed offers powerful applications in the framework of tissue engineering. Cell sensitivity to forces and to fluctuations in substrate stiffness can be harnessed to construct biomaterials that effectively guide stem and resident cells towards generating functional replacement tissue in medical patients. Moreover, fundamental understanding of the mechanical behavior of the matrix can be employed in the programming of stem cell differentiation within organ-on-chip and regenerative medicine scenarios. As an example, in a recent in vitro study [82], flexible fiber networks of different stiffness were developed by taking advantage of the swelling properties of GelMA (gelatin methacryloyl). The activity of mesenchymal stem cells (hMSCs), embedded in these matrices, was then assessed. The authors found that, even for high fiber stiffness, flexible fiber structures substantially improved cell mechanosensing, managing to establish the coexistence of fiber recruitment and a high elastic modulus of the fibers, which strengthen cellular adhesion. Mechanobiology constitutes a crucial platform for the quantitative research of mechanical scenarios with regards to their implications towards health and disease. First, it establishes the necessary tools for the mechanical characterization of the biological media in which the investigated “main” — 2025/2/25 — 15:02 — page 15 — #53 Chapter 1. Introduction 15 processess take place. Second, it aims at extracting fundamental knowledge about the mechanical variables implicated in the diverse molecular processess associated to pathological conditions, known to be induced by the perturbation of stable and healthy cellular mechanical states. Finally, from the interpretation of the obtained results, it allows for the elaboration of novel therapeutic treatments for the studied medical conditions. The work presented in this Ph.D. Thesis belongs to the study of novel computational methodologies for 3D Traction Force Microscopy, a powerful technique extensively used in mechanobiology for the indirect measuring of cellular tractions. It will be thoroughly described, regarding its theoretical foundations and state-of-the-art implementations, in the following sections. 1.2. Traction Force Microscopy Assessing the impact of mechanics on cell behavior requires the precise measurement and quantification of the involved mechanical variables within the cell-matrix interface, for which a sufficiently high spatio-temporal measurement resolution is needed. Generally, forces do not constitute an accessible magnitude within the context of an experiment, and its measurement is rather challenging. In consequence, several methodologies have been devised over time to achieve this goal [83, 84]. 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. Within these, Traction Force Microscopy (TFM) is emphasised by its effectiveness and development potential, and which has been the main choice for calculating forces exerted by cells in the ECM for the last decades [85–94]. Conceptually, TFM aims at inferring cellular forces from some kind of traced motion, e.g., the deformations resulting from the cell-matrix interaction. The technique enables the reconstruction of traction fields associated to specific cell-ECM interaction scenarios, deriving from displacement data captured via microscopy imaging techniques applied to different stages of those scenarios in the laboratory. In brief, an in vitro situation in which a cell is cultivated on top or embedded in a substrate is prepared, and monitoring of its activity is conducted. After the associated displacements are measured during the execution of the experiment, cellular tractions result from the output of a computational algorithm that relates the mechanical properties of the substrate, its deformation state and the cell morphological features. “main” — 2025/2/25 — 15:02 — page 16 — #54 Chapter 1. Introduction 16 Traditionally, TFM has been extensively performed employing 2D in vitro cultures where cells are seeded on top of a substrate, with tractions being reduced to act on a direction contained in the cell-ECM contact plane. 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 [95–99]. However, the dimensionality of the surrounding microenvironment is crucial for characterizing cell behavior [3, 100]. 3D ECM-mimicking media serves as a platform for a more accurate traction characterization, in which all possible spatial directions for cell-ECM interactions are considered, yielding a more phyisiologically relevant methodology for the assessment of the geometry and contractile behavior of biological cells [101, 102]. Synthetic materials with 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 [3, 100, 103]. More importantly, in the recent years, the availability of synthetic hydrogel environments for cell culture has helped the development of experimental protocols allowing the execution of more robust 3D TFM in vitro experiments [104–112]. The general 3D TFM workflow is summarised as follows: cells are fully embedded in the ECM-mimicking material, which contains spatially distributed fiducial markers (often fluorescent particles) known as "beads". Once cell activity is initiated, a set of microscopy images are generated from which information about cell morphology and the motion of the beads is retrieved. In particular, two sets of images are generated, one corresponding to the cell stressed state (reference configuration), and the other corresponding to the cell relaxed state attained after cell removal or lysis. By means of cross-correlation techniques, a measured displacement field is obtained by the tracking of the beads motion between the two recorded activity states. Finally, the measured displacements can be interpolated and mapped to a prescribed mesh grid of points and fed to a computational algorithm that, through a certain mathematical formulation that incorporates the constitutive behavior selected for the matrix, retrieves the tractions exerted at the cell-ECM interface. Fig. 1.1 shows a detailed illustration of the different elements that are part of the general 3D TFM workflow. An alternative approach concerning the 3D nature of the problem is the so-called 2.5D TFM, which constitutes an intermediate step between 2D and 3D TFM [16, 113, 114]. This methodology is analogous to 2D TFM in regards to the experimental setup (cell seeded on top of “main” — 2025/2/25 — 15:02 — page 17 — #55 Chapter 1. Introduction 17 a substrate of certain thickness) but, on the contrary, 2.5D TFM computes 3D tractions along the zdirection. These components are obtained from the tracking of bead displacements in 3D image stacks. In contrast to 2.5D TFM, 3D TFM requires the cell under study to be completely immersed in a three-dimensional domain throughout the course of the experiment. Compared to 2D TFM, (full) 3D TFM presents substantially higher complexity of execution with respect to both its experimental and computational components. First, imaging 3D cell-induced matrix deformations requires acquisition of large image volumes, which can take several minutes with the risk of phototoxicity of cells and photobleaching of the fluorescent markers. Moreover, the field of view needs to be large enough to contain the entire cell of interest (including any large protrusions) and sufficient space around it to visualize bead displacements. The literature provides guidelines with protocols and good practices to acquire 3D TFM data using confocal microscopy [115, 116], second harmonic generation (SHG) [117], high resolution techniques such as stimulated emission depletion (STED) microscopy [118, 119], and fast imaging techniques such as optical coherence microscopy (OCM) [120, 121]. Furthermore, acquiring the image corresponding to the relaxed state of the cell in 3D also differs from 2D. In 2D TFM it is common to remove the cells by simply detaching them from the substrate with sodium dodecyl sulphate (SDS), trypsin, proteinase K, Triton X-100 or via a micromanipulator [90, 95, 122]. In 3D, obtaining the relaxed state of the cell involves adding drugs like Cytochalasin D [105, 112, 116] or sodium azide [9] to the cell culture medium. However, it is crucial to calibrate the drug concentration and the waiting time after drug addition to verify that cells are fully relaxed when acquiring the relaxed state [15]. Alternatively, some studies have bypassed the use of drugs obtaining the relaxed state right after hydrogel polymerization and prior to cell attachment [123, 124]. Concerning the collection of the required displacement data, to measure matrix deformations from the acquired image stacks, multiple algorithms compatible with 3D exist in the literature, namely, particle tracking [125, 126], block-matching based algorithms [127] and free form deformation-based image registration [117, 128]. Finally, to calculate cell tractions, two ingredients are needed: the mechanical properties of the hydrogel and a method that calculates tractions from the measured displacements. In 2D TFM, linear elastic polyacrylamide (PAA) hydrogels are commonly used, requiring only the measurement of the elastic modulus and Poisson’s ratio to define their mechanical behavior. Such hydrogels are highly robust with low variability, and the “main” — 2025/2/25 — 15:02 — page 18 — #56 Chapter 1. Introduction 18 literature has extensively reported on the relationship between PAA concentration and its mechanical properties [91, 95, 129–131]. However, in 3D TFM fibrillar nonlinear elastic collagen hydrogels are often used, requiring time-consuming characterizations by means of rheometry to measure the storage and the loss modulus via time sweeps and the nonlinear elastic response via strain/stress sweeps [105, 132]. With regards to the computational method to recover tractions, while elegant closed-form solutions are available in 2D TFM (Green’s functions method, for example), these are generally not applicable in 3D TFM. Instead, numerical methods in the context of the Finite Element Method have been proposed in the last decade with increasing complexity, from linear elastic approaches at small strains [104, 112], to those that take into consideration nonlinear elasticity and geometric nonlinearities combined with ECM degradation models [107, 133–136]. Regarding its computational part, the field of 3D TFM has led to the development of complex and advanced methodologies during the last decade. Traction reconstruction in TFM is approached in two fundamental ways: the forward and inverse methods. On the one hand, forward methods aim to obtain the associated strain field from which to retrieve stresses and tractions in a direct manner, that is, taking the measured displacements as input data, prescribing all or a subset of them. [91, 92, 137]. The conventional approach for forward TFM found in the literature assumes prescribing the whole measured displacement field at the corresponding nodes within the ECM domain. Some authors reconstruct the strain tensor directly from the measured displacements (which are interpolated and transformed into a continuous field) via numerically differentiating them in the process of obtaining the displacement gradient [110, 138, 139]. The latter is introduced in the constitutive equation selected for the matrix to compute the corresponding stresses and tractions. Other implementations feature the modification of the measured displacements aiming at reducing the impact of measurement noise in the process of differentiating the measured displacements by, for example, employing a displacement-gradient technique [140, 141]. Normally, for 3D cases, a Finite Element Analysis software is used to solve the problem and retrieve the field variables. This methodology poses relatively low computational complexity and therefore its implementation and execution can be fairly straightforward, shining due to its high computational efficiency. 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 deriva- “main” — 2025/2/25 — 15:02 — page 25 — #63 Chapter 2. Nonlinear continuum mechanics and hyperelasticity 25 particles after the deformation process in terms of their relative material position prior deformation, F=∂x ∂XFiJ =∂xi XJ ∀i, J = 1,2,3(2.3) In this sense, Fallows the transformation of vectors expressed in the material (reference) configuration into vectors expressed in the spatial (current) configuration. Because of this, it is known as a two-point tensor. This means that Fis not expressed in either of the two configurations, but "in the middle" of both. As a consequence, indices pointing to components associated to different configurations are distinguished from one another by means of lowercase (spatial) and uppercase (material) print. Analogously, other magnitudes which operate exclusively on material or spatial elements can be defined, Material configuration Right Cauchy–Green deformation tensor,C=FTF Lagrangian or Green strain tensor,E= [C−1]/2 Spatial configuration Left Cauchy–Green or Finger tensor,b=F F T Eulerian or Almansi strain tensor,e= [1−b−1]/2 Regarding stress measures, although these are described in greater detail together with the introduction to the concept of strain energy density function in the next section, it is anticipated herein which ones are of material, two-point and spatial nature, Material.Second Piola-Kirchhoff stress tensor (S). Two-point.First Piola-Kirchhoff stress tensor (P). Spatial.Cauchy stress tensor (σ) and Kirchhoff stress tensor (τ). It is convenient to define an operator by means of which relations between the different stress measures described by the tensors presented above can be expressed. The following operators are thereby defined: the push-forward operator is defined as that which enables the transformation “main” — 2025/2/25 — 15:02 — page 26 — #64 Chapter 2. Nonlinear continuum mechanics and hyperelasticity 26 of magnitudes from the material configuration into the spatial configuration, ϕ∗[]; and the pull-back operator as that which transforms magnitudes from the spatial configuration into the material configuration, ϕ∗[]. Both operators are mutually inverse of one another, and the specific operation designated by them depends on the magnitude to which they are applied. In this way, the relations existing between the different stress measures through the push-forward operator are the following, τ=Jσ=ϕ∗[P] = P F T(2.4) τ=Jσ=ϕ∗[S] = F SF T(2.5) 2.1.2. Hyperelasticity The feature of elasticity is defined as that of a material whose constitutive equation is uniquely dependent on its current state of deformation. In this respect, hyperelasticity is more generally defined as the characteristic of elastic materials whose constitutive equation derives directly from a scalar function that relates the strain energy with the deformation gradient tensor. The latter is known as strain energy density function or elastic potential. This is possible due to the work carried out by stresses during the deformation process depending exlusively on the initial and final times, and therefore not depending on the intermediate path between the two. This function is defined as, Ψ(F(X),X) = Zt to P(F(X),X):˙ Fdt −→ ˙ Ψ = P:˙ F(2.6) where notation ˙ [] indicates deriative with respect to time. Every strain energy density function must satisfy the objectivity condition, that is, it has to remain invariant under rigid body motions. Taking into account that the deformation gradient tensor can be decomposed into the combination of a rotation and a stretch (polar decomposition), F=RU, the elastic potential needs to depend exclusively on the stretch component U. Conveniently, the elastic potential is often expressed in terms of the Right Cauchy–Green deformation tensor as U2=C, hence allowing for the reinterpretation of the expression given by (2.6) as, Ψ(C(X),X) = Zt to S(C(X),X):˙ Edt −→ ˙ Ψ = S:˙ E(2.7) “main” — 2025/2/25 — 15:02 — page 27 — #65 Chapter 2. Nonlinear continuum mechanics and hyperelasticity 27 Pairs Pand ˙ Fand Sand ˙ Eare work conjugate, enabling the definition of the constitutive equation in the mixed (two-point) and material configurations, respectively, in the following way, P(F(X),X) = ∂Ψ(F(X),X) ∂FPiJ =∂Ψ ∂FiJ (2.8) S(C(X),X) = ∂Ψ(C(X),X) ∂ESIJ =∂Ψ ∂EIJ (2.9) Moreover, in the context of an isotropic material, the elastic potential must be independent of any particular direction in the material configuration, that is, it must depend only on the three first invariants of the Ctensor, Ψ(C(X),X) = Ψ(I1, I2, I3,X)(2.10) where the invariants of Care defined as, I1=trace[C] = C: 1 (2.11) I2=trace[C2] = C:C(2.12) I3= det[C] = det[F]2=J2(2.13) Generally, the relations given by (2.8) and (2.9) are nonlinear. In the framework of a numerical iterative scheme, such as a Newton-Raphson procedure, they appear linearized with respect to an increment of the displacement field uin the corresponding configuration. In light of this, said relations are presented in the following fashion, DP[u] = A:DF[u]DS[u] = C:DE[u](2.14) where Df[u]designates the directional derivative of ftaken in the direction of u.Aand Care fourth-order tensors known as First (mixed configuration) and Second (material configuration) Elasticity Tensor, respectively. These are described by the following expressions, A=∂P ∂FAiJkL =∂PiJ ∂FkL (2.15) C=∂S ∂ECIJKL =∂SIJ ∂EKL (2.16) “main” — 2025/2/25 — 15:02 — page 28 — #66 Chapter 2. Nonlinear continuum mechanics and hyperelasticity 28 Both tensors, as they are the result of the differentiation of a function (elastic potential) with respect to the same quantity twice (deformation gradient tensor), exhibit major symmetry, AiJkL =AkLiJ CIJKL =CKLIJ (2.17) However, only the second elasticity tensor exhibits minor symmetries, as a consequence of differentiating a symmetric quantity (second PiolaKirchhoff stress tensor) with respect to another symmetric quantity (lagrangian or Green strain tensor), CIJKL =CJIKL CIJKL =CIJLK (2.18) 2.2. Theoretical description and computational implementation of the SAEN model The aim of this section is to present a thorough description of the nonlinear hyperelastic model employed throughout all the work presented in this thesis, as well as the derivation process of the corresponding fundamental mechanical relations, to be numerically implemented into the Finite Element Analysis software ABAQUS via user subroutine UMAT, according to the numerical scheme used in all herein presented numerical implementations. The model concerned was devised by the authors in [105]. 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. 2.2.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 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, “main” — 2025/2/25 — 15:02 — page 29 — #67 Chapter 2. Nonlinear continuum mechanics and hyperelasticity 29 δW := Zv σ:δεdv −Zv ρb·δudv −Z∂v t·δuda = 0 (2.19) where δudenotes the virtual displacement, vthe current volume, ρthe current density, bthe volumetric forces vector, δεthe linear virtual strain tensor, tthe traction vector resulting from the product between the Cauchy stress tensor σand the vector normal to the surface ∂v,n. A linearized form of Eq. (2.19), 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 ∆u. Being ϕka test function, δW(ϕk, δu) + DδW(ϕk, δu)[∆u]=0 (2.20) ABAQUS employs the Updated Lagrangian configuration, by means of which the Principle of Virtual Work is expressed as, δW := ZV τ:δεdV −ZV ρob·δudV −Z∂v t·δuda = 0, δε=1 2h∇δu+ [∇δu]Ti (2.21) 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. (2.21), 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, DZV τ:δεdV [∆u] = ZV D[τ:δε][∆u]dV = =ZVh∇δu:[J c J− H ]:∇δu+τ:h[∇δu]T∇∆uiidV (2.22) where c Jis the Jaumann elasticity tensor, which can be computed as, J c J=J c +H, Hijkl =1 2[δjkτil +δikτjl +δjlτik +δilτjk](2.23) being c the so called Fourth Cauchy Elasticity Tensor, c =1 Jϕ∗[C], cijkl =1 JFiIFjJ FkKFlLCIJKL (2.24) where Einstein’s summation convention has been used. “main” — 2025/2/25 — 15:02 — page 30 — #68 Chapter 2. Nonlinear continuum mechanics and hyperelasticity 30 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. J c J. 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          DDSDDE =          D1111 D1122 D1133 D1112 D1113 D1123 D2211 D2222 D2233 D2212 D2213 D2223 D3311 D3322 D3333 D3312 D3313 D3323 D1211 D1222 D1233 D1212 D1213 D1223 D1311 D1322 D1333 D1312 D1313 D1323 D2311 D2322 D2333 D2312 D2313 D2323          (2.25) 2.2.2. Computational implementation of SAEN model Steinwachs et al. [105] 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 than 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 undergoes 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 λ, Ψ′′(λ) = κo     exp[λ/do]∀λ < 0 1∀0≤λ<λs exp[[λ−λs]/ds]∀λ≥λs λ=∥F·eΩ∥ − 1 (2.26) “main” — 2025/2/25 — 15:02 — page 31 — #69 Chapter 2. Nonlinear continuum mechanics and hyperelasticity 31 where κois a material parameter with stress dimensions, dois the buckling dimensionless parameter, dsis the strain-stiffening dimensionless parameter and λsis 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, eΩis the unit vector that represents the fiber’s orientation and Fis the deformation gradient. Fig. 2.1 shows a conceptual description of the nonlinear model under consideration, illustrating the progression of the deformation process undergone by a fiber, as it transitions through the three different states a fiber can exist in: buckling, straightening and stretching. The mechanical stress of the fiber can be then defined by the expression: Ψ′(λ) = Zλ 0 Ψ′′(λ)dλ =κo     do·[exp[λ/do]−1] ∀λ < 0 λ∀0≤λ<λs λs+ds·[exp[[λ−λs]/ds]−1] ∀λ≥λs (2.27) in which integration from zero indicates that the material is not prestressed. The constitutive relation in mixed configuration can be obtained afterwards as, P=∂Ψ(λ(F)) ∂FΩ =∂Ψ ∂λ ∂λ ∂FΩ =Ψ′(λ)[F eΩ]⊗eΩ ∥F eΩ∥Ω (2.28) where Pdenotes the First Piola-Kirchhoff stress tensor and brackets indicate averaging over all directions Ωin the unit sphere. From Eq. (2.28), the Cauchy stress tensor is obtained by applying a push-forward operation as follows, σ=1 Jϕ∗[P] = 1 JP F T=1 J1 4πZ2π 0Zπ 0 Ψ′(λ)[F eΩ]⊗eΩ ∥F eΩ∥sin θdθdϕ·FT (2.29) where the integral over the unit sphere is approximated by means of the Gauss-Legendre quadrature. Contrary to other elasticity tensors, the spatial elasticity tensor cannot be directly obtained from the elastic potential. Therefore, the definition of the tangent stiffness matrix DDSDDE depends ultimately on the definition of the second elasticity tensor C, which can be transformed using Eqs. “main” — 2025/2/25 — 15:02 — page 32 — #70 Chapter 2. Nonlinear continuum mechanics and hyperelasticity 32 (2.23) and (2.24) to accommodate to the ABAQUS format. The fourth invariant of the Right Cauchy-Green deformation tensor is described as, I4=eT Ω·C·eΩ=eT Ω·FTF·eΩ=∥F·eΩ∥2= (λ+ 1)2(2.30) association which permits the rewriting of λin terms of the Right CauchyGreen deformation tensor, λ(F) = ∥F·eΩ∥ − 1−→ λ(C)=(I4)1/2−1(2.31) In this way, the constitutive relation in material configuration is obtained as, S= 2∂Ψ(λ(C)) ∂CΩ = 2∂Ψ ∂λ ∂λ ∂CΩ =D[λ+ 1]−1·Ψ′(λ)·[eΩ⊗eΩ]EΩ (2.32) where ⊗denotes tensor product. The Second elasticity tensor can then be retrieved from, C= 2∂S∗(λ(C)) ∂CΩ = 2∂S∗ ∂λ ⊗∂λ ∂CΩ (2.33) in which S∗denotes the Second Piola-Kirchhoff stress tensor of a single fiber (non-homogenized) and, ∂S∗ ∂λ =h−[λ+ 1]−2·Ψ′(λ)+[λ+ 1]−1·Ψ′′(λ)i·[eΩ⊗eΩ](2.34) ∂λ ∂C=1 2[λ+ 1]−1·[eΩ⊗eΩ](2.35) 2.2.3. Alternative formulation of the elasticity tensor It is described herein an alternative formulation of the second elasticity tensor, this time obtained from the first elasticity tensor in mixed configuration. As indicated in the previous section, the spatial elasticity tensor cannot be directly obtained from the elastic potential. Alternatively, said tensor must be obtained through the relation existing between the material elasticity tensor (2.24). Revisiting the expression (2.15), the first elasticity tensor in mixed configuration can be retrieved as, “main” — 2025/2/25 — 15:02 — page 33 — #71 Chapter 2. Nonlinear continuum mechanics and hyperelasticity 33 Figure 2.1: Conceptual description of the SAEN nonlinear model: the mechanical stiffness (Ψ′′(λ)) of a single fiber is represented in terms of λ. The plot features three distinct regions corresponding to the three different states a fiber can exist in: buckling (compression), straightening, and stretching. A=∂P(λ(F),F) ∂F=∂P ∂λ ⊗∂λ ∂FΩ +∂P ∂F:41(2.36) which with index notation becomes, AiJkL =∂PiJ ∂FkL =∂PiJ ∂λ ∂λ ∂FkL Ω +∂PiJ ∂Fmn δmkδnL = =eΩJeΩL∥F eΩ∥Ψ′′(λ)−Ψ′(λ) ∥F eΩ∥2·[F eΩ]i[F eΩ]k ∥F eΩ∥+ Ψ′(λ)1 ∥F eΩ∥δikΩ (2.37) It can be proven ([161]) that the mixed and material elasticity tensors are related by the following equation, AiJkL =δikSJL +FiN CNJMLFkM (2.38) “main” — 2025/2/25 — 15:02 — page 34 — #72 Chapter 2. Nonlinear continuum mechanics and hyperelasticity 34 Making use of one of the minor symmetries of C, it follows that, FiN CNJMLFkM =FiN CNJLM FkM −→ F·C·FT(2.39) Notwithstanding, the tensor quantity F·C·FTdoes not possess, in general, minor symmetries, so in the context given by (2.38) it is necessary to reflect such a transposition between the indices, A=Λ+hF·C·FTik↔L,AiJkL = ΛiJkL+hF·C·FTiiJLk ,ΛiJkL =δikSJL (2.40) where the notation k↔Lindicates transposition between the pair of indices kand L. Hence, an explicit expression for Ccan be derived through rearrangement of the previous expression, C=F−1·[A−Λ]K↔l·[FT]−1,CIJKL =F−1 Ii [A−Λ]iJlK[FT]−1 lL (2.41) 2.3. Mechanical characterization process with the SAEN model In all the studies presented in this thesis, the parameters that characterize the SAEN model were fitted to experimental rheological data corresponding to a simple shear experiment performed on collagen hydrogels. In this section, a description of this process is given, and a fitting with respect to real shear rheology data is exposed in order to illustrate that the model is able to reproduce the real material behavior. 2.3.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, F=   1γ0 010 001  (2.42) “main” — 2025/2/25 — 15:02 — page 41 — #79 Chapter 3. Inverse methods in Traction Force Microscopy 41 Figure 3.1: PBNIM’s rationale: conceptual 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. constrained minimization of the mismatch between the computed and measured displacements, that may or may not include regularization. In the case of constrained methods, the constraint is similar to the one defined by other authors (nullity of internal body forces within the ECM domain), which here is imposed by introducing the weak form of the Principle of Virtual Work as a Lagrangian field through a penalty into the corresponding functional to minimize. Both approaches consider the possibility of including a term that regularizes surface and volume forces. The generalized optimization scheme can be described in a simplified way as: Search the displacement field that, “main” — 2025/2/25 — 15:02 — page 42 — #80 Chapter 3. Inverse methods in Traction Force Microscopy 42 minimize(Data mismatch term +Equilibrium constraint +Regularization term) in which the cost function is minimized aiming at reducing the data mismatch term regarding new and measured displacement fields while at the same time enforcing equilibrium within the ECM domain and regularizing surface and body forces. Non-regularized approaches are performed by dropping the regularization term from the formulation. Unconstrained approaches result from dropping the equilibrium constraint term from the formulation. The PBNIM corresponds to the constrained implementation with no regularization (α= 0). 3.2.1. Unconstrained methods This formulation is based on searching for the optimal solution to the following problem: compute a displacement field that is as similar as possible to the measured one and to which regularization is applied. The objective (cost) function of the problem is, min u1 2||u−u∗||2 2+α 2ZΓ ||Tt·t||2 2ds +ZΩ ||Tb·b||2 2dv (3.1) where udenotes the searched displacements, u∗are the measured displacements, αis the regularization term parameter, tand bare the surface and body forces acting on surfaces Γ(cell surface and ECM walls) and interior domain Ω, respectively; and Tt(x, y, z)and Tb(x, y, z)are scalar fields bounded in the range [0,1] that weight the contributions to the objective function of the corresponding forces in specific regions of the domain. This problem has a minimum stationary solution of the objective function with respect to the unknown displacements. By means of the Gateaux derivative, the following expression is obtained, δu=0→u+αϕ(u) = u∗(3.2) in which ϕ(u)is defined in Eq. (3.2) as a nonlinear function that depends on uimplicitly through tand bas, ϕ(u) = 1 2 ∂ ∂uZΓ ||Tt·t||2 2ds +ZΩ ||Tb·b||2 2dv(3.3) The nonlinear nature of Eq. (3.3) makes the finding of a close-form solution for ϕ(u)impossible. Therefore, Eq. (3.2) is linearized in the proximity of u, “main” — 2025/2/25 — 15:02 — page 43 — #81 Chapter 3. Inverse methods in Traction Force Microscopy 43 ϕ(u+ ∆u)∼ =ϕ(u) + Dϕ(u)[∆u]=[ϕ]u+∂ϕ ∂uu ·∆u(3.4) with Dϕ(u)[∆u]being the directional derivative of ϕ(u)in the direction of a sufficiently small increment ∆u. The resulting linearized form of the equation is thus given by, u+α[ϕ(u) + Dϕ(u)[∆u]] = u∗(3.5) which will be, upon discretization, solved iteratively via a Newton-Raphson scheme. A convergence criterion under which ∆uis considered negligible is established, after which the iterative process is stopped and the resulting the desired solution uis obtained. 3.2.2. Constrained methods In the case of constrained methods, the objective function that represents the inverse problem is the same as the one of the unconstrained case, but the solution is forced to fulfill the equilibrium of internal and external forces within the hydrogel domain. The problem is then formulated as follows, min u1 2||u−u∗||2 2+α 2ZΓ ||Tt·t||2 2ds +ZΩ ||Tb·b||2 2dv s.t. Θ=0 (3.6) where Θdenotes the equilibrium condition equations. The constraint can be introduced in the objective function by means of a Lagrange multiplier η, min u1 2||u−u∗||2 2+η·Θ+α 2ZΓ ||Tt·t||2 2ds +ZΩ ||Tb·b||2 2dv (3.7) Similarly to the unconstrained case, the minimum stationary solution to the problem is found by differentiating the objective function with respect to variables uand ηas, “main” — 2025/2/25 — 15:02 — page 44 — #82 Chapter 3. Inverse methods in Traction Force Microscopy 44 δu=0→u+∂Θ ∂u·η+α·ϕ(u) = u∗(3.8a) δη=0→Θ=0(3.8b) The equilibrium condition has been chosen to be represented by the weak form of the Principle of Virtual Work (PVW), since it is intrinsically defined as an equilibrium equation such that its fulfillment ensures the equilibrium of internal (stress) forces with real and known acting external forces. Hence, η·Θ⇒δW(u, δu)(3.9) where δuis a kinematically admissible displacement field, typically called virtual displacement field. This virtual field is interpreted in our formulation as the lagrangian multiplier field η. The FEM implementation for solving the PDE-based equations is based on the software Simulia ABAQUS, which employs an Updated Lagrangian formulation of the PVW. This version of the PVW is expressed as, δW(u, δu) := ZΩ τ:δεdv −ZΓ t·δuds, δε=1 2h∇δu+ [∇δu]Ti (3.10) where τdenotes the Kirchhoff stress tensor, δuis the virtual displacement, vthe reference volume (hydrogel), δε the linear virtual strain tensor, tthe traction vector resulting from the product between the Cauchy stress tensor σ= (1/J)τand the vector normal to the surface Γ,n. By this definition, one can obtain the following nonlinear set of equations, u+∂δW ∂u+α·ϕ(u) = u∗(3.11a) δW(u, δu) = ZΩ τ:δεdv −ZΓ t·δuds = 0 (3.11b) in which ϕ(u)is defined equivalently that the unconstrained formulation. Analogously, the linearization of Eq. (3.11b) is performed via a first order Taylor expansion around u, δW(u+ ∆u, δu) = δW(u, δu) + DδW(u, δu)[∆u](3.12) This results into the following set of linearized equations, “main” — 2025/2/25 — 15:02 — page 45 — #83 Chapter 3. Inverse methods in Traction Force Microscopy 45 u+∂δW ∂u+α[ϕ(u) + Dϕ(u)[∆u]] = u∗(3.13a) δW(u, δu) + DδW(u, δu)[∆u] = 0 (3.13b) where ϕ(u)is linearized as in the unconstrained case. The directional derivative of the internal work DδW(u, δu)[∆u]can be separated into two components [160] corresponding to the work associated to internal and external forces, respectively, DδW(u, δu)[∆u] = DδWint(u, δu)[∆u]−DδWext(u, δu)[∆u](3.14) with, DδWint(u, δu)[∆u] = DZΩ τ:δεdv[∆u] = ZΩ D(τ:δε)[∆u]dv (3.15a) DδWext(u, δu)[∆u] = DZΓ t·δuds[∆u] = ZΓ D(t·δu)[∆u]ds (3.15b) In the Updated Lagrangian configuration, the term corresponding to the internal forces is developed by means of the Jaumann or corotational derivative of the Kirchhoff stress tensor: DδWint(u, δu)[∆u] = =ZVh∇δu : [J c J− H ]:∇δu+τ:h[∇δu]T∇∆uiidV (3.16) where c Jis the Jaumann elasticity tensor, which can be computed as, J c J=J c +H, Hijkl =1 2[δjkτil +δikτjl +δjlτik +δilτjk](3.17) being c the so called Fourth Cauchy Elasticity Tensor, c =1 Jϕ∗[C], cijkl =1 JFiIFjJ FkKFlLCIJKL (3.18) “main” — 2025/2/25 — 15:02 — page 46 — #84 Chapter 3. Inverse methods in Traction Force Microscopy 46 where Einstein’s summation convention has been used. Cis the Second Elasticiy Tensor (material configuration), given by the formula, C= 4 ∂2Ψ ∂C∂C(3.19) where Ψdenotes the strain energy density function, and Cis the deformation gradient tensor in the material configuration. Regarding Eq. (3.15b), the integration of the surface forces is treated by means of a parametrization of the surface specific to the geometry of the problem. For a more detailed description, the reader is referred to Ref. [160]. 3.2.3. Finite element discretization In this section, the numerical implementation of Eqs. (3.5) and Eqs. (3.13), based on their discretization within a FE framework, is explained in detail. First, the field variables are discretized such that u≈Ne·ui,∆u≈ Ne·∆ui,δu≈Ne·ηi,u∗≈Ne·u∗,i,t≈Ne·tiand b≈Ne·bi, where Neis the matrix that contains the shape functions in a certain element (e). On the other hand, ui,∆ui,u∗,i,ηi,tiand biare vectors whose components correspond to the discrete values of u,∆u,u∗,η,tand bat node positions i. Then, the problem domain is discretized into finite elements, such that RΩ•=PNel e=1 RΩel •, being Ωel the domain corresponding to finite element (e)and Nel the total number of elements that compose the mesh. After discretization, surface forces tand body forces bare represented as nodal forces with values contained in vector F, displayed according to the degrees of freedom of the problem. As a consequence, scalar weight fields Tt(x, y, z)and Tb(x, y, z)take the form of a square diagonal matrix Twith size the number of degrees of freedom of the problem, displaying values contained in the range [0,1] in the positions of the diagonal corresponding to each specific degrees of freedom. In this case, Ttakes values of 0 or 1 only, depending on which nodes are set to contribute to the objective function, i.e., which nodal forces positions are regularized. Due to this, the regularization term in Eqs. (3.1) and (3.6) takes the expression, α 2ZΓ ||Tt·t||2 2ds +ZΩ ||Tb·b||2 2dv−→ α 2||T·F||2 2(3.20) and ϕ(u)is rewritten as, “main” — 2025/2/25 — 15:02 — page 47 — #85 Chapter 3. Inverse methods in Traction Force Microscopy 47 ϕ(u) = 1 2 ∂ ∂u||T·F||2 2=∂F ∂u⊺ ·T⊺·T·F(3.21) which is linearized in the proximity of uin the following fashion, ϕ(u+ ∆u) = ϕ(u) + Dϕ(u)[∆u] = ∂F ∂u⊺ ·T⊺·T·Fu + +∂F ∂u⊺ ·T⊺·T·∂F ∂u+∂ ∂u∂F ∂u⊺ ·T⊺·T·Fu ·∆u(3.22) Both unconstrained and constrained approaches are also distinguished during the discretization process. Unconstrained The finite element discretization of domain and field variables yield the following expression for the linearized Eq. (3.5), ui k+1 +α Nel X e=1 hϕ(e)(ui k) + Dϕ(e)(ui k)[∆ui]i=u∗,i (3.23) in which ∆ui=ui k+1 −ui k, and ui kdenotes the value for the solution at the k-th iteration within the quasi-Newton-Raphson procedure. The different sum terms are developed as follows (see additional details in Ref. [160]), Nel X e=1 ϕ(e)(ui k) = Nel X e=1 ∂F ∂u⊺ ·T⊺·T·F(e) uk =ANel e=1 hK⊺ k·T⊺·T·Fint ki(e) (3.24) Nel X e=1 Dϕ(e)(ui k)[∆ui] = = Nel X e=1 ∂F ∂u⊺ ·T⊺·T·∂F ∂u+∂ ∂u∂F ∂u⊺ ·T⊺·T·F(e) uk ·∆ui= =ANel e=1 K⊺ k·T⊺·T·Kk+∂ ∂u[K]⊺k ·T⊺·T·Fint k(e) ·∆ui(3.25) where ANel e=1 is the finite element assembly operator, and Fint kand Kkare, respectively, the element vector of internal forces and the element tangent stiffness matrix, both evaluated at the k-th iteration. “main” — 2025/2/25 — 15:02 — page 48 — #86 Chapter 3. Inverse methods in Traction Force Microscopy 48 Since Kimplicitly depends on u, a closed-form for the expression of its derivative is not available. Due to this, an approximation is made by neglecting the contribution of second order derivative terms, Nel X e=1 Dϕ(e)(ui k)[∆ui]≈ANel e=1 K⊺ k·T⊺·T·Kk(e)·∆ui(3.26) Note that, upon convergence of the solution, this approximation does not impact the validity of the result, just the rate of convergence of the iterative procedure. Finally, after assembly, Eqs. (3.26) and (3.24) are introduced into (3.23): [1+α·K⊺ k·T⊺·T·Kk]·ui k+1 =u∗,i +α·K⊺ k·T⊺·T·[Kk·ui k−Fint,i k](3.27) where 1is the identity matrix. Constrained Upon discretization, the system of equations is expressed as, ui k+1 + Nel X e=1 ∂δW ∂u(e) uk ·ηi k+1 +α Nel X e=1 hϕ(e)(ui k) + Dϕ(e)(ui k)[∆ui]i=u∗,i (3.28a) Nel X e=1 δW(e)(ui k, δui) + Nel X e=1 DδW(e)(ui k, δui)[∆ui]=0 (3.28b) Similarly to the unconstrained case, assembly of the previous equations is performed as in Eqs. (3.24) and (3.26), ui k+1 +ANel e=1K(e) k·ηi k+1 +α·ANel e=1 hK⊺ k·T⊺·T·Fint ki(e)+ +α·ANel e=1 K⊺ k·T⊺·T·Kk(e)·∆ui=u∗,i (3.29) ANel e=1(F(e),int k−F(e),ext k) + ANel e=1K(e) k·∆ui= 0 (3.30) which yields the following system of equations, “main” — 2025/2/25 — 15:02 — page 49 — #87 Chapter 3. Inverse methods in Traction Force Microscopy 49 [1+α·K⊺ k·T⊺·T·Kk]·ui k+1 +Kk·ηi k+1 =u∗,i+ +α·K⊺ k·T⊺·T·[Kk·ui k−Fint,i k](3.31) Kk·ui k+1 =Fext,i k−Fint,i k+Kk·ui k(3.32) Upon defining ξ≡K⊺ kT⊺TKkand χ≡K⊺ kT⊺T, Eqs. (3.31)-(3.32) can be written in matrix form as, "1+α·ξKk Kk0#·"ui k+1 ηi k+1#="u∗,i +α·χ·[Kk·ui k−Fint,i k] Fext,i k−Fint,i k+Kk·ui k#(3.33) Now, the system is rearranged according to cell boundary nodes Aand hydrogel internal nodes Bas follows,      1+α·ξ(A, A)α·ξ(A, B)Kk(A, A)Kk(A, B) α·ξ(B, A)1+α·ξ(B, B)Kk(B, A)Kk(B, B) Kk(A, A)Kk(A, B)0 0 Kk(B, A)Kk(B, B)0 0      ·     ui k+1(A) ui k+1(B) ηi k+1(A) ηi k+1(B)      = =      "1 0 0 1#·"u∗,i(A) u∗,i(B)# Fext,i k(A)−Fint,i k(A) Fext,i k(B)−Fint,i k(B)       +     0 0 "Kk(A, A)Kk(A, B) Kk(B, A)Kk(B, B)#·"ui k(A) ui k(B)#     + +     "α·ξ(A, A)α·ξ(A, B) α·ξ(B, A)α·ξ(B, B)#·"ui k(A) ui k(B)# 0 0      − −     "α·χ(A, A)α·χ(A, B) α·χ(B, A)α·χ(B, B)#·"Fint,i k(A) Fint,i k(B)# 0 0      (3.34) The term of external forces Fext,i k(A)in the third row of the right hand side of Eq. (3.34) represent the nodal reaction forces exerted by the cell at the cell-hydrogel boundary nodes and it constitutes an unknown variable of “main” — 2025/2/25 — 15:02 — page 50 — #88 Chapter 3. Inverse methods in Traction Force Microscopy 50 the system. The value of the associated Lagrange multiplier is prescribed as ηi(A) = 0. This is equivalent to imposing that the equilibrium of internal forces with known external forces must be fulfilled. In fact, nodal reaction forces in the interior of the hydrogel must be zero Fext,i(B) = 0 in the absence of any body or active forces acting on the hydrogel. Eq. (3.34) can be thus rewritten as,    1+α·ξ(A, A)α·ξ(A, B)Kk(A, B) α·ξ(B, A)1+α·ξ(B, B)Kk(B, B) Kk(B, A)Kk(B, B)0  ·   ui k+1(A) ui k+1(B) ηi k+1(B)  = =   u∗,i(A) + α·ξ(A, A)·ui k(A) + α·ξ(A, B)·ui k(B) u∗,i(B) + α·ξ(B, A)·ui k(A) + α·ξ(B, B)·ui k(B) −Fint,i k(B) + Kk(B, A)·ui k(A) + Kk(B, B)·ui k(B)  − −   α·χ(A, A)·Fint,i k(A) + α·χ(A, B)·Fint,i k(B) α·χ(B, A)·Fint,i k(A) + α·χ(B, B)·Fint,i k(B) 0   (3.35) The systems defined in Eqs. (3.27) and (3.35) are numerically implemented following a Newton-Raphson iterative scheme. “main” — 2025/2/25 — 15:02 — page 57 — #95 Chapter 4. Holistic 3D TFM study with real cell morphologies 57 complex (closer to a spherical shape) cell geometries. The solidity value of each cell is plotted in Fig. 4.3. Figure 4.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). 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. 4.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 prin- “main” — 2025/2/25 — 15:02 — page 58 — #96 Chapter 4. Holistic 3D TFM study with real cell morphologies 58 Figure 4.3: Solidity index for selected cell morphologies: starshaped, spread, protrusive and spherical morphologies, in order of increasing solidity. “main” — 2025/2/25 — 15:02 — page 59 — #97 Chapter 4. Holistic 3D TFM study with real cell morphologies 59 cipal strains vector) ranged from 48% to 54% (see Table 4.1, nonlinear case) using the nonlinear model defined in Section 2 of Chapter 2. In Table 4.1, the values for the large cellular pulling force magnitude prescribed for each cell lie within the range of nN found in the literature [166]. Moreover, the resulting displacements are of the order of microns, as typically seen in the laboratory in TFM [15]. 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 4.1 also shows maximum strains and averaged maximum strains achieved in the different protrusions of the cells. Fig. 4.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 4.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. 4.2). Note that the range of values shown in Fig. 4.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. Protrusion # Force [pN] Average Maximum Mean average Mean max. Maximum max. Spread 1 2.25e3 0.1870 0.5631 0.1592 0.4884 0.653120.1215 0.2490 30.1690 0.6531 Protrusive 1 0.375e3 0.1218 0.1838 0.1064 0.2412 0.6476 20.0795 0.1048 30.1103 0.1803 40.1125 0.2103 50.0814 0.1187 60.1331 0.6476 Spherical 10.75e30.2824 0.5607 0.2158 0.3781 0.6570 20.1491 0.1955 Star-shaped 1 0.67e3 0.1303 0.2065 0.1334 0.2602 0.4863 20.1457 0.2099 30.1390 0.4863 40.1086 0.1638 50.1660 0.3173 60.1109 0.1772 Table 4.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). Hydrogel’s constitutive behavior “main” — 2025/2/25 — 15:02 — page 60 — #98 Chapter 4. Holistic 3D TFM study with real cell morphologies 60 Figure 4.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. “main” — 2025/2/25 — 15:02 — page 61 — #99 Chapter 4. Holistic 3D TFM study with real cell morphologies 61 A fibered nonlinear model was selected to describe the nonlinear mechanical behavior of the hydrogel (ECM of the in vitro setups). This model assumes a purely hyperelastic behavior and has been proven as an excellent fitting for collagen matrices in shear rheology tests (see [105] and Section 3 of Chapter 2). While previous studies have addressed the importance of viscoelasticity in cell behavior [79, 167, 168] and implemented viscoelastic models for TFM [131], the selected model neglects viscoelastic effects, as it was verified in previous studies for similar collagen hydrogels to the ones used in this chapter [105], 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 Section 2 of Chapter 2. Collagen matrices of 1.2 mg/ml concentration were prepared and tested in a shear rheology device following the protocol described in [115]. Fig. 4.5 shows the experimental results of shear rheology tests, as well as the best fit for the selected nonlinear model (see Section 3 of Chapter 2 for details about the fitting process). It can be observed, that the model accurately represents the real behavior of the hydrogel up to 50% of shear strain. Table 4.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 modulus through the expression 2G(1+ν). In the literature, one can find Poisson’s ratio values for collagen hydrogels ranging from 0.2 to 0.48 (see [106, 139, 148, 158, 169, 170]). 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. Model Parameter Value Nonlinear κo[Pa] 156.6217 do[–] 0.375 λs[–] 0.0111 ds[–] 0.0836 Linear E[Pa] 29.5076 ν[–] 0.340 Table 4.2: Fitted nonlinear and linear model parameters. Parameters E and νof the nonlinear model stand for the Young’s modulus and Poisson’s ratio, respectively (see Section 2 in Chapter 2 for the definition of the nonlinear model parameters). The assumed linear and nonlinear models in the different cases of anal- “main” — 2025/2/25 — 15:02 — page 62 — #100 Chapter 4. Holistic 3D TFM study with real cell morphologies 62 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 Engineering shear strain ( ) 0 5 10 15 20 25 30 First Piola-Kirchhoff stress (P12) [Pa] Experimental curve Nonlinear model Linear model Figure 4.5: Determination of model parameters with experimental data: experimental curve of 1.2 mg/ml collagen hydrogel corresponding to shear rheology tests is fitted with both the selected nonlinear and linear elastic constitutive models. The shaded blue area represents the variability of the different tested specimens. The solid blue line curve is the average of all specimens. ysis consider the linear and nonlinear fittings shown in Fig. 4.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 Sections 2.1, 2.2, Section 2 of Chapter 2; and Chapter 3. On the other hand, the case of the linear model uses the linearized (infinitesimal strains) versions of the forward and inverse formulations [112]. 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. 4.1, using the software ABAQUS. These displacement fields, and associated strain and traction variables, are referred to as the ground truth (GT) solutions. “main” — 2025/2/25 — 15:02 — page 63 — #101 Chapter 4. Holistic 3D TFM study with real cell morphologies 63 4.3.2. Displacement reconstruction Similarly to the TFM experimental methodology, GT hydrogel displacements were sampled for each case at discrete random (bead) positions. For all the cases of analysis, an experimentally feasible bead density of approximately 0.03 beads/µm3was assumed, which lies within the order of magnitude of experimental procedures [112]. Table 4.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. 4.6 shows a schematic representation of this process. The measured displacement field is referred to as the corrupted input displacement field for the forward and inverse methods simulations. LX[µm]LY[µm]LZ[µm]Nbeads 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 Table 4.3: Size of the hydrogel domain, number of beads (Nbeads), interbead distance for each cell, and size of the FE mesh. The inter-bead distance was calculated as the mean nearest neighbor Euclidean distance. 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. 4.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 dis- “main” — 2025/2/25 — 15:02 — page 64 — #102 Chapter 4. Holistic 3D TFM study with real cell morphologies 64 Figure 4.6: Generation of the corrupted input displacement field: conceptual description of the process to obtain corrupted displacement fields from ground truth simulations. Errors between ground truth and corrupted displacements were amplified for representation purposes. placement 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 4.3. The convergence of the results with the selected meshes was previously checked. 4.3.4. Error indicators In this study the focus is put 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, the forward and inverse traction error indicators Tfwd e,Tinv eare defined as follows, “main” — 2025/2/25 — 15:02 — page 65 — #103 Chapter 4. Holistic 3D TFM study with real cell morphologies 65 Tinv e= 100 ·1−corr[tinv pdf, tGT pdfi (4.10a) Tfwd e= 100 ·1−corr[tfwd pdf , tGT pdfi (4.10b) where tinv pdf and tfwd pdf are the probability density functions of the magnitudes of the tractions distribution tn i(ibeing the boundary nodes of the cell) for the inverse and forward solutions, respectively. tGT pdf has the same definition but for the ground truth solution. Finally, corr[f, g]is the correlation coefficients of functions fand g. 4.4. Results The results referred to the different analyses performed in this study (see Fig. 4.1) are summarized in this section. The magnitude of tractions on the boundary of the selected cells are represented in Figs. 4.7–4.9 for the low cellular pulling force case, including the ground truth solution (either assuming a linear or nonlinear behavior of the hydrogel) and the reconstructed tractions solution for each considered matrix behavior. Specifically, Figs 4.7 and 4.8 show the results of reconstructing tractions using the correct (linear and nonlinear, respectively) constitutive behavior (i.e. the behavior used to generate the ground truth, see Fig. 4.5) by means of the forward and the inverse methods. Finally, Fig. 4.9 shows the ground truth tractions for a nonlinear matrix, compared to the forward and inverse reconstructed solutions assuming a linear one. Qualitatively, a better performance of the inverse methodology versus the forward one, can be seen in Figs. 4.7–4.9 when compared with their corresponding ground truth solutions. Fig. 4.10 summarizes the errors in the reconstruction of tractions, according to the error indicators defined in Eq. (11), for the different cell morphologies, cellular pulling force levels and methodology (forward and inverse). It quantifies the superior performance of the inverse method for all the analyzed scenarios. The mean values of the traction errors, averaged for the considered pulling force levels, are represented in Fig. 4.11 versus the solidity index (as a measure of the complexity of the cell geometry) for the forward and inverse methods. Interestingly, it is shown that traction reconstruction becomes more challenging in complex geometries (far from a sphere shape and hence low solidity). A good correlation between the defined traction error indicator and the solidity index was found. Moreover, “main” — 2025/2/25 — 15:02 — page 66 — #104 Chapter 4. Holistic 3D TFM study with real cell morphologies 66 the outperformance of the inverse versus the forward method is also shown in Fig. 4.11. Mean values of traction errors, CPU time and numerical iterations during the analysis, are represented in Table 4.4. All the simulations were performed in a standard laptop (AMD Ryzen 7 4800H 2.90 GHz, 16GB RAM). All these values were averaged for the considered pulling force levels and are given for each analyzed cell morphology, hydrogel behavior and methodology (forward and inverse). Moreover, the CPU time is plotted versus traction error (mean values) in Fig. 4.12 for each considered case. It is observed that the inverse method shows a better accuracy than the forward one, for all the analyzed scenarios, but with a higher CPU time cost. On the other hand, the fact of assuming a linear behavior of the matrix for a nonlinear ground truth behavior, provides a faster but less accurate solution. Forward Inverse Model Forward time (s) Forward error (%) Iters (avg) Inverse time (s) Inverse error (%) Spread Linear 165 55.9 1 559 18.7 Nonlinear 233 31.2 3.1 1849 7.7 Protrusive Linear 112 55.8 1 234 15.1 Nonlinear 156 20.7 2.7 721 4.8 Spherical Linear 72 15.7 1 116 1.4 Nonlinear 96 1.2 4.3 594 0.48 Star-shaped Linear 124 74.7 1 308 30.5 Nonlinear 172 36.8 3.0 993 9.3 Table 4.4: Mean CPU time and mean traction error indicators (averaged over the three analyzed cellular pulling force level cases) of recovered traction solution assuming either linear or nonlinear behaviors of the matrix. Traction errors were calculated with respect to the ground truth solution for the case of the considered nonlinear matrix. 4.5. Discussion This section elaborates on the accuracy of traction reconstruction from the different perspectives analyzed in this study: the magnitude of the cellular pulling force, the effect of the complexity of the cell morphology, accuracy and efficiency of the forward vs the inverse method, and the effect of reconstructing tractions with a correct/incorrect behavior of the hydrogel (linear/nonlinear). “main” — 2025/2/25 — 15:02 — page 73 — #111 Chapter 4. Holistic 3D TFM study with real cell morphologies 73 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. [112] 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. 4.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. 4.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 4.11b and 4.11c). 4.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. 4.10 for all the analyzed cases. Also, these data are plotted in Fig. 4.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. 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 [135]. The inverse method requires the setup of a matrix system, and an iterative process in a nonlinear case [135]. Table 4.4 and Fig. 4.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 233s 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 “main” — 2025/2/25 — 15:02 — page 74 — #112 Chapter 4. Holistic 3D TFM study with real cell morphologies 74 4.4 and Fig. 4.12). 4.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. 4.10a and 4.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, [163] found that although the use of a strain-stiffening material provides higher quality traction reconstructions, errors are of similar magnitude in both materials. In this case, the explanation for the last point is found in Fig. 4.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 Generalized inverse formulation), 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. 4.10b shows. These regions are particularly represented in the spherical cell, as illustrated in Fig. 4.13. Fig. 4.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. 4.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.4 and Fig. 4.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 spehere [105]. 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.4 and Fig. 4.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.4), CPU times increases to ratios up to 3.3 (worst case), according to “main” — 2025/2/25 — 15:02 — page 75 — #113 Chapter 4. Holistic 3D TFM study with real cell morphologies 75 Table 4.4 and Fig. 4.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. 4.10c shows the traction error for this case, and Fig. 4.11c the error averaged for the considered pulling force levels. As expected, it is observed in Fig. 4.11c that the linear assumption for the hydrogel performs worse in traction reconstruction. Specifically, for the toughest case i.e. the star-shaped 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. 4.11a and 4.11c). Focusing on the small cellular pulling force case, such that geometrical nonlinear effects may be neglected, the differences found in Figs 4.10a and 4.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. 4.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. 4.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 “main” — 2025/2/25 — 15:02 — page 76 — #114 Chapter 4. Holistic 3D TFM study with real cell morphologies 76 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 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. “main” — 2025/2/25 — 15:02 — page 77 — #115 Chapter 4. Holistic 3D TFM study with real cell morphologies 77 Figure 4.11: Mean traction error indicator (for analyzed cellular pulling force level cases) of recovered tractions solution: (a) Assuming a nonlinear behavior of the matrix versus ground truth solution of the associated nonlinear matrix. Forward fitting: R2= 0.974. Inverse fitting: R2= 0.988. (b) Assuming a linear behavior of the matrix versus ground truth solution of the associated linear matrix. Forward fitting: R2= 0.986. Inverse fitting: R2= 0.969. (c) Assuming a linear behavior of the matrix versus ground truth solution of the considered nonlinear matrix. Forward fitting: R2= 0.949. Inverse fitting: R2= 0.973. “main” — 2025/2/25 — 15:02 — page 78 — #116 Chapter 4. Holistic 3D TFM study with real cell morphologies 78 Figure 4.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. “main” — 2025/2/25 — 15:02 — page 79 — #117 Chapter 4. Holistic 3D TFM study with real cell morphologies 79 Figure 4.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. “main” — 2025/2/25 — 15:02 — page 80 — #118 Chapter 4. Holistic 3D TFM study with real cell morphologies 80 “main” — 2025/2/25 — 15:02 — page 81 — #119 Chapter 5 Multiscale framework with heterogeneous matrices 5.1. Introduction The previous chapter emphasized the great importance lying in the correct characterization of the ECM-mimicking material employed in 3D TFM experiments. As it was described in the introductory chapter, there are few TFM studies in which heterogeneity of the ECM is considered. In [171] the authors implement a methodology for solving 3D TFM cases while considering ECM heterogeneity, although not due to ECM degradation. Nevertheless, the applicability of the methodology is limited, since they consider small strains and assume the material to behave as linear elastic. Also, the authors utilize a forward method to perform traction reconstruction, which, as it was concluded in the study presented in Chapter 4 (see also [133, 135]), exhibits a significantly reduced efficacy when compared to inverse methods. Song et al. [107] proposed an in silico study in which the authors perform 3D TFM while accounting for heterogeneity of matrix mechanical properties, while considering nonlinearity of the ECM mechanical behavior and conducting traction reconstruction through and inverse method. However, the law that correlates ECM degradation with the material parameters is not defined by the authors, but obtained through the inverse methodology simultaneously to tractions (this is the context of the work shown in Chapter 7 of this thesis), which represents an additional layer of computational complexity and may render it not applicable to certain data sets. Additionally, the degraded region is imposed by the authors and limited to the vicinity of the cell surface, which may not represent a “main” — 2025/2/25 — 15:02 — page 82 — #120 Chapter 5. Multiscale framework with heterogeneous matrices 82 realistic case. In this chapter, a complete 3D TFM framework which considers ECM degradation retrieved from the solution of a set of equations belonging to a previously described MMP-induced degradation mathematical model [157] is presented (the referred paper was briefly commented on in the introduction to this thesis). This model was adapted for its application to individual complex cell geometries in 3D. Degradation is then correlated with a nonlinear hyperelastic model, that was calibrated with synthetic shear rheology tests corresponding to different levels of ECM degradation. The shear rheology curves were obtained from computational (virtual) simple shear tests applied to synthetically generated fibered matrices. Shear rheology tests corresponding to real collagen hydrogels tested in the laboratory, were used to fit and calibrate such a virtual tool. As a result, TFM accuracy is assessed in these heterogeneous matrices referred to traction reconstruction versus a (controlled) synthetic ground truth solution in which degradation of the matrix is included. The analysis of cellular forces in several scenarios is then conducted: under the presence and absence of ECM degradation, and forward / inverse reconstructions for each case of study. In contrast to the methodology presented in [171], the one described herein accounts for finite strains and nonlinear constitutive behavior of the ECM; and tractions are recovered by means of both forward and inverse methods. Regarding [107], the approach concerned in this chapter considers the simulation of realistic time-dependent ECM density maps by implementing a thorough mathematical model based on enzymatic reactions in 3D TFM within the framework of MMP-induced ECM degradation. Moreover, the methodology explicitly takes into consideration the ECM microstructure to establish the relationship between ECM degradation and corresponding material parameters. The next sections include: (i) a brief description of the study, (ii) the mathematical model for MMP-induced degradation of the ECM and its computational implementation, (iii) the virtual shear rheology tests and their corresponding synthetic data using an in-house fibered matrix generator tool, and (iv) the multiscale assembly of these methodologies in the forward/inverse TFM simulations. Finally, the obtained results are discussed and the main conclusions of the study are highlighted at the end of the chapter. “main” — 2025/2/25 — 15:02 — page 89 — #127 Chapter 5. Multiscale framework with heterogeneous matrices 89 applied (and equivalently, the previously described boundary conditions) are represented in Figure 5.4. Figure 5.4: Different steps comprising the computational implementation of the model in COMSOL Multiphysics: geometry, mesh, MMP-emitting regions (labelled in blue), and definition of prescribed field variable cancer density (c(x, y, z)), which triggers the emission of MMPs in the selected MMP-emitting regions. Results A degradation process throughout 28 dimensionless time steps (from t=0 to t=27) was simulated. The values of the parameters that govern the proposed degradation model (Eqs. (5.1a)-(5.1f)) were fitted in order to “main” — 2025/2/25 — 15:02 — page 90 — #128 Chapter 5. Multiscale framework with heterogeneous matrices 90 obtain a maximum ECM degradation level of 50% (ECM density level of 50%) at the end of the degradation process at t= 27. Since the scope of this study is computational, the calibrated parameters do not necessarily coincide with values found in the literature, although some of them remain unchanged with respect to the ones used in [157, 172], which are either extracted from the literature or estimated by the authors. Table 5.1 contains the fitted nondimensional values for the parameters. Parameters δ1,δ2,µvand αv3were calibrated by means of a trial and error procedure. Dv2and Dv4were modified only with respect to their non-dimensional values according to the dimensions of this particular problem, that are different to those presented in [157, 172]. The resultant ECM density maps throughout the simulated degradation process are shown in Figure 5.5. Only four representative time snapshots are shown: t= 0,9,18 and 27 (dimensionless time). These degradation levels were calibrated with the experimental measurements presented in [173]. In this work, the authors provide data corresponding to the volume of DQ collagen (which becomes fluorescent while being degraded) per cell within a period of 24 hours. By computing the volume of each tetrahedron of the finite element mesh of the problem, the total degraded volume at t= 9 (equivalent to approximately 1 day, according to the dimensionless definition of this variable [157]), is obtained performing the sum of the degraded volume fractions of the tetrahedra for the corresponding ECM density map. The computed quantity was 8456 µm3, which lies within the range given in the referred study (1000 −10000 µm3/cell). The referred experimental study was established at similar circumstances than the ones presented in this study, although not strictly the same. Therefore, this calibration should be understood as a qualitative estimation to show the applicability of the degradation model, without compromising its validity. Figure 5.6 plots variables v2,v3,v4and v5of the degradation model at t= 27. 5.4. Mechanical characterization of the hydrogel including degradation This section describes the multiscale approach to obtain the (virtual) macroscopic behavior of fibered hydrogels, as well as the tool developed for producing randomly generated interconnected fiber network microstructures. This methodology is used to generate synthetic shear rheology tests corresponding to different levels of ECM density (variable v1in the degra- “main” — 2025/2/25 — 15:02 — page 91 — #129 Chapter 5. Multiscale framework with heterogeneous matrices 91 Figure 5.5: Results of the simulated degradation process: ECM (dimensionless) density maps (XY cross midsection) for 4 different timepoints (t= 0,9,18 and 27). “main” — 2025/2/25 — 15:02 — page 92 — #130 Chapter 5. Multiscale framework with heterogeneous matrices 92 Figure 5.6: Results of the simulated degradation process: variables v2,v3,v4and v5(XY cross midsection) at t= 27. These variables correspond to the concentrations of MMP-2, MT1-MMP, TIMP2 and the generic intermediate complex, respectively. “main” — 2025/2/25 — 15:02 — page 93 — #131 Chapter 5. Multiscale framework with heterogeneous matrices 93 Parameter Dimensionless value δ10.02 δ20.1 µv0.002 Dv2129 k21 5 k22 0.195 βv20.1 αv31.55 k31 27.4 k32 2 βv30.1 Dv4129 αv44 k41 5 k42 27.4 k43 2 k51 27.4 k52 0.195 k53 2 Table 5.1: Calibrated nondimensional values for the parameters used in the simulations. dation model previously described, corresponding v1= 1 and v1= 0 to 100% and 0% of the initial ECM density, respectively). The explanation given here is contextual to the study concerned. A more detailed description of the model is given in the Appendix to this chapter. Only heterogeneity of the mechanical properties of the medium due to matrix degradation is considered, and therefore all the synthetically generated fiber microstructures are considered to behave homogeneously, at the initial time of analysis. 5.4.1. Fibered matrix generator and shear rheology tests The generator creates randomly distributed and interconnected fiber networks and is aimed at replicating the microstructure of fibered materials, such as collagen hydrogels. Fibers are constructed from sequences of nonlinear beam finite elements. The main input parameters of the generator “main” — 2025/2/25 — 15:02 — page 94 — #132 Chapter 5. Multiscale framework with heterogeneous matrices 94 are the fiber volumetric density, length of a single fiber and the curvature of a single fiber (see the Appendix to this chapter). These parameters were calibrated by inspecting Second-Harmonic Generation (SHG) images of real collagen hydrogels obtained in the laboratory. Values representative of the average were retrieved from the segmentation of a set of fibers, ensuring that they remained within the ranges found in the literature. Figure 5.7 shows a comparison between a real collagen hydrogel microstructure and the synthetically generated one. Figure 5.7 (left) shows a SHG image of a real collagen hydrogel, which correspond to a cross section of a 3D sample of the hydrogel. Equivalently, Figure 5.7 (right) shows a cross section of the reproduced 3D microstructure using the in-house fiber matrix generator. Table 5.2 includes the values employed in the matrix generation along with those that are typically present in the literature. Figure 5.7: Calibration of the fibered matrix generator parameters: real and generated (synthetic) microstructures of 1.2 mg/ml collagen hydrogel. The generated matrices are then subjected to computational (virtual) simple shear tests, by using the referred multiscale approach (see the Appendix at the end of this chapter). 1.2 mg/ml collagen type I hydrogels were prepared following the protocol described in [115]. Briefly, a mix is prepared containing 400 µL rat tail collagen (collagen R, 1.2 mg/mL; “main” — 2025/2/25 — 15:02 — page 95 — #133 Chapter 5. Multiscale framework with heterogeneous matrices 95 Symbol Definition Value Range found in the literature Efiber Elastic modulus of a single fiber 4.3[MPa] Tens-hundreds of MPa [174, 175] Dfiber Diameter of a single fiber 600 [nm] Hundreds of nm [174, 176] Lfiber Length of a single fiber 10 [µm] Tens of µm [177] Table 5.2: Calibrated parameters for matrix generation and virtual shear rheology tests. Matrix Bioscience), 400 µL bovine skin collagen (collagen G, 4 mg/mL; Matrix Bioscience), 191 µL NaHCO3 (23 mg/mL; Supplier), 191 µL 10x DMEM (Biochrom), 802 µL distilled water, and 14.37 µL NaOH (1 M) to adjust the pH to 10. The resultant 2 mL final collagen solution is pipetted into a 35-mm petri dish and polymerized in an incubator at 37 °C, 95% relative humidity, and 5% CO2 for 1 h. After polymerization, 2 mL of PBS is added to prevent dehydration of collagen gels during microscopy imaging. The parameters of the fiber model were calibrated from real shear rheology tests applied to these collagen hydrogels. Figure 5.8 shows a comparison between the experimental data and the calibrated synthetic one. Once the values of the parameters are set, the degraded matrices are generated by reducing the fiber volumetric density to the considered values (Figure 5.9), and applying the virtual simple shear tests over these degraded microstructures. 5.4.2. Nonlinear hyperelastic model Once the synthetic shear rheology data have been generated for the different levels of ECM density, the resulting curves are used to calibrate the parameters of a macroscopic nonlinear hyperelastic model [105]. This model assumes a purely hyperelastic behavior and takes into account the fibered nature of collagen hydrogels. This model is the one used in all the works presented in this thesis for representing the hyperelastic behavior of the concerned collagen hydrogels at the macroscopic scale. For a thorough description of its formulation and computational implementation, the reader is referred to Sections 2 and 3 of Chapter 2 of this document (or alternatively [133]). The parameters of the model are first fitted to the synthetically generated simple shear curve corresponding to a virgin (non-degraded) ECM (100% initial density, or v1= 1, Table 5.3). Then, just the parameter κois considered as a fitting parameters to calibrate the model to the synthetic curves corresponding different levels of ECM density. The values contained “main” — 2025/2/25 — 15:02 — page 96 — #134 Chapter 5. Multiscale framework with heterogeneous matrices 96 0 0.05 0.1 0.15 0.2 0.25 0.3 Engineering shear strain ( ) 0 5 10 15 20 25 30 35 40 45 First Piola-Kirchhoff stress (P12) [Pa] 1.2 mg/ml Synthetic data Real data Figure 5.8: Validation of the fibered matrix generator with experimental data: experimental and fitted synthetic shear rheology data for a 1.2 mg/ml collagen hydrogel with 100% ECM density (v1= 1). “main” — 2025/2/25 — 15:02 — page 97 — #135 Chapter 5. Multiscale framework with heterogeneous matrices 97 Figure 5.9: The progressive effect of degradation on the ECM’s microstructure illustrated: different degraded ECM microstructures, synthetically obtained using the developed fiber generator (see the Appendix for more details). “main” — 2025/2/25 — 15:02 — page 98 — #136 Chapter 5. Multiscale framework with heterogeneous matrices 98 in the resulting dataset (κoand ECM density) are interpolated in order to extract the corresponding fiber stiffness to the desired ECM density values in the associated maps. This strategy is followed to avoid expensive multiscale FE2simulations for each specific density level of the ECM. Figure 5.10 shows the synthetically generated shear rheology curves along with their counterparts produced with the hyperelastic model, where one can appreciate the degree of accuracy with which the model approximates well the mechanical behavior of collagen hydrogels. Moreover, Table 5.4 includes the considered values of ECM density (percentages) and single fiber stiffness. Parameter Value κo[Pa] 420.16 do[–] 0.375 λs[–] 0.0219 ds[–] 0.0323 Table 5.3: Fitted values of the model parameters for the macroscopic hyperelastic nonlinear model in the absence of degradation (non-degraded ECM). ECM density κo[Pa] 100% 420.16 89.9% 344.60 79.8% 304.98 69.6% 275.77 59.4% 249.51 49.4% 209.81 38.9% 163.54 28.5% 130.22 17.8% 69.50 7.71% 32.08 Table 5.4: Parameter κo(single fiber stiffness) of the macroscopic hyperelastic nonlinear model, for different levels of ECM density (percentages). “main” — 2025/2/25 — 15:02 — page 105 — #143 Chapter 5. Multiscale framework with heterogeneous matrices 105 Figure 5.13: In silico 3D TFM results (DEG-DEG): ground truth and both forward and inverse TFM reconstruction results (magnitude of displacements, Frobenius norm of the logarithmic strain tensor, and magnitude of tractions) considering ECM degradation (t= 27). Insets in the figures represent a detail of the variables localized in the tips of the cell. “main” — 2025/2/25 — 15:02 — page 106 — #144 Chapter 5. Multiscale framework with heterogeneous matrices 106 Figure 5.14: In silico 3D TFM results (DEG-NODEG): ground truth and both forward and inverse conventional (homogeneous) TFM reconstruction results (magnitude of displacements, Frobenius norm of the logarithmic strain tensor, and magnitude of tractions) neglecting ECM degradation (t= 27). Insets in the figures represent a detail of the variables localized in the tips of the cell. “main” — 2025/2/25 — 15:02 — page 107 — #145 Chapter 5. Multiscale framework with heterogeneous matrices 107 t = 0 t = 9 t = 18 t = 27 10 15 20 25 30 35 40 Traction error [%] DEG-DEG Inverse method Forward method (a) Ground truth cases considering degradation - Traction reconstruction considering ECM degraded mechanical properties. t = 0 t = 9 t = 18 t = 27 10 20 30 40 50 60 70 Traction error [%] DEG-NODEG Inverse method Forward method (b) Ground truth cases considering degradation - Conventional traction reconstruction with homogeneous ECM properties. Figure 5.15: In silico 3D TFM results: traction error metric magnitude corresponding to the considered timepoints of the study, for both forward (red boxes) and inverse (blue boxes) formulations. “main” — 2025/2/25 — 15:02 — page 108 — #146 Chapter 5. Multiscale framework with heterogeneous matrices 108 in the laboratory. With a sufficient amount of data and validation, it is assumed that experimental ECM degradation measurements will be greatly reduced, complemented and accurately reproduced by the model. Indeed, experimentally, degradation due to proteinase activity can be visualized and characterized in both space and time by means of live-imaging with DQ-collagen [178]. DQ-collagen is initially non-fluorescent, but as it starts to become degraded, it releases fluorescent particles that can be imaged by means of confocal microscopy. This can be used to characterize the rate (speed) of degradation and, furthermore, to tune the degradation model parameters. In addition, only MT1-MMP and MMP-2 -induced degradation was considered in this work. This constitutes a limitation since, in reality, cancer cells are able to secrete a great variety of additional substances that influence the ECM remodeling process (other MMPs, as well as proteases pertaining to other families, such as ADAMs/ADAMTSs, among others). The relation between the degradation model outcomes (i.e. ECM density) and corresponding ECM mechanical properties was established through a multiscale methodology that provides virtual rheology tests for different ECM densities. In this context, the efficacy of the in-house fibered matrix generator tool is proven in regards to providing synthetic rheology data corresponding to generated collagen microstructures that mimic real collagen hydrogels. The shear rheology curves obtained were accurately fitted by the chosen hyperelastic nonlinear model (used to reconstruct macroscopic tractions), which allowed for the retrieval of the relationship between ECM density and material model parameter κo. Nonetheless, a homogeneous behavior of the hydrogel was considered at the initial time of analysis, as commonly considered in TFM. This limitation can be addressed in future works as an application of the proposed computational framework. In this case, the heterogeneous mechanical behavior of each material macroscopic point of the domain can be obtained through the initial matrix density distribution of fibers, similarly to ECM degradation. For that purpose, an accurate reconstruction of the fibered microstructure of the hydrogel is required for the whole domain of analysis (including the embedded cell), and can be obtained, for instance, from SHG microscopy or fiber staining images. Moreover, at the time of reproducing the collagen hydrogel microstructure, constant geometric parameters (length and curvature of fibers) representative of the average values found in the examined images were considered for all fibers in the synthetic microstructures. In the future application of the methodology to a real in vitro experiment, “main” — 2025/2/25 — 15:02 — page 109 — #147 Chapter 5. Multiscale framework with heterogeneous matrices 109 a reconstruction of the exact material microstructure resulting from the segmentation of the 3D images and subsequent meshing will be employed, for a more realistic representation of the assessed cell microenvironment. Regarding the TFM results, there are several interesting conclusions that can be drawn from this study. First, from a simulation point of view, displacements and strains are underestimated at the same time that tractions become higher when neglecting degradation. This is the result of collagen degradation in the vicinity of the boundary of the cell, which induces a local reduction of the stiffness (softening) of the ECM in the location in which the cell exerts forces. Second, with respect to performance of both forward and inverse methodologies at traction reconstruction, the clear outperformance of the inverse method is concluded, as it was previously found in Chapter 4. Finally, it is demonstrated that neglecting ECM degradation (as in conventional TFM) results in a significant overestimation of tractions and non negligible errors in all analyzed cases, but specially for advanced degradation stages. The study of the effect of degradation in traction reconstruction may therefore be more relevant for intermediate and long time scales, and not much for short time scales. However, this also suggests it may also be relevant to study the mechanical behavior of cells under the effect of matrix degradation, to analyze the associated transient dynamics during the whole process and its effect on the mechanical behavior of cells. This suggests that the presented approach may be an available methodology to accurately reconstruct tractions in TFM in the presence of metalloproteinase-mediated ECM degradation. Appendix: Details about the fibered matrix generator tool The extracellular matrix (ECM) is composed by a complex tapestry-like structure of collagen fibers. These fibers are connected to each other within a non-fibrillar supporting matrix. Fibers are subdivided into fibrils, fibrils are composed of collagen molecules, and collagen molecules are in turn formed by ordered groupings of alpha and amino acid chains. This hierarchical structure makes the ECM’s macroscopic behavior to be largely dependent on its microstructural organization. Insight into the mechanical behavior of fibered structures is thus of great interest in mechanobiology, for it will enable tissue engineering of biomimicking materials. Reproducing structures with similar properties to the ECM facilitates in vitro studies to resemble in vivo material behavior [179, 180] that provides a suitable en- “main” — 2025/2/25 — 15:02 — page 110 — #148 Chapter 5. Multiscale framework with heterogeneous matrices 110 vironment for the organization and regulation of cells, potentially aiding in developing new therapeutic procedures [181]. Collagen fibers are found in a variety of individual curved shape configurations, such as twisted, crimped or non-regular curved fibers [179, 182–184]. Uniaxial deformations applied to a curved fiber straightens it, resulting in a mechanical response characterized by small stresses. As the distance between the fiber’s ends approach its natural length, the fiber gets more rigid, the resulting stresses becoming larger in order for the fiber to maintain strain rate. Eventually, the fiber shows linear behavior once it is completely straightened [184–186]. In the context of matrices of interconnected fibers, fibers’ axial direction becomes in time aligned with the stretch direction, prompting the stiffening of the matrix [187]. Thus, the variably curved shape of the fibers at the fiber scale (micro) introduces nonlinearity into the stress-strain response of the matrix at the tissue scale (macro) due to the stiffening resulting from fibers stretching and aligning with the loading direction. In this appendix, the multiscale formulation underlying the in-house fibered matrix generator tool employed in this study is briefly described. The way in which the tool concerned was applied to characterize the non-linear macroscopic behavior of the employed randomly curved fibered microstructures is described in Section 5.4.1. As a disclaimer, none of the symbols utilized in the definition of the different variables herein presented have any correspondence with those employed all through the rest of the thesis. This appendix is self-contained, and is part of a bigger study carried out by the author and the supervisor of this thesis; to which the reader is referred for the consultation of more details about the formulation and results of the model concerned [188]. A.1. Single fiber generation procedure The mechanical characterization of a curved single fiber is described hereunder, in regards to the methodology adopted for the generation of randomly curved single fibers, as well as to the mathematical background vinculated. Validation of the approach concerned is given at the end of the segment by comparing the solution provided by the model to a known analytical scenario. Methodology This segment describes the methodology to build different 3D randomly curved collagen fibers (see Fig. A.1). The input parameter Lis defined as “main” — 2025/2/25 — 15:02 — page 111 — #149 Chapter 5. Multiscale framework with heterogeneous matrices 111 the distance between both ends of the fiber. This distance is discretized into npoints, in which an auxiliary polar system (r, θ)is defined. The polar angle is randomly selected from a uniform distribution such that θ∈[0,2π). Analogously, the fiber eccentricity from the main axis, ris selected from a uniform distribution r∈[0,2R](see Fig. A.1a). Therefore, the coordinates of the different discretized points along the fiber are given as: xi=rcosθ , r sin θ , Li n, i = 1...n. (A.1) The parameter nrepresents in Eq. (A.1), the number of twists of the curved fiber along its length. On the other hand, Ris defined from Fig. A.1a as follows: R=s(1 + P)L n2 −L n2 ,(A.2) with P(P≥0) being a relative length parameter of the fiber length versus the initial length of the end points. Figure A.1: Schematics of the fiber generator: (a) Definition of the discrete points xi(for n= 3). (b) Different individual randomly-curved fiber in 3D representation, for input parameters n= 3, 5 and 10; and Lf/L = 1.11, 1.21 and 1.4 for blue, green and magenta fibers, respectively. Once the 3-D spatial points xihave been created (coarse discretization of the fiber), the curved fiber is smoothed by defining a spline through “main” — 2025/2/25 — 15:02 — page 112 — #150 Chapter 5. Multiscale framework with heterogeneous matrices 112 points xi. This curve is then discretized (fine discretization of the fiber) into NE elements. The final length of the fiber Lfis then defined as PNE i=1 li, being lithe element length. Different generated individual fibers can be seen in Fig. A.1b. Mathematical formulation The generated curved fiber is discretized into NE Euler-Bernouilli beam finite elements. Each element contains 2 nodes and 6 degrees of freedom (3 translations and 3 rotations) per node. The updated lagrangian formulation is adopted, such that the fiber geometry is updated at each current time tof analysis: xi,t+∆t=xi,t + ∆ui(A.3) with xi,t+∆tand xi,t the position vectors of node iat configurations t+∆t and t, respectively, and ∆uithe vector of displacements of node ifrom configuration tto t+∆t. This vector ∆uiis obtained after the discretization and assembly of global vector and matrix following a structural matrix finite element analysis [189]: ANE e=1∆Fe=ANE e=1 nKe(xe i,t)·∆ueo(A.4) being Athe assembly operator; and ∆Feand ∆uethe incremental nodal vector of structural forces and displacements/rotations at element e, respectively. Ke(xe i,t)is the matrix of element ecomputed in the configuration twith element nodal coordinates xe i,t. After assembly, the (global) system in Eq. (A.4) yields, ∆F=Kt·∆u(A.5) The total nodal force vector is then updated as: Ft+∆t=Ft+ ∆F(A.6) The finite element implementation of Eq. (A.5) was developed in Matlab®. The computer implementation was validated in Fig. A.2 versus the exact solution of a nonlinear cantilever rod subjected to a point bending moment. “main” — 2025/2/25 — 15:02 — page 113 — #151 Chapter 5. Multiscale framework with heterogeneous matrices 113 Figure A.2: Validation of the single fiber computational implementation: Exact solution of a nonlinear cantilever rod subjected to a point dimensionless bending moment ˆ M=ML EI versus its finite element implementation. The analytical solution is an arc with dimensionless radius ˆ R= 1/ˆ M[1]. The rod was discretized into 150 elements, and the moment step was ∆ˆ M=π/500. A.2. Fibered matrix generation procedure This section describes the multiscale mechanical analysis of a fibered matrix, composed of randomly curved fibers, from the mechanical interaction of fibers in a representative volume element (RVE). Matrix generator The generation of the RVE is established under the basis of the definition of RVE, that is a selected volume of the fiber microstructure that statistically represents the heterogeneity of the matrix ||[190]. The RVE is composed of isolated fibers Nfiand crosslink fibers Nfx, such that the total number of fibers of the RVE is Nf =Nfi+Nfx. Therefore, the volume concentration “main” — 2025/2/25 — 15:02 — page 114 — #152 Chapter 5. Multiscale framework with heterogeneous matrices 114 of fibers is defined as Vc=Nf/VRV E, with VRV E being the volume of the RVE. The algorithm to generate fibered matrices proceed as follows: Box 1: Algorithm to generate fibered matrices. 0. Set the volume concentration of fibers Vcand the volume of the RVE, and consequently, Nf,Nfiand Nfx. 1. Generate Nfiisolated fibers (see Section A.1: Methodology). 2. FOR m= 1..Nfi 2.1 Randomly select Nmx crosslinking nodes of fiber m. 2.2 FOR j= 1..Nmx i. Create a set with potential connecting candidate crosslinking nodes k, for fibers n= 1..Nfi(n=m) such that L·(1 −ϵ)≤Lj−k≤L·(1 + ϵ). ii. Randomly select node jfrom the set. iii. Generate a crosslinking fiber from nodes j−k(with length L) as in Section A.1: Methodology. END FOR. END FOR. Nmx = 1 is set in Box 1, meaning that a crosslinking fiber is generated per isolated fiber. Therefore, Nmx represents the degree of crosslinking of the matrix. On the other hand, the endpoints length of the crosslinking fiber is set to the same length Lof the isolated fiber, according to item i in Box 1, up to a tolerance ϵwhich is set to 0.01 in the associated code. Multiscale formulation The overall macroscopic mechanical behavior of a fibered matrix and its evolution, is obtained from the micromechanical interaction of fibers within the RVE in a multiscale fashion (see Fig. A.3). In the macroscale, the (Cauchy) stress tensor and the (logarithmic) strain tensor are defined as σM I,t(XI,t),εM I,t(XI,t), respectively, in a macroscopic (Gauss) material point Iof the macroscale (with coordinates XI,t), for load time t(see Fig. A.3). On the other hand, in the microscale, the variables Fm t(xi,t),um t(xi,t), associated to an Euler-Bernoulli beam in a material point (node) iof the mi-