Validation of contact mechanics models for atomic force microscopy via finite elements analysis and nanoindentation experiments

Abstract

In this work, we have validated the application of Hertzian contact mechanics models and corrections for the analysis of force vs. indentation curves, acquired using spherical indenters on linearly elastic samples, by means of finite elements simulations and AFM nanomechanical measurements of polyacrylamide gels possessing a thickness gradient. We have systematically investigated the impact of both large indentations and vertical spatial confinement (bottom effect) on the accuracy of the nanomechanical analysis performed with the Hertz model for the parabolic indenter compared to the Sneddon model for the spherical indenter. We demonstrated the accuracy of the combined corrections of large indentations and bottom effect for the Hertz model proposed in the literature in the framework of linearized force vs. indentation curves acquired using spherical indenters, as well as a validation of a new linearized form of the Sneddon model. Our results show theoretically and numerically, as well as through nanomechanical measurements on gels, that the corrected Hertz model allows to accurately quantify the Young's modulus of elasticity of linearly elastic samples using spherical tips in a wide range of indentation and thickness values that are relevant for the study of cellular systems.

---

1. Introduction

The micro- and nanoscale characterisation of the mechanical properties of systems and device components is of increasing importance in several fields like biology and biophysics, where the elastic properties of cells, tissues, and extracellular matrix (ECM) can affect the behaviour and fate of an organ,^1–6 or medical engineering, where microdevices with specific elastic properties are employed.^7,8

Atomic Force Microscopy (AFM) is an ideal tool for the quantitative and non-destructive characterisation of the mechanical properties of biological and non-biological samples at the sub-micrometre scale.^9–11 AFM provides high spatial and force resolution and is very versatile, including the ability to work in physiological solution and controlled environments, the freedom of choosing the best tip dimensions and geometry to match the typical length scales of the system under investigation.^12–14

While sharp AFM tips are mandatory when high spatial resolution is necessary, spherical probes (colloidal probes, CPs) possess several characteristics that make them suitable for the investigation of mechanical properties of soft or biological samples.¹⁴ Indeed, a well-defined interaction geometry (sphere on flat, sphere on sphere) and a reduced stress and strain in mechanical tests, make the application of contact mechanics models more reliable, and the interpretation of results less ambiguous, and the selection of the tip radius to match the characteristic length scales of the system under investigation provides better averaging of the mechanical signals together with a system-adapted spatial resolution.^13,15

A simple contact mechanics model describing the deformation of elastic solids was proposed by Heinrich Hertz in 1882.^16,17 The Hertz model can be used to describe the indentation of a purely elastic, semi-infinite half space under the pressure exerted by a paraboloidal indenter. It is well known that the application of the Hertz model to real systems is based on a series of assumptions and, in some cases, gross approximations: the sample must be uniform and isotropic, interfacial adhesion must be absent, with no friction between the indenter and the sample, the strain and stress must be sufficiently small to ensure a linear elastic response and the sample must not exhibit a constrained mechanical response due to its spatial confinement and/or finite dimensions (the bottom effect problem^18–21). In addition, to apply the Hertz model to data collected using spherical indenters, the indentation δ of the probe must be small compared to the radius R of the tip (δ ≪ R), otherwise the parabolic approximation of the spherical profile will be inaccurate (the large indentation problem^22,23).

Microscopic systems like single cells or thin tissue slices for histological analysis make the assumptions behind the use of the Hertz model to fit nanoindentation data hard to be respected. Cell height usually varies between 5–15 µm and its internal structure is extremely heterogeneous in all three dimensions, forcing to indent up to a few microns to characterise the overall mechanical response, and not just the elastic contribution of the cell membrane coupled to the actin cortex.^15,24

Achieving large indentions up to a few microns is even more critical for tissue or tissue-derived samples. Indeed, these samples are usually far more heterogeneous than cellular systems, across a broader range of length scales (from 50 nm to 100 μm), which requires probing a larger volume with the indenter. In addition, their surface may present micron-scale roughness due to the above-mentioned structural complexity and to the slicing process, which produces irregular interfaces, due to cell detachment and disentangled ECM fibres.

Working with cells or other finite-thickness samples (like tissue or ECM slices) exposes therefore to both large indentation and bottom effect issues, even if nonlinear effects are excluded; these effects impact on the accuracy and precision of mechanical parameters obtained using the Hertzian contact model.

To mitigate the large indentation issue in the framework of the Hertz model, it is possible to increase the radius R of the tip to reduce the δ/R ratio; this measure however can limit the maximum indentation achievable or mandate the use of very rigid cantilevers, with consequent loss of force sensitivity, beside causing a severe loss of spatial resolution. Alternatively, the model developed by Sneddon to describe the indentation of an elastic half-space by a spherical indenter²⁵ can be used; this model (hereafter simply called the Sneddon model) does not suffer from the constraint δ ≪ R. Unfortunately, Sneddon's equation for the contact radius of the spherical probe cannot be cast in an analytic close form but requires numerical methods to be solved (see section 2.2).

A typical mitigation of the bottom effect typically consists of limiting the maximum indentation to a small fraction (well below 1/10) of the sample thickness, or height, which may keep from sensing the elastic contributions of the deeper layers of the system under investigation. To overcome the above-mentioned limitations of the Hertz theory, the Bottom Effect Correction (BEC)^18–20,26 and the Large Indentation Correction (LIC)^22,23 have been proposed. The aim of such corrections is to extend the applicability of the Hertz model for the paraboloidal indenter to those cases where the estimation of the Young's modulus E is biased due to the bottom effect and the large indentations, respectively.

Despite the increasing occurrence of application of BEC in published reports (LIC is by far less considered), a systematic validation and a suitable integration of both correction methods in a single experimental and data analysis approach are still missing. This is also due to the difficulty in producing suitable test samples for the nanomechanical investigation allowing to probe the influence of the varying thickness, spanning the range from the strong spatial confinement to the bulk conditions. The scarce availability of reference samples for nanomechanical tests limits the standardisation of the experimental and data analysis protocols.^27,28

In this work, we have validated the application of Hertzian contact mechanics models and corrections in the framework of linear elasticity for the analysis of force vs. indentation curves (FCs), acquired using spherical indenters, by means of simulations using finite element analysis (FEA) and AFM nanomechanical measurements of polyacrylamide gels possessing a thickness gradient. FEA is a suitable tool for the simulation of the mechanical response of materials subject to the stress applied by an indenter under controlled conditions.^26,29–32 We have used FEA to simulate FCs on ideal elastic samples.^29,30

We here discuss the accuracy of the Hertz model^33,34 compared to the Sneddon model in bulk systems and in conditions of vertical spatial confinement. We assess the accuracy of the existing large indentation and bottom effect corrections for the Hertz model. We present a simplified linearised version of the Sneddon model, which can be implemented efficiently in the same data analysis framework we developed for the use of the linearised Hertz model. Eventually, we demonstrate that LIC and BEC can be coupled in a single correction function for the Hertz model and provide an experimental validation of the proposed methods based on AFM nanoindentations measurements with colloidal probes on polyacrylamide test samples, fabricated on purpose to possess a gradient of thickness, in a wide range of indentation and thickness values that are relevant for the study of cellular systems.

2. Theoretical models and corrections

We present here a concise summary of the theoretical models used in this work, under the assumption of linear elasticity. For a deeper presentation and discussion of contact mechanics models, the reader is referred to some recent review papers, see ref. 17 and 35.

2.1 Hertz model

The Hertz model describes the non-adhesive, frictionless contact between two uniform, isotropic elastic bodies, within the framework of linear elasticity (small strain and stress), thus considering the deformations small compared to the dimensions of the bodies.³⁶ Originally developed for describing the elastic contact between two spherical bodies,¹⁶ the Hertz model strictly applies to the case of a paraboloidal indenter. It can be used to describe indentation by a spherical indenter provided the indentation is small compared to the tip radius (δ ≪ R); in this case the parabolic profile represents a fairly accurate approximation of the circular one. In the case of indentation by an infinitely rigid paraboloid of an elastic half space with Young's modulus E and Poisson's ratio ν, the relation between the applied force F_Hertz and the indentation δ is:where R is the radius of curvature of the indenter, identified in apical region of the tip. The contact region has circular cross section, and a simple relation links the tip radius R, the indentation δ and the contact radius a (Fig. S1a, blue circles):

From eqn (2) it follows:

2.2 Sneddon model for the spherical indenter

Sneddon developed a solution³⁷ for the case of the non-adhesive, frictionless indentation of an elastic half-space by a rigid spherical indenter, represented by the system of eqn (3) and (4):

The Sneddon model is not subject to the limitation δ ≪ R. The drawback of this model is the lack of an analytical solution linking the force to the indentation, since eqn (4) cannot be inverted analytically to obtain the relation a = a(δ); hence, a numerical solution is required. It is worth noting that eqn (2), valid for the Hertzian contact, represents the limit of eqn (4) when a/R → 0, which justifies the use of the simpler Hertz model when using a spherical indenter as long as the indentation δ is sufficiently small compared to the radius R (see Fig. S1a, red line).

2.3 Large indentation corrections for the Hertz model

