The mechanical properties of plasma membranes have been shown in recent years to play an important part in physiology and cell biology. An archetypal system for study is the red blood cell (RBC) membrane, the elastic character of which determines many of the flow properties of blood and enables the red blood cells to transverse the capillaries. The RBC membrane consists of a mixed lipid-protein bilayer, an attached intracellular network called the membrane cytoskeleton, and an extracellular glycocalyx (a layer of glycoproteins covering the outer surface of the cell), thought to have little or no affect on the mechanical properties. The main component of the cytoskeletal network is the protein spectrin. Spectrin cross-links with very short actin filaments, resulting in junctional nodes that are approximately five- or sixfold coordinated in spectrin. Although the red cell is unique, in that it contains no cytoplasmic cytoskeleton or organelles, the complexes of the proteins that make up its associated membrane skeletal network and give it its elastic properties also occur in the plasma membranes of other cell types.

The RBC at rest assumes a biconcave discoid shape with a diameter of ~8 μm and it is capable of passing through capillaries with less than half the diameter. To allow the RBC to deform inside a capillary, the resistance to deformation of the RBC membrane must not be too large. The resistance, on the other hand, must not be too small, otherwise the cell's integrity would not be preserved during normal flow in the circulatory system. Bending deformability of the RBC membrane is also shown by its mechanical out-of-plane and in-plane fluctuations with amplitudes in the range of 10–400 nm (1–3).

A long-standing problem in the study of RBC structure is that the observation of thermal fluctuations of nanometer-size seems to be consistent only with a vanishingly small shear modulus (2), whereas static deformation experiments such as micropipette aspiration (4) report a large shear modulus. Strey et al. (2) studied the flickering of RBC to derive the mean-square amplitudes of the first three azimuthal eigenmodes and they obtained a shear modulus, μ, 100-times smaller than the static value measured by micropipette aspiration, μ = 6 × 10⁻⁶ N/m (4), with a fluorescent multiparticle tracking system μ = 1–10 × 10⁻⁶ N/m (3) or by optical tweezers, μ = 2.5 × 10⁻⁶ N/m (5) and μ = 11.1–17.7 × 10⁻⁶ N/m (6). They concluded that the shear resistance for small deformations vanishes and therefore it does not contribute to thermal shape fluctuations.

We have chosen to analyze the problem considering both shear and bending moduli in the attempt to reconcile a finite value of the shear modulus with the presence of the nanometer-scale fluctuations observed experimentally. The essential assumptions of the model are that the membrane preserves a constant global area under thermal fluctuations and its mechanical behavior can be described by two material parameters, an out-of-plane bending modulus, B, which describes the resistance to bending, and an in-plane shear modulus, μ, which describes the resistance to shape changes of the membrane surface. It is usually assumed that the bending modulus, B, is dominated by the microscopic properties of the lipid bilayer. Equivalently the shear modulus, μ, is thought to be dominated by the microscopic properties of the cytoskeleton, although protein-protein interactions or domains of lipids in solid-state form can also contribute. Throughout the article, we shall use “membrane” to refer to the composite structure formed by the lipid bilayer, its integral proteins, and associated membrane skeleton.

Our model envisages the RBC membrane as a two-dimensional network consisting of N particles arranged in a regular triangulation. We emphasize that these are virtual particles that represent discrete areas of the membrane. For small deformation, a two-dimensional network with sixfold symmetry can be described by two elastic moduli—the area compression modulus K_A and the shear modulus, μ (7). In our model we keep the global surface area constant, so we can assume that K_A ~ ∞. The two-dimensional network is embedded in a three-dimensional closed surface where it is allowed to fluctuate in all directions (see Fig. 1 a). Equations of motion of the virtual particles are solved in a molecular-dynamics fashion (8,9) to derive the time-evolution of the system and calculate the time-dependent properties.

FIGURE 1 Particle-network representing the whole-cell model is reported in a. In b we illustrated a portion of the network and the dihedral angle ijml between the adjacent triangles Δijm and Δjml.

Each particle is tethered to its six nearest-neighbors (see Fig. 1 b) by linear springs, giving a potential energy of

(1)

A dihedral angle potential is applied between adjacent triangles. This potential is used in molecular dynamics simulations to describe interactions arising from torsional forces in molecules (8,9). The expression of the dihedral angle potential we used is (10)

(2)

We stress here that the bending and shear moduli are not independent (16), so in general it is not possible to change k (or D) and independently change the value of μ (or B). Details are given below about the derivation of μ and B that we adopted in our work.

