Introduction

In recent years interest has increased in the measurement of the viscoelastic properties of soft biological samples motivated by their correlation with disease, differentiation, or cellular transformation1,2,3,4. A large variety of methods have been introduced for measuring cellular mechanical properties including micropipette aspiration5, 6, stretching or compression between two microplates7,8,9, optical tweezers10, 11, and magnetic twisting cytometry12, 13. However, indentation with the atomic force microscope (AFM) remains one of the most popular methods for probing the nanoscale properties of soft samples like cells, tissues and hydrogels14,15,16,17.

In the AFM, the elastic properties of live cells are usually evaluated from force versus displacement (F-Z) curves. Then the Hertz model or its modifications are applied to the approach part of the F-Z curve to extract Young’s modulus (E Hertz ), the elastic parameter used for characterization of the sample’s mechanical properties18, 19. However, such models assume a purely elastic nature of the sample, while in reality most biological samples are viscoelastic. Viscoelasticity is revealed in a clear hysteresis between the approach and retraction parts of curves20; the indentation speed dependence of E Hertz values extracted from force curves with the Hertz model11; the observations of force relaxation at constant indentation depth and the creep at constant loading force21.

Approaches other than the standard F-Z curves are usually used to obtain the viscoelastic properties of samples with AFM in both the time22,23,24 and frequency domains25,26,27,28,29,30. These generally require modifications in the AFM apparatus and/or in the data acquisition protocol. Equally importantly, each approach has its own sets of measurement uncertainties.

If a standard F-Z curve could also be used to quantify viscoelastic properties, it would allow one standard method with well quantified uncertainties31 to be used for both viscoelasticity and elasticity measurements. This has not been possible to date, we believe, due to the lack of a mathematical/computational framework that allows the post-processing of force-displacement data to extract the relevant viscoelastic constitutive parameters.

Here we propose a new method to extract viscoelastic properties of soft samples like cells and hydrogels directly from standard AFM F-Z curves. This method is based on the theoretical model developed by Ting32 for the problem of the indentation of a linear viscoelastic half-space by a rigid axisymmetric indenter for arbitrary load history. We combine Ting’s model together with related numerical procedures to process the data from both approach and retraction phases of the experimental F-Z curves. The method is validated with finite element (FE) simulations, with experiments performed on polyacrylamide hydrogel (PAAm) samples, and by comparison with the existed AFM viscoelasticity measurement techniques. The method is then applied to measure the viscoelastic properties of several benign and cancerous cell lines (NIH 3T3, NMuMG, MDA-MB-231, MCF-7), and to study the sensitive changes in viscoelastic properties related to tumorigenesis including TGF-β-induced epithelial-to-mesenchymal transition on NMuMG cells and Syk expression-induced phenotype changes in MDA-MB-231 cells. The method also provides a way to test which constitutive viscoelastic relations best fit the acquired data. We tested applicability of two commonly used relaxation models, standard linear solid (SLS) and power-law rheology (PLR), on data collected on hydrogels and the above cell lines.

We show that the developed method is highly robust, does not require changes in experimental design and is compatible with findings from conventional AFM microrheology. By bringing rigorous viscoelastic analysis within the reach of standard, automated AFM force curves, the method is expected to greatly enhance the wider use and adoption of AFM for living cell and soft biomaterials visco-mechanics assays.

Results

Experimental AFM force curves on live cells display viscoelastic behaviour that simply cannot be captured or explained by elastic contact mechanics models. The schematic setup for the method (Fig. 1a) is a commonly used indentation protocol in which the piezo scanner moves the cantilever base relative to the sample with predefined vertical speed profile, which is identical for the approach and retraction phases. The data from such an experiment is the so called force curve (F-Z curve) representing force vs scanner displacement or, after processing, force vs indentation (F-δ curve) dependency. It can be seen from the typical force-indentation curve obtained on a living cell (and also polyacrylamide (PAAm) gels, Fig. S9) that Hertz’s model fits the approach phase fairly well. Not surprisingly, the framework of Hertz contact mechanics19, 33 is commonly used for the description of AFM force curves, at least for shallow indentations and low indentation rates1, 18, 34,35,36,37. However, even at low piezo speeds such as 2 μm/s, viscoelastic behaviour can be observed in the hysteresis between approach and retraction phases, which is not captured by the Hertz’s model (Fig. 1a).

Figure 1
figure 1

Flowchart of the method and used viscoelastic models. (a) The flowchart. At the first step, a force-displacement curve obtained with a standard indentation protocol is pre-processed to determine zero force level and corrected for hydrodynamic drag (if needed). At the second step, the contact point is located and F-Z curve is transformed to the force-indentation coordinates. Apparent Young’s modulus E Hertz is calculated with the Hertz’s model. Inset: scheme of a spherical probe indenting a half-space. At the third step, the indentation time history is used to model the force curve with the Ting’s model and prescribed viscoelastic model. The difference between the modelled and experimental force curves (e) is minimized with the fitting algorithm. Viscoelastic parameters providing the lowest e value are acquired as the output. At the final step, contact point position could be adjusted and the step 3 repeated to obtain the best fit. Details of the processing are given in the text. (b) Normalized relaxation modulus (E/E 0) for standard linear solid (SLS) and power-law rheology (PLR) relaxation models. Insets: schematic spring-dashpot representation and equation for the relaxation modulus. The SLS model is a spring in parallel with a spring-dashpot combination, leading to an exponential relaxation with a single relaxation time τ. The PLR model can be imagined as an infinite number of spring-dashpot combinations in parallel, leading to a continuous relaxation spectrum and power-law decay (α is the power law exponent).

