A resonant structure designed for probing the elastic properties of suspension and adherent cells in liquid environments

Structures exploiting resonance phenomena can outperform equivalent structures operated in a static mode [19]. Despite this potential advantage, resonant cantilevers are rarely used for determining Young's modulus of living cells. Indeed, a major issue is that dynamic performances of AFM cantilevers are dramatically deteriorated as soon as they are surrounded by a viscous medium. Thereby, the quality factor of typical resonant AFM cantilevers immersed in liquids rarely exceeds 10–30 [20, 21].

To extract Young's modulus of a single living cell while taking full benefit of a dynamic mode, even in the presence of a liquid environment, we propose the suspended structure shown in figure 1(a). It consists of a planar structure clamped at its extremities. It also incorporates two rectangular apertures. Finite Element Analysis (FEA) conducted with COMSOL V.4 shows that if the structure is mechanically excited at the proper frequency, it can provide an antisymmetrical vibration mode where the two outer beams (OBs) oscillate in antiphase. In the mean time, the third central beam (CB) remains immovable, as depicted in figure 1(b).

Zoom In Zoom Out Reset image size

Measuring the force applied to a cell is a prerequisite to extract its Young's modulus. Conceptually, if a static force normal to the structure's plane is applied to the half-span of the CB, the whole structure deflects. Therefore, the OBs are forced to deflect as well. When the OBs oscillate as illustrated in figure 1(b), interactions between static and dynamic behaviors occur. Thereby, the static deflection imposed by the normal force affects their initial resonance frequency. In fact, even slight static pre-deflection of beam resonators may significantly impact their dynamics [22]. Accordingly, it is possible to recover the magnitude of the normal force applied upon the CB by monitoring frequency shifts of the OBs, as it is demonstrated in section 3.

The balanced mode of vibration shown in figure 1(b) is of particular interest for determining Young's modulus of individual cells. Indeed, it enables high sensitivity rates offered by resonators without exposing cells to vibrations. Thereby, the area located at the half-span of the CB can be adapted to place suspension or adherent cells in liquids (see figure 2). A main advantage offered by such a configuration is that, unlike conventional AFM cantilevers, the oscillating OBs do not need to be immersed in the liquid in order to measure a force applied to a cell. Therefore, major energy losses are avoided and high quality factor can be guaranteed. Because the OBs always oscillate in air, potential difficulties related to laser alignments in liquids are eliminated and measurements based on vision techniques are not required. It is also worth noting that since the CB is inherently force sensitive, low-cost and commercially available microindenters can be used to apply forces on cells during experiments. Moreover, artifact measurements due to capillary effects are minimized since the microindenter is not directly used as a force sensor.

Zoom In Zoom Out Reset image size

A static analysis including nonlinear effects is first conducted in order to detail how the OBs exactly behave during large deflections of the whole structure. Due to symmetry, and without loss of generality, our analysis can be limited to one-fourth of the structure, as sketched in figure 3.

Zoom In Zoom Out Reset image size

This quarter structure, similar in shape to a tuning fork, is composed of three segments: the overhang (black), the central tine (dark gray) and the outer tine (light gray). All segments are assumed to satisfy Bernoulli's beam theory. One-dimensional coordinate functions are used to approximate the displacement field of each segment. To guarantee enough degrees of freedom, the flexural displacement w_{1, 2, 3}(x) of each segment is modeled by third-order polynomial expressions:

$$\begin{equation} \left\lbrace \begin{array}{@{}c}w_1 (x) = a_0+a_1\, x+a_2\, x^2 + a_3\, x^3 \\ w_2(x) = a_4+a_5\, x+a_6\, x^2 + a_7\, x^3 \\ w_3(x) = a_8+a_9\, x+a_{10}\, x^2 + a_{11}\, x^3.\end{array}\right. \end{equation} \tag{ 1 }$$

The three polynomials in (1) are then used to calculate the total potential energy stored by the tuning fork during deflection:

$$\begin{equation} U_{tf} = U_b + U_s, \end{equation} \tag{ 2 }$$

where U_b is the sum of strain energies developed by each segment during bending:

$$\begin{equation} U_b = \sum _{i=1}^3 \, \frac{1}{2}\, E I_i \int _{c_i}^{d_i} \bigg(\frac{{\rm d}^{2}w_i}{{\rm d}x^2}\bigg)^2 \, {\rm d}x. \end{equation} \tag{ 3 }$$

