A Continuum-Tensegrity Computational Model for Chondrocyte Biomechanics in AFM Indentation and Micropipette Aspiration

Mechanical stimuli are fundamental in the development of organs and tissues, their growth, regeneration or disease. They influence the biochemical signals produced by the cells, and, consequently, the development and spreading of a disease. Moreover, tumour cells are usually characterized by a decrease in the cell mechanical properties that may be directly linked to their metastatic potential. Thus, recently, the experimental and computational study of cell biomechanics is facing a growing interest. Various experimental approaches have been implemented to describe the passive response of cells; however, cell variability and complex experimental procedures may affect the obtained mechanical properties. For this reason, in-silico computational models have been developed through the years, to overcome such limitations, while proposing valuable tools to understand cell mechanical behaviour. This being the case, we propose a combined continuous-tensegrity finite element (FE) model to analyse the mechanical response of a cell and its subcomponents, observing how every part contributes to the overall mechanical behaviour. We modelled both Atomic Force Microscopy (AFM) indentation and micropipette aspiration techniques, as common mechanical tests for cells and elucidated also the role of cell cytoplasm and cytoskeleton in the global cell mechanical response.


INTRODUCTION
Apart from genes and chemical factors, mechanical stimuli are important regulators for the development of organs and tissues, their growth and connected processes such as remodelling, regeneration or disease (References 20, 21, 29 to cite a few). This strong influence is due to the fact that cells can sense and respond to mechanical signals, by converting them into biochemical responses (this process is called mechanotransduction, firstly coined in 1998 6 ). Moreover, external mechanical forces control the shape, type and function of a living cell by altering the internal balance of the cell itself, 7 thus, modifying the internal prestress of cell main subcomponents (the cytoskeleton, the membrane, the cytoplasm and the nucleus) and consequently influencing the biochemical signals produced by the cell. 21,31,40 The evidence that biochemistry and signals transmission can be modified by cell mechanics variations is also strictly connected to apparently unrelated diseases manifestation, where an abnormal mechanotransduction, sometimes combined with alterations of the Extra Cellular Matrix (ECM) mechanical properties, influences a disease development and spreading. 21,29 In particular, for oncological pathologies, tumour cells are usually characterized by a decrease in the mechanical properties that provides the cell with a higher deformation and mobility, thus suggesting that cell mechanics can be directly linked also to its metastatic potential. 9,16,20,21,42 For this reason, the study of cell mechanics and related properties is facing a growing interest in the last years, from both an experimental and a computational point of view.
Various experimental approaches have been used to describe the passive response of cells, 18 as the AFM indentation, [8][9][10][11] which is able to locally stress the cell with compression forces, thanks to a microscale cantilever with a spherical tip, or the Micropipette Aspiration (MPA), 12,15,17,27,34,38,39 by means of microscale pipette tips that can apply a negative pressure to the cell surface, thus exciting the cell with tensile stresses. Other techniques, as the Magnetic Tweezers (MT), are more interested in determining the mechanical properties of cell subcomponents such as the cytoskeleton. 4 However, experimental cell response measurements are usually characterized by a huge variability due to cell phenotypes and types, 8,11,12,26,28,33,34,38 shape, 8 source 9,26,33 and aging. 30 This being the case, it is quite difficult to fully describe the mechanical behaviour of a cell by means of in in-vivo tests, and even more difficult to extract the role of each subcomponents.
For this reason, in-silico computational models have been developed through the years, to overcome such limitations and uncertainties, while proposing a valuable tool to understand cell mechanical behaviour. [1][2][3]5,24,28,32,41 Some of these models describe the cell as an homogeneous viscoelastic material, [1][2][3]41 others specify the main subcomponents adopting a tensegrity structure to model the cytoskeleton, 5,24,25,28 as stated and assessed in past works. 22,23,40 Computational models usually aim at describing cells undergoing AFM indentation, accounting for cell mechanics variability 3 or the effect of an altered cytoskeleton to predict the mechanics of cancer cells, 24 as well as MPA, to understand the influence of the pipette shape and the material properties. 1,15,41 For this purpose, thanks to a computational approach, we used finite elements to analyse the mechanical response of a cell experiencing both AFM indentation and MPA techniques, in order to: (i) develop a homogeneous model that is able to catch the global viscoelastic response of a cell subjected to both compression and tensile mechanical stimuli; (ii) propose a combined continuous-tensegrity model to observe how each subcomponent contributes to the overall mechanical behaviour; (iii) compare the two models to elucidate the influence of the cytoplasm and the cytoskeleton, by varying their mechanical properties. We firstly set both the homogeneous and the continuum-tensegrity models for the description of the AFM indentation tests, then we validated their applicability with the MPA simulations, highlighting from these two types of tests some useful insights on the subcellular components roles.