The Hertz model is simple (analytic and linearizable, see below), and it is therefore desirable to use it to fit indentation data obtained with spherical indenters and arbitrarily large indentations, typically up to δ = R. To this purpose, in this work we have selected two different corrections, presented in the papers published by Kontomaris et al.²² and by Müller et al.²³

Both correction functions 1/Ω_K/M depend on the nondimensional ratio γ_R(δ) = δ/R and transform the experimentally measured force–indentation curve F(δ) into an apparent Hertzian curve F_Hertz(δ) (eqn (1)), which can be fitted by the Hertz model across the whole range of indentations:

The function Ω_K presented in Kontomaris et al.²² (from now on called LIC_K) is a power series in γ:

with coefficients: c₁ = 1.0100000, c₂ = −0.07303003, c₃ = −0.1357000, c₄ = 0.0359800, c₅ = −0.0040240 and c₆ = 0.0001653.

The function Ω_M presented in Müller et al.²³ (from now on called LIC_M) is a polynomial in γ:

Both correction functions act pointwise on the force curve, since their values depend on the indentation through the ratio γ (see Fig. S2a). The contact radii corrected for the large indentation effect are in excellent agreement with the Sneddon's radius (see Fig. S1a and Note S1).

2.4 Bottom effect corrections for the Hertz model

The presence of a rigid substrate underneath the sample represents a boundary condition, which, in the case of a large indentation to thickness ratio, can cause strong perturbation of strain and stress fields and therefore influence the obtained elastic modulus.^17,28 The bottom effect is a special case of the more general case of three-dimensional spatial confinement of the elastic body, where the strain and stress fields are constrained also laterally.³⁸ This could be relevant, for example, when cells are at confluence, tightly arranged within a nearly two-dimensional monolayer, and firmly connected through cadherin bonds.

We have selected two correction functions for the bottom effect, developed for the paraboloidal indenter, in the case of linear elastic response: the first one (1/Δ_D) was developed by Chadwick, and published in Dimitriadis et al.;¹⁸ the other correction function (1/Δ_G) was developed by Garcia and published in Garcia et al.¹⁹ (see Fig. S2b). The correction presented in the paper by Long et al.,²⁶ which considers the vertical confinement for a non-linear neo-hookean material, as well as others, have not been investigated in this work. Nevertheless, such polynomial corrections, as well as others LICs, could be readily incorporated within the presented protocol. These corrections can be applied, similarly to the large indentation case, to transform the experimentally measured force–indentation curve F(δ) into an apparent Hertzian curve F_Hertz(δ), describing the case of an infinitely extended elastic half space across the whole range of indentations:

The functions presented by Dimitriadis et al.¹⁸ (from now on called BEC_D) for a paraboloidal indenter refers to the following boundary conditions (BCs): (i) a sample bonded to the rigid substrate (Δ^bonded_D); (ii) a sample that is allowed to slide over it (Δ^unbonded_D); (iii) a sample that is partially bonded to the substrates (Δ^mixed\_BCs_D, like for single cells¹³):

Δ^bonded_D = 1 + 1.133χ + 1.283χ² + 0.769χ³ − 0.0975χ⁴ (9)

Δ^unbonded_D = 1 + 0.884χ + 0.781χ² + 0.386χ³ + 0.0048χ⁴ (10)

Δ^mixed\_BCs_D = 1 + 1.009χ + 1.032χ² + 0.578χ³ − 0.0464χ⁴ (11)

Following Garcia et al.,¹⁹ The sign of the 4^th order coefficient in eqn (9) has been changed with respect to the original work of Dimitriadis et al. (the value of the 4^th order coefficient in eqn (11) has been also changed, accordingly); the same corrections should be applied to eqn (2)–(4) of our previous work.¹³ In these equations, the driving nondimensional parameter is , where h is the height (or thickness) of the sample measured with respect to the hard substrate. The parameter χ represents the ratio between the Hertzian contact radius a (see eqn (2)) and the sample height h, which suggests that bottom effect is negligible in the limit a ≪ h (like for sharp probes, even though a higher stress would be produced). Noticeably, the bottom effect depends on the ratio of horizontal (a) to vertical (h) dimensions, and not simply on δ/h, which in turn reminds us that the elastic deformation of a body is a truly three-dimensional process. Moreover, this also suggests that large spherical tips are more prone to bottom effect than sharp pyramidal ones.

The function Δ_G presented by Garcia et al.¹⁹ (from now on called BEC_G) for a paraboloidal indenter on a sample bonded to the rigid substrate is:

Δ_G = 1 + 1.133χ + 1.497χ² + 1.469χ³ + 0.755χ⁴ (12)

It is worth noting that the corrective factors of the BECs reported in eqn (9)–(12) are computed based on the assumption of incompressibility of the sample, thus imposing a Poisson's ratio ν = 0.5. Since the FEA simulations are performed for an elastic medium with Poisson's ratio ν = 0.49, we report the corrective factors Δ_D and Δ_G for the bonded case with ν = 0.49, truncated to the fourth order (eqn (13) and (14)). While Dimitriadis et al.¹⁸ report the equations to calculate the BEC coefficients for an arbitrary Poisson's ratio, in the work of Garcia et al.¹⁹ explicit formulae are not present; from the equations reported by Garcia et al., we derived approximate expressions for the BEC coefficients for an arbitrary Poisson's ratio (see Note S2 for details).

Δ_D,_ν=0.49 = 1 + 1.112χ + 1.237χ² + 0.712χ³ − 0.133χ⁴ (13)

Δ_G,_ν=0.49 = 1 + 1.112χ + 1.444χ² + 1.374χ³ + 0.645χ⁴ (14)

Based on Garcia's formulae, we also derived the expression of the contact radius in confined geometry for a spherical tip up to the 5^th order in χ, extending Garcia's 1^st order result (see Note S2 and Fig. S1b–d). The corrected contact radius is in good agreement with the predictions of finite elements simulations.

3. Materials and methods

3.1 Finite elements analysis

Finite elements simulations were performed using ANSYS Mechanical (ANSYS Student, the free version of ANSYS software), to produce ideal FCs describing the mechanical response of elastic films in different regimes of spatial confinement, i.e., with different δ/R and ratios (Fig. 1a). The large indentations option, which retains the quadratic terms in the deformation tensor,³⁹ was used in the simulations, to secure convergence. The resulting material shows, within the mesh-related error (Note S6, and Fig. S6), ideal elastic behaviour across the whole tested indentation range, for both bulk and thin film configurations, as witnessed by the linearity of the effective strain (ε*) vs. effective stress (σ*) figures of merit (according to Tabor's definition:³⁰ ε* = 0.2a/R, σ* = F/πa² as well as to the similar Kalidindi's definition,^40–42 see Fig. S3, S4 and Note S3). The linearity of the simulated systems was further checked through the equivalent, local von Mises stress and strain (Fig. S5 and Note S4).

Fig. 1 (a) Representation of the simulated system. R is the tip radius, h and L the thickness and the width of the simulated elastic medium, respectively. The medium is defined solely by its nominal Young's modulus E and the Poisson's ratio ν. As the only imposed boundary condition, the bottom face of the media is bonded to a rigid substrate. (b) Comparison of FCs simulated using FEA for the different systems and the different contact mechanics models (the theoretical curve for the TF-S system includes already the bottom effect through the corrective function Δ_D).

We took advantage of the axial symmetry of the system, which allowed to create two-dimensional meshes and take full advantage of the reduced number of nodes available (N_max = 1.28 × 10⁵). The simulated systems consist of a rigid spherical or paraboloidal tip with a wide range of tip radii R indenting an elastic medium of different thickness h to simulate both the bulk and the thin film regimes. The lateral dimensions of the simulated systems L were chosen to be large enough to neglect spatial confinement-related effects (see Table 1);³⁸ the validity of this hypothesis was confirmed for each simulation by observing that the vertical and lateral displacements at the boundaries were negligible.

Table 1 Summary of the different systems simulated for this study. In the last three rows, which correspond to the three Thin-Film systems, each of the three parameters, indentation δ, height h and tip radius R, were varied, while fixing the other two

Model name	Tip shape	Tip radius R [μm]	Sample thickness h [μm]	Sample width L [μm]	Young's modulus E [MPa]	Poisson's ratio ν	Corrections to the Hertz model	Indentation δ [μm]
B-S (Bulk-Sphere)	Sphere	5	500	120	0.5	0.49	LIC	5
B-P (Bulk-Paraboloid)	Paraboloid	5	500	120	0.5	0.49	None	5
TF-P (Thin Film-Paraboloid)	Paraboloid	5	10	120	0.5	0.49	BEC	2.5
TF-S (Thin Film-Sphere)	Sphere	5	10	120	0.5	0.49	BEC, LIC, BEC + LIC	0.05–3.1
TF-S (Thin Film-Sphere)	Sphere	10	8–60	200	0.5	0.49	BEC, LIC, BEC + LIC	1
TF-S (Thin Film-Sphere)	Sphere	1–65	20	50–400	0.5	0.49	BEC, LIC, BEC + LIC	1

---