Full size image

Therefore, to account for the viscoelastic nature of the sample, more advanced tip-sample interaction models are needed. Ting’s model32 provides a solution for the problem of indentation of a linear viscoelastic half-space by a rigid axisymmetric indenter for any load history, including approach and retraction phases of F-Z curves. The following equations describe Ting’s solution for indentation of a viscoelastic sample with a rigid spherical indenter:

$$F(t,\delta (t))=\{\begin{array}{cc}\frac{4\sqrt{R}}{3(1-{\nu }^{2})}{\int }_{0}^{t}E(t-\xi )\frac{{\rm{\partial }}{\delta }^{\frac{3}{2}}}{{\rm{\partial }}\xi }d\xi , & 0\le t\le {t}_{m}\\ \frac{4\sqrt{R}}{3(1-{\nu }^{2})}{\int }_{0}^{{t}_{1}(t)}E(t-\xi )\frac{{\rm{\partial }}{\delta }^{\frac{3}{2}}}{{\rm{\partial }}\xi }d\xi , & {t}_{m} < t\le {t}_{ind}\end{array},$$
(1)
$${\int }_{{t}_{1}(t)}^{t}E(t-\xi )\frac{\partial \delta }{\partial \xi }d\xi =0,$$
(2)

where F is the force acting on the cantilever tip; δ is the indentation depth; t is the time initiated at the time of initial contact (t m is the duration of approach phase; t ind is the duration of complete indentation cycle); t 1 is the auxiliary function determined by the equation (2); ξ is the dummy time variable required for the integration; E(t) is the Young’s relaxation modulus; v is the Poisson’s ratio of the sample (assumed to be time-independent); and R is the radius of the indenter. A time-dependent relaxation modulus E(t) allows for arbitrary linear viscoelastic constitutive equations to be used with Ting’s model. The relaxation modulus is generally a decaying function that can be described by rheological models, like the standard linear solid (SLS) and power-law rheology (PLR) models used here. The corresponding E(t) functions are:

$$E(t)={E}_{\infty }+({E}_{0}-{E}_{\infty }){e}^{-\frac{t}{\tau }},$$
(3)
$$E(t)={E}_{0}{(1+\frac{t}{t^{\prime} })}^{-\alpha }.$$
(4)

The normalized SLS and PLR relaxation functions are presented in Fig. 1b. The SLS model is characterized by three parameters, E 0, E ∞, and τ representing the instantaneous and long term elastic modulus and relaxation time, respectively. In the SLS model the extent of relaxation can be determined by the E ∞/E 0 and τ/t ind (Deborah’s number) ratios. The PLR model is characterized by two parameters, the instantaneous elastic modulus E 0 and the power law exponent α; t′ is a small time offset which does not affect relaxation behaviour at experimental times). A larger α value means larger relaxation; materials exhibit a solid-like behaviour at α = 0, and a fluid-like behaviour at α = 1.

We have developed numerical procedures to: first, simulate the F-Z curves with Ting’s equations (1) and (2) and prescribed viscoelastic relaxation model; and, second, to fit the experimental F-Z curves to the Ting’s equations to extract the viscoelastic constitutive parameters of the sample for a chosen viscoelastic relaxation model (Fig. 1, Fig. S2). All algorithms were implemented in MATLAB and are described in the Materials and Methods section.

To test the validity of the derived algorithms, we used finite element (FE) analysis. The FE model represented the indentation experiment where the spherical probe indents the sample with predefined viscoelastic constitutive relation (SLS or PLR). The model-predicted (numerically simulated) and FE-simulated curves were markedly similar for the same input parameters (Fig. 2). Some discrepancy appeared at higher depths and may be associated with finite-size effects. Then the fitting procedure was applied to the FE-simulated curves. An excellent agreement between FE-simulated curve and the results of the fitting procedure (R2 > 0.99) was observed for both SLS and PLR models; the difference between fitted parameters and parameters used in the FE-model was below 4%. Additional FE simulations to check the effects of finite sample thickness and nonconstant tip velocity that arises from the cantilever deflection are described in the Supplementary Information, Section B.

Figure 2
figure 2

Validation of the developed algorithms with FE simulations. (a) Axisymmetric FE model, consisting of the rigid spherical indenter with 2 μm radius and viscoelastic sample with 15 μm height and radius. (b) Magnified view of the contact area, there the sample was meshed more finely. (c) Result of FE simulation and the F-δ curve modelled using Ting’s model with the same input parameters (SLS model with E 0 = 2 kPa, E ∞ = 1 kPa, and τ = 0.1 s).

Full size image

Then, we applied the developed method directly to the force curves obtained on living NIH 3T3 fibroblasts and found that both viscoelastic models, SLS and PLR, fit the experimental F-δ curve fairly well (Fig. 3a). Indeed, both relaxation functions are quite similar when data spans only a single time decade, but differences between them are expected to become apparent on a logarithmic scale over several time decades (Fig. 3b). The SLS model has plateau regions at both short and long times, corresponding to E 0 and E ∞ modulus, while PLR model does not.