Finite Element Models of the Cell and Its Subcomponents
Three-dimensional (3D) models of a cell were realized with the finite element software Abaqus Explicit 2019 (Dassault Systemes Simulia Corp., Providence, RI).
The model simulated a 16 lm diameter cell (with reference to Ref. 28), composed of all those features that mainly contribute to the cell mechanics, such as the cytoskeleton (microfilaments and microtubules), the cytoplasm, the cell membrane and the nucleus (Fig. 1). A homogeneous continuous model (CM) and a combined continuous-tensegrity model (CMT) were developed to analyse the mechanical response of a cell undergoing AFM indentation and aspiration through micropipette.
In the CM the cell was represented with a 3D solidsphere composed only by the cytoplasm, discretized by means of about 100.000 linear hexahedral elements, 110.000 nodes and 360.000 degrees of freedom. Within the CTM, the nucleus was represented, with reference to the literature 28, as an ellipsoid (major axis 8 lm and minor axis 5 lm) while the membrane was modelled with a shell part with a constant thickness of 6 nm. In this case the cytoskeleton was described with a tensegrity structure of six compression bearing struts and twenty-four tensional cables that are able to mimic the real behaviour of the microtubules and the microfilaments, respectively. 5,22,23,28 The connection between each strut and cable created twelve common nodes, representing the 'receptor' sites where actin filaments clustered at adhesion complexes. Regarding their geometry, the microtubules had a cross-sectional area of 190 nm 2 , while the microfilaments were thinner, with an area of 18 nm 2 . 28 The CTM was discretized with a number of linear hexahedral elements between 70.000 and 400.000 depending on the model complexity, about 35.000 linear quadrilateral elements for cell membrane and 3000-linear truss elements for the cytoskeleton. In the model, the nodes of the membrane were coincident with the underlying cytoplasm and with cytoskeletal receptor sites, and tie constraints were assigned to cell subcomponents. The total number of nodes was within the range 105.000-440.000 resulting in 370.000-1.500.000 degrees of freedom.
Even if the model has the potential to mimic all kinds of cells, this work primarily focused on chondrocytes (specialized cells present in the cartilage) and could be addressed to chondrosarcoma cells (malignant tumour cells that origin from chondrocytes), thanks to a larger availability of data in the literature with respect to other cells.

Finite Element Model of Atomic Force Microscope (AFM) Indentation
The Atomic Force Microscope indentation was performed on a rounded cell adherent to a rigid substrate, where cell height and contact radius with the substrate were assumed with reference to Ref. 28, thus about 14 lm and 6 lm, respectively. A spherical rigid body simulated the cantilever tip of the AFM. Different analyses were performed with 5 lm and 10 lm tip diameter sizes, to observe possible changes induced by the indenter size in the load-displacement response of both the CM and CTM, such as non-linear contribution of the cell subcomponents. The contact between the cell and the indenter was assumed to be frictionless, for the tangential behaviour, and hard contact type, for the normal behaviour. The bottom nodes, at the cell-substrate interface, were constrained in all three translational degrees of freedom. These constrained points mimicked the focal adhesion sites in cells adherent to a substrate, thus implies that the substrate is adopted as rigid. In some recent works the influence of a solid substrate has been studied, especially when a large probe indents a spread cell, thus some bottomeffects corrections have been proposed. 13,14 However, thanks to the comparison with the Hertz model of the contact between two spheres and the models outputs, we noticed no undesired effects due to the applied boundary conditions, probably also because we are adopting an almost spherical configuration of the cell, instead of a spread one.
The analyses were performed with a two-step dynamic explicit simulation. Similarly to experimental protocols reported in literature, 10 during the first step, a 1500 nm displacement was applied by the spherical indenter to the top of the cell; then this loading phase was followed by a relaxation one in which the cantilever was held in place for up to 60 s (Fig. 2a). In this way, stress relaxation tests were performed on the central region of the cell using a 9.5 lm/s approach velocity. 8