Four models were studied (see Table 1): Bulk-Sphere (B-S), Bulk-Paraboloid (B-P), Thin Film-Paraboloid (TF-P) and Thin Film-Sphere (TF-S); the first two models mimic the indentation of an infinite linearly elastic half space by a spherical and a paraboloidal indenter, respectively, whereas the latter two models mimic the indentation of a linearly elastic thin film by paraboloidal and spherical indenters, respectively. The material was simulated with nominal Young's modulus E_nom = 0.5 MPa and Poisson's ratio ν = 0.49. The choice of ν = 0.49 was made to simulate a nearly incompressible material⁴³ (the choice ν = 0.5 keeps finite elements simulations from converging). In all models, the elastic layer is bonded to a rigid substrate, and the layer-tip contact is frictionless and non-adhesive. In the range of applied forces, the tip deformation is negligible since its Young's modulus is six orders of magnitudes greater than that of the elastic material. The studied contact mechanics models were fitted to each single FC generated by a single finite element simulation from 0.5% up to 99.5% of the maximum indentation.

The B-S model consists of a rigid spherical tip of radius R = 5 μm indenting a medium with thickness h_bulk = 100R. This model was used to investigate the large indentation effect; the bottom effect is not expected.

The B-P model consists of a rigid paraboloidal tip with radius of curvature R = 5 μm indenting a medium with thickness h_bulk = 100R. This model was used as reference for the large indentation effect. Neither large indentation nor bottom effects are expected for a parabolic indenter.

The TF-P model consists of a rigid paraboloidal tip with radius of curvature R = 5 μm indenting a medium with thickness h_film = 2R. This model was used to investigate the bottom effect. No large indentations effect is expected for a parabolic indenter.

In the TF-S model indentation, thickness and tip radius are varied to investigate both large indentations and bottom effect in different regimes. The details of the different thin film systems simulated are reported in Table 1.

Representative simulated FCs are shown in Fig. 1b. Simulated FCs were fitted using the Hertz model with LIC and BEC and the Sneddon model. The maximum indentation in all the cases shown in Fig. 1 was δ = 3 μm. It is possible to observe the good agreement with the theoretical models in all three configurations. Moreover, the impact of the indenter geometry and of the spatial confinement is negligible if the δ/R ratio is below 0.2.

For larger δ/R values, irrespective to the indenter geometry, the bottom effect becomes dominant; it is evident that it is necessary to apply much higher forces to indent of the same amount a spatially confined material (compare TF-S with B-S and B-P). The Sneddon and the Hertz models represent the reference contact mechanics models to describe the indentation of an elastic half-space by spherical and parabolic indenters (B-S and B-P systems, respectively), irrespective to the value of the maximum indentation.

Comparing B-S and B-P systems, it can be noticed that when indentation increases, the force needed to obtain the same indentation with the parabolic indenter increases; this is due the difference in the contact area (see Fig. S1a). Indeed, the parabolic contact area grows faster than the spherical one. The discrepancy between the simulated TF-S FC and the Hertzian FC corrected for the bottom effect is caused by the fact that the Hertz model is exact for the parabolic indenter, while the simulated curve assumed a spherical indenter. A suitable application of the large indentation correction is supposed to correct further this discrepancy (we will discuss this in section 4.1.2).

3.2 Data analysis

3.2.1 Hertz model. Linearizing the experimental FCs and the Hertz model (eqn (1)) would allow to apply a simple linear regression to the data to extract the Young's modulus. However, both measured force and indentation values (F′,δ′) are typically shifted with respect to the true values (F,δ) in eqn (1). In terms of the measured quantities, the Hertz equation is:

In the case of the force, as long as adhesion and friction are negligible (in principle, the Hertz model assumes no adhesion and friction at all; when measuring in liquid, adhesion is typically very small, as the van der Waals interactions are effectively screened by water), the offset F₀ is due in general to a misalignment of the optical beam deflection system and can be easily subtracted during the pre-processing of the raw force curves. Without loss of generality, we will assume in the following that the offset on the force axis is not present, i.e., F′ = F.

On the other hand, an apparent indentation δ′ = δ + δ₀ is computed when rescaling the distance axis in the raw force curve, where δ₀ represents the unknown location of the contact point; δ₀ must be therefore considered as a free parameter of the fit, as well as the Young's modulus. The presence of the offset δ₀ in eqn (15) does not allow to linearise the equation by taking the logarithm of both sides. Following ref. 13 and 44, we adopt a convenient way to linearise eqn (15) through the variable transformation F* = F^2/3, which leads to the equation:

δ′ = αF* + δ₀ (16)

where the parameter α depends on the Young's modulus E:and eqn (16) holds for δ > δ₀.

The contact point δ₀ represents the intercept of the δ′ vs. F* curve, while the Young's modulus E can be calculated from the slope α, since both radius R and Poisson's ratio ν are supposed to be known. The free parameters of the fit, δ₀ and E, are not correlated in eqn (16), and can therefore be determined independently; this is not the case in the original nonlinear Hertz model (eqn (1)). Details on how to identify the contact point are described in Note S5.

It is straightforward to generalise this approach to include LIC and BEC, or both. The correction functions 1/Ω (eqn (6) and (7)) and 1/Δ (eqn (9)–(12), (13) and (14)), can be integrated into the pseudo force F* as follows:¹³

where eqn (18) corrects for the large indentation effect, while eqn (19) corrects for the bottom effect.

In this work, we propose a way to combine both corrections into a single large indentations and bottom effect correction by defining the pseudo-force F* as:

The resulting rescaled pseudo-Hertzian vs. δ′ curve is expected to be largely unaffected by both the large indentation issue and can then be fitted by the linearised Hertz model (eqn (16) and (17)) across the whole indentation range achieved experimentally. We demonstrate in this work that it is possible to retrieve the correct/intrinsic Young's modulus with an accuracy better than 3%.

As originally developed, the bottom effect compensation factor does not include any large indentation effect, and vice versa, the large indentation correction factors do not include the bottom effect. Due to the multiplicative nature of such corrections, we propose in this work that they can be applied sequentially (eqn (20)), implicitly assuming that the bottom effect can be compensated even in the presence of a spherical tip rather than a paraboloidal one, so that the subsequent application of the large indentation correction can transform the force curve into an equivalent Hertzian curve for a parabolic indenter. This assumption is supported a posteriori by the results of finite element analysis presented in this work, which shows that the Young's modulus is restored with great accuracy upon the simultaneous application of the two corrections as proposed in this work. Our proposed model is an approximation, which consider negligible the coupling between the two effects; this coupling could be represented by an additional multiplicative factor, in the form 1 + f(χ·γ), with the function f depending to the leading order on the product of the two nondimensional factors χ and γ. As it will be shown in section 4.1.3 (Fig. 3f) and in Note S8, neglecting the coupling term may give rise to unexpected (yet small) discrepancies in the large χ·γ regime where, applied separately, BEC and LIC perform reliably.

3.2.2 Sneddon model. If a spherical indenter is used and the bottom effect is negligible, the model developed by Sneddon can be used. This model accounts for arbitrarily large indentations, with the drawback of the lack of analytical solution. We assume here that both the measured force F′ and indentation δ′ contain unknown offset: F′ = F + F₀; δ′ = δ + δ₀. Eqn (3) and (4) are then replaced by eqn (21) and (22):

To obtain the Young's modulus from the experimental F′ vs. δ′ curve, the system of eqn (21) and (22) must be numerically solved, treating E, δ₀ and F₀ as free parameters. This can be done by first solving numerically eqn (22) using an initial guess for δ₀, to obtain the relation a(δ′), then substituting a in eqn (21) and evaluating the quadratic distance between the obtained curve and the experimental data. The best values of the free parameters E, F₀ and δ₀ are those that minimise this distance. This method is time consuming, given that typically several hundreds of force curves must be processed.

3.2.3 Linearised Sneddon model. We propose a faster method that once again takes advantage of linearisation. As for the Hertz model case, we assume that the vertical offset F₀ can be easily determined and subtracted. We then assume that the non-corrected, linearised Hertz model (eqn (16)) provides a rather accurate estimation of δ₀ when applied to the first 10% of the indentation range, i.e. far from the large indentation limit. As for the force offset, also the indentation offset can be subtracted from the indentation axis, which allows to solve numerically only once eqn (4), to obtain a(δ).

We then observe that the force in eqn (21) depends linearly on the variable s* (with dimensions of an area) defined as:

in terms of which eqn (21) can be rewritten as:where we have left for convenience a residual force offset ΔF₀. A linear regression of eqn (24) provides the Young's modulus value.

While the standard approach, treating E, δ₀ and F₀ as free parameters in eqn (21) and (22), requires a complete minimisation procedure for each FCs (the numerical inversion of eqn (22) representing the bottleneck), the new approach requires the numerical inversion of eqn (4), to be done only once, and then a series of linear regressions, which can be parallelised easily due to their algebraic nature, leading to a tremendous cut of computation time: the processing of a sets of FCs, also called force volumes (FVs), consisting of hundreds of FCs, takes several minutes on a standard personal computer when following the standard procedure, and just a few seconds with the linearised one. The free parameters ΔF₀ and E in eqn (24) are independent, therefore a residual offset of the force does not affect the accuracy of the evaluation of the Young's modulus, assuming that δ₀ can be determined with good accuracy.

3.3 AFM nanomechanical measurements