Figure 3
figure 3

Comparison of SLS and PLR models in experiments with PAAm hydrogels and NIH 3T3 fibroblasts. (a) Experimental F-δ curve obtained on a fibroblast with both PLR and SLS model fits. (b) Normalized relaxation functions for PRL and SLS models with adjusted parameters on a logarithmic scale, inset – same functions on a linear scale. (c,d) Modelled F-δ curves with different indentation times. Left – SLS (τ = 0.1 s), right – PLR (α = 0.15). (e,f) Experimental F-δ curves obtained on PAAm hydrogel (left) and fibroblast (right) with different indentation times. The offset is added to the force for clarity. Black lines are SLS (e) or PLR (f) model fits. (g,h) Viscoelastic parameters α and τ for PAAm hydrogels (left, combined data for experiments on 3 gels, mean ± s.d.) and NIH 3T3 fibroblasts (right, combined data for experiments on 12 cells, mean ± s.d.) as a function of piezo displacement speed.

Full size image

In an effort to distinguish which viscoelastic model best describes the experimental data, numerical simulations using Ting’s model were performed by varying the indentation time t ind (Fig. 3c and d). As can be seen from the simulated curves using SLS constitutive relations, at both low (≪1) and high (≫1) τ/t ind ratios the approach-retraction hysteresis is small and a complete curve might be described by single modulus E ∞ or E 0 (Fig. 3c). Numerical simulations using the PLR constitutive relation, on the other hand, show that the amount of relaxation varied only slightly (Fig. 3d), the hysteresis area was always present and demonstrated slight growth with t ind . These numerical studies suggest that experimental data should be acquired with different values of t ind to better distinguish the applicability of a specific viscoelastic constitutive relation.

By varying the piezo displacement speed on NIH 3T3 cells and PAAm hydrogels we found that the PLR constitutive relation better describes the viscoelastic data over a wider range of timescales for live cells while the SLS better described the data for the PAAm hydrogels. Specifically, the piezo scanner displacement speed was varied in a range of 64–16,000 nm/s, which covers the majority of published AFM indentation experiments, thus leading to indentation time in a range of 25 ms–6 s (for a ~500 nm indentation depth). For NIH 3T3 cells, the amount of hysteresis (hysteresis area) in F-δ curves was always significant and increased slightly with scanner velocity (more greatly for the highest speed point), indicating prevalence of the PLR model (Fig. 3f and Supplementary Fig. S11). Indeed, after processing with the PLR model, we obtained approximately the same values of power law exponent α for all used scanner displacement speeds except the highest (16 μm/s), for which the α value increased ~30%. For the same data, the SLS model provided dramatically decreasing (by two orders) relaxation time τ values proportionally to the piezo speed, meaning that the results were highly dependent on experimental conditions (Fig. 3h). These results are consistent with the previous observation of power-law decay in step-hold force relaxation experiments38, 39 and in microrheology2, 26, 40 experiments on cells. In case of the PAAm hydrogels, the hysteresis area almost vanished at indentation times larger than 1 s (piezo speed below 0.5 μm/s); and such force curves might be used for acquisition of E ∞ modulus (Fig. 3e,g and Fig. S11C). At shorter times, the hysteresis area gradually increased. The SLS model provided consistent results for indentation times greater than 0.1 s with τ about 0.04 s. At shorter times (higher piezo speeds), τ values started to decrease, indicating the presence of several relaxation mechanisms. Application of the PLR model, however, resulted in gradually increased α values, indicating its poor applicability for describing PAAm hydrogel relaxation. Since it better describes our live cell data, we used the PLR constitutive relation for processing the rest of the data obtained on living cells.

To compare the presented method with the previous techniques (step-hold stress relaxation and microrheology) the experiments were conducted on PAAm hydrogels (3 samples) and NIH 3T3 cells (40 different cells) at the same locations in the random order. The comparison showed a reasonably good agreement between all three methods, the details are given in the Supplementary Information, Section C.

With the described Ting’s model-based method and the PLR relaxation model, we measured the viscoelastic properties of several cell lines: benign murine NIH 3T3 fibroblasts and epithelial NMuMG cells; human breast cancer cell lines MDA-MB-231 and MCF-7. The representative examples of the experimental F-δ curves for each cell line are shown in Fig. 4. In agreement with previous work41, E 0 values exhibited a log-normal distribution and α values were distributed normally (Supplementary Fig. S12). The benign cell lines NIH 3T3 and NMuMG had higher values of E 0 and lower values of α than did the cancer cell lines MDA-MB-231 and MCF-7 (p < 0.01 between all parameters except E 0 for pair of cancer cells and α for the pair of benign cells), which were both softer and more viscous (Table 1, Fig. 5). The results are in line with the data obtained in previous work2, where AFM microrheology experiments were performed in the frequency domain on the same cell lines.

Figure 4
figure 4

Force curves obtained on different cell lines: NIH 3T3, NMuMG, MDA-MB-231, MCF-7. The fit with PLR model; E Hertz , short-term modulus E 0 and power-law exponent α values are shown.

Full size image
Table 1 Viscoelastic parameters of studied cell lines. Mean ± s.d. (median ± m.a.d.).
Full size table
Figure 5
figure 5

Viscoelastic parameters of studied cell lines. Box plots of apparent Young’s modulus E Hertz , short-term modulus E 0 and power-law exponent α. Difference between all distribution except those marked is significant at the p < 0.01 level. All cell lines were at 60–80% confluence during the experiment.