In general the free energy of a system is given by two contributions, F = E–TS. In the case of a spring-network, for example, the energetic contribution E is given by the sum of over all the particles. The entropic contribution TS depends in general on the number of configurations available to the system at a particular state. Both contributions influence the system elasticity. A clear sign that the elasticity is entropy-dominated is that the shear modulus increases when the temperature increases. It is not clear yet which of the two contributions is the dominant one in determining the RBC membrane elasticity (7). In the present model the elasticity of the system is dominated by the energetic contribution. This should not restrict our analysis since we study membrane fluctuations at a given temperature and with a given shear modulus. Furthermore, the constraint is readily overcome since removing bonds between neighboring particles within our model allows simulation of a system where the entropy can be the dominant contribution to the elasticity (17).

As in any standard molecular dynamics simulation, we define the temperature, T, of the system as

(3)

Since the lipid bilayer is relatively incompressible, we require that the global surface area of the membrane, defined as the sum of the areas of all triangles by which the surface is tessellated, remains constant. To do this we apply Lagrange's method of undetermined multipliers. We chose the global area, A_tot, to be the value given by triangles whose sides are all equal to the equilibrium length, r₀:

(4)

The system used in the simulations consists of N particles (vertices) embedded in a closed surface. For the initial particle-configuration we followed Discher et al. (20), who used two parallel sheets in the form of a hexagon. Each particle on the perimeter of one sheet is connected to two particles on the perimeter of the other sheet. There are six “corners” on each hexagonal sheet, each of which has only fivefold coordination. The minimum number of fivefold defects required by topology for the triangulation of a spherical surface is 12 and so this configuration represents one possible triangulation of a spherical surface. All particles move in space but the connectivity is fixed, in that each particle has a fixed set of neighbors: 12 particles have five neighbors and all the others have six neighbors.

The molecular dynamics simulations were performed using the DL_POLY simulation package version 2.12 developed in Daresbury Laboratory, Cheshire, UK (10). The equations of motion of the particles were integrated using a Verlet leapfrog integration algorithm in conjunction with a Nosé-Hoover thermostat. We developed a subroutine that applies Lagrange's method of undetermined multipliers to keep the total surface area constant.

Because of computational time-demand considerations we used N = 5582 for number of particles. In most simulations we used an equilibrium length r₀ = 100 nm. The area per particle, A_particle, is given by

(5)

We used two values of the spring constant, k₁ = 5 × 10⁻⁶ N/m and k₂ = 25 × 10⁻⁶ N/m, and two values of the dihedral-angle potential constant, D₁ = 6 × 10⁻²⁰ J and D₂ = 60 × 10⁻²⁰ J. We chose these values because the derived shear and bending moduli are in the range of the experimental measurements. We observed that using higher values of the k and D gives an increasingly rigid behavior of the model-network. On the other hand, lower values of k and D generate crumpled configurations, which were not worth studying with the present model that lacks a repulsive interaction between particles and a constraint on the cell volume. Therefore, we performed four main simulations combining the two chosen values of k and D. As we report in the following section, during the simulations we were able to keep track of the thermal fluctuations and to estimate the shear and bending moduli of the membrane.

We observed a common feature of the evolution of the shape of the cell models. After an equilibration phase (~10⁷ time steps), the cell model attains an oblate shape for any of the values of k and D used (see Fig. 2). The oblate shape persists up to the maximum simulation time we reached (4 × 10⁷ time steps). This is also evident in the behavior of the correlation function (reported in Fig. 3). The correlation function is defined as where is the z-coordinate unit vector, …_t is the time-average, and is the instantaneous average of the unit vectors perpendicular to the triangles that fall within r and r+Δr from the cap-center. In any flat region of the membrane, f(r) has a value close to unity.

FIGURE 2 Overall shape of the cell: top and side views for four different cases. Case LS, k = 25 × 10−6 N/m, B = 6 × 10−20 J. Case LL, k = 25 × 10−6 N/m, B = 60 × 10−20 (more ...)

FIGURE 3 Correlation function, versus r/r0. Case LS (), k = 25 × 10−6 N/m, B = 6 × 10−20 J; Case LL (□), k = 25 × 10−6 N/m, B = 60 × 10−20 J; Case (more ...)

It is not the purpose of our current work to investigate the overall shape of the red blood cell. Further work is needed to assess how the shape is affected by the presence of the 12 defects at the rim of the cell. Moreover a mechanism to keep the total volume constant needs to be introduced in the model to study the overall shape of the cell. However, we observed that once the system has reached equilibrium, the volume fluctuations account for <1% of the total volume. This feature corresponds well to the cell behavior observed in the experiments, where the volume remains essentially constant. This allows us to infer that our findings are not significantly affected by volume fluctuations.