Finite Element Model of Micropipette Aspiration (MPA)
The Micropipette Aspiration tests simulated the interaction with a rounded visco-hyperelastic cell and a round rigid cylindrical pipette with a smooth-edged mouth and 900 nm-round fillet. In this case, several simulations were realized, considering different ratios between the cell and the micropipette radii, more precisely with ratio values of 1.5, 2, 3, 4 and 5.5. The contact between the cell and the micropipette was assumed to be frictionless (tangential behaviour), and hard contact type (for the normal behaviour).
The Micropipette Aspiration tests were performed with a two-time step dynamic explicit simulation, similarly to the AFM indentation. During the first step, lasting 1 s, the cell was exposed to a fluid negative pressure variation ÀDP, within the micropipette. After this phase, the creep response of the cell was analysed by keeping a constant ÀDP and measuring the projection length L P of the cell inside the pipette for up to 60 s (Fig. 2b).
Since there was a concentration of the deformation gradient near the pipette fillet, a dense mesh in the upper part of the cell was adopted, for both the cytoplasm and the plasma membrane.

Mechanical Properties of Cell Subcomponents
Both the CM and the CTM were tested with viscoelastic and visco-hyperelastic material properties. The Neo-Hookean formulation was used for the hyperelastic material model. The mechanical properties adopted in this study are reported in Tables 1 and  2, with reference to the several analysed configurations.
Both AFM indentation and MPA analyses lasted between 1 and 72 h, depending on the model complexity and the number of required steps, running contemporary on 20 threads of a High-Performance Computing Server Fujitsu Primergy RX4770 equipped with four Intel Xeon E7 8890 v4 processors, 512 GB RAM and SSD HD.

AFM Indentation: Loading Phase and Stress Relaxation Behaviour
A typical representation of the AFM indentation experiments is usually performed by reporting the force-displacements curves obtained at the cantilever tip. A good analytical model usually employed in this case is the Hertz contact model, 10 reported in Eq. (1), which describes the interaction between two spheres, one of which is assumed to be infinitely rigid. Given that the cantilever tip can be indeed described as an infinitely rigid body, it is possible to use this equation to describe the force-displacement relationship during the loading phase of an indentation experiment: d is the indentation, E el is the Young's modulus of the cytoplasm, m is the Poisson's ratio and R is the equivalent radius calculated using Eq. (2): where R cell is the radius of the cell and R tip is the radius of the cantilever spherical tip. When introducing the viscoelastic behaviour of the cell (i.e. the stress relaxation phenomenon) the relationship between the force and the displacement could be described with a Solid Linear Standard model (SLS) with a deformation given by the one of the Hertz model, as described in Eq. (3) 9,10 : where s r and s e are the relaxation times under constant load and constant deformation respectively, while E R corresponds to the stiffness k 1 of the Maxwell's standard linear solid model (Fig. 3b). By fitting Eq. (3) to a force-displacement curve, it is possible to obtain a standard linear solid representation of the cell's viscoelastic response, where: where E 0 and E 1 are the instantaneous and long-term elastic response respectively and m is the Poisson's coefficient of the cell. 9 Firstly, the AFM experiment was simulated with a homogeneous computational model, similarly to past works. 1,35,41 Both a linear elastic and a Neo-Hookean hyperelastic formulations were adopted and compared with respect to the Hertz analytical model, to confirm the almost equal response. In Fig. 3 we reported the comparison between the linear elastic Hertz model and different model configurations all with the Neo-Hookean material formulation.
The force-displacements curves obtained using the two different indenter sizes (namely R = 2.5 lm and 5 lm) revealed a similar behaviour of both material formulations, close to the Hertz predictions, thus they can be used as a reference for the comparison with the complex continuum-tensegrity model (Fig. 3a).
In addition, stress relaxation data were obtained from the second step of the simulations, and reported with reference to the Hertz viscoelastic model (Eq. 1) in Fig. 3b.
In order to qualitatively and quantitatively describe the contribution of the cytoplasm versus the cytoskeleton, the mechanical properties of each subcomponents were altered using various parameters combinations, as reported in Table 3. The parameter Q = E el1 /E el2 = 12.78 is the ratio between the elastic moduli of the cytoplasm of cell type 1 (E el1 ) and cell E el is the Young's modulus assuming a linear elastic formulation for the loading phase; E R is the relaxed modulus (with reference to 8 ), which corresponds to k 1 of the Standard Linear Solid model (Fig. 3b)  Parameters notation is the same as reported in Table 1.