Full size image

It is interesting to check if the correlation exists between extracted viscoelastic parameters. Although cell lines with higher E 0 had lower values of α, inside each cell line the parameters E 0 and α did not correlate with each other. The apparent Young’s modulus measured with the Hertz’s model E Hertz demonstrated a weak negative correlation with α, which was larger for cancer cells (Pearson’s r from −0.1 for NIH 3T3 to −0.6 for MDA-MB-231, p < 0.05). A strong, significant negative correlation was observed between α and the E Hertz /E 0 ratio (≈−0.9, p < 0.001), meaning that E Hertz values depend highly on the amount of relaxation. The example of cancer cells MDA-MB-231 and MCF-7 shows the same effect: both have close E 0 values, but E Hertz is two times lower for MCF-7 cells due to larger α value.

Then we applied the method to study changes in cell mechanical properties in two processes of relevance to cancer tumorigenesis. In the first case study, we compared the difference between control and induced Syk-expressing MDA-MB-231 cells. The Syk protein-tyrosine kinase plays roles in tumour progression, acting as an inhibitor of cellular motility and metastasis in highly invasive cancer cells. Here, both an increase in E 0 and decrease in α in Syk-expressing cells was observed (Table 2), which together led to an even higher increase in E Hertz , consistent with prior work42 (see the Supplementary Information Section E for details).

Table 2 Viscoelastic parameters of MDA-MB-231 cells without (control) and with doxycycline induction (Syk-expressing).
Full size table

In the second case study, we investigated alterations in the mechanical properties of normal mouse mammary epithelial cells (NMuMG) during epithelial-to-mesenchymal transition (EMT). NMuMG cells are extremely sensitive to TGF-β (known driver of EMT), undergoing morphological changes upon TGF-β treatment. Previous research with an AFM microrheology technique showed that α values decrease after TGF-β-induced EMT43. The same trend was observed here (Table 3). Moreover, the E 0 value was approximately preserved during the transition, while α values decreased almost two-fold (see the Supplementary Information Section E for details). Such solidification of cells is probably associated with the assembly of stress fibers and increased pre-stress stored in the membrane and actin cortex43.

Table 3 Viscoelastic parameters of NMuMG cells before and during EMT.
Full size table

Discussion

Viscoelastic properties of cells have been investigated by AFM using several different experimental setups. In the time domain, creep and stress relaxation experiments have been conducted in a large number of studies23, 44, 45. Generally, a Heaviside step-loading is assumed before the hold phase (step-hold experiments)46,47,48 which greatly simplifies the analysis. However, in an actual experiment, the maximum loading speed is limited by hydrodynamic and inertial effects and step-loading conditions are impossible to implement. The fit of the total ramp-hold curve49 or introduction of the ramp correction factor should be implemented for proper analysis50, 51. Difficulties in interpretation may also arise at long relaxation times due to thermal drift in the instrument52 and from active cell responses to the applied forces53.

AFM experiments in the frequency domain are also popular for measuring cell viscoelastic properties. Here the cantilever is sinusoidally oscillated with fixed small amplitude and several frequencies during the period of indentation2, 25,26,27, 43, 54, 55, or during the scanning process28, 30. Frequency-dependent complex Young’s or shear modulus is measured from the amplitude and phase shift of the cantilever displacement. Such experiments generally require more experienced users and have their own set of uncertainties. Thus the extraction of cell viscoelastic properties directly from force curves remains by far the most easily implemented and attractive option for users. Several authors attempted to assess viscoelastic properties directly from force curves, by analysis of approach-retraction hysteresis11, 20, 35, 56, 57 or indentation-rate dependency of apparent elastic modulus11, 36, 58. However, the obtained parameters are not based on a complete viscoelastic solution like the Ting’s model and are thus defined as “apparent”. The main advantage of the method described here is that it has a firm theoretical foundation for viscoelastic analysis – Ting’s model – and can be applied directly to conventional AFM force curves. Therefore, it shares the limitations and uncertainties of the usual AFM experiments, including cantilever stiffness and optical lever sensitivity calibration5. Other limitations are absence of strong adhesion forces during the retraction; drifts and active cell response at low indentation speeds; relatively high computational costs. The method could be applied to the most indentation experiments (ramp-hold, sine function with a single peak), where indentation and force histories are available, conducted not only with AFM. All indentation data points are used for processing and no lag times before measurements are required. Measurement time may be adjusted by varying the piezo speed; fast force curves might be analysed with appropriate sampling rate. Different viscoelastic relaxation models can be used and compared. The method may further benefit from instruments with direct measurement of cantilever deflection and velocity as in laser Doppler vibrometry59.

This work clearly demonstrates that elastic contact mechanics is insufficient to describe force curves on live cells. What, if any, are the drawbacks to using E Hertz to distinguish cell lines and phenotypes? Because of the ongoing relaxation, the Hertz’s model provides some effective, or apparent, value of Young’s modulus, E 0 < E Hertz  < E ∞, which depends on the indentation rate. Since E Hertz is determined by the combination of E 0 and α values, different cells with the same E Hertz value can have different combinations of E 0 and α values. Another example, in the case of cancer cells, low E Hertz values might be caused by both low E 0 and high α values. Thus, E Hertz cannot provide a unique characterization of the viscoelastic phenotype of cells.

