Modeling the Local Viscoelastic Behavior of Living Cells Under Nanoindentation Tests

Aiming at multiscale numerical investigations of fibrous connective tissues, the present work proposes a finite strain viscoelastic model suitable to represent the local mechanical response of living cells in conjunction with finite element simulations of nanoindentation tests. The material model is formulated in a thermodynamically consistent framework based on a variational constitutive approach. As a case of study, a numerical investigation of the local compress-ible response of fibroblast cells is addressed. Moreover, a set of cyclic, stress relaxation and creep experiments are simulated in order to investigate the capability of the model to predict rate dependence of cells, showing sound agreement with experimental data. The proposed model may be used as a suitable tool for a better understanding of stiffness and energetic dissipation of cells.


INTRODUCTION
Withstanding mechanical loadings is the primary function of many biological connective tissues such as tendons, ligaments and cartilage.The mechanical response to external loadings are highly dependent on their hierarchical microstructure, cellular organization and interactions between one another (Lavagnino et al., 2015).
On its hierarchy inside the tissue, cells are subjected to a complex mechanical environment (Qi et al., 2006).These mechanical stimuli are transduced by cells into biochemical signals leading to tissue adaptation.The particular mechano-chemical pathways that induce changes in tissue mor-Latin American Journal of Solids and Structures 14 (2017) 844-860 phology are usually known as mechanotransduction mechanisms (Lavagnino et al., 2015;Wang, 2006).
Regardless their physiological functions, living cells are structural units holding particular nonlinear behaviors, including viscous and plastic effects (Bonakdar et al., 2016;Haase and Pelling, 2015;Nawaz et al., 2012).Accordingly, cells may present important mechanical contributions in the macroscopic response of tissues, where multiscale analysis can be employed to investigate these micromechanical behaviors.
In multiscale numerical analyses based on representative volume elements (RVE), the homogenized behavior relies strongly on how accurate the behavior of each material phase is represented by corresponding specific material models.In addition, these models must show sound agreement with experimental data at that scale, specifically in cases where these phases present nonlinear behavior and are subjected to finite strains.
Several measurement techniques have been employed to assess the mechanical response of living cells, such as micropipette aspiration, magnetic and optical tweezers and atomic force microscopy (AFM) indentation, also known as nanoindentation test (Haase and Pelling, 2015;Janmey and McCulloch, 2007;Kollmannsberger and Fabry, 2011).Among these experimental methods, the later provides a convenient structure to finite element simulations.For instance, whereas the cell sample can be considered locally as a homogenous medium, the tip of the indenter is modeled as a rigid body.
With these motivations in mind and aiming future numerical multiscale investigations of fibrous soft tissues, the present work proposes a finite strain viscoelastic model suitable to represent the local mechanical response of living cells in conjunction with finite element simulations of nanoindentation tests.The proposed model is formulated in a thermodynamically consistent framework based on the variational formalism addressed in Ortiz and Stainier (1999) and Radovitzky and Ortiz (1999), particularized in Fancello et al. (2006) and Vassoler et al. (2012) to viscoelastic materials.In addition, an alternative incremental solution strategy is proposed.
This manuscript is organized as follows.Section 2 presents the theoretical background employed to model the viscoelastic behavior of cells.Aiming finite element simulations, the continuum constitutive equations are discretized in time, leading to a variational constitutive update algorithm, shown in section 3.In order to verify if the material model is able to represent the available experimental data, section 4 addresses a set of numerical simulations of nanoindentation experiments of fibroblast cells.In addition, the compressible and viscoelastic responses of the model are studied.Particularities of the finite element simulations and further discussions about the mechanical behavior of cells are highlighted in section 5. Finally, appendices provide operational details and the analytical expression for the consistent material tangent modulus.