The center of the cell cap is almost flat and parallel to the x-y plane. We exploited this feature and decided to monitor a small portion of the membrane in this region to study the thermal fluctuations. We consider 19 particles located at the center of the cell cap (the lighter patch in Fig. 4), which correspond to an area of ~1.6 × 10⁵ nm². We define the average of the z-displacement of the 19 particles as the out-of-plane instantaneous fluctuations of the patch. We also keep track of the positions of single particles, which correspond to fluctuations of the smallest area accessible to our model (~9740 nm²). More specifically, we monitor the particle at the center of the cap and four other particles at orthogonal positions around the cell (see dashed arrows in Fig. 4) and monitor the displacements of the particles in the membrane plane (in-plane fluctuations). Our resolution is comparable with the resolution used in flicker measurements. Tuvia et al. (19), for example, measured membrane fluctuations over a cell patch with an area of 2.5 × 10⁵ nm² using point-dark-field microscopy.

FIGURE 4 View of the cell from the top. The lighter spot represents the 19-particle patch. The four dashed-arrows indicate the positions where the rim fluctuations are taken. The solid arrows indicate the directions along which the amplitudes, U, are taken.

To derive an estimate of the shear modulus we consider the 19-particle patch and we use the average of the x- and y-displacements of the 19 particles as the in-plane fluctuations of the patch. The estimate of the shear modulus is obtained using the expression reported by Lee and Disher (3),

(7)

(8)

To derive the bending modulus we used the formulae reported by Helfrich and Servuss (14),

(9)

(10)

From each simulation we can easily calculate and ΔA_patch/A_patch, so we can derive the relative values of B and σ using Eqs. 9 and 10. In the case in which the lateral tension is zero, Eqs. 9 and 10 simplify to

(11)

(12)

In Table 1 we report the main results of our simulations. The first two rows show the values of k and D we used as inputs in four different simulations, the results of which are displayed in the columns LS, LL, SS, and SL. The following six rows report the root-mean-squared (RMS) amplitude, U, of the fluctuations:

• Cap out-of-plane fluctuations,
• Cap in-plane fluctuations, and in the x- and y-directions, respectively.
• Rim out-of-plane fluctuations,
• Rim in-plane fluctuations, and along the smaller and bigger curvatures, respectively.

TABLE 1 Main results from four simulations

The first notable result in Table 1 is that the fluctuations, in particular the out-of-plane fluctuations, are in the nanometer scale (~20–60 nm) as observed in experimental measurements (1,2). This is significant because these fluctuations are observed in conjunction with values of the shear and bending moduli that are of the same order of magnitude of those measured experimentally. Case LS in Table 1 is particularly interesting: is 55 nm; the calculated bending modulus, B = 3.7 × 10⁻²⁰ J, agrees well with experimental measurements, 1.3–3 × 10⁻²⁰ J (18), and 2.3 × 10⁻²⁰ J (24), as does the calculated shear modulus, μ = 6.7 × 10⁻⁶ N/m, with the relative experimental values, 6 × 10⁻⁶ N/m (4), 1–10 × 10⁻⁶ N/m (3), and 2.5 × 10⁻⁶ N/m (5). Moreover, and calculated with our simulation (14 and 15 nm), are quite close to the experimental values, 17–19 nm (3). Fig. 5 shows our results for μ and B together with an indication of the experimental measurements obtained with different techniques.

FIGURE 5 Values of the shear modulus μ and bending modulus B calculated as described in the text compared with experimental values measured by different methods as reported in the literature. The calculated values, LS, SS, LL, and SL (·), refer (more ...)

For each of the four cases, and coincide within the variance. This means that the flat portion of the membranes behaves like an isotropic two-dimensional network. Furthermore, in cases LL, SS, and SL, and coincide within the variance, showing an isotropic behavior. In case LS, on the other hand, where the differences between the two principal curvatures are more pronounced (see Fig. 2), and do not coincide. The fluctuations in that region are not isotropic. Lee and Discher (3) studied the thermal fluctuation of the red blood cell by measuring the displacement of a fluorescent nanoparticle attached to the membrane. They found isotropic fluctuations when the particle was centrally attached in a spherical contour of the cell, but anisotropic fluctuations when the particle was located in a highly strained region where the principal curvatures were different. We observed the same phenomenon with the out-of-plane fluctuations. The different values of and reported in Table 1 are also due to the different curvature in the two regions.

In Fig. 6 we show the histogram of the out-of-plane fluctuation amplitude for the case LS (k = 25 × 10⁻⁶ N/m, D = 6 × 10⁻²⁰ J) and for the case LL (k = 25 × 10⁻⁶ N/m, D = 6 × 10⁻¹⁹ J). Thermal fluctuations follow a Gaussian distribution. We observe in both cases shown in Fig. 6 that the amplitude can be approximated by a Gaussian distribution and the amplitude for the case LL is narrower due to the higher value of the dihedral constant D. This is expected since a higher value of D gives a higher value of the bending modulus, which induces a lower RMS amplitude in the out-of-plane fluctuations. Equivalently, the RMS amplitude of the in-plane fluctuations is lowered by higher values of the shear modulus, which are obtained using a higher value of the spring constant k.