Here we used consistency between output results as a criterion for the applicability of a relaxation model: the model describes data well if it yields the same viscoelastic parameters (τ or α) for different experimental conditions (different piezo speeds). The SLS model provided consistent results for PAAm hydrogels for indentation times > 0.1 s (indentation rates < 4 μm/s), but failed at shorter times, where other relaxation mechanisms are probably manifest (Fig. 3g). The PLR model described data for cells better than the SLS model for indentation times from 50 ms–6 s, corresponding to piezo extension rates of 64–8000 nm/s, which cover most of the general AFM indentation experiments1, 34, 57. However, this finding might not hold outside of this range of experimental conditions. Here, an increase in power law exponent was observed at the lowest used indentation time of ~25 ms. Indeed, it was shown that the PLR model with a single exponent does not cover complete time/frequency ranges in rheological experiments on cells, where two or three regimes with different power law exponents were observed45, 60,61,62. In this study, at the highest used piezo speeds, we likely see a transition to the regime with the higher power law exponent (0.75), associated with lateral bending fluctuations of semiflexible cytoskeleton filaments found both in cells and in F-actin solutions62,63,64.

The basis of the PLR model is the so called soft glassy rheology (SGR) theory, developed for soft glassy materials65. This theory assumes that the observed general scale-free behaviour is a natural consequence of disorder and metastability of material internal structures. The cytoskeleton may represents such a structure in cells66. It consists of many disordered elements, which are held together by weak attractive forces and, as a result, trapped in energy wells. Scale-free (power-law) rheology arises from a wide distribution of energy well depths and element lifetimes. The power-law exponent α is related to the effective temperature of the material (amount of the agitation energy in the system), which determines the probability of elements to jump between the energy wells and reflects the system’s dynamics. The jumps of elements between wells are the origin of fluid-like behaviour, and the higher effective temperature (power law exponent) leads to more pronounced fluid-like features of the material. Thus, in the limit of α = 0 soft glassy materials behave like elastic solids and in the limit of α = 1 behave like viscous liquids. The measured exponent value of live cells here and before66 suggests that their rheology resembles that of a soft glassy material close to the glass transition. SGR theory implies that changes in the level of internal disorder and the effective temperature associated with cytoskeleton contraction or remodelling modulate cell rheological behaviour, and thus can provide a conceptual framework for such processes like cell migration, wound healing, invasion, metastasis and embryonic development66. Past studies also found relationships between PLR parameters in different cell groups2, 54, 66, 67. The stiffest group displayed the lowest power-law exponent, whereas the softest group displayed the highest. Here we also can observe such correlations between modulus E 0 and α for groups of benign and cancerous cells, however, not inside the single cell line population. Differences in cytoskeleton structure between benign and cancerous cells are well known, but the relationships between cytoskeleton structure and E 0 and α parameters requires additional research.

Other relaxation models also could be used with the present method in a straightforward manner using corresponding E(t) functions. Generalized Maxwell model with two or more relaxation times increased the quality of the fit relative to the single relaxation time like in SLS model44, 47, 50. But it is not completely clear if it is really capturing some physical phenomena or quality of the fit increased just due to the introduction of additional fitting parameters. Poroelasticity based viscoelastic models could also be used with an approximate solution obtained by finite-element simulations48, 68, 69. However, poroelastic effects seem to be relatively small at the shallow indentation and timescales used here. Fractional variation of the SLS model was recently used for description of cell rheology70. Mathematical analysis shows that it produce relaxation spectra similar to that of the PLR model71, and the choice of the model is more a question of convenience and simplification of the analysis.

It should be noted that adhesion forces are not accounted for in the presented Ting’s model. However, in experiments conducted here adhesion forces were low (the ratio of the maximum adhesive force to the maximum loading force was below 5%) and presumably could be safely ignored. Larger adhesive forces, however, may be responsible for the prominent part of the approach-retraction hysteresis and more complicated models should be used for the description of F-Z curves72,73,74,75,76. Yet, sometimes it is possible to diminish the adhesive force between the probe and the cell surface by proper cleaning and hydrophobic modification of the probe38.

In summary, we presented here general method bringing rigorous viscoelastic analysis within the reach of conventional AFM force curves. The method is highly robust, does not require changes in experimental design and is compatible with findings from conventional AFM microrheology. We believe it greatly enhances the potential of wider use and adoption of AFM for cell mechanics assays.

Material and Methods

Cells

NIH 3T3, MDA-MB-231, and MCF-7 cell lines were cultured in Dulbecco’s modified Eagle’s medium (DMEM) containing 10% FBS and 1% antibiotic/antimycotic solution (Invitrogen, Carlsbad, CA) in a humidified 5% CO2 atmosphere at 37 °C. For the NMuMG cell line, DMEM containing 10% FBS, 2 mM L-glutamine, 100 U/ml penicillin–streptomycin, and 10 μg/ml insulin (Sigma-Aldrich, USA) was used. Experiments were also conducted with the line of MDA-MB-231 cells expressing Syk-EGFP upon incubation with tetracycline analog77; and the expression of Syk-EGFP was induced where indicated by the addition of 1 μg/mL doxycycline for 24 h. Epithelial-to-mesenchymal transition (EMT) was induced in NMuMG cells by addition of TGF-β1 (R&D Systems, USA) in a final concentration of 10 ng/mL for 48 and 96 h. Prior to AFM experiments, cells were plated on glass-bottom cell culture dishes (FluoroDish, WPI, FL) and grown for a period of 1–2 days to a final confluence of ~70%. NMuMG cells were also prepared at higher confluence (monolayer) for EMT experiments.