The nanomechanical measurements were performed in liquid (MilliQ water) using a Bioscope Catalyst AFM (Bruker) mounted on top of an inverted optical microscope (Olympus IX71). The whole system is isolated from the ambient noise by an active antivibration base (DVIA-T45, Daeil Systems) located inside an acoustic enclosure (Schaefer, Italy). We used custom colloidal probes, fabricated and calibrated as described in ref. 14, produced attaching spherical borosilicate glass beads with different nominal radii R of 5 and 10 μm to tipless cantilevers (Nanosensors; TL-FM-50, TL-TM-50). The spring constants of the AFM probes, with nominal values of 3–5 N m⁻¹ (FM) and 45 N m⁻¹ (TM), were calibrated using the thermal noise method^45–47 and corrected for the contribution of the added mass of the sphere and for the back shift of the loading point.^48,49 The measured and corrected spring constant were k₁ = 3.52 ± 0.12 N m⁻¹, k₂ = 3.24 ± 0.07 N m⁻¹, k₃ = 159.7 ± 3.26 N m⁻¹ and the corresponding radius of the fabricated colloidal probes after calibration turned out to be R₁ = 4989 ± 35 nm, R₂ = 7309 ± 13 nm, and R₃ = 8898 ± 54 nm.

Both thick (bulk) and variable-thickness polyacrylamide samples, produced as described in Note S7, were studied by collecting FVs in different macroscopically separated locations. Each FV typically covered an area of 150 × 150 µm² and consisted of 15 × 15 or 20 × 20 FCs (Fig. S9a). In each experiment, the maximum applied force was set to achieve a maximum indentation δ ≈ R. Raw curves consist of the recording of deflection signal from the photodetector (in V units) vs. the distance travelled by the z-piezo (in nm); these curves are rescaled into force vs. indentation curves (shortly force curves, FCs) according to standard procedures, as described in ref. 13. To this purpose, the deflection sensitivity of the optical beam deflection apparatus (in units of nm V⁻¹) was calculated using the contactless SNAP procedure,²⁷ using the previously calibrated spring constant as reference. The ramp rate was set at 1 Hz in all experiments.

From each FV, the local topographic map was reconstructed as the map of the contact points from each FC.⁵⁰ The mean-subtracted topographic map represents the map of local height variations around the mean thickness value of the sample in the region where the FV has been collected.

Knowledge of the mean thickness of the reference sample in the region of the FV is crucial to implement the BEC. We measured the local thickness of the PA sample (Fig. S8a) by repeatedly engaging the AFM tip on the film and on the substrate and recording the distance Z travelled by z-microtranslator of the AFM head (used for the coarse approach of the tip to the sample), together with the XY coordinates of the precision motorised stage. The resolution of the stepper motors is as good as 100 nm, therefore allowing to reconstruct a “low-resolution” topographic map Z(X,Y) of the PA film (Fig. S8b). The points measured on the glass slide were used to fit and subtract an eventual baseline of the low-resolution topographic map.

Having characterised the thickness gradient of the PA film, it was then possible to associate to each FV the mean thickness value of the measured portion of the sample by recording the XY position of each FV and looking at the corresponding Z value in the PA film topographic map (Fig. S8b); the local height corresponding to each FC in the FV was then precisely determined by adding to the mean Z position of the FV the local height variation from the mean-subtracted, unflattened topographic map reconstructed from the FV (Fig. S8c).

3.3.1 Fitting procedure of force curves and error analysis. The different contact mechanics models were fitted to each FC in the FV, from 5% up to 99% of the maximum indentation, to get the Young's modulus values. The median values of each FV were obtained by fitting a Gaussian curve to the Young's modulus distributions in semi-log₁₀ scale (Fig. S9c, E_med = 10log₁₀(E/P_a)|_best, where log₁₀()|_best is the centre of the Gaussian curve).¹³

When multiple FV corresponded to the same condition (e.g. Fig. 5a, b and 7b), median Young's modulus values from different FVs were averaged, to obtain a mean median value; the corresponding error was calculated as the standard deviation of the mean^13,51 Individual median Young's modulus values are collected into violin plots.

3.4 Error analysis in FEA

Since the Hertz model is mathematically exact for a paraboloidal indenter on a flat elastic medium for arbitrarily large indentations, we assume that the discrepancy between the nominal Young's modulus (E_nom) and the one obtained fitting the Hertz model to the simulated force curve (E_Hertz) depends only on simulation-related artifacts (such as the finite mesh size, Note S6) and can therefore be assumed as the smallest error that can be obtained. The error associated to the fit of the Sneddon model for the spherical indenter to the simulated curves for arbitrary indentations can be estimated in a similar way. Fig. 2a shows that these errors are well below 1% across the widest range of indentations considered in this work (0 < δ/R ≤ 1).

Fig. 2 Relative discrepancies of the obtained Young's modulus E compared to the nominal one (E_nom = 0.5 MPa) for the bulk simulations for different ranges of indentation. The legends report the contact mechanics models used for the fit and the simulated systems (B-S and B-P) (see also Table 1). (a) Results for the appropriate matching of contact mechanics model and indenter geometry (Hertz model for the parabolic indenter, Sneddon and linearised Sneddon models for the spherical indenter). (b) Results for the Hertz model applied to the B-S system, with and without LIC.

4. Results

4.1 Validation of contact mechanics models and corrections by FEA

We have used FEA to study the range of applicability of the Hertz and Sneddon model, as well as LIC and BEC (the models used in the simulations are listed in Table 1). The Hertz model is exact for the paraboloidal indenter, while the Sneddon model is exact for the spherical indenter; both models, when used for the appropriate indenter geometry, are valid in principle for arbitrarily large indentations, if the material response is linearly elastic.

4.1.1 Hertz and Sneddon models. Artificial FCs produced simulating the B-S and B-P systems were fitted using the Hertz, Sneddon, and Sneddon linearised models up to increasing δ/R ratios. Fig. 2a shows that on B-S systems the result of the fit performed using the Sneddon model agrees within 0.7% with the nominal Young's modulus; a slightly larger maximum discrepancy (1%) is observed when the B-P system is analysed using the simple Hertz model.

Noticeably, the linearised Sneddon (Fig. 2a, purple crosses) model performs similarly, if not better, than the standard Sneddon model, across the whole range of indentations, which represents a validation of our linearised Sneddon approach. Eventually, our results confirm that when the appropriate contact mechanics model is applied with respect to the indenter geometry (i.e., Hertz model over a parabolic indenter and Sneddon model over a spherical indenter), at least for bulk films, the result is independent on the range of indentation selected for the fit, up to δ/R = 1, as far as the system is linearly elastic.

4.1.2 Hertz model and large indentations. The Hertz (section 2.1) and Hertz + LIC (section 2.3) models have been used to fit the FCs for the B-S system. The results are shown in Fig. 2b. When the Hertz model is applied for indentations up to the tip radius, it underestimates the true Young's modulus of the media, with an error of approximately 10% when δ = R; the error decreases almost linearly as the maximum indentation decreases, following the same trend of the LIC corrections (see Fig. S2a). The Hertz model, which is exact for the paraboloidal indenter, does not accurately reproduce the indentation by a spherical indenter when the condition δ ≪ R is not met. The Hertz + LIC models recover the true Young's modulus within less than 1% (by default) for all indentation ranges up to δ = R (the correction proposed by Müller provides slightly more accurate results).

4.1.3 The Hertz model and the simultaneous correction of bottom and large indentation effects. The effect of the spatial (vertical) confinement of the sample was explored by fitting Hertz, Hertz + BEC and Hertz + BEC + LIC models to the FCs simulated on the TF-S systems (see Table 1). Different values of the maximum indentation δ, sample thickness h and tip radius R were used to cover a wide range of χ values, up to χ = 0.4. We have extended the simulation work of Garcia et al.¹⁹ by considering different routes in the parameter space leading to the same maximum χ: variable indentation δ at fixed R, h; variable medium thickness h at fixed δ, R; variable tip radius R at fixed δ, h.

Fig. 3 shows the relative discrepancies of the obtained Young's modulus E from the nominal one used in the simulations of the TF-S systems. Fig. 3a, b and c show the results of the application of BEC and LIC separately, whereas Fig. 3d, e and f show the results of the combination of LIC and BEC.

Fig. 3 Relative discrepancies of the obtained Young's modulus E compared to the nominal one (E_nom = 0.5 MPa) for the TF-S systems for different maximum values of the parameter In the first column of plots (a–c), LIC and BEC were applied separately. In the second column (d–f), the combined correction LIC + BEC was applied (notice the dramatic reduction of the ΔE/E_nom range). (a and d) The maximum indentation was varied while fixing the medium thickness and the tip radius (third row of Table 1). (b and e) The medium thickness was varied while fixing the indentation and the tip radius (fourth row of Table 1). (c and f) The tip radius was varied while fixing the indentation and the medium thickness (fifth row of Table 1).