In (3), c₁ = 0, d₁ = l/2, c_{2, 3} = l/2 and d_{2, 3} = l. I_i are the moments of inertia of the three segments with I₁ = bh³/4 and I₂ = I₃ = bh³/12, E represents Young's modulus of the structure's material, the width b and the length l are defined in figure 3 and h is the structure's thickness.

In addition, U_s is the sum of energy contributions due to the midplane stretching of the segments that occurs during deflection:

$$\begin{equation} U_{s} =\sum _{j=1}^3 \frac{E A_j }{4l}\Bigg[ \int _{c_j}^{d_j} \bigg(\frac{{\rm d}w_j}{{\rm d}x}\bigg)^2 {\rm d}x \Bigg]^2, \end{equation} \tag{ 4 }$$

where c₁ = 0, d₁ = l/2, c_{2, 3} = l/2 and d_{2, 3} = l. A_j in (4) represents the cross section areas of the three segments with A₁ = 3bh and A_{2, 3} = bh. The total potential energy function can hence be written as:

$$\begin{equation} \Phi = U_{tf}-W, \end{equation} \tag{ 5 }$$

where W = F w₂(l) is the work done by the punctual force F.

Applying the principle of minimum potential energy, one sets for each unknown coefficient a_i:

$$\begin{equation} \frac{\partial \, \Phi }{\partial \, a_i}=0. \end{equation} \tag{ 6 }$$

The problem is then augmented with Lagrange multipliers by considering a set of constraints. To obtain satisfactory results, mechanical constraints do not need to be considered. Enforcing geometrical constraints is enough to provide an accurate solution. Assuming that the central tine and the outer tine are rigidly linked to the overhang at x = l/2, and that the central and outer tines are terminated by sliding conditions (see figure 3), the following set of conditions can be enforced:

$$\begin{equation} w_1(0)=\frac{{\rm d}\,w_1}{{\rm d}x}(0)=0 \end{equation} \tag{ 7 }$$

$$\begin{equation} w_1(l/2) = w_2(l/2) = w_3(l/2) \end{equation} \tag{ 8 }$$

$$\begin{equation} \frac{{\rm d}w_1}{{\rm d}x}(l/2) = \frac{{\rm d}w_2}{{\rm d}x}(l/2) = \frac{{\rm d}w_3}{{\rm d}x}(l/2) \end{equation} \tag{ 9 }$$

$$\begin{equation} \frac{{\rm d}w_2}{{\rm d}x}(l) = \frac{{\rm d}w_3}{{\rm d}x}(l) = 0. \end{equation} \tag{ 10 }$$

The augmented system can then be numerically solved for different values of F. To compute such a system, an algorithm such as the one reported in [26] was used with the following numerical values: l = 12.5 mm, b = 0.25 mm, h = 0.1 mm and E = 212 GPa.

To validate this modeling approach, FEA was used as a reference tool. COMSOL simulations were conducted in the nonlinear deflection mode with the values of l, b, h and E aforementioned. We also specified a density of 8030 kg m⁻³ and a Poisson's ratio of 0.29 for the structure's material. These values were selected in accordance with the material utilized for the prototype of section 4. For a point forceF = 260 mN applied at the half-span of the CB, figure 4 shows the deflection profile of the structure. Analytical results are in good agreement with FEA simulations and permit to predict the deflection of the structure accurately. It can also be seen that the deflection amplitude of the OB is smaller than the deflection of the CB.

Zoom In Zoom Out Reset image size

This intuitive result is also confirmed by figure 5 which compares the deflection amplitude of both beams as a function of the force F applied. Analytical and FEA results are again in accordance. Figure 5 also clearly illustrates that the stretching (i.e. nonlinear) effect progressively dominates as the beams deflection increases. Thereby, for deflections above ∼50% of the structure's thickness, the structure becomes clearly stiffer.

Zoom In Zoom Out Reset image size

We now demonstrate that, if the force F is applied while the OBs oscillate as shown in figure 1(b), the variations of deflection imposed by F actually impact the initial resonance frequency of the OBs. To demonstrate such a coupling between static and dynamic behaviors, a symmetric structure again simplifies the analysis. Indeed, valuable insights into the structure's dynamics can be grasped by restricting the analysis to a single OB (see figure 6).