Gel preparation

Solutions of acrylamide and bis-acrylamide in PBS were polymerized with 10% ammonium persulfate and TEMED (Sigma-Aldrich, USA). Concentrations of acrylamide (10%) and bis-acrylamide (0.01%, w/v) were chosen to prepare a PAAm gel with approximately 1–2 kPa Young’s modulus, which is close to the eukaryotic cell stiffness. For AFM experiments, a total solution volume of 75 μL was allowed to polymerize between circular 18 mm (top) and square 24 mm (bottom) glass coverslips (VWR, USA). The bottom coverslip was aminosilanized with 3-aminopropyltriethoxysilane (APTES, Sigma-Aldrich, USA) and further coated with 0.5 (v/v)% glutaraldehyde solution in PBS (Fischer Scientific, USA) for at least 30 minutes. The upper coverslip was silanized with dichlorodimethylsilane (DCDMS, Sigma-Aldrich, USA). After 30 minutes of polymerization, the top coverslip was removed and the sample was extensively washed with PBS. Prepared PAAm gels had a ~200 μm thickness. The PAAm gels were kept in PBS for several days before the experiments to reach equilibrium. AFM measurements were performed in PBS containing 0.1% Triton X-100 detergent (Sigma-Aldrich, USA) to decrease probe-gel adhesion.

AFM

AFM measurements were performed using a commercial MFP-3D-Bio atomic force microscope (Asylum Research, Oxford Instruments, USA) mounted on an IX-71 inverted optical microscope (Olympus, Japan). The AFM is equipped with a heated stage and in experiments with cells the sample temperature was kept constant at 37 °C. Tipless AFM cantilevers CSC38 (rectangular, Micromash Inc., Estonia) or BL-TR400PB (triangular, Olympus/Asylum Research, USA) were modified with 5 μm diameter silicon dioxide beads (Microspheres-Nanospheres, Corpuscular, NY, USA). The bead was glued to the end of the cantilever using UV-curable glue (Optical Adhesive No. 71, Norland Products, USA) under control of the inverted optical microscope. The typical spring constant of both cantilevers is 0.02–0.05 N/m. The accurate value was determined using the laser Doppler vibrometry system (Polytec MSA-400 Micro System Analyzer from Polytec GmbH, Waldbronn, Germany) using the equipartition theorem as described in refs 78 and 79. The radius of the probe was calculated after scanning the test grating TGT01 (Micromash Inc., Estonia). Before and after measurements, the relationship between the photodiode signal and cantilever deflection (sensitivity factor S) was calibrated by recording several force curves at a bare region of the glass coverslip and averaging its slope. Typically, F-Z curves were taken at 2 μm/s piezo displacement speed along the Z axis using the closed loop feedback circuit, so the piezo movement was corrected with the capacitive sensors implemented into the scanner of the AFM. Data sampling rate was 2 kHz in most experiments, but higher for higher piezo speeds. The piezo displacement range was 2–5 μm (typically 3 μm) to obtain a non-contact region that is long enough for the baseline determination. The force set point was chosen individually for all samples to obtain maximal indentation depth around 500 nm. At least 60–70 cells per cell line (3 F-Z curves per cell above its central part, where the cell height is large enough to diminish substrate effects), and 3 PAAm gel samples (10 F-Z curves at 3 random locations for each sample) were analysed. In experiments with varied scanner displacement speeds the force curves were acquired in a random order, 3 curves per each speed.

The protocols for stress relaxation and microrheology experiments are described in the Supplementary Information Section C.

Adaptation of Ting’s solution for the AFM force curve processing

Raw AFM F-Z curves are presented as raw photodiode signal (ΔV) measured in Volts versus Z scanner displacement in nm (displacement of the cantilever base). First steps in the F-Z curves preprocessing are generally the same and independent of the mechanical model used, including conversion of the ΔV signal to force units (nN), determination of the zero force level and the contact point position. The force acting on the cantilever is calculated by Hooke’s law as F = k * q, where k is a cantilever spring constant and q is the cantilever deflection. The latter one is determined from the raw photodiode signal (ΔV) with the sensitivity factor S [nm/V]: d = S * ΔV.

The zero force level (baseline) was determined from the precontact region of the approach curve with the linear fit, then subtracted from the whole curve. The contact point was determined by the bi-domain fitting procedure36. Part of the approach curve, including both the precontact and contact regions, was fitted with a two-regime regression model (linear for noncontact and Hertz model for contact region), both Young’s modulus in Hertz’s theory framework (E Hertz ) and the location of the contact point (Z 0) were simultaneously estimated, and indentation depth is calculated as δ = Z − Z0 − d. After that step, the resulted curves represent force vs indentation dependencies (F-δ curves).