As long as one considers the effect of large indentation and finite thickness separately (Fig. 3a–c), the same trend is observed: as χ increases, the relative discrepancy increases, in absolute terms. Fitting the Hertz model without considering any correction leads up to an 80% overestimation of the Young's modulus for χ = 0.4 (red circles). Noticeably, the bottom effect is mostly responsible of this large discrepancy, the large indentations effect contributing for maximum 10% as δ/R approaches 1, confirming the results obtained on bulk systems (Fig. 2b and 3a, red vs. black markers). The bottom effect dominates the inaccuracy of the Hertz model already at moderate values of χ; a discrepancy of 5–8% is observed already at χ = 0.05. We also notice (Fig. 3a, red vs. grey markers) that the large indentation effect is opposite to the bottom effect (Fig. 3a, red vs. black markers), which leads to a marked overestimation of the Young's modulus. As long as a maximum 10% underestimation of the Young's modulus is acceptable, the application of BEC only represents a valid solution.

When both BEC and LIC are applied, the maximum relative discrepancy of the obtained Young's modulus from the nominal value can be further reduced to ±3% when Garcia's BEC is used (Fig. 3d–f), in the interval 0 < χ < 0.4. When coupled to a BEC, LIC_M and LIC_K perform similarly; however, when δ/R is small (down to δ/R = 0.01), LIC_M associated to BEC_G leads to a smaller Young's modulus discrepancy (Fig. 3d and f); indeed, in this small δ/R regime, LIC_M is more accurate (see Fig. S2a).

For similar maximum χ values, BEC_D may perform worse than BEC_G, in combination with a LIC (Fig. 3e and f). This result may seem at odd with the fact that BECs depend only on χ, therefore we could expect, for the same maximum χ, a similar difference between the accuracy of the two BECs (i.e., a similar relative discrepancy of the Young's modulus); instead, Fig. 3d–f show that, while BEC_G overall performance is rather independent on which parameter is varied, BEC_D is not. An explanation of the observed behaviour is that BEC_G uses a more general expression for the contact radius, which takes explicitly into account the vertical confinement (Garcia's contact radius a_G can be written as a_G = a_H(1 + p₁χ + …), where a_H is the Hertzian contact radius and p₁ is a constant, see eqn (S23) in ref. 19 and Note S2), while BEC_D assumes a purely Hertzian contact radius. It follows that Garcia's BEC, through the χ-dependent contact radius, can account for the different combinations of R, h, δ parameters, therefore providing a similar performance irrespective to how a specific maximum χ value is obtained; on the contrary, BEC_D, based on the Hertz radius, which is accurate only in the limit of the semi-infinite elastic layer, is more sensitive to specific combination of those parameters. Considering that , we can expect a stronger residual dependence of the relative Young's modulus discrepancy on the two R/h and δ/h ratios for BEC_D. Fig. 3d–f show indeed that the difference between BEC_D and BEC_G is stronger when the R/h is larger, while the δ/h ratio is far less impacting.

Fig. 3f shows an anomalous and systematic downward deviation when χ < 0.1 and δ/R > 0.2. This can be due neither to a poor performance of LIC, since we have validated its goodness up to δ = R (see Fig. 2b), nor to a poor performance of BEC, as Fig. 3d and e clearly show for χ < 0.1. As we previously discussed at the end of section 3.2.1, the combined approach of BEC + LIC in a multiplicative fashion basically neglects the coupling term linking BEC to LIC. Reasonably, being χ and γ the two relevant parameters for BEC and LIC respectively, the coupling term is expected to depend, at least to the leading order, on their product χγ; we can therefore expect a larger discrepancy as the value of this product increases. Indeed, when χ < 0.1 in Fig. 3d and e, the product χγ remains well below 0.01, whereas in Fig. 3f, χγ = 0.025 for χ = 0.1 and reaches χγ = 0.05 when χ = 0.05. To further elucidate this aspect, we have performed additional simulations at low χ regimes and moderate χγ values, see Note S8 and Fig. S11.

Our simulations confirm that due to their multiplicative nature, LIC and BEC, the large indentation and bottom effect correction functions, can be applied together to the FCs obtained using a spherical indenter to obtain sets of FCs that can be fitted using the standard Hertz model. Depending on the geometry of the system, the most accurate correction, according to our simulations, is generally obtained combining BEC_G and LIC_M, though coupling with LIC_K provides only negligible worsening of the result in the limit of very small δ/R ratios. It is worth noting that, despite BEC only depends on the parameters , the accuracy of the correction depends also on the way a specific maximum χ is obtained (i.e., on the specific combination of R, h, δ).

4.2 Experimental validation of contact mechanics models and corrections

AFM nanoindentation experiments were performed on polyacrylamide (PA) gels, both bulk and exhibiting a gradient of thickness (from the thin film to the bulk regime), produced according to the protocol described in Note S7 and shown in Fig. S7. For the analysis of the experimental force curves, we selected LIC_M and BEC_G, respectively, since this combination performed better according to the FEA study presented in the previous sections.

4.2.1 Hertz and Sneddon models in the large indentation regime. Fig. 4 shows a representative set of FCs acquired using a colloidal probe with R = 8.9 μm in a region of a PA gel with thickness of approximately 400 μm, which we can consider as bulk (χ = 0.018). Fig. 4a shows the force vs. indentation curves, Fig. 4b shows the linearised FCs according to eqn (16); a deviation from linearity of the non-corrected rescaled FCs is observed in Fig. 4b for δ > 5 μm, corresponding to δ/R ≃ 0.56. Fig. 4c shows the same FCs linearised and corrected for the large indentations effect, according to eqn (18). The linearity of the corrected FCs is remarkable, especially considering that this trend is conserved up to a maximum indentation of approximately 9 μm (δ/R ≃ 1).

Fig. 4 The force vs. indentation curves acquired on a 700 μm thick PA gel, (a) original FCs; (b) FCs rescaled according to eqn (16); (c) FCs rescaled and corrected according to eqn (18) (LIC). In (d) the effective stress σ* vs. strain ε* curve (Tabor; see Note S3) after the application of LIC is shown.

That the linearity of the rescaled, corrected FCs shown in Fig. 4c is a sign of a linear elastic behaviour of the PA sample, is further suggested by the observed scaling of the effective stress vs. effective strain curves calculated using the Tabor's relations:³⁰ ε* = 0.2a/R, σ* = F/(πa²), after the application of LIC (Fig. 4d). The stress–strain scaling for the PA film is linear up to a strain of 20%; the maximum stress is approximately 70 kPa, about half the Young's modulus of the PA gel, as measured in the bulk region. The result reported here is due to a combination of factors: the excellent linear elastic behaviour of PA gels and the use of large colloidal probes, which are effective in keeping stresses and strains low, distributing them across a large volume in the sample. Incidentally, the large deformation field is what makes the bottom effect important in nanomechanical measurements performed with colloidal probes. Sharp tips are far less affected by the bottom effect, though they can easily stimulate the nonlinear elastic response of the material.

Fig. 5a shows the results of the experiments performed on two bulk samples made from the same working solution and mechanically tested using two different probes: probe [R₁, k₁] and probe [R₂, k₂] (described in section 3.3), yet achieving the same maximum indentation (δ ≃ 4.2 μm), corresponding to maximum δ/R values of 0.83 and 0.58, respectively. The Young's moduli of each bulk sample obtained using the different contact mechanics models (Sneddon, linearised Sneddon, Hertz and Hertz + LIC_M) are shown. These experiments confirmed that the application of the Hertz model without LIC leads to the underestimation of the obtained Young's modulus in a δ/R-dependent manner (see red violin plots in Fig. 5a).