Constitutive Modeling of Dissipative Materials in a Variational Framework
The constitutive modeling of dissipative materials can be cast in a thermodynamically consistent framework by defining a potential function, in present case the Helmholtz free energy ( ) , y F a , and a so-called dissipation pseudo-potential function ( ) , ; , ¡ F F a a   (de Souza Neto et al., 2009;Jirásek Latin American Journal of Solids and Structures 14 (2017) 844-860 and Bazant, 2002).The argument of the free energy, in the case of an isothermal process, is the thermodynamic state defined by the deformation gradient F and a set of internal variables a .The dissipation potential ¡ , on the other hand, depends also on the corresponding rates F  and a  , where the semicolon in the argument list indicates the state variables play a role of parameters.For the sake of clarity, only the dependence on the rates is kept henceforward in the notation of this potential.
It is hypothesized that , which allows the following decoupling: Based on these assumptions, it is shown in Ortiz and Stainier (1999) and Radovitzky and Ortiz (1999) that a potential  of the form may be defined such that the optimality condition of the minimization problem ( ) provides a thermodynamically consistent equation for the evolution of the internal variables: The first Piola-Kirchhoff stress tensor is computed from the first partial derivative of the rate potential with respect to F  .Moreover, the rate potential evaluated at the minimizer argument of (2.2) defines a reduced potential red  , whose total derivative with respect to F  also results in the first Piola-Kirchhoff stress tensor: plays an important role in present framework, since its incremental counterpart presents a compelling mathematical structure for the calculation of the incremental updates of the internal variables.

Thermodynamic Potentials for a Class of Viscoelastic Models
The classic multiplicative decomposition of the deformation gradient F is used herein to separate elastic e F and viscous v F contributions: , : det 0, : det 0 Similarly, total and viscous rate of deformations are formally defined as : sym , , : sym , skew , where the assumption of null inelastic spin is used (Anand and Gurtin, 2003;Gurtin and Anand, 2005).
The proposed model follows the rheological assembly illustrated in Figure 1.According to this, the set of internal variables is reduced to , and the Helmholtz free energy and dissipation potential take the form: (2.9) The superscript ¥ represents the time-independent response.As can be noted, the dissipation potential j of equation (2.1) is null.Moreover, the arguments of the potentials (2.9) are objective entities (Gurtin et al., 2010).In view of (2.9), the rate potential (2.2) reduces to the expression: (2.10) with v d depending on v F  by means of (2.8).Equation (2.10) provides a continuum background for the variational approach introduced in section 2.1.The incremental counterpart of the constitutive equations and numerical issues are addressed in the next section.

Incremental Potential and Updating Rule for the Internal Variable
An incremental counterpart of the rate potential (2.10) is used to directly compute the incremental internal variable updates needed in finite element schemes.Considering a time increment where the variables and v F  are discrete approximations for the rates v d and v F  .One can note that v f is evaluated at an intermediate time n t J + inside the time increment, where the parameter 0,1 J é ù Î ë û is closely related to the discretization rule employed for v F  .The consistency of the incremental form (3.1) can be verified if The most common option to obtain an update expression for the internal variable , is the well-known exponential mapping (Weber and Anand, 1990): ( ) This choice has been extensively used in combination with the spectral parameterization of v d  to handle arbitrary isotropic hyperelastic potentials or in combination with Hencky's isotropic model to obtain some numerical advantages (Fancello et al., 2006;Fancello et al., 2008;Vassoler et al., 2012).Moreover, under an incompressible viscous flow assumption, the incremental relation is automatically fulfilled.( )

Another alternative expression to obtain the updated value
(3.4) When 1 J = , a fully implicit integration scheme is recovered.In this case, substituting equation (3.4) into (3.3)leads to the following updating rule: ( ) where I is the second order identity tensor.
In view of the parameterization (3.5), the resulting variational problem is formulated in a fully tensorial format.While this implies the simultaneous evaluation of a set of independent components, the formulation remains able to handle anisotropic behaviors.Moreover, neither the spectral decomposition of v d  , nor the first and second derivatives of these spectral quantities are needed.In addition, if the viscous flow is considered incompressible, it is enforced into the variational problem by means of a Lagrangian functional (further details in next section).
To end this section it is worth mentioning that, even both approaches are quite different, is not the goal of the present work to compare and quantify its numerical advantages or drawbacks.Solids and Structures 14 (2017) 844-860