BIOMEDICAL ENGINEERING SOCIETY
type 2 (E el2 ), and it was used to increase or decrease the mechanical properties between these cell types, obtaining other four sets of mechanical properties ( Table 4). The results of these simulations are reported in Fig. 4a while normalized values with respect to the maximum force reached in each group is shown in Fig. 4b.
It is possible to observe how the overall mechanical response is more affected by a decrease in the mechanical properties of the cytoskeleton rather than by its increase.

Micropipette Aspiration (MPA)
In order to represent the mechanical behaviour of a cell during micropipette aspiration, it is customary to display the displacement of the cell during the application of a constant negative pressure inside the chamber of the micropipette. The half-space viscoelastic model frequently used to represent the evolution of the displacement of the cell has been developed by Sato et al. 34 which represents an exten-sion of the previous half space elastic model by Theret et al. 37 This viscoelastic model assumes an incompressible behaviour; thus, the aspiration length is obtained from the following equation: where L p is the projection length of the cell inside the pipette, R p is the micropipette radius, DP is the pressure applied to the cell's surface, E 2 and E 1 are the elastic constants of the springs of the corresponding Maxwell standard linear solid representation of the cell, s is the characteristic time of the creep response and finally / is the ''punch coefficient'' usually reported to be / % 2:1. 19,41 However, the half-space analytical model presents some limitations when applied to the cell biomechanics, as also already reported in some works, 1,41 since cell behaviour is compressible, with an average Poisson's coefficient of 0.35-0.4 and strains are not infinitesimal.
In Fig. 5a it is possible to observe the comparison between compressible and incompressible models of  These values were used to evaluate the role of the tensegrity structure with respect to the overall mechanical response of the computational model.
the micropipette aspiration during the loading phase using different values of the ratio between the micropipette diameter and the cell diameter. Some curves are obtained from a work of Baaijens et al. 1 and are compared with our simulations and the half-space model developed by Theret et al. 37 In Fig. 5b it t is possible to observe the influence of the Poisson's ratio during the loading phase, comparing our model with different results by Baaijens et al. 1 and the half-space model, highlighting the effects of adopting a compressible model instead of an incompressible one.
Since the ratio D c /D p has shown to play a significant role on the resulting stimulus-response curves, different ratios have been further investigated (Fig. 5c) it and compared with previous results obtained from literature. 1 To account for possible coupling and scale effects, both the micropipette and the cell radii were alternatively changed, resulting in different combinations of D c /D p values, as reported in Table 4.
The computational model has been tested for its viscoelastic response as well. Two of the simulations obtained using two values of the D c /D p ratio, as well as the experimental curve from the work of Baaijens et al. 1 are reported in Fig. 5d. Then, once the CM has been validated by comparison with both computational and experimental data obtained from literature, the CTM has been included within the finite element simulation of micropipette aspiration, and compared with respect to the former model. Different orientations of the cytoskeleton could alter the results; thus, two possible tensegrity icosahedron dispositions were used. Figure 6 reported the two analysed configurations and summarizes the results obtained during the loading phase in these simulations.