Here, we also included correction for the effect of hydrodynamic drag forces for the F-Z curves obtained at piezo speeds higher than 2 μm/s. This effect manifested as the separation between precontact regions of approach and retraction curves, which increases with the piezo speed. We adapted the procedure from80 to account for the hydrodynamic forces. As shown in previous work80, 81, the hydrodynamics forces are proportional to the probe velocity and the probe-sample separation. First, the baseline was determined as a median between precontact regions of approach and retraction curves, then these regions were independently fitted with the polynomial function (second order), normalized per the probe velocity, and then calculated hydrodynamic force were subtracted from both approach and retraction curves (Fig. S1). Otherwise, the hydrodynamic drag will make some contribution into the approach-retraction hysteresis of the contact part.

The Hertz’s model, which is traditionally used for description of the indentation with a spherical indenter, has the following form:

$$F(\delta )=\frac{4\sqrt{R}}{3(1-{\nu }^{2})}{E}_{Hertz}{\delta }^{\frac{3}{2}},$$
(5)

where R is the effective radius of curvature of the probe-sample system, 1/R = 1/R probe  + 1/R sample , or, in an assumption of the flat sample surface, it is just a radius of the spherical probe. E Hertz is the Young’s modulus and ν is the Poisson’s ratio of the sample (the probe is assumed to be infinitely rigid in comparison to the soft samples like cells). Models for other probe geometries were introduced later33, 37 and can be used in corresponding experimental setup (conical, pyramidal, cylindrical indenter and others).

Lee and Radok utilized the elastic–viscoelastic correspondence principle to obtain the solution for indentation of a viscoelastic sample when contact area increases with time82. According to this principle, if the solution of the elastic problem is known, then the solution to the corresponding viscoelastic problem may be found by replacing elastic modulus after introduction of appropriate hereditary integral operator:

$$F(t,\,\delta (t))=\,\frac{4\sqrt{R}}{3(1-{\nu }^{2})}{\int }_{0}^{t}\psi (t-\xi )\frac{\partial {\delta }^{\frac{3}{2}}}{\partial \xi }d\xi ,$$
(6)

where ψ(t) is a sample relaxation function and ξ is the dummy time variable required for the integration. Henceforth, we will use Young’s modulus relaxation function, or, in other words, time-dependent Young’s relaxation modulus, E(t)83 to preserve consistency with the initial Hertz’s model equation. Young’s (E) and shear modulus (G) are related through G = E/[2(1 + v)] and, since for cells Poisson’s ratio v = 0.5 is widely accepted, G(t) = E(t)/3. Accordingly, a different prefactor should be used in the Eq. (6) for the shear relaxation modulus. As mentioned earlier, the Lee and Radok solution holds true only for the case when contact area increases with time, for the approach and holding phases in the AFM indentation experiments. But as we showed here, usage of only the approach phase data for extraction of the viscoelastic parameters might lead to erroneous results (see the Supplementary Information Section F for details).

Ting presented a more general approach to solving viscoelastic contact problem32. This approach can be applied to an arbitrary history of the contact area. For the general approach-retraction indentation cycle with a spherical indenter, Ting’s solution of the considered viscoelastic problem has the form presented in equation 1. The first part of the equation is valid for the approach part of F-δ curve (0 ≤ t ≤ t m , t m is a time then approach phase ends; contact area a(t) increases during this phase) and it is equal to the solution of Lee and Radok. The second part of the equation is valid for the retraction part (t > t m ). Here, the auxiliary function t 1 is introduced and defined by a(t 1) = a(t), t 1(t) < t m , meaning that the contact area at time t during retraction is equal to the contact area at time t 1 during approach. t 1 values should be found that satisfy the formula in the equation (2). Using the t 1 (t) function, the contact area a(t 1 (t)) and the effective indentation δ 1 (t 1 (t)) during retraction were calculated. The effective indentation could be imagined as the indentation relative to the relaxing sample surface, as opposite to the indentation δ relative to the initial sample surface before relaxation. The modelled and experimental F-δ curves (Figs 2c and 4) indicate that the force during retraction becomes zero before indentation δ returns to the zero value, but zero force means zero contact area and, therefore, zero effective indentation (see also Supplementary Movie 1). This occurs when sample relaxation rate is slower than the cantilever retraction speed, also leading to the approach-retraction hysteresis.

The equations (1) and (2) provide a relation between the force, F(t), the indentation, δ(t), which are both known experimentally, and the relaxation function E(t). The rate of the indentation raised to the power 3/2 \(\frac{\partial {\delta }^{\frac{3}{2}}}{\partial \xi }\) was calculated by numerical differentiation of the \({\delta }^{\frac{3}{2}}(t)\) data and then smoothed with the moving average filter (5 points). It should be noted, that during general AFM indentation experiments, the vertical piezo (cantilever base displacement) speed is constant, so both loading and indentation speed vary nonlinearly. The force changes nonlinearly with time because of the nonlinear change of the contact area, and the indentation speed during indentation is less than the piezo speed because of the increasing cantilever deflection (also nonlinear) during the indentation.