Constitutive Algorithm
The local constitutive algorithm is composed by two main stages: the solution of the variational principle and the stress evaluation.Accordingly, the general concepts introduced in section 2 are particularized to the present class of viscoelastic models.

Solution of the Variational Principle
In view of the parameterization (3.5), the viscous rate of deformation becomes the primal variable of the variational problem.Moreover, in the present modeling, the viscous flow is considered incompressible.Therefore, the minimum principle (2.3) can be presented as where represents the isochoric viscous space and ym  the space of symmetric second order tensors.The kinematic constraint ( ) where 1 n g + is a Lagrange multiplier.Consequently, the minimum principle (3.6) can be rewritten as an unconstrained problem, such that, ( ) ( ) It is important to note that the symmetric nature of v d  is not included as a constraint into the problem (3.8).However, its symmetry is considered throughout the directional derivatives required by the Newton's procedure, resulting in suitable Fréchet operators (see Itskov (2002) and Jog (2007) for further details about directional derivatives in symmetric spaces).In this particular case, the solutions of the principle (3.8) always return symmetric tensors

Stress Evaluation and Consistent Material Tangent Modulus
Once the solution of (3.8) is obtained, the stress can be updated as follows.Taking into account the incremental potential (3.1), the incremental first Piola-Kirchhoff stress tensor is given by Latin American Journal of Solids and Structures 14 (2017) 844-860 The partial derivatives of the Helmholtz strain energies in relation to , : where 1 n+ S is the second Piola-Kirchhoff stress tensor.In addition, the time-independent and elastic second Piola stresses are defined as It should be noted that, due to the directional derivative of equation (3.9), the elastic second Piola-Kirchhoff stress tensor in (3.10) is consistently pulled back to the referential configuration.
Within the framework of a conventional nonlinear finite element code, the consistent tangent modulus must be provided (Simo and Taylor, 1985).Therefore, a closed form for the consistent material tangent modulus is provided in appendix B.

Choice of Specific Potentials
A common assumption in the biomechanical modeling of cells is considering them as an incompressible material.However, under physiological conditions, cells can experience large volume changes (Hamann et al., 2002;Zlotek-Zlotkiewicz et al., 2015).Based on these observations, the following Neo-Hookean and quadratic potentials are chosen to represent the viscoelastic behavior of living cells: : 2 is the set of constitutive parameters to be identified according to experimental data.In the expression y ¥ , the parameter k ¥ controls changes in volume, while m ¥ accounts for distortional effects.To ensure an incompressible behavior, k ¥ typically takes values of ( ) 10 10 m ¥ - (Bonet and Wood, 2008).