Fig. 5 (a) Obtained Young's modulus on two different bulk samples using two probes with different radii (the maximum indentation was the same, corresponding to maximum δ/R values of 0.83 and 0.58). The experimental FVs have been analysed using the five different contact mechanical models described in the Methods and reported in the label of the horizontal axis of the graph; each single point in a violin plot corresponds to the median YM value of a single FV, analysed with the described models. (b) Relative discrepancy between the Young's modulus obtained using the different models and the value obtained using the Sneddon model. The experimental results are shown as violin plots (single Young's modulus median values have been omitted for sake of clarity) together with the results of the simulations (single dots), described in the previous sections. The results obtained by applying the simple Hertz model without corrections are shown in red, while the results obtained applying the linearised Sneddon model are shown in purple. The grey violin plots represent the results obtained using the Hertz model corrected with the Müller polynomial (LIC_M); the dashed violin-plots represent the results obtained by fitting the contact mechanics models to a smaller portion of the indentation axis.

The agreement between Sneddon and linearised Sneddon models (Young's modulus values agrees with less 1%) demonstrates that the assumption of negligible adhesion in experiments was reasonable; indeed, the proposed linearisation of the implicit Sneddon system works only under the hypothesis that any adhesion offset can be neglected. Using a linearised form of the Sneddon model represents a great advantage in data analysis. However, given that the results obtained by Hertz + LIC deviate from those of the Sneddon model by only 0.6%, the corrected Hertz model represents an even better tool, due to its simplicity, for the fit of the FCs, when spherical tips are used.

Plotting the relative discrepancy ΔE/E_bulk = (E − E_bulk)/E_bulk highlights the large indentation artifacts and the efficacy of the correction procedure. We normalized the data with respect to the Young's modulus E_bulk obtained through the Sneddon model and showed in Fig. 5b. Given that FEA showed that on bulk samples this model is very accurate, implicitly including the large indentation effects; it represents the experimental equivalent of the nominal Young's modulus in the simulations; results from both experiments and simulations are shown together (in the case of simulation, E_bulk corresponds to the nominal E_nom). Fig. 5b confirms that the large indentation effect causes an underestimation of the true Young's modulus. When a spherical tip is used, neglecting the correction of the large indentation effect, i.e., applying the simple Hertz model, causes, in the worst case (δ/R = 1 ) an underestimation of the true Young's modulus value as large as 10%. Furthermore, if a smaller portion of the same experimental indentation curves is fitted with the Hertz model (the dashed violin-plots in Fig. 5b) the underestimation of the Young's modulus decreases in agreement with the predictions of finite elements analysis. The linearised Sneddon model for the spherical indenter and the Hertz + LIC_M provides the same level of accuracy of the standard Sneddon model for all δ/R ratios up to unity.

These results show that using spherical tips it is possible to exploit the whole indentation range (at least up to δ = R) for the fit, regardless of the tip radius, provided the correct model is used: Sneddon (also linearised) or Hertz + LIC. The use of parabolic tips in combination with the Hertz model would provide, in principle, no error at all, since for this geometry the Hertz model is exact. In practice, however, parabolic tips are not present on the market. Sharp pyramidal tips with blunted apexes⁵² mimic the parabolic profile for small indentations; however, sharp tips would likely enhance nonlinear elastic effects already at relatively small δ/R values, which is supported by the evidence that experimentally obtained Young's modulus values are typically largely overestimated.¹⁵ The availability on the market of nearly paraboloidal tips with larger radii of curvature (>100 nm) would allow in principle to apply the simple Hertz model with no need of LIC across a broad range of δ/R values.

4.2.2 The Hertz model and the simultaneous correction of bottom and large indentation effects. Given the negligible differences between the results obtained using the Sneddon model and the Hertz + LIC model and considered the computational advantages of the linearised corrected Hertz models (as discussed in section 3.2.1), we have used the Hertz + LIC_M model to further investigate the BEC correction. Indeed, to our knowledge, models implementing BEC on the original implicit Sneddon system are not available.

To experimentally test the efficacy of combined BEC and LIC, we acquired several FVs on samples possessing a thickness gradient probe [R₃, k₃] and a maximum indentation of δ = [5–8] μm, corresponding to δ/R = 0.56–0.89, then we applied the simple Hertz model and the two corrections coupled together, resulting in LIC_M + BEC_G.

Fig. 6 shows a representative set of FCs acquired in a thin region (approximately 47 μm) of a PA gel possessing a gradient of thickness. Fig. 6a shows the force vs. indentation curves, Fig. 6b shows FCs linearised according to eqn (16). No significant deviations from linearity are observed in the rescaled FCs acquired in the presence of a bottom effect, at odd with the case of the large indentation effect; the bottom effect mainly affects the overall slope of the rescaled force–indentation curve rather than distorting it. Fig. 6c shows the same FCs linearised and corrected for both large indentation and bottom effect, according to eqn (20). The linearity of the corrected FCs is very good, up to a maximum indentation of approximately 5 μm (δ/R ≃ 0.56, χ ≈ 0.15). It can be noticed that the slope of the corrected linearised FCs is smaller than that of the non-corrected FCs, which is consistent with the expected stiffening of the film caused by the bottom effect.

Fig. 6 The force vs. indentation curves acquired on a 47 μm thick PA gel. (a) Original FCs; (b) FCs rescaled according to eqn (16); (c) FCs rescaled and corrected according to eqn (20) (LIC_M + BEC_G). In (d) the effective stress σ* vs. strain ε* curve (Tabor; see Note S3) is shown.

Fig. 6d shows that also in the case of the thin PA gel the effective stress vs. effective strain curves, after the application of BEC_G and LIC_M, are linear up to the maximum induced strain of 15% and a maximum stress which is less than half the Young's modulus of the film; this confirms that also the PA sample with a gradient of thickness behaves like a linear elastic material.

Fig. 7a shows the Young's modulus vs. thickness experimental curves, obtained using the Hertz and Hertz + LIC_M + BEC_G models. First, it can be noticed that the Young's modulus obtained using the standard Hertz model has a marked dependence on the thickness (red symbols), for thickness up to 250 μm; for larger values of the thickness, the obtained Young's modulus value tends to converge. However, the bulk Young's modulus value obtained by the Hertz model is not accurate, since it is affected by the large indentations effect (at such large thicknesses, the bottom effect is negligible). The reference bulk Young's modulus value, E_bulk (green line), has been calculated as the mean of the YM values obtained from those FVs where the differences between the values obtained using Hertz + LIC_M (not shown) and Hertz + LIC_M + BEC_G, respectively, is less than 1%; indeed, the two models are expected to provide the same correct bulk Young's modulus value in the bulk region. In agreement with the predictions of the LIC model and the results of FEA in the limit δ/R → 1, the difference between E_bulk and the Young's modulus obtained by the simple Hertz model, is approximately 10%.

Fig. 7 (a) The Young's modulus obtained on PA samples with a thickness gradient, as a function of the local thickness. The results obtained fitting the Hertz and Hertz + LIC_M + BEC_G models, respectively, to the experimental data are shown. The red dashed line represents the apparent bulk Young's modulus, as obtained by the Hertz model at large thicknesses (χ < 0.01). The green dashed line represents the corrected bulk Young's modulus, E_bulk, as obtained by the Hertz + LIC_M + BEC_G model in the regions where the bottom effect is negligible (χ < 0.01). (b) The relative discrepancy between the obtained Young's modulus and its bulk value as a function of the parameter χ. The experimental data (violin plots) are shown along with the results of FEA simulations (symbols) performed with the same experimental geometry.

Fig. 7a shows that bottom effect is clearly impacting already at relatively small values of χ (χ > 0.03), which for the δ/R ratios used in our experiments (large indentations and tip radius) corresponds to a thickness of approximately 200–250 μm. It is common habit to consider the bottom effect important when the δ/h ratio exceeds 0.1; our results confirm that the parameter , rather than the δ/h ratio, should be considered to assess the opportunity of applying the bottom effect correction (especially when the tip radius is large). A partial mitigation of this caveat valid for large indentation experiments comes from the fact that in the intermediate 0.15 > χ > 0.03 region, bottom and large indentation effects are both present, with opposite effects.

Applying both LIC_M and BEC_G (Fig 7a, green symbols) allowed us to obtain an experimental Young's modulus vs. thickness curve that is flat across the widest range of thickness values, converging to the bulk Young's modulus value for large thicknesses. BEC_G very effectively removed the dependence of the Young's modulus on thickness in the 60–250 μm thickness range.

The relative discrepancy against χ allows to better characterise the scaling of the obtained Young's modulus and to compare experimental results with the predictions of the finite elements simulations. In Fig. 7b the observed experimental discrepancies before and after the application of the Hertz + LIC_M + BEC_G model are shown, along with the results of simulations performed with the same experimental geometry. When both BEC and LIC are applied, the bulk Young's modulus is recovered within 1% up to χ ≈ 0.15 above which the discrepancy increases, for reasons that will be discussed later. The results of the nanomechanical measurements with both large indentation and bottom effect corrections applied on the PA gel are in very good agreement with those of FEA, following the same trend (for χ < 0.15). This result, and the evidence of the linear elastic response (Fig. 4c and 6c), confirm the accuracy of the proposed combined correction of bottom and large indentation effects for the Hertz model when a spherical indenter is used.

For thickness below 60 µm, corresponding to χ ≳ 0.15, the Hertz model corrections are not able to recover the bulk Young's modulus. At present we do not have a conclusive explanation for this unexpected behaviour, yet we can consider several evidences and make some hypotheses. Nonlinear elastic effects must be considered first, as potentially responsible for the residual nonlinearity observed in the corrected Young's modulus vs. thickness curve. However, we tend to exclude that nonlinear elasticity plays a dominant role here because the Tabor's effective stress–strain curves provide clear evidence of linear elastic behaviour of our films, in both bulk and constrained configurations (Fig. 4d and 6d). This is consistent with the fact that PA is known as a highly linear elastic material. Moreover, the corrected FCs show excellent linearity across the entire range of indentation explored (Fig. 4c and 6c), as predicted by the Hertz model for a linear elastic material. We are prone to consider the unexpected stiffening of the PA gel as a real phenomenon; in fact, our presumed reference elastic PA samples are likely to possess a gradient of elasticity, and in particular are stiffer in the thinner regions. This is further supported by the fact that the FEA simulations shown in Fig. 7b, performed with the same experimental geometry (probe radius, sample's thickness and maximum indentation) of the indentation experiments, provide the nominal Young's modulus within 3% for all the range of χ explored.

Our failure in producing a reference elastic sample with uniform elastic properties, irrespective to the local thickness, could be due to the specific synthesis process that we used (see Note S7). The monomer and the crosslinker were mixed with a photoinitiator, and the solution were exposed to UV radiation for a relatively long time. Indeed, exposure time was increased to 30 minutes after we moved the UV lamp away from the sample to avoid overheating. It is not unlikely that a combination of nonuniform illumination and light absorbance across the sample thickness, consequent differential heating and diffusion of different species in the polymer solution may alter the crosslinking ratio and the material density in the thinner regions, increasing both, leading to a gradient of rigidity and a stiffening. Alternatively, chemical initiation of the crosslinking process could be exploited.^53–55 Thickness dependence of the Young's modulus of polymers has been observed and attributed to different organisation of polymeric chains during the synthesis process of thin films (for example as an effect of shear stress for spin-coated materials⁵⁶). As the thickness of the substrate-bound film decreases, a strongly nonuniform, thickness-dependent swelling of the spatially constrained hydrogel (all measurements have been carried out in MilliQ water) can lead to reorganisation of the polymeric chains and to an internal pre-stress state (increasing with decreasing thickness), balanced by the surface tension of the PA gel, which could impact on the results of the nanomechanical measurements.⁵⁷ Moreover, the presence of a continuous gradient of thickness leads to an asymmetric redistribution of stresses inside the sample as the tip approaches the thin end and weakens the assumption of an axisymmetric contact. However, due to the smallness of the tilt angle, we consider the error deriving from the non-axisymmetric loading of the probe negligible; in any case, this error is of systematic nature and common to all the measurements, and we expect that its impact is minimised by the normalisation adopted with respect to the bulk value of the Young's modulus. To mitigate the potential impact of nonuniform swelling and stress distribution inside the PA films, it would be desirable to produce samples with a discrete, rather than continuous, variation of thickness, designing the mould to obtain flat terraces of different height rather than a uniform slope across the sample.

We further investigated whether the anomalous stiffening observed experimentally could be associated to thickness-dependent poroelastic effects, related to the migration of water through the polymeric network as a result of the applied load. As discussed in detail in Note S9, we cannot rule out a contribution of poroelasticity, likely due to a reduced diffusivity of water in the thinner region, where, as we argue, both crosslinking ratio and material density increase.

Non-water-based, non-swelling elastic materials could represent a valid alternative to hydrogels, although potential issues related to higher Young's modulus values, viscoelasticity and elastic nonlinearity should be considered. Eventually, by increasing the radius of the tip and the maximum indentation, it would be possible to explore a similar range of γ and χ values while measuring thicker regions, where the above-mentioned effects are expected to be less important.

The residual nonlinearity observed in our experiments deserves further investigations, which in turn stresses the importance of developing suitable test samples for the validation of mechanical models. It would be highly desirable to produce a highly elastic sample across a broad range of Young's modulus values [0.5–1000] kPa, with controlled nonlinearity. The lack of reliable, uniform and well characterizable elastic reference samples for nanomechanical analysis is at present still an open issue.

In consideration of the importance of applying the proposed corrections to biologically relevant samples, like single cells and tissue slices, it must be considered that the relevant parameter for bottom effect is , rather than the absolute thickness of the sample. Even restricting to χ = 0.15 the validity of our approach as a consequence of the observed residual discrepancy, by the proper choice of the tip radius (between 2 and 10 μm) and indentation (between 0.5 and 3 μm), it would be possible to measure and correct reliably the Young's modulus of a sample down to thicknesses of a few microns with χ remaining below 0.15 (see Fig. S10). In fact, according to the literature, the range of application of the proposed corrections extends far beyond χ = 0.15, as demonstrated for example in ref. 15 and 28.

5. Conclusions

In this work, we have validated the application of Hertzian contact mechanics models and corrections in the framework of linear elasticity for the analysis of force vs. indentation curves acquired using spherical probes by means of finite elements simulations and nanoindentation measurements on polyacrylamide gels with controlled thickness. We have systematically investigated the impact of both large indentations (δ/R → 1) and vertical spatial confinement on the accuracy of the nanomechanical analysis performed with the Hertz model.

Our results demonstrate that, on bulk systems (thick samples, for which ), the standard Hertz model can be used to fit force vs. indentation curves acquired using colloidal probes up to indentation δ ≤ 0.2R with an error below 2%, while both standard and linearised versions of the Sneddon model and the Hertz + LIC model can be used to fit force curves up to δ = R with approximately 1% accuracy. The linearised Sneddon and the Hertz + LIC models turned out to be as accurate as the standard Sneddon model, which is computationally more demanding.

We also demonstrated that the bottom effect leads to significant (>5%) deviations from the bulk Young's modulus value for χ as small as 0.05, in the large indentation regime. We demonstrated that it is possible to combine BEC and LIC into a single correction function to apply the Hertz model to the FCs, greatly expanding its range of use and reliability, thus allowing to characterise samples possessing ample variations of thickness (height) on the scale of the tip dimensions. Such combined corrections are relevant when using micrometer-sized spherical tip in the large indentation regime (both δ/R and χ are large) and allow to take full advantage of the use of colloidal probes in the nanomechanical analysis.

The results of the nanoindentation measurements carried out on PA gels agree with the predictions of FEA and contribute to validate the proposed combined corrections of the Hertz model for force acquired using a spherical indenter. The agreement between FEA and experimental results is due to a combination of factors: the excellent linear elastic behaviour of PA gels and the use of large colloidal probes, which are effective in keeping stresses and strains low, distributing them across a large volume in the sample. Incidentally, the large deformation field produced by large spherical tips makes the bottom effect important in nanomechanical measurements performed with colloidal probes compared to sharp tips at similar indentations; however, sharp tips, besides possessing a not well characterizable geometry, can easily stimulate the nonlinear elastic response of the material. Test of the linearity of the elastic response should be always carried out, for example by looking at the scaling of the Tabor's effective stress vs. strain curves.

The data analysis framework we presented, which is based on the linearisation of the Hertz equation, allows to include and test several different corrections among the one already published and those that will be published in the future.

The knowledge of the local sample height, or thickness, is crucial to implement the bottom effect correction. AFM offers the unique opportunity of obtaining both the topographic and the elastic maps from the same set of force curves; in those cases where the mean thickness of the sample cannot be obtained from the FCs (for example, because the underlying substrate is not accessible), an independent measure of the sample height must be obtained (for example by fluorescence confocal microscopy,²⁸ profilometry, or exploiting the raw approach of the AFM head until engage, as described in this work).

The identification of reliable contact mechanics models and the standardisation of efficient and accurate data analysis methods are still open fields of investigation within the AFM and the nanoindentation community. Given the typical dimensions of the systems and indentations parameters, both large indentations and bottom effect have a strong impact on the nanomechanical characterisation of living cells as well as thin tissue and ECM slices and should be routinely implemented in the experimental and data analysis procedures. Overall, the accurate correction of both large indentations and bottom effect can reduce the variability between results coming from different laboratories and increase the reproducibility of the nanomechanical experiments, supporting the standardisation effort.

Moreover, our study demonstrated that an essential part of this standardisation effort must be the development of reliable, uniform and well characterizable elastic reference samples with negligible elastic nonlinearity and controlled Young's modulus in physiologically relevant range [0.5–1000] kPa.

Author contributions

Conceptualisation: MC, AP; methodology – probe fabrication and characterisation: HH, MC, LDF; methodology – finite element simulations: LDF; data curation and analysis: HH, MC, LDF; original draft writing and editing: MC, LDF; draft revision: all authors; supervision: MC, AP; resources, funding, and project administration: AP. Authors’ contributions were allocated adopting the terminology of CRediT Contributor Roles Taxonomy.

Conflicts of interest

There are no conflicts to declare.

Data availability

The original data for this article (the raw force–distance curves, saved in the proprietary file format of Bruker, as well as the simulated force–distance curves, saved in ascii format), are available upon reasonable request from the authors.

The Matlab scripts developed by the authors to analyse the data are also available, upon reasonable request.

Matlab (Mathworks, https://it.mathworks.com/products/matlab.html) is a commercial software. Versions>R2020 of Matlab have been used to analyse the data using our custom scripts.

Supplementary information (SI) is available. See DOI: https://doi.org/10.1039/d5nr05344g.

Acknowledgements

This research was partially funded by the European Union Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement no. 812772, project Phys2BioMed. We acknowledge financial support under the National Recovery and Resilience Plan (NRRP), Mission 4, Component 2, Investment 1.1, Call for tender no. 104 published on 2.2.2022 by the Italian Ministry of University and Research (MUR), funded by the European Union – NextGenerationEU– project title Impact of chromatin organization on nuclear stiffness and cell migration – CUP G53D23002540006 – Grant Assignment Decree No. 2022KRWA7Y_002 adopted on 15-06-2023 by the Italian Ministry of Ministry of University and Research (MUR). We thank Andrei Shvarts and Lukasz Kaczmarczyk from the University of Glasgow for support and discussions on Finite Element Simulations; we thank Ricardo Garcia from the Institute of Science and Materials of Madrid, Andreas Stylianou from the European University of Cyprus and Stylianos Vasileios Kontomaris from the Metropolitan College of Athens for useful discussions.

References

1. K. Takahashi, Y. Kakimoto, K. Toda and K. Naruse, J. Cell. Mol. Med., 2013, 17, 225.
2. M. L. Jackson, A. R. Bond and S. J. George, Cardiovasc. Drugs Ther., 2023, 37, 997.
3. S. W. Verbruggen and L. M. McNamara, in Mechanobiology in Health and Disease, Elsevier, 2018, pp. 157–214.
4. D. E. Ingber, Ann. Med., 2003, 35, 564.
5. Mechanics of Cells and Tissues in Diseases, ed. M. Lekka, D. Navajas, M. Radmacher and A. Podestà, De Gruyter, 2023.
6. M. Chighizola, T. Dini, C. Lenardi, P. Milani, A. Podestà and C. Schulte, Biophys. Rev., 2019, 11, 701.
7. H. Liu, L. A. MacQueen, J. F. Usprech, H. Maleki, K. L. Sider, M. G. Doyle, Y. Sun and C. A. Simmons, Biomaterials, 2018, 172, 30.
8. F. Iberite, M. Piazzoni, D. Guarnera, F. Iacoponi, S. Locarno, L. Vannozzi, G. Bolchi, F. Boselli, I. Gerges, C. Lenardi and L. Ricotti, ACS Appl. Bio Mater., 2023, 6, 2712.
9. Y. F. Dufrêne, T. Ando, R. Garcia, D. Alsteens, D. Martinez-Martin, A. Engel, C. Gerber and D. J. Müller, Nat. Nanotechnol., 2017, 12, 295.
10. D. J. Müller and Y. F. Dufrêne, in Nanoscience and Technology, Co-Published With Macmillan Publishers Ltd, UK, 2009, pp. 269–277.
11. H. Holuigue, E. Lorenc, M. Chighizola, C. Schulte, L. Varinelli, M. Deraco, M. Guaglio, M. Gariboldi and A. Podestà, Sensors, 2022, 22, 2197.
12. E. Lorenc, H. Holuigue, F. Rico and A. Podestà, in Mechanics of Cells and Tissues in Diseases, ed. M. Lekka, D. Navajas, M. Radmacher and A. Podestá, Berlin, Boston: De Gruyter, 2023, vol. 1,  DOI:10.1515/9783110640632-006.
13. L. Puricelli, M. Galluzzi, C. Schulte, A. Podestà and P. Milani, Rev. Sci. Instrum., 2015, 86, 33705.
14. M. Indrieri, A. Podestà, G. Bongiorno, D. Marchesi and P. Milani, Rev. Sci. Instrum., 2011, 82, 023708.
15. A. Kubiak, M. Chighizola, C. Schulte, N. Bryniarska, J. Wesolowska, M. Pudelek, M. Lasota, D. Ryszawy, A. Basta-Kaim, P. Laidler, A. Podestà and M. Lekka, Nanoscale, 2021, 13, 6212.
16. H. Hertz, Journal fur die Reine und Angewandte Mathematik, 1881, 92, 156–171.
17. L. Lacaria, A. Podestà, M. Radmacher and F. Rico, in Mechanics of Cells and Tissues in Diseases, ed. M. Lekka, D. Navajas, M. Radmacher and A. Podestá, Berlin, Boston: De Gruyter, 2023, vol. 1,  DOI:10.1515/9783110640632-003.
18. E. K. Dimitriadis, F. Horkay, J. Maresca, B. Kachar and R. S. Chadwick, Biophys. J., 2002, 82, 2798.
19. P. D. Garcia and R. Garcia, Biophys. J., 2018, 114, 2923.
20. N. Gavara and R. S. Chadwick, Nat. Nanotechnol., 2012, 7, 733.
21. P. D. Garcia, C. R. Guerrero and R. Garcia, Nanoscale, 2020, 12, 9133.
22. S. V. Kontomaris and A. Malamou, Eur. J. Phys., 2021, 42, 025010.
23. P. Müller, S. Abuhattum, S. Möllmert, E. Ulbricht, A. V. Taubenberger and J. Guck, BMC Bioinf., 2019, 20, 465.
24. Á. dos Santos, A. W. Cook, R. E. Gough, M. Schilling, N. A. Olszok, I. Brown, L. Wang, J. Aaron, M. L. Martin-Fernandez, F. Rehfeldt and C. P. Toseland, Nucleic Acids Res., 2021, 49, 340.
25. I. N. Sneddon, Proceedings of the Glasgow Mathematical Association, 1965, vol. 7, p. 48.
26. R. Long, M. S. Hall, M. Wu and C.-Y. Hui, Biophys. J., 2011, 101, 643.
27. H. Schillers, C. Rianna, J. Schäpe, T. Luque, H. Doschke, M. Wälte, J. J. Uriarte, N. Campillo, G. P. A. Michanetzis, J. Bobrowska, A. Dumitru, E. T. Herruzo, S. Bovio, P. Parot, M. Galluzzi, A. Podestà, L. Puricelli, S. Scheuring, Y. Missirlis, R. Garcia, M. Odorico, J.-M. Teulon, F. Lafont, M. Lekka, F. Rico, A. Rigato, J.-L. Pellequer, H. Oberleithner, D. Navajas and M. Radmacher, Sci. Rep., 2017, 7, 5117.
28. S. Pérez-Domínguez, S. G. Kulkarni, J. Pabijan, K. Gnanachandran, H. Holuigue, M. Eroles, E. Lorenc, M. Berardi, N. Antonovaite, M. L. Marini, J. Lopez Alonso, L. Redonto-Morata, V. Dupres, S. Janel, S. Acharya, J. Otero, D. Navajas, K. Bielawski, H. Schillers, F. Lafont, F. Rico, A. Podestà, M. Radmacher and M. Lekka, Nanoscale, 2023, 15, 16371.
29. C. Valero, B. Navarro, D. Navajas and J. M. García-Aznar, J. Mech. Behav. Biomed. Mater., 2016, 62, 222.
30. D. C. Lin, D. I. Shreiber, E. K. Dimitriadis and F. Horkay, Biomech. Model. Mechanobiol., 2009, 8, 345.
31. C. E. Wu, K. H. Lin and J. Y. Juang, Tribol. Int., 2016, 97, 71.
32. K. D. Costa and F. C. P. Yin, J. Biomech. Eng., 1999, 121, 462.
33. S. V. Kontomaris, A. Stylianou, A. Georgakopoulos and A. Malamou, Micron, 2023, 164, 103384.
34. S. V. Kontomaris and A. Malamou, Mater. Res. Express, 2020, 7, 033001.
35. S.-V. Kontomaris and A. Malamou, Eur. J. Phys., 2022, 43, 015010.
36. A. C. Fischer-Cripps, The Hertzian Contact Surface, 1999.
37. I. N. Sneddon, Int. J. Eng. Sci., 1965, 3, 47.
38. K. B. Park, M. S. Kim, J. H. Kim, S. K. Kim and J. M. Lee, J. Polym. Eng., 2019, 39, 432.
39. L. D. Landau, E. M. Lifshitz, A. M. Kosevich and L. P. Pitaevskii, Theory of Elasticity, 3rd edn, Butterworth-Heinemann (Elsevier), Oxford, 1986, vol. 7,  DOI:10.1016/C2009-0-25521-8.
40. S. R. Kalidindi and S. Pathak, Acta Mater., 2008, 56, 3523.
41. D. K. Patel and S. R. Kalidindi, Acta Mater., 2016, 112, 295.
42. B. R. Donohue, A. Ambrus and S. R. Kalidindi, Acta Mater., 2012, 60, 3943.
43. C.-E. Wu, K.-H. Lin and J.-Y. Juang, Tribol. Int., 2016, 97, 71.
44. P. Carl and H. Schillers, Pfluegers Arch., 2008, 457, 551.
45. H. J. Butt and M. Jaschke, Nanotechnology, 1995, 6, 1.
46. H. J. Butt, B. Cappella and M. Kappl, Surf. Sci. Rep., 2005, 59, 1.
47. M. Chighizola, J. Rodriguez-Ramos, F. Rico, M. Radmacher and A. Podestà, in Biomedical Methods, ed. M. Lekka, D. Navajas, M. Radmacher and A. Podestà, De Gruyter, 2023, pp. 105–128.
48. M. Chighizola, L. Puricelli, L. Bellon and A. Podestà, J. Mol. Recognit., 2021, 34, e2879.
49. J. Laurent, A. Steinberger and L. Bellon, Nanotechnology, 2013, 24, 225504.
50. J. Domke and M. Radmacher, Langmuir, 1998, 14, 3320.
51. E. Lorenc, L. Varinelli, M. Chighizola, S. Brich, F. Pisati, M. Guaglio, D. Baratti, M. Deraco, M. Gariboldi and A. Podestà, Sci. Rep., 2023, 13, 12175.
52. F. Rico, P. Roca-Cusachs, N. Gavara, R. Farré, M. Rotger and D. Navajas, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2005, 72, 1.
53. J. R. Tse and A. J. Engler, Curr. Protoc. Cell Biol., 2010, 47, 10.16.1–10.16.16.
54. A. K. Denisin and B. L. Pruitt, ACS Appl. Mater. Interfaces, 2016, 8, 21893.
55. S. Sheth, E. Jain, A. Karadaghy, S. Syed, H. Stevenson and S. P. Zustiak, Int. J. Polym. Sci., 2017, 2017, 5147482.
56. M. Liu, J. Sun, Y. Sun, C. Bock and Q. Chen, J. Micromech. Microeng., 2009, 19, 035028.
57. J. M. Long and G. F. Wang, Mech. Mater., 2013, 56, 65.

---

Footnotes

† Present affiliation: Mechanobiology Institute (MBI), National University of Singapore, 5A Engineering Drive 1, 117411, Singapore.

‡ Present affiliation: Institut national de la santé et de la recherche médicale (INSERM), Dynamo Group, 163 avenue de Luminy, Marseille, France.

---