Abstract
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 compressible 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.
Keywords:
Living cells; Viscoelasticity; Constitutive model; AFM nanoindentation; Multiscale
1 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., 2015Lavagnino, M., Wall, M.E., Little, D., Banes, A.J., Guilak, F., Arnoczky, S.P., (2015). Tendon mechanobiology: Current knowledge and future research opportunities. Journal of Orthopaedic Research 33, 813822.).
On its hierarchy inside the tissue, cells are subjected to a complex mechanical environment (Qi et al., 2006Qi, J., Fox, A.M., Alexopoulos, L.G., Chi, L., Bynum, D., Guilak, F., Banes, A.J., (2006). IL1beta decreases the elastic modulus of human tenocytes. Journal of applied physiology 101, 189195.). These mechanical stimuli are transduced by cells into biochemical signals leading to tissue adaptation. The particular mechanochemical pathways that induce changes in tissue morphology are usually known as mechanotransduction mechanisms (Lavagnino et al., 2015Lavagnino, M., Wall, M.E., Little, D., Banes, A.J., Guilak, F., Arnoczky, S.P., (2015). Tendon mechanobiology: Current knowledge and future research opportunities. Journal of Orthopaedic Research 33, 813822.; Wang, 2006Wang, J.H.C., (2006). Mechanobiology of tendon. Journal of Biomechanics 39, 15631582.).
Regardless their physiological functions, living cells are structural units holding particular nonlinear behaviors, including viscous and plastic effects (Bonakdar et al., 2016Bonakdar, N., Gerum, R., Kuhn, M., Spörrer, M., Lippert, A., Schneider, W., Aifantis, K.E., Fabry, B., (2016). Mechanical plasticity of cells. Nature Materials 16.; Haase and Pelling, 2015Haase, K., Pelling, A.E., (2015). Investigating cell mechanics with atomic force microscopy. Journal of The Royal Society Interface 12.; Nawaz et al., 2012Nawaz, S., Sánchez, P., Bodensiek, K., Li, S., Simons, M., Schaap, I.A.T., (2012). Cell ViscoElasticity Measured with AFM and Optical Trapping at SubMicrometer Deformations. PLoS ONE 7.). 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, 2015Haase, K., Pelling, A.E., (2015). Investigating cell mechanics with atomic force microscopy. Journal of The Royal Society Interface 12.; Janmey and McCulloch, 2007Janmey, P.A., McCulloch, C.A., (2007). Cell mechanics: integrating cell responses to mechanical stimuli. Annual review of biomedical engineering 9, 134.; Kollmannsberger and Fabry, 2011Kollmannsberger, P., Fabry, B., (2011). Linear and Nonlinear Rheology of Living Cells. Annual Review of Materials Research 41, 7597.). 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 (1999Ortiz, M., Stainier, L., (1999). The variational formulation of viscoplastic constitutive updates. Computer Methods in Applied Mechanics and Engineering 7825, 419444.) and Radovitzky and Ortiz (1999Radovitzky, R., Ortiz, M., (1999). Error estimation and adaptive meshing in strongly nonlinear dynamic problems. Computer Methods in Applied Mechanics and Engineering 172, 203240.), particularized in Fancello et al. (2006Fancello, E., Ponthot, J.P., Stainier, L., (2006). A variational formulation of constitutive models and updates in nonlinear finite viscoelasticity. International Journal for Numerical Methods in Engineering 65, 18311864.) and Vassoler et al. (2012Vassoler, J.M., Reips, L., Fancello, E.A., (2012). A variational framework for fiberreinforced viscoelastic soft tissues. International Journal for Numerical Methods in Engineering 89, 16911706.) 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.
2 CONTINUUM FORMULATION
2.1 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 ψ(F, α), and a socalled dissipation pseudopotential function ϒ(de Souza Neto et al., 2009de Souza Neto, E.A., Peric, D., Owen, D.R.J., (2009). Computational Methods for Plasticity: Theory and Applications, Engineering.; Jirásek and Bazant, 2002Jirásek, M., Bazant, Z., (2002). Inelastic analysis of structures. John Wiley & Sons, Ltd.). 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 α. The dissipation potential ϒ, on the other hand, depends also on the corresponding rates , and , 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., F, α) (
It is hypothesized that = (F, α), which allows the following decoupling:
Based on these assumptions, it is shown in Ortiz and Stainier (1999Ortiz, M., Stainier, L., (1999). The variational formulation of viscoplastic constitutive updates. Computer Methods in Applied Mechanics and Engineering 7825, 419444.) and Radovitzky and Ortiz (1999Radovitzky, R., Ortiz, M., (1999). Error estimation and adaptive meshing in strongly nonlinear dynamic problems. Computer Methods in Applied Mechanics and Engineering 172, 203240.) 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 PiolaKirchhoff stress tensor is computed from the first partial derivative of the rate potential with respect to also results in the first PiolaKirchhoff stress tensor:. Moreover, the rate potential evaluated at the minimizer argument of (2.2) defines a reduced potential 𝒫red, whose total derivative with respect to
Potential 𝒫^{red} 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.
2.2 Thermodynamic Potentials for a Class of Viscoelastic Models
The classic multiplicative decomposition of the deformation gradient F is used herein to separate elastic F ^{e} and viscous F ^{v} contributions:
Based on this, total, elastic and viscous right CauchyGreen strain tensors are given, respectively, by
Similarly, total and viscous rate of deformations are formally defined as
where the assumption of null inelastic spin is used (Anand and Gurtin, 2003Anand, L., Gurtin, M.E., (2003). A theory of amorphous solids undergoing large deformations with applications to polymers and metallic glasses, International Journal of Solids and Structures 40, 129.; Gurtin and Anand, 2005Gurtin, M., Fried, E., Anand, L., (2010). The mechanics and thermodynamics of continua. Cambridge University Press.).
The proposed model follows the rheological assembly illustrated in Figure 1. According to this, the set of internal variables is reduced to α: = {F ^{v}}, and the Helmholtz free energy and dissipation potential take the form:
The superscript ∞ represents the timeindependent response. As can be noted, the dissipation potential φ of equation (2.1) is null. Moreover, the arguments of the potentials (2.9) are objective entities (Gurtin et al., 2010Gurtin, M., Fried, E., Anand, L., (2010). The mechanics and thermodynamics of continua. Cambridge University Press.).
In view of (2.9), the rate potential (2.2) reduces to the expression:
with d ^{v} depending on 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.
3 INCREMENTAL FORMULATION
3.1 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 ∆t = _{tn} _{+1}  _{tn} , a possible general expression for this incremental potential is
where the variables d ^{v} and are discrete approximations for the rates . One can note that ϕv is evaluated at an intermediate time tn +ϑ inside the time increment, where the parameter ϑ ∈[0, 1] is closely related to the discretization rule employed for . The consistency of the incremental form (3.1) can be verified if ∆t → 0 ⇒𝒫inc = 𝒫. and
The most common option to obtain an update expression for the internal variable Fvn+1, in function of the discrete tensor dv, is the wellknown exponential mapping (Weber and Anand, 1990Weber, G., Anand, L., (1990). Finite deformation constitutive equations and a time integration procedure for isotropic, hyperelasticviscoplastic solids. Computer Methods in Applied Mechanics and Engineering 79, 173202.):
This choice has been extensively used in combination with the spectral parameterization of Fancello et al., 2006Fancello, E., Ponthot, J.P., Stainier, L., (2006). A variational formulation of constitutive models and updates in nonlinear finite viscoelasticity. International Journal for Numerical Methods in Engineering 65, 18311864.; Fancello et al., 2008Fancello, E., Vassoler, J.M., Stainier, L., (2008). A variational constitutive update algorithm for a set of isotropic hyperelasticviscoplastic material models. Computer Methods in Applied Mechanics and Engineering 197, 41324148.; Vassoler et al., 2012Vassoler, J.M., Reips, L., Fancello, E.A., (2012). A variational framework for fiberreinforced viscoelastic soft tissues. International Journal for Numerical Methods in Engineering 89, 16911706.). Moreover, under an incompressible viscous flow assumption, the incremental relation tr(F ^{v} _{n} _{+1}) =1 is automatically fulfilled.) = 0 ⇔ det ( to handle arbitrary isotropic hyperelastic potentials or in combination with Hencky’s isotropic model to obtain some numerical advantages (
Another alternative expression to obtain the updated value F ^{v} _{n} _{+1} is based on the classical Euler integration rule:
When ϑ = 1, 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 , 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.
3.2 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.
3.2.1 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 ℐso: = {ym det[F ^{v} _{n} _{+1}(ym the space of symmetric second order tensors. The kinematic constraint J ^{v} _{n} _{+1}:= det(F ^{v} _{n} _{+1})=1 is taken into account by means of the Lagrangian functional)] =1} represents the isochoric viscous space and 𝒮 ∈ 𝒮
where ϒ_{n} _{+1} is a Lagrange multiplier. Consequently, the minimum principle (3.6) can be rewritten as an unconstrained problem, such that,
The solution of (3.8) defines the internal variable updates algorithm. Once the optimal solution F ^{v} _{n} _{+1} is updated by equation (3.5). The proposed solution strategy is based on the full NewtonRaphson procedure, where further operational details are addressed in appendix A. It is important to note that the symmetric nature of Itskov (2002Itskov, M., (2002). The derivative with respect to a tensor: some theoretical aspects and applications. ZAMM  Journal of Applied Mathematics and Mechanics 82, 535544.) and Jog (2007Jog, C.S., (2007). Foundations and Applications of Mechnics Volume I: Continuum Mechanics, 2 edition. ed. Alpha Science Intl Ltd.) for further details about directional derivatives in symmetric spaces). In this particular case, the solutions of the principle (3.8) always return symmetric tensors is obtained from (3.8), the internal variable 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 .
3.2.2 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 PiolaKirchhoff stress tensor is given by
The partial derivatives of the Helmholtz strain energies in relation to F _{n} _{+1} result in
where S _{n} _{+1} is the second PiolaKirchhoff stress tensor. In addition, the timeindependent 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 PiolaKirchhoff 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, 1985Simo, J.C., Taylor, R.L., (1985). Consistent tangent operators for rateindependent elastoplasticity. Computer Methods in Applied Mechanics and Engineering 48, 101118.). Therefore, a closed form for the consistent material tangent modulus is provided in appendix B.
3.3 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., 2002Hamann, S., Kiilgaard, J.F., Litman, T., AlvarezLeefmans, F.J., Winther, B.R., Zeuthen, T., (2002). Measurement of Cell Volume Changes by Fluorescence SelfQuenching. Journal of Fluorescence 12, 139145.; ZlotekZlotkiewicz et al., 2015ZlotekZlotkiewicz, E., Monnier, S., Cappello, G., Le Berre, M., Piel, M., (2015). Optical volume and mass measurements show that mammalian cells swell during mitosis. The Journal of Cell Biology 211, 765774.). Based on these observations, the following NeoHookean and quadratic potentials are chosen to represent the viscoelastic behavior of living cells:
where {µ ^{∞} , κ ^{∞} , µ ^{e} , η ^{v}} is the set of constitutive parameters to be identified according to experimental data. In the expression ψ ^{∞}, the parameter κ ^{∞} controls changes in volume, while µ ^{∞} accounts for distortional effects. To ensure an incompressible behavior, κ ^{∞} typically takes values of ((10^{3}  10^{4}) µ ^{∞} (Bonet and Wood, 2008Bonet, J., Wood, R.D., (2008). Nonlinear continuum mechanics for finite element analysis, 2nd Editio. ed, Zhurnal Eksperimental’noi i Teoreticheskoi Fiziki. Cambridge University Press.).
4 NUMERICAL SIMULATIONS AND RESULTS
4.1 Compressibility of Fibroblast Cells
AFM indentation experiments were carried out by Nawaz et al. (2012) in fibroblast cells. Authors reported cyclic experiments under indentation depth of (0.38 ± 0.02 μm at 0.8 μm/s indentation rate, retrieving 72 ± 5 pN indentation forces. The relative energy dissipated by the hysteresis loop reached a value of (30±5% (Figure 3a).
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.
(a) Comparison between experimental and numerical results obtained from the nanoindentation tests. The experimental points were plotted based on data provided by Nawaz et al. (2012). (b) Sensitivity of the numerical results in relation to the incompressible response of the model. The related constitutive parameter sets are shown in Table 1.
Numerical simulations of these experiments were run with a rigid 1.98 μm diameter spherical indenter and a numerical sample with diameter of 10 μm and height of 4 μm. 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 usersubroutine 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 κ ^{∞}, the parameter µ ^{∞} was kept fixed and the viscoelastic parameters µ ^{e} and η ^{v} were adjusted in order to fit the numerical curves to the range of the experimental data (see Table 1).
Figure 3a compares the forcedisplacement curves of the indenter obtained from the experiments and Simulation 1. Figure 3b show the corresponding results obtained with all simulations. Finally, Figure 4 depicts the displacement and the vonMises (Cauchy) stress fields at maximum penetration of the indenter for Simulation 4.
Displacement (a) and vonMises stress (b) fields at the instant of the maximum penetration of the indenter for Simulation 4.
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. Figure 5a display the uniaxial first Piolastretch curves, while the corresponding time history of the Jacobians are shown in Figure 5b. All simulations were run under an engineering strain rate of 0.05s^{1}.
Stressstretch 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 Table 1.
4.2 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 depth at indentation rates of 0.1, 0.5, 1.0 and 10 μm/s, was run. The corresponding forcedisplacement curves are displayed in Figure 6a. Different hystereses are seen, depending on the indentation rate. The relative percentage of the energy dissipated at each hysteresis loop is shown in Figure 6b. Moreover, a quadratic trend curve is fitted in order to better illustrate how the dissipated energy changes within the range of indentation rates.
(a) Cyclic nanoindentation tests at different rates. (b) Relative loss of energy by the hysteresis loops computed from cyclic curves.
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; MorenoFlores et al., 2010). Based on this motivation, a set of relaxation and creep simulations under different load conditions were tested, as shown in Figure 7. Concerning the relaxation tests, three indentation depths of 0.1, 0.3 and 0.5 μm were applied at a 10 μm/s rate, and then the indenter was kept fixed for 10 seconds. The corresponding relaxation curves are displayed in Figure 7a. 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 7b.
Mechanical responses predicted by the viscoelastic model when submitted to nanoindentation tests of stress relaxation (a) and creep (b) under different loading conditions.
5 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 3b, it can be seen that when incompressibility is enforced, the forceindentation response changes its curvature within the firsts 0.1 μ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 5b). However, further numerical and experimental investigations are required to support this hypothesis.
Taking into account the cyclic numerical experiments displayed in Figure 6a 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 6b. 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 MorenoFlores et al. (2010MorenoFlores, S., Benitez, R., Vivanco, M. dM, TocaHerrera, J.L., (2010). Stress relaxation and creep on living cells with the atomic force microscope: a means to calculate elastic moduli and viscosities of cell components. Nanotechnology 21.).
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 (1998Vinckier, A., Semenza, G., (1998). Measuring elasticity of biological materials by atomic force microscopy. FEBS Letters 430, 1216.)). 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.
Acknowledgments
The authors would like to thank the financial support provided by the Brazilian funding agencies CAPES  (Coordination for the Improvement of Higher Education Personnel) and CNPq  (National Council for Scientific and Technological Development).
References
 Anand, L., Gurtin, M.E., (2003). A theory of amorphous solids undergoing large deformations with applications to polymers and metallic glasses, International Journal of Solids and Structures 40, 129.
 Belytschko, T., Liu, W., Moran, B., (2000). Nonlinear Finite Elements for Continua and Structures, Chichester, New York, John Wiley. John Wiley & Sons, Ltd.
 Bonakdar, N., Gerum, R., Kuhn, M., Spörrer, M., Lippert, A., Schneider, W., Aifantis, K.E., Fabry, B., (2016). Mechanical plasticity of cells. Nature Materials 16.
 Bonet, J., Wood, R.D., (2008). Nonlinear continuum mechanics for finite element analysis, 2nd Editio. ed, Zhurnal Eksperimental’noi i Teoreticheskoi Fiziki. Cambridge University Press.
 de Souza Neto, E.A., Peric, D., Owen, D.R.J., (2009). Computational Methods for Plasticity: Theory and Applications, Engineering.
 Fancello, E., Ponthot, J.P., Stainier, L., (2006). A variational formulation of constitutive models and updates in nonlinear finite viscoelasticity. International Journal for Numerical Methods in Engineering 65, 18311864.
 Fancello, E., Vassoler, J.M., Stainier, L., (2008). A variational constitutive update algorithm for a set of isotropic hyperelasticviscoplastic material models. Computer Methods in Applied Mechanics and Engineering 197, 41324148.
 Gurtin, M., Fried, E., Anand, L., (2010). The mechanics and thermodynamics of continua. Cambridge University Press.
 Gurtin, M.E., Anand, L., (2005). The decomposition F = FeFp, material symmetry, and plastic irrotationality for solids that are isotropicviscoplastic or amorphous. International Journal of Plasticity 21, 16861719.
 Haase, K., Pelling, A.E., (2015). Investigating cell mechanics with atomic force microscopy. Journal of The Royal Society Interface 12.
 Hamann, S., Kiilgaard, J.F., Litman, T., AlvarezLeefmans, F.J., Winther, B.R., Zeuthen, T., (2002). Measurement of Cell Volume Changes by Fluorescence SelfQuenching. Journal of Fluorescence 12, 139145.
 Itskov, M., (2002). The derivative with respect to a tensor: some theoretical aspects and applications. ZAMM  Journal of Applied Mathematics and Mechanics 82, 535544.
 Janmey, P.A., McCulloch, C.A., (2007). Cell mechanics: integrating cell responses to mechanical stimuli. Annual review of biomedical engineering 9, 134.
 Jirásek, M., Bazant, Z., (2002). Inelastic analysis of structures. John Wiley & Sons, Ltd.
 Jog, C.S., (2007). Foundations and Applications of Mechnics Volume I: Continuum Mechanics, 2 edition. ed. Alpha Science Intl Ltd.
 Kollmannsberger, P., Fabry, B., (2011). Linear and Nonlinear Rheology of Living Cells. Annual Review of Materials Research 41, 7597.
 Lavagnino, M., Wall, M.E., Little, D., Banes, A.J., Guilak, F., Arnoczky, S.P., (2015). Tendon mechanobiology: Current knowledge and future research opportunities. Journal of Orthopaedic Research 33, 813822.
 MorenoFlores, S., Benitez, R., Vivanco, M. dM, TocaHerrera, J.L., (2010). Stress relaxation and creep on living cells with the atomic force microscope: a means to calculate elastic moduli and viscosities of cell components. Nanotechnology 21.
 Nawaz, S., Sánchez, P., Bodensiek, K., Li, S., Simons, M., Schaap, I.A.T., (2012). Cell ViscoElasticity Measured with AFM and Optical Trapping at SubMicrometer Deformations. PLoS ONE 7.
 Ortiz, M., Stainier, L., (1999). The variational formulation of viscoplastic constitutive updates. Computer Methods in Applied Mechanics and Engineering 7825, 419444.
 Qi, J., Fox, A.M., Alexopoulos, L.G., Chi, L., Bynum, D., Guilak, F., Banes, A.J., (2006). IL1beta decreases the elastic modulus of human tenocytes. Journal of applied physiology 101, 189195.
 Radovitzky, R., Ortiz, M., (1999). Error estimation and adaptive meshing in strongly nonlinear dynamic problems. Computer Methods in Applied Mechanics and Engineering 172, 203240.
 Simo, J.C., Taylor, R.L., (1985). Consistent tangent operators for rateindependent elastoplasticity. Computer Methods in Applied Mechanics and Engineering 48, 101118.
 Vassoler, J.M., Reips, L., Fancello, E.A., (2012). A variational framework for fiberreinforced viscoelastic soft tissues. International Journal for Numerical Methods in Engineering 89, 16911706.
 Vinckier, A., Semenza, G., (1998). Measuring elasticity of biological materials by atomic force microscopy. FEBS Letters 430, 1216.
 Wang, J.H.C., (2006). Mechanobiology of tendon. Journal of Biomechanics 39, 15631582.
 Weber, G., Anand, L., (1990). Finite deformation constitutive equations and a time integration procedure for isotropic, hyperelasticviscoplastic solids. Computer Methods in Applied Mechanics and Engineering 79, 173202.
 ZlotekZlotkiewicz, E., Monnier, S., Cappello, G., Le Berre, M., Piel, M., (2015). Optical volume and mass measurements show that mammalian cells swell during mitosis. The Journal of Cell Biology 211, 765774.
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 NewtonRaphson procedure.
A.1 Stationarity Condition
In view of the variational principle (3.8), the stationarity condition is defined through the variation dℒ = 0, which results in the following equations:
The partial derivatives of (A.1) result in,
and,
where M _{en+} _{1} := C_{e n+} _{1}S_{e n+} _{1} are the elastic Mandel stress tensor.
A.2 NewtonRaphson Procedure
Concerning to the full NewtonRaphson procedure, one can define from the nonlinear equations (A.1) the quantities,
where r _{n+} _{1} is a residual and x_{n+} _{1} represents the set of unknown variables. Therefore, an increment δ related to a current iteration of the Newton’s algorithm, is computed from the linear system of equation,
which can be represented by the compact arrangement J _{n+} _{1} ⋆ δx_{n+} _{1}= r _{n+} _{1}. The symbol ⋆ represents the appropriate product between the Newton’s increments δx_{n+} _{1} 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,
In equations (A.8), 𝕀_{sym} is the fourth order symmetric identity tensor. In addition, are introduced the tensor products (A ⊗ B)_{ijkl}: =(A)_{ij} (B)_{kl} and (AB)_{ijkl}: =(A)_{ik} (B)_{jl} , for any second orders tensors A and B.
APPENDIX B  CONSISTENT MATERIAL TANGENT MODULUS
Taking into account a total Lagrangian formulation (Belytschko et al., 2000Belytschko, T., Liu, W., Moran, B., (2000). Nonlinear Finite Elements for Continua and Structures, Chichester, New York, John Wiley. John Wiley & Sons, Ltd.), the linearization of the equilibrium equations results in the material tangent modulus
where one can note that if 𝒫^{red} _{inc} is a convex potential of C _{n} _{+1}, the tensor ℂ_{n} _{+1} will present major symmetry. The total derivative of the second PiolaKirchhoff stress tensor is given by
The partial derivatives of (B.2) result in
where,
in which the fourth order operators ∂S^{e} _{n+1}/∂^{e} _{n} _{+1} were introduced in equations (A.8). The total derivative dC _{n+1} of equation (B.2) is one of the solutions of the following linear system of equation:/d and ℂ
resulting from the total derivative of the residual (A.4). The linear system (B.5) can be presented in a convenient compact form:
One can note that from equation (B.6), the derivative ∂r _{n+1}/∂ x_{n+1} is the Jacobian matrix previously defined in equation (A.6). Finally, the partial derivatives of the set ∂r _{n+1}/∂C _{n+1}are computed as
where,
Publication Dates

Publication in this collection
June 2017
History

Received
07 Feb 2017 
Reviewed
07 Mar 2017 
Accepted
13 Mar 2017