Compressibility of Fibroblast Cells
AFM indentation experiments were carried out by Nawaz et al. (2012)  As can be seen from Figure 2, the cell sample is considered locally homogeneous and it was discretized by 8,130 quadratic hexahedral elements with reduced integration.Due to symmetry, only a quarter of the sample was modeled.Moreover, frictionless contact was considered in the interface between sample and indenter.
The finite element simulations presented herein were run in the software Abaqus, where the proposed viscoelastic model was implemented into the user-subroutine UMAT.Using this finite element setup, a set of numerical simulations evaluated the capability of present material model to represent the experimental data, considering different compressibility levels.
To this aim, for increasing values of k ¥ , the parameter m ¥ was kept fixed and the viscoelastic parameters e m and v h were adjusted in order to fit the numerical curves to the range of the experimental data (see   The related constitutive parameter sets are shown in Table 1.
In order to gain a better understanding of how the parameters of Table 1 affect the incompressible response of the model, simulations of homogeneous uniaxial tensile tests were performed.1.

Cyclic, Stress Relaxation and Creep Tests
Identification procedures often recover material parameters that successfully reproduce the used set of experimental data.However, it is not unusual these parameters fail to predict the behavior of the same material under a mechanical condition different from (but close to) that identified.In order to gain confidence on the behavior of the proposed model using the identified parameters, a brief sensitivity analysis considering different mechanical tests was carried out.
The same finite element model of the previous section with constitutive parameters of Simulation 1 (Table 1) was used.
A first set of numerical tests considering cyclic indentations with 0.5 m m depth at indentation rates of 0.1 , 0.5 , 1.0 and 10 m/s m , was run.The corresponding force-displacement curves are displayed in Figure 6-a.Different hystereses are seen, depending on the indentation rate.The relative percentage of the energy dissipated at each hysteresis loop is shown in Figure 6-b.Moreover, a quadratic trend curve is fitted in order to better illustrate how the dissipated energy changes within the range of indentation rates.Besides cyclic tests, stress relaxation and creep AFM indentation tests are also found in the literature as efficient experiments to get information on the viscoelastic properties of cells (Haase and Pelling, 2015;Moreno-Flores et al., 2010).Based on this motivation, a set of relaxation and creep simulations under different load conditions were tested, as shown in

DISCUSSIONS AND FINAL REMARKS
As can be seen from Figure 3, all numerical simulations were able to predict the local viscoelastic response of fibroblast cells under nanoindentation (in view of the standard deviations) using spherical indenter.In addition, the displacement and stress fields depicted in Figure 4 emphasize the numerical results seem to present minimum influence of the boundaries of the sample.Concerning to the numerical curves plotted in Figure 3-b, it can be seen that when incompressibility is enforced, the response changes its curvature within the firsts 0.1 m m displacement, diverging from the usual expected behavior.Based on these numerical results, one can hypothesize that cells seem to respond locally, under indentation tests, as a compressible material (see time histories of the volumetric jacobians plotted in Figure 5-b).However, further numerical and experimental investigations are required to support this hypothesis.
Taking into account the cyclic numerical experiments displayed in Figure 6-a one can observe that, for the identified material parameters, the dissipation decreases when the indentation rate increases.The energetic loss is quantified by the hysteresis loop as shown in Figure 6-b.Moreover, one can see the percentage of hysteresis presents a nearly quadratic decay over rates.
According to the relaxation and creep simulations shown in Figure 7, the material reaches equilibrium states beyond 5 seconds in stress relaxation and 10 seconds in creep.These time scales predicted by the model show sound agreement with experimental data reported in Moreno-Flores et al. (2010).
It is worth mentioning the use of analytical expressions based on Hertz contact model and indentation tests to estimate the cell stiffness is common in the literature (technical details may be found elsewhere, Vinckier and Semenza (1998)).This approach, however, is quite sensitive to the indentation depth, since the Hertz expressions were tailored for linear elastic materials and infinitesimal strains.In addition, it is unable to consider rate dependent effects.In contrast to these analytical approaches, the finite element simulation procedure presented in this work seems to provide a more accurate tool for a better understanding of stiffness and energetic dissipation of cells.
Within a multiscale numerical framework of biological tissues, the present constitutive model may be employed to represent the microscopic material phases composed of cells and/or other isotropic viscoelastic components as well, a necessary step for future developments in the mechanobiology field, where the knowledge of the complex mechanical responses of cells are issues of major concern.

APPENDIX A -SOLUTION OF THE LOCAL CONSTITUTIVE PROBLEM
Herein are present the operational details related to the solution strategy employed to solve the variational constitutive principle (3.8).Furthermore, the nonlinear equations resulting from the stationarity condition are solved by the full Newton-Raphson procedure.

A.1 Stationarity Condition
In view of the variational principle (3.8), the stationarity condition is defined through the variation , which results in the following equations: The partial derivatives of (A.1) result in, ( ) and, where , : Latin American Journal of Solids and Structures 14 (2017) 844-860 where 1 n + r is a residual and represents the set of unknown variables.Therefore, an increment d related to a current iteration of the Newton's algorithm, is computed the linear system of equation, which can be represented by the compact arrangement . The symbol  repre- sents the appropriate product between the Newton's increments and the Jacobian matrix In order to achieve the quadratic convergence rate provided by the Newton's method, the derivatives of equation (A.6) must be computed precisely.Accordingly, the second partial derivatives of (A.6) result in where it is introduced the fourth order operators, Latin American Journal of Solids and Structures 14 (2017) 844-860 In equations (A.8), sym  is the fourth order symmetric identity tensor.In addition, are introduced the tensor products ( ) ( ) ( ) and ( ) ( ) ( , for any second orders tensors A and B .

APPENDIX B -CONSISTENT MATERIAL TANGENT MODULUS
Taking into account a total Lagrangian formulation (Belytschko et al., 2000), the linearization of the equilibrium equations results in the material tangent modulus


will present major symmetry.The total derivative of the second Piola-Kirchhoff stress tensor is given by The partial derivatives of (B.2) result in ( ) where, ( ) ( ) : 4 , : : : 2) is one of the solutions of the following linear system of equation: resulting from the total derivative of the residual (A.4).The linear system (B.5) can be presented in a convenient compact form:

Figure 1 :
Figure 1: Schematic representation of the rheological model.

F
solution of (3.8) defines the internal variable updates algorithm.Once the optimal solution opt v d  is obtained from (3.8), the internal variable 1 v n+ is updated by equation (3.5).The proposed solution strategy is based on the full Newton-Raphson procedure, where further operational details are addressed in appendix A.
in fibroblast cells.Authors reported cyclic experiments under indentation depth of 0.38 0.02 m m   at 0.8 m/s m indentation rate, retrieving 72 5 pN  indentation forces.The relative energy dissipated by the hysteresis loop reached a value of 30 5 %   (Figure 3-a).Latin American Journal of Solids and Structures 14 (2017) 844-860 Numerical simulations of these experiments were run with a rigid 1.98 m m diameter spherical indenter and a numerical sample with diameter of 10 m m and height of 4 m m .

Figure 2 :
Figure 2: Finite element model employed for the simulations of nanoindentation tests.(a) Perspective view of the model showing the spherical indenter and the cell sample in a quarter symmetry.(b) Top view of the numerical sample emphasizing the mesh refinement.

Figure 3 -
Figure3-a compares the force-displacement curves of the indenter obtained from the experiments and Simulation 1.Figure 3-b show the corresponding results obtained with all simulations.Finally,

Figure 4 Figure 3 :
Figure3-a compares the force-displacement curves of the indenter obtained from the experiments and Simulation 1.Figure 3-b show the corresponding results obtained with all simulations.Finally,

Figure 4 :
Figure 4: Displacement (a) and von-Mises stress (b) fields at the instant of the maximum penetration of the indenter for Simulation 4.

Figure 5 :
Figure 5: Stress-stretch curves (a) and time history of the volumetric jacobians (b) obtained from simulations of homogeneous uniaxial tensile tests.The related constitutive parameter sets are shown in Table1.

Figure 7 .
Concerning the relaxation tests, three indentation depths of 0.1, 0.3 and 0.5 m m were applied at a 10 m/s m rate, and then the indenter was kept fixed for 10 seconds.The corresponding relaxation curves are displayed in Figure 7-a.Concerning the creep tests, three indentation forces of 10, 25 and 50 pN were progressive and linearly applied in 0.05 seconds, and then kept constant for 16 seconds.The corresponding creep curves are shown in Figure 7-b.

Figure 6 :Figure 7 :
Figure 6: (a) Cyclic nanoindentation tests at different rates.(b) Relative loss of energy by the hysteresis loops computed from cyclic curves.

Table 1 :
Table1).Constitutive parameters related to the numerical curves shown in Figure3.