In principle, it is possible to extract E(t) function by replacement of the integral with the Riemann sum and then by calculation of E(t) for each time point. However, due to the fact that the time derivative of the indentation data is noisy and not quite accurate in the proximity of contact point, the produced error is difficult to evaluate. Furthermore, even after acquisition of a result for E(t), it should be compared with known viscoelastic models. An alternative model-based method we propose is to directly begin with the chosen viscoelastic model for E(t) and to see if this model can represent well the experimental force curve after implementation of the fitting procedure. Here, we compared two relaxation models, which are quite simple but yet frequently used to describe viscoelasticity of cells, namely, standard linear solid (SLS, also known as the three-element model, the Zener model and Prony series with one element) and power-law rheology (PLR) models, described by the equations (3) and (4) respectively. For the SLS model E 0 is the instantaneous modulus [E 0 = E(t → 0)], E ∞ is the infinite (long-term, equilibrium) modulus [E ∞ = E(t → ∞)] and τ is the relaxation time of the material. For the PLR model one can find different representations in the literature from as simple as E(t) = E 1 t −α (E 1 is modulus at t = 1 s)39. Here we used the following expression for the modified power law71, 84, which removes zero time singularity and allows more direct comparison with the SLS model:

$$E(t)={E}_{\infty }+\frac{{E}_{0}-{E}_{\infty }}{{(1+\frac{t}{t^{\prime} })}^{\alpha }},$$
(7)

where E 0 and E ∞ have the same meaning as in SLS model, t′ is a small time offset (here was set equal to the sampling time 5*10−4 s), and α is the power law exponent. E ∞ is often considered to be 0 for cells especially in the rheology experiments, here we also adopted this assumption. We have also tried to leave it as a fitting parameter, and values close to 0 were obtained in most cases (data not shown). So the final expression looks simpler and only two parameters E 0 and α determine the relaxation behaviour (Equation 4). Since t′ parameter is small relative to the time of the experiment, this equation is almost equal to the initial power law equation E(t) = E 1 t −α, and E 1 could be calculated as E 1 ≈ E 0 t′α.

To extract the viscoelastic parameters from the experimental data, we used a fitting procedure with adaptations as described below. In the fitting procedure, the least squares error between the AFM indentation data and the models described above is minimized,

$$e=\sum _{i=1}^{n}{({F}_{i}^{expt}-{F}_{i}^{Ting})}^{2},$$
(8)

where e is an error function (norm of residuals), which together with coefficient of determination R2 provides the information about goodness of fit with chosen relaxation model. The benefit of such method is that all data points in curve are used, and if some of them are affected by noise, the contribution of the noise will be diminished proportionally to the total number of data points. The non-linear Levenberg–Marquardt or the trust-region-reflective least squares algorithms is often used for such nonlinear optimization problem. The latter one showed better results in our preliminary convergence analysis and therefore was used in this study. The algorithm was implemented in the MATLAB (MathWorks, USA), and all integrals were evaluated numerically in the same software. Implementation of numerical analysis might require high sampling rates (small time steps Δt), but here we did not noticed effects of sampling rates on obtained results in the range 100–4,000 Hz (Δt = 0,01–0,00025 s) for 2 μm/s piezo velocity, so 2 kHz sampling frequency was used. It gives approximately 1 point per 1 nm of piezo displacement, the same ratio was preserved for the higher piezo velocities.

We noticed that the best fit was obtained by adjusting the contact point position slightly. The initial estimate for the contact point was obtained with the elastic assumptions of the Hertz’s model and might not be equal to the optimal contact point position for the viscoelastic model. Henceforth, the contact point position was varied in the region close to the initial estimate, then the solution providing the lowest norm of residuals was chosen as the final position85 (see the Supplementary Information Section A for details). Such behaviour is probably associated with the uncertainties in the beginning of the approach phase caused by long-range electrostatic forces, steric forces from the brush layer on the surface of cell86, surface roughness, or contamination of the probe surface – all factors that mostly affect the beginning of the approach curve, but whose contributions decrease with indentation depth. Indeed, the resulting fit shows that larger residuals are observed in the beginning of the curve (Fig. S3). In experiments with varied indentation depths (Figs S10 and S16), strong variation of viscoelastic parameters were also found at small depths.

The algorithm presented here could be combined with bottom-effect correction models for the sample of a finite thickness in a straightforward manner (see the Supplementary Information Section A for details). The custom MATLAB codes used in the analyses are available from the corresponding author upon request.

Finite element analysis

Finite element (FE) simulations were conducted in Abaqus CAE (version 14, Simulia Corp., Providence, RI). The axisymmetric system consisted of the rigid spherical indenter with 2 μm radius and the viscoelastic cylindrical sample with 15 μm height and radius (Fig. 2). The indentation depth and speed were selected to be 500 nm and 1 μm/s respectively, both for approach and retraction. The sample was meshed more finely in the contact region. The bottom surface of the sample was constrained, and the contact between the indenter and the sample was considered to be frictionless.

E 0 of the sample was set to 2 kPa. For SLS model τ was set to 0.1 s and E ∞ to 1 kPa. PLR model could not be directly prescribed in this FE software, so it was approximated as the Prony series expansion including six terms, with coefficients adjusted in MATLAB. The effectiveness of such approach was showed previously87, 88. Power law exponent α was chosen to be 0.2. More details could be found in the Supplementary Information Section B.

Statistical Analysis

All statistical analyses were performed using OriginPro 2016 software (OriginLab Corporation, MA). A 2-tailed Pearson’s correlation coefficient was used to characterize the correlation between measured viscoelastic parameters. A non-parametric Mann–Whitney U test was used to determine the statistically significant differences between the groups. Since most of the data were distributed not normally, results are presented both as mean ± standard deviation and median ± median absolute deviation. The percentiles in the box-and-whisker plots are 10%, 25%, 50%, 75% and 90%, the dots are the outliers.