The software functionalities were split into different self-contained libraries to help organize the code while making it more readable, maintainable, testable and reusable. Similarly to the AFM data analysis PyJibe developed by Paul Muller (Müller et al., 2019), the functionalities of PyFMLab are split on three main Python libraries organized in a hierarchical manner (Figure 1):
PyFMRheo
The source code for PyFMRheo is available (from zenodo López-Alonso, 2023a) (also available from GitHub)
As the main building block of PyFMLab, this library includes a wide variety of models and predefined routines for the nanomechanical characterization of biological samples, correction for hydrodynamic drag, calibration of the z-piezo for oscillatory measurements and non-contact calibration of AFM probes. As PyFMReader and PyFMRheo have been developed with the intention of being used together, the predefined routines available in PyFMRheo expect input FDCs and thermals as objects defined by PyFMReader. Nevertheless, this library can be used on its own, with the user being able to decide which library to use to load their data, which allows it to be used for the analysis of any data acquired on lab made or other commercial instruments. All models and tools required for processing AFM force curves and thermal tunes are fully accessible, allowing the users to compose their own custom data processing scripts. Examples on how to use each of the predefined routines are provided as interactive Jupyter notebooks. In addition, a tutorial on how to implement a contact model is provided in the dedicated GitHub site.
Elastic models.
Contact models describing the force (F) versus indentation (δ) relationship allowing to obtain the Young’s modulus (E) from the approach segment of FDC for several AFM tip geometries have been included (see chapter 2.2 Contact Mechanics Lacaria et al. in Lekka et al., 2023a). By default, the fitting routine considers the point of contact and force offset as fitting parameters.
Additionally, the library includes a model considering the finite thickness of the sample (Table 2) developed by Garcia and Garcia (Garcia & Garcia, 2018), where the thickness of the sample is considered when determining the Young’s modulus.
Due to the complexity of automatically determining the thickness of the sample at each pixel, in particular for force maps, the models incorporating Garcia and Garcia bottom effect correction are implemented in the PyFMRheo analysis library but are not accessible through the graphical user interface (GUI).
Default routine implementation
1. The raw data of the curve is pre-processed to obtain the tip position (δ) by subtracting the cantilever vertical deflection (d) to the piezo position (z).
2. An estimate of the point of contact (PoC) is computed using the ratio of variances (RoV) method based on Gavara, 2016. This is computed in two small windows before and after each point of the deflection signal (d):
Where RoVi is the ratio of variances at each data point i, var is the variance, d is the deflection signal and N is the window size inputted by the user.
3. The spring constant and initial estimation of the PoC are used to obtain the force vs indentation curve.
4. Using the parameters defined by the user, the corresponding model is fitted to the force vs indentation curve to obtain the refined PoC and the Young’s modulus.
5. The default fitted parameters are the Young’s modulus (E), the PoC and the force offset.
Viscoelastic models. Several models have been proposed to determine the viscoelastic force-indentation relationship of a sample. These models are numerical (Efremov et al., 2017) or analytical (Brückner et al., 2017; Garcia et al., 2020; Eroles et al., 2023) solutions and are based on the seminal works by Lee and Radock (Lee & Radok, 1960), Graham (Graham, 1967) and Ting (Ting, 1966; Ting, 1968). The selection and implementation of the models in PyFMLab were performed focusing on the determination of the mechanical properties of living mammalian cells and extracellular matrices. Such samples have been shown to exhibit a power law response both at low and high frequencies (Balland et al., 2006; Fabry et al., 2001; Jorba et al., 2017; Jorba et al., 2019; Rigato et al., 2017). PyFMRheo includes both numerical (Table 3) and analytical (Table 4) models for a power law response material.
In the models, F is the force acting on the cantilever tip, δ is the indentation which is a function of the time, t. Following (Brückner et al.), the time t1 is obtained from the auxiliary function determined by:
where ξ is the dummy time variable required for the integration, E(t) is the Young’s or elastic relaxation modulus described by a power-law rheology (PLR) model:
where E0 is the scaling factor or instantaneous elastic modulus, t0 an arbitrary time scale assumed to be 1s, and β is the power-law exponent or fluidity (β = 0 for material with solid-like behaviour, β = 1 for material with fluid-like behaviour).
Assuming that the indentation is proportional to the different approach and retract velocities, t1 is found to be:
where va is the constant approach velocity, vr is the constant retract velocity (Brückner et al., 2017).
C~ is the geometry dependent coefficient. The definition of this parameter based on the tip geometry is included in Table 5. In addition to the models listed in Table 1–Table 5, other models have been implemented in pyFMRheo, including a simple linear model to estimate the sample stiffness, the Derjaguin-Muller-Toporov (DMT) model and the Kontomaris approximation of a spherical indenter (Kontomaris & Malamou, 2021).
Table 1. Implemented contact elastic relationships between force (F) and indentation (δ).
R corresponds to the tip radius, v to the Poisson’s ratio, E to the Young’s modulus, θ to the semi-opening angle of the cone or the pyramid face, FHertz (z–d0) is the force obtained from applying the Hertz model and FAdhesion (d0) is the adhesion force according to DMT theory (Derjaguin et al., 1975).
| Tip geometry | Formula | Reference |
|---|
| Paraboloidal | | (Hertz, 1882) |
| Conical | | (Love, 1939; Sneddon, 1965) |
4-sided regular
pyramid | | (Barber & Billings, 1990;
Bilodeau, 1992) |
| Paraboloidal | | (Derjaguin et al., 1975) |
Table 2. Implemented bottom effect corrected elastic relationships between force (F) and indentation (δ).
R corresponds to the tip radius, v to the Poisson’s ratio, E to the Young’s modulus, h to the sample thickness and θ to the semi-opening angle of the cone or the pyramid face.
| Tip geometry | Formula | Reference |
|---|
| Paraboloidal | | (Garcia & Garcia, 2018) |
| Conical | |
4-sided regular
pyramid | |
| Paraboloidal | | (Kontomaris & Malamou, 2021) |
Table 3. Implemented numerical viscoelastic relationships between force (F) and indentation (δ).
R corresponds to the tip radius, v to the Poisson’s ratio, E to the elastic relaxation modulus, t1 to the integration auxiliary function and θ to the semi-opening angle of the cone or the pyramid face.
| Tip geometry | Force curve
segment | Equation | Reference |
|---|
| Paraboloidal | Approach |
| (Efremov et al., 2017) |
| Retract | |
| Conical | Approach | |
| Retract | |
4-sided regular
pyramid | Approach | |
| Retract | |
Table 4. Implemented viscoelastic analytical relationships between force (F) and time (t) assuming power law rheology (E(t) = E0(t/t0)–β).
Γ corresponds to the gamma function, R to the tip radius, θ to the semi-opening angle of the cone or the pyramid face, v to the Poisson’s ratio, E0 to the scaling factor or instantaneous elastic modulus, β to the power law exponent or fluidity, t0 to the scaling time, t1 to the integration auxiliary function, va to the approach velocity and C~ to the tip geometry coefficients are defined in Table 5.
| Tip geometry | Force curve
segment | Formula | Reference |
|---|
| Paraboloidal | Approach | | Based on
(Brückner et al., 2017) |
| Retract | |
Conical / 4-sided
regular pyramid | Approach | |
| Retract | |
Table 5. Geometry coefficients used in the implemented viscoelastic analytical models (Table 4).
R corresponds to the tip radius and θ to the semi-opening angle of the cone or the pyramid face.
| Tip geometry | Coefficient |
|---|
| Paraboloid | |
| Cone | |
4-sided regular
pyramid | |
Default viscoelastic routine implementation
1. The raw data of the curve is pre-processed to obtain the tip position by subtracting the cantilever vertical deflection to the piezo position.
2. An estimate of the point of contact (PoC) is computed using the ratio of variances (RoV) method based on Gavara, 2016.
3. The spring constant and initial PoC estimation are used to obtain the force vs. indentation curve.
4. The elastic model corresponding to the AFM tip geometry is fitted to the approach segment of the force vs. indentation curve to obtain an estimation of the Young’s modulus and a better estimation of the PoC.
5. The corresponding viscoelastic model is then fitted to the full approach-retract cycle of the force vs indentation curve to obtain the power-law fluidity exponent (β), the instantaneous elastic modulus (E0) and the PoC. The viscous drag force (force felt by the cantilever due to the viscosity of the medium) can be corrected, from the correction factor specified by the user and the respective approach and retract velocities.
Determination of the complex shear modulus from active microrheology measurements. Dynamic mechanical analysis (DMA), or microrheology measurements, consist of applying deformations in a cyclical manner to a material in the following manner:
1. The cantilever approaches and indents the sample.
2. The cantilever oscillates at a single frequency. This step is repeated in order to measure all the desired frequencies (see AFM measurements section below for details).
3. The cantilever is retracted from the sample.
These cyclic forces result in variable deformations. As the stress (force applied to an area of the material) or strain (material deformation) applied is oscillatory (sinusoidal), the response of the sample is also sinusoidal. If a material is purely elastic, the strain and stress will be perfectly synchronous, at all times, meaning that, when the force is applied, the material will deform, and when the force stops being applied, the sample will recover its original shape. If the material is purely viscous, an applied stress will induce a strain response that is both delayed in time and not recovered once the application of stress ceases. Generally, a material will exhibit both elastic and viscous behaviors, resulting in a strain response that will lag (phase angle 0 – 90°) behind the applied stress. From DMA, the following frequency dependent parameters can be obtained:
• Complex modulus (G*): Overall resistance of the material to deformation including a real part (elastic, G′) and an imaginary part (viscous, G″).
• Storage modulus (G′): Elasticity of the materials and its ability to store energy.
• Loss modulus (G): Viscosity of the material and its ability to dissipate energy.
• Loss tangent (η): Degree of solid- or liquid-like mechanical behavior of the material.
Table 6 shows the values for G′, G″ and η for purely elastic and purely viscous materials. Mammalian cells and extracellular matrices (ECM) are viscoelastic materials, with both an elastic and viscous component. In 2003 Alcaraz et al. defined a methodology for the microrheological characterization of biological samples with AFM. When a sample is cyclically indented in AFM, a certain stress is applied to the sample, causing it to deform. Both the stress and the deformation of the sample can be measured through force-indentation curves, where the force can be used to monitor the stress applied, and the indentation signal can be used to monitor the corresponding deformation of the sample.
Table 6. Values of the storage modulus (G′), loss modulus (G″) and loss tangent (η) for purely elastic and viscous materials.
E is the Young’s modulus, v is the Poisson’s ratio and ω is the angular frequency.
Purely elastic
material | Purely viscous
material |
|---|
| G′ = 0 |
| G″(ω) = 0 | G″(ω) = μ |
| η(ω) = 0 | η(ω)→∞ |
Models for two AFM tip geometries have been implemented and are shown in Table 7.
Table 7. Implemented frequency domain rheological models for the complex shear modulus (G*) as a function of frequency (ω).
Where F(ω) and δ(ω) are the Fourier transform of the force and indentation, respectively, δ0 corresponds to the working indentation at which small amplitude oscillations are applied, R corresponds to the tip radius and θ to the semi-opening angle of the cone or the pyramid face.
| Geometry | Function | Reference |
|---|
| Paraboloidal | | (Mahaffy et al., 2000;
Rico et al., 2005) |
| Pyramidal | | (Alcaraz et al., 2003) |
The default routine expects data which has been acquired by indenting the sample and applying single frequency oscillations at consecutive segments, as shown in Figure 5.
Default rheology routine implementation
1. The raw data of the curve is pre-processed to obtain the tip position by subtracting the cantilever vertical deflection to the piezo position.
2. The force curve is pre-processed to obtain the piezo position vs tip position.
3. An estimate of the point of contact (PoC) is computed from the approach segment using the ratio of variances (RoV) method based on (Gavara, 2016). This allows determining the value of the working indentation (δ0).
4. If the user provides piezo calibration data and selects the option, the amplitude and phase of the z-piezo signal are corrected.
5. Both signals are then detrended by subtracting the rolling average of the signal with a window size equal to the number of points of each period.
6. The transfer function (Hd) is computed considering the ratio of the Fourier transforms of the small amplitude indentation (δ(ω)) as the input signal and the measured force (F(ω)) as the output signal.
7. If the user provides a value for the drag force at contact, the transfer function values are corrected accordingly.
8. G*(ω) is computed using the corresponding frequency domain rheological contact model. For a paraboloid geometry: . From the values of G*(ω), η(ω) is computed.
Z-Piezo characterization. When characterizing the microrheological properties of cells or ECM, any delay between the force and indentation signals inherent to the instrument will result in an over- or underestimation of the viscoelastic properties of the sample. By acquiring oscillatory measurements in a non-deformable substrate, like glass, the inherent phase angle (ϕPZT) between the deflection and z movement signals of the instrument can be measured and used to correct the measurements acquired on the sample.
In addition to the previous correction, in this routine, the ratio between the vertical deflection and the z-piezo signals is computed. This value can then later be used to correct the z-piezo signal acquired on the sample.
Default z-piezo characterization routine implementation
1. Both the z-piezo displacement (z) and the vertical deflection of the cantilever are detrended (d) by subtracting the rolling average of the signal with a window size equal to the number of points of each period.
2. The transfer function is computed considering zPZT as the input signal and dv as the output signal.
3. From the transfer function, the phase between z and d signals at each applied frequency are obtained.
4. Finally, the signal amplitude quotient between z and d signals is computed.
Determination of viscous drag coefficient. To minimize the noise in the measurements and to provide a suitable environment when working with living cells, most AFM microrheological measurements are performed in liquid. Therefore, when the cantilever oscillates, the surrounding fluid applies force on the cantilever by adding resistance to the motion (hydrodynamic drag). In the same manner as the piezo lag, the hydrodynamic drag results in an overestimation of the viscous properties of the measured sample. To account for this effect, Alcaraz and collaborators (Alcaraz et al., 2002) proposed a correction applicable to AFM measurements on soft samples in liquid at low Reynold numbers (Re < 1). They approached this problem by modelling the AFM cantilever as a sphere close to a rigid wall. The characterization of the hydrodynamic on a spherical body close to a rigid wall has been addressed by multiple authors (Brenner, 1961; Cox & Brenner, 1967). Based on these works, Alcaraz and collaborators estimate the drag force b(h) when approaching the sample to be:
Where η is the viscosity of the fluid, aeff is the effective radius of the cantilever, heff is the effective height of the tip and h is the distance between the tip and the sample.
Using this scaled spherical model, the drag force at contact (b(0)) can be estimated by measuring the drag force b(h) at different cantilever-sample separations and extrapolating the data to h = 0.
To compute b(h) the default routine expects data acquired with a specific protocol:
1. The AFM tip approaches and indents the sample (cell, tissue, gel, etc.).
2. The cantilever retracts a certain distance defined by the user.
3. The cantilever is oscillated at a single frequency defined by the user.
4. Steps 2 and 3 are repeated until the maximum desired distance from the sample has been reached.
An example of data collected using this protocol is shown on Figure 6.
Default viscous drag routine implementation
1. If the user provides piezo calibration data and selects the option to correct the amplitude of the z-piezo, the signal is corrected.
2. The transfer function (Hd) is computed considering the indentation (δ) as the input signal and the measured force (F) as the output signal.
3. If the user provides piezo calibration data, the transfer function values are corrected using the piezo amplitude and lag.
4. b(h) is then computed using the following expression:
Non-contact cantilever calibration. Many indirect spring constant calibration methods have been developed, with the Sader method being the most widely used, thanks to its ease of use and applicability to different cantilever geometries (rectangular, V-shaped) (Sader et al., 1999; Sader et al., 2005; Sader et al., 2012). This method relies on measuring the quality factor (Q) and resonant frequency (ωR) from the thermal spectrum of the cantilever to estimate the spring constant (k) by relying on hydrodynamic theory. When the cantilever oscillates, the surrounding fluid (air, water, etc.) applies force on the cantilever in two ways: by adding resistance to the motion (hydrodynamic drag) and by adding mass to the cantilever. By relating the measured ωR and Q factor of the cantilever oscillating in a fluid to the values of the cantilever oscillating in vacuum, an analytical model to determine the spring constant of the cantilever for the first eigenmode (k1) can be derived, by taking into consideration the geometry of the cantilever together with the density and viscosity of the fluid (Sader et al., 2005). The spring constant of a rectangular cantilever can be determined applying the following formula:
Where ρ is the density of the medium, b is the thickness of the cantilever, L is the length of the cantilever, Γi are the imaginary components of the hydrodynamic function, ωR is the resonance frequency of the cantilever and Q is the quality factor of the cantilever. Although this expression is limited to rectangular cantilevers, the Sader method was later extended to cantilevers of arbitrary shape (Sader et al., 2012)
With the goal of allowing to calibrate the spring constant of cantilevers with a reference standard, Sader et al. developed a virtual instrument (that allows the user to obtain the spring constant of a cantilever regardless of its shape by measuring ωR and Q. This method is referred to as the global calibration initiative (GCI) (Sader et al., 2012; Sader et al., 2016).
In 2006, Higgins et al. proposed a non-contact method to obtain the optical lever sensitivity (or deflection sensitivity) from the thermal spectrum of the cantilever, assuming that the spring constant remains unchanged in air and liquid. This approach requires a previous calibration of the spring constant in air with the Sader or other method (Higgins et al., 2006). As this method yields a deflection sensitivity different to the one obtained by contact-based methods, the result is then multiplied by a correction factor to obtain the correct value.
Although the Sader and GCI methods enhance the precision of nanomechanical measurements when compared to conventional contact-based approaches (Sader et al., 2012; Sader et al., 2016), they both assume cantilevers with a high Q value, which is not always the case for most cantilevers in liquid. For this reason, in 2020, Sumbul and collaborators (Sumbul et al., 2020) assessed the accuracy and precision of these methods. In this study, they observed that using the single harmonic oscillator (SHO) model to fit the thermal spectra together with the GCI method was less prone to systematic uncertainties and provided higher accuracy in the determination of k and the deflection sensitivity (Sumbul et al., 2020). The implementation of the calibration routine in this package has been based on these works.
Default spring constant calibration routine implementation
1. The SHO model is fit to the first eigenmode peak of the thermal data to obtain the white noise of the signal (Awhite, the amplitude of the peak (B), the resonance frequency (fR and Q factor of the cantilever.
2. From these values the spring constant can be computed using the general Sader method. This step is only valid for rectangular and V-shaped cantilevers.
3. If the user provides thermal data acquired in air, their username and password for the GCI web application, the software will get the GCI spring constant through the application programming interface of the GCI web application.
4. The deflection sensitivity (or inverse of the optical lever sensitivity, invOLS) is computed using the method proposed by Higgins et al. and Sumbul et al. The user can specify the factors to correct for the static versus dynamic deflection (sometimes referred as χ) and the fundamental oscillation mode (sometimes referred as β) of the cantilever (Chighizola et al., 2023). If no correction factor is specified, the default values () of 0.90 (rectangular) and 0.87 (V-shaped) are used (Pirzer & Hugel, 2009; Stark et al., 2001; Sumbul et al., 2020).