DISCUSSION
It has been assessed that a tumour cell exhibits a variation in its mechanical behaviour with respect to a healthy cell, which implies also a change in the mechanotransduction of signals that regulate its normal routine processes. For this reason, in-silico com- Four parameters combinations to stress both the role of the cytoplasm and the cytoskeleton to the overall mechanical cell response.  Table 3. (b) normalized results with respect to the maximum force obtained for both cell type 1 and 2.
putational models have been developed throughout the years to deepen the knowledge about cell mechanical behaviour without experimental uncertainties. With this aim, in this work we developed a finite element model of a cell and its subcomponents in order to evaluate and quantify the influence of each part, when the same cell is subjected to compression and tensile mechanical stimuli. The first approach consisted in adopting a continuum model (CM) with a homogeneous material to simulate AFM indentation on a single cell and the related stress relaxation behaviour. Both the loading and stress-relaxation phases were confirmed by Hertz model (Fig. 3), where an optimal correspondence was reached. When the tensegrity structure and the other components were added to the computational model (CTM, Fig. 1), a similar trend was noticed, but characterized by larger values, especially when increasing the indenter size. In particular, the load-displacement curve of the CTM slightly deviated from both the linear Hertz model and the CM when the indenter tip was greater. This aspect highlighted the contribution of the cell subcomponents with respect to the CM, when the behaviour became non-linear (greater displacements). Moreover, these effects became more evident if other parameters combinations were considered, as shown in Fig. 4, where the Young's modulus of the cytoplasm was kept constant and the other subcomponents' elastic moduli were variated by a factor Q. Computational results showed that a stiffer cytoskeleton, membrane and nucleus contribute in enhancing the response of the cell subjected to a compression force, but this influence is even more noticeable when these subcomponents are characterized by a softer behaviour, which strongly affects the overall cell mechanics (Fig. 4a). Similar trends were found by changing cell type from 1 to 2, which consists of a softer cytoplasm (one order of magnitude lower). When considering the normalized results in order to analyse only the influence of cell subcomponents ( (Fig. 4b), it is possible to infer that they significantly modify a cell mechanical behaviour, regardless of the  25 , they observed that the cytoskeleton is the most involved part to carry the reaction force during the AFM tip indentation. From our insights, we can also state that cytoplasm is able to influence the global mechanical response of a cell undergoing AFM indentation, since its variation led to significant different behaviours of the cell (Fig. 4), but cell subcomponents (especially the cytoskeleton) are the ones able to tune the overall cell behaviour in a not uniform way. This is evident in Figs. 4a and 4b, where different mechanical parameters combinations were used: by reducing the mechanical properties of the subcomponents of one order of magnitude, the mechanical response of the cell significantly decreases, while increasing them of the same amount does not influence the global response in the same way.
When dealing with MPA, other useful aspects emerged from the computational modelling. Theret's elastic model 37 has been adopted through the years to obtain the elastic properties of the cell (i.e., the elastic modulus) by describing a linear dependence between the normalized aspiration length and the applied pressure (e.g., Refs. 10,31,34,36), as reported by Eq. (6). Theret's model assumes the cell as a homogeneous incompressible half space and considers only infinitesimal strains. However, it has been observed 1 that this common approach is not suitable for cells characterized by a Poisson's coefficient m < 0.5 and subjected to large elongations. Moreover, in MPA, cells and micropipettes sizes are comparable, thus contrasting the hypothesis of the half space model.
Theret's solution is reported in Figs. 5a, 5b with a dashed blue line, in contrast with our continuum model and previous studies by Baaijens et al. 1 . When a compressible material is adopted, it emerges the need of a computational approach to catch the cell biomechanics during MPA. In agreement with Baaijens' predictions, our results for a single cell, modelled with only cytoplasm as a homogeneous compressible material, revealed a non-linear variation of the aspiration length (L p ) when increasing the applied negative pressure.
The influence of cell and micropipette diameters (D c and D p respectively) was also investigated, by varying one or both in the range D c /D p from 1.5 to 5.5 (Fig. 5c). As reported in previous works, 1,41 differences in problem size (i.e., the ratio D c /D p ) appear to strongly modified L p ; results also highlighted that the parameter which governs this behaviour is not only D c or D p , but the ratio D c /D p . Figure 5c reports the aspiration lengths of different combinations of D c /D p , normalized by the cell radius instead of the micropipette one, since it usually may vary between cells. This visualization allows to clearly identify the differences by varying D c /D p , in agreement with other proposed models in literature.
When the CTM was used also to model MPA, in order to evaluate possible variations in the final L p, two cytoskeleton orientations were analysed, namely configuration 1 and configuration 2 (Fig. 6a). Slight oscillations with respect to the CM were observed between the two configurations, but some more evident changes arose when the ratio D c /D p decreases, meaning that a larger portion of the cell is subjected to the negative pressure applied by the pipette. These oscillations highlighted that also in the MPA the cytoskeleton plays a role and some configurations (e.g., configuration 2) entail a more significant tensegrity structure involvement (dashed blue lines in Fig. 6b), thus supporting the previous insights obtained in the AFM simulations. Moreover, in MPA both CT and CTM were modelled with the same material parameters adopted in the AFM, thus validated, thanks also to the experimental comparison reported in Fig. 5d. These two applications of the CTM and the obtained results confirmed the great importance and advantages in considering the main cell subcomponents when modelling a single cell subjected to various mechanical stimuli. Even if the here presented model presents some limitations (such as the lack in connections between the nucleus and the other parts or the missing interactions of the cell with the substrate), it represents a step forward in understanding the differences between cell types from a mechanical point of view. In particular, with our model it is possible to explore how a variation in the mechanical properties of a single subcomponent may affect the overall cell behaviour. Moreover, the main benefit in adopting a tensegrity structure is its potentiality in shape adaptation, a key feature to include when dealing with living cells.
In the present work we proposed a computational model to mimic the biomechanical behaviour of living cells, starting from a continuum model, and then adopting a combined continuum-tensegrity approach in order to elucidate how a variation in the parameters of cell subcomponents may alter the global cell response. Being able to model this intrinsic variability is a key factor in cell mechanics, since it is also associated with a different behaviour between healthy and tumour cells. The analysis of the here proposed CM with reference to other past results from the literature confirmed its applicability for these future purposes, while the outcomes of the CTM with respect to the CM underlined the non-negligible mechanical contributes of the cell subcomponents. In particular, the results highlighted that, when a single cell is subjected to AFM indentation and micro-pipette aspiration (MPA) the cytoskeleton may strongly alter the overall cell biomechanics.
Indeed, these FE models represent a useful tool for the mechanical investigation of both living and cancer cells, revealing to be also a valuable steppingstone in the studies of the mechanical processes that undergo during the different stages of tumour cells and to overcame the complexity in studying the neoplasms, in particular when referring to the inter and intra tumour variability.

FUNDING
Open access funding provided by Universita`degli Studi di Padova within the CRUI-CARE Agreement.

Open Access
This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.or g/licenses/by/4.0/.

ACKNOWLEDGMENTS
This work has been supported by the Italian Ministry of Education, University and Research (MIUR), project PRIN 20177TTP3S (InteMA-PreMed: Integrated mechanobiology approaches for a precise medicine in cancer treatment), project PRIN 2017HFPKZY (Modelling of constitutive laws for traditional and innovative building materials) and project PRIN 20209F3A37 (Sustainable modelling of materials, structures and urban spaces including economic-legal implications).