Zoom In Zoom Out Reset image size

To a first approximation, one can choose to model the OB as a hinged–hinged beam terminated by two rotational springs of stiffness k_r1 (see figure 6). The presence of rotational springs is justified by the facts that extremities of the OB are attached to the overhangs. Intuitively, it can be foreseen that, although the overhangs can be considered as axially immovable ends, they should not act as ideal clamping supports. In the literature, elastic rotational springs are often used to model such a flexibility [27–29].

Again, energy approaches can be exploited to conveniently handle the configuration of figure 6. Because the OB is terminated by rotational springs attached to pinned ends supports, it is conventional to assume a displacement function in the form of a sine:

$$\begin{equation} w(x)=D_s \, \sin \Big(\frac{\pi x}{l}\Big), \end{equation} \tag{ 11 }$$

where D_s is the midspan deflection of the OB.

D_s is related to the force F applied upon the central steady beam, and can be estimated from the previous static analysis (see figure 3):

$$\begin{equation} D_s= w_3(l)-w_1(l/2). \end{equation} \tag{ 12 }$$

Conceptually, it is however more convenient to consider that D_s is rather due to an unknown equivalent virtual force F_v (see figure 6). With an energy approach, a relationship between D_s and F_v can be derived. Considering that F_v is applied at the beam midspan, the potential energy stored by the beam is the sum of the following contributions:

$$\begin{equation} U_T= U_b+U_s+U_{{\rm rs}}\,. \end{equation} \tag{ 13 }$$

In (13), U_b is the bending energy of the beam:

$$\begin{equation} U_b=\frac{E I }{2}\, \int _0^{l} \bigg(\frac{{\rm d}^2w}{{\rm d}x^2}\bigg)^2 \, {\rm d}x, \end{equation} \tag{ 14 }$$

U_s is the energy developed during the midplane stretching of the beam:

$$\begin{equation} U_s=\frac{E A}{8\,l}\, \Bigg[ \int _0^{l} \bigg(\frac{{\rm d}w}{{\rm d}x}\bigg)^2 \, {\rm d}x \Bigg]^2 \end{equation} \tag{ 15 }$$

and U_rs is the energy stored by the two rotational springs of stiffness k_r1:

$$\begin{equation} U_{{\rm rs}}= \frac{k_{r_1}}{2} \, \theta _0^2 + \frac{k_{r_1}}{2} \, \theta _l^2, \end{equation} \tag{ 16 }$$

where θ₀ and θ_l are the slopes $\frac{{\rm d}w}{{\rm d}x}$ evaluated at x = 0 and x = l, respectively.

Considering that the work done by F_v is W₁ = F_v w(l/2), the total potential energy function Φ₁ is:

$$\begin{equation} \Phi _1= -F_v\, D_s + \Bigg[ \frac{\pi ^2 k_{r_1}}{l^2} + \frac{Ebh^3\pi ^4}{48l^3} \Bigg] D_s^2+ \frac{Ebh\pi ^4}{32l^3} D_s^4\,. \end{equation} \tag{ 17 }$$

Minimizing (17) with respect to D_s yields a cubic force-centered-deflection law:

$$\begin{equation} F_v= k_1\,D_s + k_3\,D_s^3, \end{equation} \tag{ 18 }$$

where

$$\begin{equation} k_1= \frac{Ebh^3\pi ^4}{24l^3} + \frac{2k_{r_1} \pi ^2}{l^2} \qquad \mathrm{and} \qquad k_3= \frac{ Ebh\pi ^4 }{8l^3}\,. \end{equation} \tag{ 19 }$$

In (19), k₁ and k₃ are linear and nonlinear spring constants, respectively.

As in [30], it can now be assumed that small deflections of the beam about a mean deflection D_s can be described approximately by a single stiffness value. An effective spring constant is found by deflecting the beam from its equilibrium position by an arbitrary amount Δ_d (see figure 6), so that:

$$\begin{equation} F_v= k_1\, (D_s + \Delta _d) +k_3 (D_s+\Delta _d)^3. \end{equation} \tag{ 20 }$$

Considering that vibration amplitudes Δ_d are sufficiently small, terms proportional to Δ²_d and Δ³_d in (20) can be neglected. Then, an equivalent stiffness of the OB may be approximated by:

$$\begin{equation} k_{{\rm eq}}= \frac{{\rm d}F_v}{{\rm d}\Delta _d} \approx k_1 + 3\,k_3\,D_s^2. \end{equation} \tag{ 21 }$$

As often in structures involving coupled beams, vibrating beams can be modeled to a first approximation as one degree of freedom oscillators (e.g., [31, 32]). If the single OB is assumed to oscillate as an undamped lumped-parameter system, the static deflection D_s affects its natural frequency as follows:

$$\begin{equation} \frac{f}{f_0}\approx \bigg[1+\frac{3\,k_3}{k_1}\,D_s^2\bigg]^{1/2}. \end{equation} \tag{ 22 }$$

As previously, the numerical values E = 212 GPa, b = 0.25 mm, l = 12.5 mm and h = 0.1 mm were used to compute (22). An order of magnitude for the rotational stiffness k_r1 was also found by using [33]:

$$\begin{equation} k_{r1}=\frac{4\,E\,b\,h^3}{3\,l}. \end{equation} \tag{ 23 }$$

Equation (23) gives k_r1 ≈ 5.6 × 10⁻³ Nm rd⁻¹. Nonetheless, (23) is originally intended to estimate the rotational spring constant of anchors for clamped–clamped beams. In our case, the OBs are attached to the overhangs (see figure 3) which do no act as ideal clamps. Therefore, the value of k_r1 was slightly adapted to fit prestressed modal analyses conducted with COMSOL. For k_r1 = 3.45 × 10⁻³ Nm rd⁻¹, figure 7 proves that (22) accurately predicts the amount of frequency change engendered by a force F applied to the CB.

Zoom In Zoom Out Reset image size

Alternatively, and because displacements of the CB can also be predicted with respect to the force F applied, the amount of frequency change can also be estimated as a function of a displacement imposed on the CB (see figure 8). It is of interest to note that in both cases, quasi-linear regions are predicted around points P. In particular, a force sensitivity of 2.91 Hz mN⁻¹ is predicted if the force applied to the CB exceeds 100 mN (see figure 7).

Zoom In Zoom Out Reset image size

Experiments were first carried out to investigate the static deflection of the prototype. To bend the structure, an indenter terminated by a metal bead with a diameter of 500 μm was placed beneath the CB. The stiffness of the metal bead was much higher than the stiffness of the structure. With a micropositioning stage (Physik Instrumente M112-1DG), incremental motions of 10 μm along the z direction were imposed to the CB. Considering the prototype symmetry, deflection data were acquired only for one OB.

Figure 10 demonstrates that the prototype deflected as expected during experiments, even though the deflection of the OBs tends to be 12% overestimated by theory when the deflection of the CB is above 30 μm. Several uncertainties sources could explain this slight estimation error (e.g., values of Young's modulus and Poisson's ratio used for calculations, uncertainties related to the fabrication process, the clamping of the prototype, etc).

Zoom In Zoom Out Reset image size

The dynamic behavior of the prototype was then explored. In the first step, the indenter was removed and no force was applied to the CB. To actuate the prototype, the PZT element was driven with a sinusoidal voltage. The peak–peak amplitude of this sine signal after the power amplifier was 9V. The excitation frequency was then swept with the function generator in order to find the resonance modes of the prototype. As illustrated by figure 11(a), the antisymmetrical mode of interest was found at 3180 Hz.

Zoom In Zoom Out Reset image size

When the PZT element was driven with a voltage supply of 9 V, the peak–peak oscillation amplitude of the OBs was 22 μm. During oscillations of the OBs, a lot of attention was paid to check if potential vibrations were transferred to the CB. As a matter of fact, they were very limited. As shown in figure 11(b), they never exceeded 500 nm, that is to say 2% of the oscillation amplitude of the OBs. Droplets of water were also deposited upon the CB with a micropipette. No particular effect on the droplets was observed. It was also confirmed that the presence of the small amount of liquid on the CB did not alter the vibration mode.

In addition, we explored the frequency response of the antisymmetrical vibration mode. With the OBs oscillating in air, a quality factor of 700 was obtained (see figure 12). This quality factor was not deteriorated when a liquid droplet was deposited on the CB. This is a high value when compared to typical quality factors of AFM cantilevers in liquids (10–30).

Zoom In Zoom Out Reset image size