FIGURE 6 Histogram of the out-plane fluctuation amplitude. (Dotted line) Case LS (k = 25 × 10−6 N/m, D = 6 × 10−20 J). (Solid line) Case LL (k = 25 × 10−6 N/m, D = 6 × 10−19 (more ...)

In Table 1, simulations using the same dihedral constant give similar bending moduli, suggesting that the bending resistance is not affected significantly by the value of the spring constant k for the range of values we used. The results for the shear modulus show a more interesting behavior. The shear modulus for cases LL and SL obtained with Eq. 7 agrees well with the prediction of Eq. 8. The shear moduli obtained with Eq. 7 for cases LS and SS, instead, are 40% and 30% lower than the prediction of Eq. 8. We believe that this deserves more detailed analysis, which we report in the Appendix.

The lateral tension reported in Table 1 ranges between 0.5 and 2.6 × 10⁻⁶ N/m. These values are three or four orders-of-magnitude smaller than the isotropic lysis tension, 3–12 × 10⁻³ N/m (25). We stress that the lateral tension proscribes the use of Eqs. 11 and 12, which are valid only for vanishingly small values of the lateral tension.

The model we implemented showed itself able to describe the elastic properties of the red blood cell membrane related to thermal fluctuations. We demonstrated that it is possible to observe nanometer-scale fluctuations in conjunction with values of the shear and bending moduli similar to those obtained experimentally. In particular, given the microscopic parameters k and D and the assumption of constant area, the model reproduces the experimental values of the fluctuation amplitudes, the shear and bending moduli. We must emphasize that the constraint of constant area is an important feature of the model, since without it the model would show, in general, different dynamics. Moreover, it is essential to use Eqs. 9 and Eq. 10 properly. We found that the approximate expressions in Eqs. 11 and 12 give incorrect and inconsistent values of the bending modulus, due to the presence of a non-vanishingly small lateral tension. This highlights the importance of measuring ΔA_patch/A_patch (or related quantities) when studying membrane fluctuations experimentally and of properly choosing the minimum wavelength appearing in Eqs. 9 and 10.

In this section we analyze the behavior of the shear modulus reported in Table 1. The shear modulus for cases LL and SL agrees well with the prediction of Eq. 8, which is valid for T = 0 K. This suggests that the shear modulus is not affected by the temperature or by the out-of-plane fluctuations of the system. In cases LS and SS, instead, the shear modulus obtained with Eq. 7 is 40%, and 30%, respectively, lower than the prediction of Eq. 8. We believe that this softening of the shear modulus can be explained using rescaling arguments reported by Aronovitz and Lubensky (13,26). These authors used renormalization group techniques to study the crumpling transition and the fluctuations around a flat phase of flexible D-dimensional crystalline membranes embedded in d-dimensional space. They reported that for fluctuations of wavelength, λ, greater than the nonlinear length, ξ_nl, the shear modulus is renormalized and behaves as μ(λ) = (1/λ)^η, the nonlinear length being defined as for D = 2, where and are the bare elastic moduli divided by the temperature. Using T = 309 K and the values of D and μ_T=0 (Eq. 8) to derive ξ_nl, the following values are obtained: and For the shear modulus to be renormalized and then to undergo softening, the critical exponent η must be positive and the fluctuation wavelengths of the system must be larger than ξ_nl. In cases LL and SL the nonlinear length is comparable with the size of the system, and for this reason the renormalization of the shear modulus is less significant. In cases LS and SS, instead, ξ_nl is much smaller than the size of the system and the conditions for the softening of the shear modulus are met.

To corroborate further the conclusion that the softening of the shear modulus can be ascribed to renormalization, we performed simulations on a hexagonal sheet (open surface) with the same characteristics as the LS case. We used hexagonal sheets with three different diameters, namely, 40r₀, 60r₀, and 100r₀. We obtained the following values for the shear modulus: μ(40r₀) = 9.9(0.1), μ(60r₀) = 9.0(0.3), and μ(100r₀) = 8.2(0.3) 10⁻⁶ N/m, which show that the shear modulus diminishes as the size of the system increases. These are preliminary results and more systematic work should be done, but they seem to show an interesting behavior. Considering that our system is realistic in terms of size and range of elastic moduli, we can infer that the renormalization group may play an important role in the understanding of red-blood cell mechanical properties.