The coupling between static and dynamic behaviors was then investigated. To that end, oscillations of the OBs were monitored while the CB was bent with the indenter. For a total displacement of the indenter of 250 μm, the oscillation frequency of the OBs evolved as shown in figure 13.

Zoom In Zoom Out Reset image size

Contrary to theoretical predictions, it is experimentally observed that for small deflections of the CB (i.e. deflections lower than 100 μm), the frequency of the OBs actually decreases. Such a behavior has been reported for buckled and deflected beams subjected to axial loads [35–39]. Nevertheless, this decrease finally appears as a transition period. Indeed, once the stretching effect dominates at large deflections, a behavior very similar to the one initially predicted is retrieved. In particular, it is worth noting that the curve in figure 13 offers a linear displacement sensitivity of 5.75 Hz μm⁻¹ around the point S. This displacement sensitivity is actually ∼17% better than the one expected initially by theory (4.8 Hz μm⁻¹).

From the analysis of figure 13, it is clear that the prototype provides a much better sensitivity when the CB is largely bent. The CB was hence maintained permanently curved (i.e. even with the indenter removed) to exploit the region around the point S. Obviously, this axial translation compressed the whole structure. As a result, the initial resonance of the antisymmetrical vibration mode decreased from 3180 to 3080 Hz. Nevertheless, no significant impact on the dynamic behavior of the prototype was noted.

To estimate the force F applied to the CB during small displacement intervals, the force developed by our indenter was characterized. By pressing the indenter against a precision electric scale (Kern 440–33), a linear relationship was found (see figure 14(a)). Then, frequency variations with respect to the force generated by the indenter were measured. With the prototype pre-deflected, only a linear increase of the frequency was observed, as expected (see figure 14(b)). This linear increase corresponded to a force sensitivity of 2.56 Hz mN⁻¹. Although a pre-deflection of the whole structure was not included in our theoretical analysis, the force sensitivity measured experimentally was only ∼12% lower than expected.

Zoom In Zoom Out Reset image size

With the suspended structure, the force applied to a cell can be measured. In addition, the cell deformation must also be determined. Conventionally, force-deformation information is then used in conjunction with Hertz theory to retrieve Young's modulus of the cell. Hertz theory, however, is subjected to a number of important assumptions. In particular, when dealing with a living cell, the condition of negligible force adhesion is not always fulfilled. Furthermore, for suspension cells, Hertz theory requires to know the radius of curvature of the cell to be probed. As an alternative, we preferred to calibrate our structure so that a sufficient estimation of Young's modulus of the cell could be rapidly obtained without the use of a descriptive model. For calibration, two gel samples T7 and T5 were bought from the company Gelmec. They had a Young's modulus of 37.5 kPa and 119 kPa, respectively.

The gel samples were manually prepared so that their size was as similar as possible to the size of the cells targeted. Thereby, efforts were made to obtain small cubic pieces whose edges measured approximately 500 μm. The gel samples were then cautiously placed on a flat indenter and compressed upon the CB of the structure. To probe only the elastic properties of the gels, the amount of compression applied did not exceed 10% of their thickness [40]. Moreover, the velocity of the translation steps was kept very slow (4 μm s⁻¹) to minimize the occurrence of viscoelastic effects [18, 41].

The same protocol was repeated with a biological entity. Because embryos or oocytes could not be obtained easily, we used a lobster egg with a diameter of 500 μm as a viable alternative (see figure 15). Linear frequency variations induced during the compression of the gel samples and the lobster egg are compared in figure 16(a).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To determine more precisely the elastic modulus of the lobster egg, slopes S_T5 and S_T7 indicated in figure 16(a) were used as single values plotted against Young's modulus of the gels (see figure 16(b)). Since only the elastic properties of the materials were probed, a linear regression equation linking these two values could be determined [42]:

$$\begin{equation} S_{\rm egg}= 0.02914 \times {E_{\rm egg}} - 0.2329, \end{equation} \tag{ 24 }$$

where S_egg is the slope of frequency variations measured for the lobster egg and E_egg is its Young's modulus. For S_egg = 2.03 Hz mN⁻¹, (24) yields a Young's modulus of 78 kPa for the lobster egg. Although no reference value has been found for lobster eggs in the literature by the authors, Young's modulus obtained is in accordance with orders of magnitude usually reported for most living cells [43–45].