Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2876725

Formats

Article sections

- Abstract
- 1 INTRODUCTION
- 2 PROBLEM DESCRIPTION
- 3 Multiscale formulations
- 4 Results
- 5 CONCLUSIONS AND DISCUSSION
- References

Authors

Related links

Phys Rev E Stat Nonlin Soft Matter Phys. Author manuscript; available in PMC 2010 May 26.

Published in final edited form as:

Phys Rev E Stat Nonlin Soft Matter Phys. 2010 March; 81(3 Pt 1): 031904.

Published online 2010 March 4. doi: 10.1103/PhysRevE.81.031904PMCID: PMC2876725

NIHMSID: NIHMS202221

Department of Structural Engineering, University of California, San Diego, La Jolla, CA 92093

See other articles in PMC that cite the published article.

To quantitatively predict the mechanical response and mechanically induced remodeling of red blood cells, we developed a multiscale method to correlate distributions of internal stress with overall cell deformation. This method consists of three models at different length scales: in the complete cell level the membrane is modeled as two distinct layers of continuum shells using finite element method (Level III), in which the skeleton-bilayer interactions are depicted as a slide in the lateral (*i.e.* in-plane) direction (caused by the mobility of the skeleton-bilayer pinning points) and a normal contact force; the constitutive laws of the inner layer (the protein skeleton) are obtained from a molecular-based model (Level II); the mechanical properties of the spectrin (Sp, a key component of the skeleton), including its folding/unfolding reactions, are obtained with a stress-strain model (Level I). Model verification is achieved through comparisons with existing numerical and experimental studies in terms of the resting shape of the cell as well as cell deformations induced by micropipettes and optical tweezers. Detailed distributions of the interaction force between the lipid bilayer and the skeleton that may cause their dissociation and lead to phenomena such as vesiculation are predicted. Specifically, our model predicts correlation between the occurrence of Sp unfolding and increase in the mechanical load upon individual skeleton-bilayer pinning points. Finally a simulation of the necking process after skeleton-bilayer dissociation, a precursor of vesiculation, is conducted.

Among all types of cells, erythrocyte (red blood cell, or RBC) possesses one of the simplest and best characterized molecular architectures. Without a nucleus, a mature erythrocyte contains a cytosol enclosed within a highly flexible yet surprisingly strong membrane. Essential to its structural integrity and stability is this composite membrane consisting of a lipid bilayer supported from inside by a protein skeleton. The connection between the skeleton and the lipid bilayer is achieved at pinning points made of transmembrane proteins.

Despite extensive investigations in the past few decades, there are still many remaining questions about the mechanics of erythrocyte. For example, it is still not fully understood what determines its resting shape (this is the first of eight mysteries about RBC proposed by Hoffman [1]). Herein a pivotal issue is the effect of the protein skeleton upon cell shape. Although a stomatocyte-discocyte-echinocyte sequence was obtained numerically based on the bilayer-coupled hypothesis [2] and the stabilizing function of the skeleton in maintaining the biconcave shape was explored [3], the relaxed reference shape of the skeleton remains controversial. Indeed, if a spherically relaxed skeleton is applied, to obtain the biconcave shape the elasticity of the skeleton must be significantly reduced [4]. Otherwise, nonspherical (*e.g.* biconcave [5] or oblate [2, 3]) relaxed shapes must be assumed. These are beyond the state-of-the-art knowledge about RBC. Moreover, much is unknown about responses of the cell in large deformations. One remaining issue is the strength of the skeleton-bilayer linkage [6]. Under sufficiently large dissociation forces this linkage may rupture, causing cell remodelings such as vesiculation. The recent understanding is based upon the adhesion energy theory [7]. Being essentially phenomenological, this theory does not offer much insight upon the molecular origin of the lipid-skeleton dissociation. In large deformations, the effects of Sp unfolding [8] and dissociation of Sp head-to-head connections [9] upon the mechanical behavior of the cell are also unexplored.

These problems are important not only because RBC serves as a model system for general cell biomechanics, but also because many diseases are related to defects of the inter-protein and protein-to-lipid linkages in the cell membrane [10]. Some of these defects will change the mechanical properties of the cell and its resting shape. Others may induce structural failures of the cell under large loading. For example, in hereditary elliptocytosis (HE), the weakening of the skeleton network reshapes the cell to be elliptical. Cells with abnormal shapes are often destroyed by the spleen, leading to anemia. Mechanically induced cell damage is more pronounced within artificially created flow fields associated with mechanical circulatory support systems [11].

To pave the way for a molecular-level understanding of mechanical responses of erythrocytes as well as the underlying conditions for mechanically triggered structural remodeling and failure, it is vital to quantitatively characterize the mechanical forces acting on the interprotein and protein-to-lipid linkages within the membrane. Toward this end there is also the need to describe the process whereby the protein skeleton, while vertically connected to the lipid bilayer, alters its lateral morphology and density as it deforms. Thus the coupled phenomena of skeletal rearrangements during deformations and skeleton-bilayer interaction are of first order importance to overall mechanical response as well as remodeling processes such as vesiculation, which involves a separation of the skeleton from the bilayer, and related protein sorting events.

In this study we explore a finite element method (FEM) to simulate the membrane as two distinct layers. Although energy minimization methods have been successfully used in membrane mechanics [2, 4, 12], FEM is still a good candidate with many advantages, *e.g.* the robustness for contact/adhesion calculation and time-dependent problems, the high computational efficiency, as well as the simplicity in coupling with other methods [5, 13, 14]. Unlike the coarse-grained models [4, 12], in our model we explicitly compute the interaction between the bilayer and the skeleton. Since we aim at molecular-detailed prediction of force distribution within the cell with maximum accuracy at large deformations, we use a multiscale representation that accounts for a molecular-level understanding of the Sp/actin junctional complex (JC) in the protein skeleton and the force-induced Sp unfolding. A specific aim of this study is to evaluate the normal association or dissociation force between the lipid bilayer and the protein skeleton – a sufficiently large dissociation force may trigger phenomena such as vesiculation and tether formation (in both cases segregation of the skeleton from the bilayer is observed) [6, 15, 16, 17, 18].

This paper is organized as follows. In the next section we summarize the development of a set of multiscale models, including a constitutive model of the Sp, a molecular-detailed model of the protein skeleton, and a double-layer FEM model of the complete cell. This set of models is then applied to predict the resting shape and simulate two experiments, micropipette aspiration and cell stretching with optical tweezers. In these studies, we not only compare our predictions with benchmark results, but also document the contact force between the bilayer and the skeleton as well as the post-dissociation behavior. Finally, conclusions and discussion are provided.

An erythrocyte possesses a thin skeletal network coupled with a lipid bilayer. The network is composed primarily of flexible Sp dimers and relatively rigid actin protofilaments, and structurally organized into approximately 33,000 JCs. As shown in Fig. 1, each JC contains a central piece of actin protofilament as well as (up to) six Sp dimers. Horizontally, the JCs are linked *via* the Sp dimer/tetramer reaction (*i.e.* the head-to-head association that connects two dimers into a tetramer). Vertically, this membrane skeleton is connected to the lipid bilayer at pinning points called suspension complexes (SC). Each SC consists of ankyrin, protein 4.2, and band 3 (a transmembrane protein). In addition, the actin protofilaments are also connected to the lipid bilayer through protein 4.1 and glycophorin C (another transmembrane protein). This linkage is usually referred to as the secondary connection. Both band 3 and glycophorin C can drift within the lipid bilayer, rendering a horizontal mobility of the skeleton-bilayer connection.

Inter- and intra-protein reactions are essential to the structural integrity and mechanical response of the cell [19]. In the erythrocyte membrane there exists a dynamic equilibrium as a result of dynamic balances between associated and dissociated states of many inter-protein and protein-to-lipid linkages. These balances are directly affected by mechanical loads. The connections that may rupture under mechanical loads include: the head-to-head association between Sp dimers that links the neighboring JCs [20, 21], the band 3–ankyrin “bridge” that is the connection between the SC and the lipid bilayer [22], the protein 4.1–band 3 connection [23], and the band 3–lipid bilayer connection [18]. These dissociations will dramatically change the mechanical properties of the membrane, causing structural instability and even permanent damage. For example, as suggested by micropipette aspirations and flow channel experiments, the rupture of the band 3–lipid linkage will cause detachment of the skeleton from the lipid bilayer, leading to creation of vesicles (vesiculation) or tethers [16, 18]. Vesiculation is associated with membrane loss, and usually with changes in the shape and physiology of the cell. To understand conditions related with these cell remodelings, numerical models are required to correlate cell deformations to the mechanical loads on these connections.

A multiscale model is illustrated in Fig. 2. At the complete cell level (Level III) (Fig. 2a), the membrane is modeled as two vertically connected layers (vertical to the mean plane of the membrane), the lipid bilayer and the protein skeleton. Both layers are modeled as continuum shells. The constitutive properties of the inner layer (the protein skeleton) are evaluated by using a molecular-detailed model of a JC (Level II). As depicted in Fig. 2b, in this JC description the skeleton is composed of a quasi-hexagonal network of roughly 6 Sp dimers (the “spokes”) joined to an actin protofilament. Furthermore, we apply a constitutive model (Level I) to calculate the tension-strain relation of each Sp including its folding/unfolding reactions (Fig. 2c).

Multiscale models of an erythrocyte: (a) the double-layer continuum shell model (Level III), (b) the molecular-detailed JC model (Level II. For clarity, the secondary connections with glycophrin C/protein 4.1 are not shown), and (c) the constitutive model **...**

In the following we provide detailed descriptions for each of these models.

Multidomain proteins such as Sp can undergo overstretching due to unfolding of domains or multiple repeats [24, 25, 26, 27, 28] (*ref.* Fig. 2c). As illustrated in AFM experiments, the transient force-extension curve of Sp displays a trademark sawtooth pattern related to unfolding of the domains [24]. Each peak in this curve corresponds to the unfolding of one or more than one domain [26]. The exact characteristic of this curve is dependent on the rate of extension. With an infinitely slow extension rate, the quasi-static (equilibrium) force-extension curve contains a strain-stiffening part before the first unfolding event. Afterwards, a distinctive plateau appears, where the extension force remains almost a constant due to the successive unfolding of domains.

Based on this description, we have developed a dynamic/quasi-static model of a Sp [8]. In this model, we consider a Sp with *N* domains stretched by a force *F*. When there is *N _{f}* folded domains and

$$\frac{x}{{\mathit{\text{N L}}}_{f}}=(1-{\varphi}_{u})\frac{{x}_{f}}{{L}_{f}}+{\varphi}_{u}\frac{{x}_{u}}{{L}_{u}}\phantom{\rule{thinmathspace}{0ex}}\left(\frac{{L}_{u}}{{L}_{f}}\right)\phantom{\rule{thinmathspace}{0ex}}.$$

(1)

In the equilibrium state, by considering the balance between the unfolding and the folding processes *via* an Arrhenius rate relation we express *ϕ _{u}* as

$${\varphi}_{u}=\frac{\text{exp}\phantom{\rule{thinmathspace}{0ex}}\left(\frac{(F-{F}_{1/2})\mathrm{\Delta}\mathrm{\Delta}x*}{{k}_{B}T}\right)}{1+\text{exp}\phantom{\rule{thinmathspace}{0ex}}\left(\frac{(F-{F}_{1/2})\mathrm{\Delta}\mathrm{\Delta}x*}{{k}_{B}T}\right)},$$

(2)

where ΔΔ*x** = Δ*x _{f→u}*−Δ

Using this model,we have quantitatively reproduced the experimentally measured force-extension relation of a Sp [8].

Our model of a JC is based upon the three-dimensional description of a single JC unit by Sung and Vera [29] (also see [30]) and implemented within the mechanical model of Zhu *et al.* [31] and Zhu and Asaro [8]. According to this depiction (Fig. 1), in a JC the protofilament functions as a mechanical axis, anchoring 3 pairs of Sp. Each Sp pair is arranged in a back-to-back fashion and the force-extension relation comes from previous Sp model (Level I). The result is a quasi-hexagonal network composed of JCs.

The dynamic and quasi-static responses of individual and multiple JCs interacting with the lipid bilayer have been studied *via* a hybrid scheme that simulated the response to thermal fluctuations and applied displacements or stresses [8, 31].

In the complete cell model, both the lipid bilayer and the protein skeleton are modeled as continuum shells (Fig. 2a). The elastic properties of the skeleton layer comes from the JC model (Level II).

Following Evans and Skalak [32], we calculate the constitutive properties of a JC undergoing finite deformations based upon the strain energy Φ stored in the Sp. In this approach, an arbitrary in-plane deformation is achieved by stretching along two orthogonal axes, axis 1 and axis 2, with stretching ratios λ_{1} and λ_{2} (λ_{1} ≥ λ_{2}). Two independent deformation parameters, α = λ_{1}λ_{2} − 1 and $\beta =({\lambda}_{1}^{2}+{\lambda}_{2}^{2})/(2{\lambda}_{1}{\lambda}_{2})-1$ are also defined. It is seen that α represents area change, and β a change of aspect ratio or eccentricity. The shear modulus of the skeleton, *µ _{s}*, is then given as

$${\mu}_{s}={\frac{1}{{A}_{0}}\frac{\partial \mathrm{\Phi}}{\partial \beta}|}_{\alpha},$$

(3)

where *A*_{0} is the projected area without deformation occupied by the JC (*i.e.* the area of the hexagon formed by the six SCs). As α is constant, it corresponds to an anisotropic deformation. The potential energy stored inside each Sp dimer is calculated by integrating its internal tension times the extension, starting from the natural state. Φ is then evaluated by summing up the total potential energy in the six dimers.

The isotropic tension in an equibiaxial deformation (β = constant) can be expressed as (see [32])

$$\overline{T}={\frac{1}{{A}_{0}}\frac{\partial \mathrm{\Phi}}{\partial \alpha}|}_{\beta}-C/{A}^{2},$$

(4)

where the second term corresponds to a steric effect in the form suggested by [12] (see also [33]), *i.e.* the compression of the network causes a repulsive force due to interactions among components of the protein skeleton. The initial isotropic tension (prestress) at natural state is denoted by ${{\overline{T}}_{0}=\overline{T}|}_{\alpha =0,\beta =0}$. The linear area modulus *K _{s}* is given as

$${K}_{s}={\frac{\partial \overline{T}}{\partial \alpha}|}_{\alpha =0}.$$

(5)

Using the single-JC model we performed simulations of the quasi-static and ensemble-averaged response of a JC in in-plane deformations. Thermal fluctuations, which may have potential impact on local dynamic responses and cause phenomena such as mode switching [31], are not included in this simulation. The shear and linear area moduli of the protein skeleton are extracted by using Eqns. 3 to 5.

The lipid bilayer consists of a nearly incompressible fluid-like membrane. It therefore possesses large area modulus and a quite small shear modulus. For numerical stability we take *µ _{b}*, the shear stiffness of the lipid bilayer, to be a finite value 10

The mechanical response of the two continuum layers and their interactions are solved using the finite element method. Detailed description of this method is in the Appendix.

In the following sections we apply our model to study problems involving mechanical responses of RBC, including the resting shapes as well as cell deformations and skeleton-bilayer interactions in two canonical tests: optical tweezer stretching and micropipette aspiration.

It is well accepted that the resting shapes of RBC are related to the mechanical properties of the composite membrane consisting of the lipid bilayer and the skeleton network. However, the roles played by each of these components are unclear. On the one hand, by washing away the lipid bilayer using non-ionic surfactants, Svoboda *et al.* [34] showed that the remaining protein skeleton was no longer biconcave. On the other hand, in some diseases (*e.g.* hereditary elliptocytosis) a weakened skeleton network changes the cell to an elliptical shape. The implication is that both the skeleton and the lipid bilayer contribute to the resting shape.

Accordingly, an important factor is the relaxed (zero shear energy) reference state of the protein skeleton. The simplest choice is to use a spherical shape as the reference state. However, it was shown that under this assumption the cell rested at a cup shape [4]. To recover the biconcave shape the elasticity of the skeleton has to be significantly reduced [4]. Although the biconcave shape can be stabilized by using the biconcave shape itself [5] or an oblate ellipsoidal shape [2, 3] as the relaxed skeleton reference, experiments are needed to explain why red cells have those nonspherical relaxed skeletons.

Using our model, we simulate the resting shapes of RBC and compare with predictions made in existing studies. We first neglect the effects of the protein skeleton and study the dependence of the cell shape upon its internal volume *V*. By assuming that the spontaneous curvature *C*_{0} = 0, we plot the cell shapes at three different values of *V/V _{sphere}* in Fig. 3 (

Resting shapes of a RBC when the lipid bilayer is considered: (a) stomatocyte (*V/V*_{sphere}=0.59), (b) oblate (*V/V*_{sphere}=0.65) and (c) prolate (*V/V*_{sphere}=0.8).

We then take into account the protein skeleton and re-simulate the case with *p _{f}*=5.625nm,

Resting shapes of the RBC when both the lipid bilayer (light color) and the skeleton network (dark color) are considered: (a) cup shape (*c*_{0} = 0, *µ*_{s} = 11 pN/µm, *V/V*_{sphere}=0.65), (b) biconcave shape (*c*_{0} = 2.6, *µ*_{s} = 0.4 pN/µm, **...**

To explain the paradox about the shear stiffness of the skeleton, Li *et al.* [4] proposed that over large time scales the skeleton might be fluidic due to remodeling and possesses a much smaller shear stiffness than measured under finite deformation rates. This appears to be reasonable. However, it fails to explain why weakened skeleton causes dramatic shape change such as observed in elliptocytosis (at *C*_{0} = 0, the biconcave shape is the ground state for *µ _{s}* → 0 at

Alternatively, we find that with some positive spontaneous curvature the biconcave shape can be obtained with a small skeleton shear stiffness. For example, if we assume that the reduced spontaneous curvature *c*_{0} = 2.6 (*c*_{0} is the reduced spontaneous curvature defined as *c*_{0} = *C*_{0} · *R _{0}* and

Complete phase diagram with respect to shear stiffness and spontaneous curvature should be explored further, which is beyond the scope of this paper.

Through comparisons with benchmark results, the aforementioned simulations confirm the validity and accuracy of our models (especially the Level III FEM model). They have also demonstrated that additional studies, both theoretical and experimental, are required to explain the resting shape of RBC.

In the following we concentrate upon mechanically induced deformations of the cell in two canonical experiments, optical tweezer stretching and micropipette aspiration. In these simulations, we consider the bilayer with a spontaneous curvature *C*_{0} = 0 for simplicity and use the initial configuration as the relaxed reference state (shear-free state) of the skeleton. Following Discher *et al.* [12], two different scenarios are considered. In the first scenario the skeleton has no initial tension inside it (stress-free case with _{0} = 0). In this case *p _{f}*=5.625 nm,

We test the model through comparisons with experimental measurements of RBC deformability through optical tweezers. In this setup, an erythrocyte is stretched by two attached beads, whose motions are optically controlled by laser beams. In our model, the initial shape of the cell is biconcave and is mathematically depicted as [39],

$$z=\pm 0.5{R}_{0}\sqrt{1-\frac{{x}^{2}+{y}^{2}}{{R}_{0}^{2}}}\phantom{\rule{thinmathspace}{0ex}}[{C}_{1}+{C}_{2}\frac{{x}^{2}+{y}^{2}}{{R}_{0}^{2}}+{C}_{3}{\phantom{\rule{thinmathspace}{0ex}}\left(\frac{{x}^{2}+{y}^{2}}{{R}_{0}^{2}}\right)}^{2}]\phantom{\rule{thinmathspace}{0ex}},$$

(6)

where *C*_{1} = 0.21, *C*_{2} = 2.03, *C*_{3} = −1.12, and *R*_{0} = 3.91 µm. (*x, y, z*) is a Cartesian coordinate system with its origin located at the centroid of the undeformed cell. As reported by Dao *et al.* [33], the stretching force is applied by two silica beads, which are attached at the opposite ends of the cell over a small oval region with a diameter of 1~2 µm (Fig. 5a). In our simulation, this diameter is chosen to be 1.5 µm.

(a) Cell deformation stretched by the optical tweezers as predicted by our multiscale model (stress-free case). (b) Axial and transverse diameters of the cell as functions of the stretching force for the stress-free and the prestress cases as compared **...**

The force versus displacement curve is shown in Fig. 5b. As we see, the model predictions match well with the experimental measurements. The prestress cases are softer than the stress-free case because the persistence length is larger. The case with _{0} = −30 pN/µm is a little bit stiffer than the case with _{0} = −15 pN/µm. It also provides best comparison with the experiment among the three cases.

We simulate the canonical micropipette aspirations and obtain the correlation between the applied pressure and the aspiration length *L* (*i.e.* the length of the cell sucked into the pipette). In this simulation, a rigid cylindrical surface is employed to represent the pipette. The interaction between the outer layer of the RBC and the pipette is simulated by using a master-slave algorithm similar to the one used to study the bilayer-skeleton interaction (see the Appendix). As indicated in experiments [40], during the aspiration the membrane is usually separated from the pipette by a small gap of fluid so that the friction between them is insignificant and thus not considered in our model. We further simplify the fluid pressure distribution inside the pipette as a uniform pressure difference Δ*P* applied on the cap region of the lipid bilayer and a linear distribution along the aspiration length (the pressure difference equals zero at the entrance).

Following Waugh and Evans [38], we study a flaccid, unswollen cell aspired from the dimple region (Fig. 6a). The initial shape of the cell is depicted by Eq. (6). The cell is aspirated into a pipette with an inner radius *R _{p}* of 0.668 µm.

In Fig. 6b, it is seen that the normalized aspiration length *L/R _{p}* depends almost linearly on the two-dimensional pressure defined as Δ

To further test the capacity of our model, we use it to study areal deformation of the skeleton and compare with experimental measurements as well as predictions by the coarse-grained model. The areal variation of the protein skeleton can be denoted as the density ratio ρ/ρ_{0}, *i.e.* the density of skeleton-attached proteins ρ normalized by its value ρ_{0} in the undeformed state. This ratio is related to λ_{1} and λ_{2} by ρ/ρ_{0} = 1/(λ_{1}λ_{2}).

In Fig. 7a we plot the density profiles predicted by the current model for stress-free and prestress cases, and compare them with the reported result using the coarse-grained model for prestress case with _{0} = −15 pN/µm [12]. In Fig. 7b we plot how the skeleton density changes with *L/R _{p}* at the cap and the entrance regions predicted by the prestress case with

Fig. 7a demonstrates that for the same prestress level (_{0} = −15 pN/µm), the skeleton density of the cap region predicted by our model is lower than that by the coarse-grained model [12]. This is partially attributed to the fact that in the coarse-grained model the deformed shape of the cell is artificially assumed (perfect semi-sphere is assumed for the cap), but in our model the deformed shape is directly computed based on continuum mechanics. Although our computed deformed shape is just slightly different from the assumed shape by Discher *et al.*, additional constraints usually make the structure stiffer. The difference is also due to another fact that at *L/R _{p}* = 8 the number of junctional complexes is relatively small in the cap region where the surface is curved so that the difference between our continuum description and the discretized description in the coarse-grained model may be pronounced. Since in reality there are approximately 33,000 JCs in the cell while 18,434 JCs are considered in the coarse-grained model [12], the real density profile should be between our result and that of the coarse-grained model.

Fig. 7a also shows that the prestress level influences the skeleton density significantly. In fact, the accurate prestress level, even the type of the prestress (precompression or pretension), is still controversial [12, 34, 42]. In our following simulations we use _{0} = −30 pN/µm, which matches the experimental data as shown in Fig. 7b (and also in Figs. 5b and and6b6b).

As mentioned in the Appendix, in this case the initial shape of the cell is spherical and it is thus impossible to simultaneously conserve area and volume. Indeed, if we assume that the cap area is semi-spherical and the part outside of the pipette remains spherical during aspiration, it can be geometrically proven that when the surface area is conserved the volume loss in this case should be 27% as *L* = 8*R _{p}*, a value consistent with our numerical simulation.

For all simulations below we use the following parameters: *p _{f}*=11.118 nm,

Hereafter we concentrate upon the vertical contact force between the two layers during the micropipette aspiration of a cell from its biconcave state.

Fig. 8a displays the distributions of the contact pressure *p _{c}* between the two layers during a micropipette aspiration with an applied pressure of 256 pN/µm

To overcome the effect of internal tension *f* inside the protein skeleton, there must exist (a) negative (disassociation) contact force between the protein skeleton and the lipid bilayer in locations with convex shape, or (b) positive (association) contact **...**

In our simulation, negative (dissociation) contact pressure with peak value of about −130 pN/µm^{2} is recorded in the area close to the cap. To relate the contact pressure *p _{c}* to mechanical loading on the molecular structure of the protein skeleton, it is convenient to define a contact force per JC as

Sp unfolding may play an important role in determining the skeleton density, and subsequently *f _{jc}*, near the cap region. Unfolding will unstiffen the skeleton, causing larger area expansion and thus increasing

In the previous section, we studied the contact force between the two layers before the dissociation. Our model is also capable of exploring the subsequent response after skeleton/bilayer dissociation. When the negative pressure in micropipette aspiration is sufficiently high, a vesicle will be separated from the cell. A vesicle is an entity formed by part of the cell membrane which is deficient in skeleton proteins such as band 3 [16], suggesting a scenario wherein its formation is caused by the skeleton separating from the lipid bilayer. It was observed in experiments that before a vesicle is created a region of reduced caliber, *i.e.* a neck, was formed in the middle section between the cap and the entrance [16]. The location of the neck may appear to be at odds with the fact the skeleton-bilayer dissociation should first occur close to the cap, where the maximum dissociation force occurs.

To shed light on this phenomenon, we simulate the post-dissociation behavior by assuming that the two layers separate when the magnitude of the dissociation force per JC reaches a critical value * _{jc}*=20 pN for the case with

FEM simulation of the necking process before vesiculation (lipid bilayer in light color and skeleton in dark color)

We also found that under the aforementioned conditions (*i.e.* external pressure of 256 pN/µm^{2} and *F*_{1/2}=7.5 pN) the skeleton/bilayer dissociation occurs before Sp unfolding only if * _{jc}* <~4.0 pN, otherwise the separation occurs after Sp unfolding. Therefore, Sp unfolding may play a critical role in skeleton/bilayer separation. Additional experiments (especially the experiments to determine key parameters such as

It is necessary to point out that this simulation is qualitative rather than quantitative – to accurately capture the dynamic separation process an accurate evaluation of the damping matrix **C** (see the Appendix) is required. This matrix summarizes the effect of the viscosity of the lipid bilayer and the protein skeleton, the friction between these two, and the effects of the surrounding fluid. Studies in this direction are underway and will be included in future reports.

By using a multiscale modeling technique we have quantitatively studied the resting shape of erythrocyte, the quasi-static mechanical response prompted by mechanical loads, as well as the interaction force between the protein skeleton and the lipid bilayer. The primary purpose of this study is to build a framework for understanding erythrocyte mechanics upon which new knowledge from future experiments can be incorporated. Such a model analysis will also help us design new experiments to illustrate the exact mechanisms of mechanically induced membrane morphological changes, and to achieve quantitative predictions of the occurrence of remodeling (such as vesiculation) of normal or defected erythrocytes.

Compared with existing models of the RBC membrane [4, 12, 44, 45], our model has the following characteristics: (1) Our multiscale approach not only delivers accurate predictions of complete cell response (due to the involvement of the detailed molecular structure and responses at different levels), but also allows us to address physical mechanisms at different length scales and to correlate mechanical loads on the cell with detailed stress distributions within the composite structure. This model has predicted phenomena that had never been found by other models (*e.g.* bifurcation, mode switching, and stress-induced unstiffening due to unfolding) [31, 8]; (2) our model explicitly incorporates the local interactions between the skeleton and the bilayer, as well as the inter- and intra-molecule interactions inside the skeleton; (3) This model is inherently dynamic and capable of studying time-dependent responses at different length scales.

One possible application of this model is the prediction of the critical contact force between the lipid bilayer and the protein skeleton that triggers dissociation, leading to physiologically important phenomena such as vesiculation and tether formation. Among the existing experiments, observations of tether formations from erythrocytes inside flow channels (see for example [46]) or induced by cantilevers [6, 17] provide the only measure to quantitatively predict this critical dissociation process. In a typical flow channel setup, erythrocytes are allowed to sediment inside a channel consisting of two parallel plates. The substrate is coated with bovine serum albumin (BSA) so that most cells do not adhere to the bottom with large attaching areas. When external flows are introduced the cells deform while one (in some cases more than one) point remains attached to the substrate. Long membrane strands (tethers) may appear when the fluidic shear surpasses a threshold value (~1.5 dyn/cm^{2}, or 0.15 pN/µm^{2}) [46]. Although it is possible to duplicate deformations of the erythrocyte inside the flow channel and extract the contact force by using our model, these simulations suffer from uncertainties about the location and the detailed configuration of the contact area between the cell membrane and the substrate. In this context, the accuracy of a continuous complete cell model is also questionable due to the small size of the contact area between the cell and the substrate. A better approach will be to simulate a set of systematic micropipette aspirations by slowly increasing the negative pressure inside the pipette (to enforce the quasi-static condition) till vesiculation occurs. On the other hand, tangential forces between the skeleton and the bilayer associated with dynamic mechanisms may play an important role in the dissociation process (see for example, the tether formation experiments reported by Borghi and Brochard-Wyart [47]). To take this into account, it is critical to simulate the time-dependent viscoelastic behavior, rather than the quasi-static response, of the cell.

Other improvements of the model are also required. Rather than a rough view of total vertical loads on each JC achievable by our current model, a detailed distribution of mechanical loads among various skeleton-bilayer pinning points, as well as mechanical load on each inter-protein or protein-to-lipid connection are needed. A possible measure to achieve these is by concurrent multiscale models. For example, by modeling in molecular detail a piece of the skeleton while representing the rest of the skeleton as continuum shell, we will be able to relate lateral tension inside the protein skeleton to forces on the connection points. Besides Sp unfolding, other processes such as dissociation of Sp head-to-head linkages may also play important roles in determining structural responses of the cell in large deformations [9]. For comprehensive modeling these mechanisms should also be considered. In addition, more sophisticated thermal fluctuation model and fluid-skeleton-bilayer interaction model are required to accurately study close-coupling between the skeleton and the bilayer [48].

This work was supported by the National Heart, Lung, and Blood Institute under award number R01HL092793 .

Numerically, we employ the finite element method to solve the mechanical response of the continuum shells in the complete cell model. The development of this finite element model is summarized below.

For nonlinear solid mechanics problem in small length scales (*i.e.* the inertia effect is negligible), the momentum equation governing the mechanical response of either one of these layers can be stated with updated Lagrangian description as:

$$c\frac{\partial \mathbf{u}}{\partial t}-\nabla \xb7\mathbf{\Theta}=\mathbf{0},$$

(7)

with boundary conditions

$$\{\phantom{\rule{thinmathspace}{0ex}}\begin{array}{ccc}\hfill \text{\hspace{1em}}\mathbf{u}=\mathbf{g}\hfill & \text{on\hspace{1em}}\hfill & {\mathrm{\Gamma}}^{D}\text{\hspace{1em}\hspace{1em}}\hfill \\ \mathbf{\Theta}\xb7\mathbf{n}={\mathbf{f}}^{\mathit{\text{ext}}}\hfill & \hfill \text{on}\hfill & \hfill {\mathrm{\Gamma}}^{\mathit{\text{ext}},N}\hfill \\ \hfill \mathbf{\Theta}\xb7\mathbf{n}={\mathbf{f}}^{\mathit{\text{cnt}}}\hfill & \hfill \text{on}\hfill & \hfill {\mathrm{\Gamma}}^{\mathit{\text{cnt}},N}\hfill \end{array}$$

(8)

where *c* is the viscous damping coefficient. **u** is the displacement vector. is the spatial gradient operator. **Θ** is the Cauchy stress. Γ^{D} is the Dirichlet part of the boundary. Γ^{ext,N} is the Neumann part of the boundary where external forces are applied. Γ^{cnt,N} is the Neumann part of the boundary where contact interaction forces are applied. **n** is the normal vector to the surface. **f**^{ext} represents the external forces (*e.g.* the pressure force inside a micropipette). **f**^{cnt} is the contact force between the two layers, and between the outer layer and boundaries (*e.g.* micropipettes). **g** is the specified displacements on the boundary. Generally speaking, *c*, **Θ**, **f**^{ext}, **f**^{cnt} and **g** are functions of **u**, **v** and *t*, where **v** = **u**/*t* is the velocity.

We model both the outer layer (the lipid bilayer) and the inner layer (the protein skeleton) as congregations of shell elements. Although the *C*^{1}-conforming thin shell elements are more accurate [49] (also see [13, 14] for an application of this method to model the lipid bilayer), for simplicity and numerically robustness, in our current study we choose the *C*^{0} explicit Hughes-Liu elements [50]. In some special cases, *e.g.* aspiration of an erythrocyte into a cylindrical micropipette from its dimple region, the configuration is axisymmetric so that axisymmetric Hughes-Liu shell elements can be employed to reduce computational effort.

After finite element discretization (detailed formulation of this method can be found in [50] and axisymmetric formulation can be found in [51]), the governing equation (7) is re-expressed in matrix form as

$${\mathbf{\text{Cv}}}_{g}+{\mathbf{f}}_{g}^{\mathit{\text{int}}}={\mathbf{f}}_{g}^{\mathit{\text{ext}}}+{\mathbf{f}}_{g}^{\mathit{\text{cnt}}},$$

(9)

where **C** is the assembled damping matrix. In this study, we concentrate upon quasi-static cases so that **C** is a numerical parameter which affects the convergence *via* dynamic re-laxation but not the final result. Specifically, we use a lumped damping matrix based on row-sum technique [52]. **v**_{g} represents the assembled velocity vector. ${\mathbf{f}}_{g}^{\mathit{\text{ext}}}$ is the global external force vector and ${\mathbf{f}}_{g}^{\mathit{\text{cnt}}}$ is the global contact force vector. ${\mathbf{f}}_{g}^{\mathit{\text{int}}}$ is the global internal force vector related to constitutive equations.

Based upon Eq. (9), the configuration of the thin layer is updated *via* a forward Euler algorithm.

To close our finite element formulations, we relate the in-plane components of stress **Θ** with the deformation *via* the constitutive law provided by Evans and Skalak [32]. For the inner layer (skeleton network), we use

$${\mathrm{\Theta}}_{\mathit{\text{ij}}}h=\overline{T}{\delta}_{\mathit{\text{ij}}}+\frac{{\mu}_{s}}{{J}^{2}}\phantom{\rule{thinmathspace}{0ex}}({b}_{\mathit{\text{ij}}}-\frac{{b}_{11}+{b}_{22}}{2}{\delta}_{\mathit{\text{ij}}})\phantom{\rule{thinmathspace}{0ex}},\text{\hspace{1em}}i,j=1,2.$$

(10)

$J=\left|\begin{array}{cc}{F}_{11}\hfill & {F}_{12}\hfill \\ {F}_{21}\hfill & {F}_{22}\hfill \end{array}\right|$ is the in-plane Jacobian, where *F _{ij}* (

For the outer layer (lipid bilayer), since it is nearly incompressible, we use

$${\mathrm{\Theta}}_{\mathit{\text{ij}}}h={K}_{b}(J-1){\delta}_{\mathit{\text{ij}}}+\frac{{\mu}_{b}}{{J}^{2}}\phantom{\rule{thinmathspace}{0ex}}({b}_{\mathit{\text{ij}}}-\frac{{b}_{11}+{b}_{22}}{2}{\delta}_{\mathit{\text{ij}}})\phantom{\rule{thinmathspace}{0ex}},\text{\hspace{1em}}i,j=1,2,$$

(11)

where *K _{b}* = 5×10

Finally, we update the transverse shear stresses Θ_{23} and Θ_{31} by using the linear elastic model expressed as

$${\dot{\mathrm{\Theta}}}_{23}h=G{\dot{\gamma}}_{23},\text{}{\dot{\mathrm{\Theta}}}_{31}h=G{\dot{\gamma}}_{31},$$

(12)

where _{23} and _{31} are the relevant strain rates related to the local deformation. *G* is the transverse shear stiffness (different for the inner layer and the outer layer). In our model, since both the inner layer and the outer layer are very thin, the transverse shear should be negligibly small. We find that as long as the transverse shear stiffness is sufficiently large, its actual value has no influence on the results.

Considering a shell with thickness *h* described by Evans and Skalak [32], its bending stiffness *k _{c}* and area modulus

$${k}_{c}={\displaystyle {\int}_{-\frac{h}{2}}^{+\frac{h}{2}}{y}^{2}\frac{K}{h}\mathit{\text{dy}}=\frac{{\mathit{\text{Kh}}}^{2}}{12}.}$$

(13)

This formula is implied in the previous finite element formulation. For the lipid bilayer, we use *h* = 2.2 nm and *K* = *K _{b}* = 5 × 10

The bending stiffness of the protein skeleton is negligibly small (experiments show that the persistence length of a Sp is only around *O*(1) nm [24], indicating that it has a small bending stiffness). It is found that if its thickness *h* is chosen to be comparable to that of the bilayer (in our simulations we choose *h* = 2 nm), this characteristics can be duplicated in our model. For example, with small deformations a typical value of the area modulus of the skeleton is around 20 pN/µm, leading to a bending stiffness of 6.67 × 10^{−24} J, which is five orders of magnitude smaller than that of the lipid bilayer.

In its free state the lipid bilayer may also possess a spontaneous curvature *C*_{0}, referring to the naturally formed curvature of a free patch of bilayer [57]. To account for it, we replace the deformation gradient tensor **F**, defined as $\left[\begin{array}{ccc}{F}_{11}\hfill & {F}_{12}\hfill & {F}_{13}\hfill \\ {F}_{21}\hfill & {F}_{22}\hfill & {F}_{23}\hfill \\ {F}_{31}\hfill & {F}_{32}\hfill & {F}_{33}\hfill \end{array}\right]$, with , where is the deformation gradient caused by external loading and is the initial deformation gradient caused by the spontaneous curvature, which is defined initially before the simulation.

The interaction between the outer layer and the inner layer in the vertical direction is modeled as uniformly distributed penalty springs. Tangentially, the two layers are allowed to slide viscously against each other (herein the exact magnitude of the viscous drag in the tangential direction rendered by the mobility of the transmembrane proteins band 3 and glycophorin C is irrelevant since we consider quasi-static cases through dynamic relaxation). A master-slave penalty contact formulation is employed [58]. The outer layer is treated as the master surface and the nodes on it are called master nodes, whereas the inner layer is considered as the slave surface and the nodes on it are called slave nodes. All the slave nodes are projected to the master surface, and the distances between the slave nodes and the master surface are calculated. The contact stress depends linearly upon these distances with a penalty stiffness. The contact force, defined as the contact stress times the area of the element, is distributed to the slave nodes and the master nodes. The penalty stiffness is tested numerically so that it is sufficiently large to enforce the contact constraint accurately, and small enough so that it does not introduce numerical instability during time integration.

The internal volume of the cell is conserved through the following penalty algorithm

$$\mathrm{\Delta}{P}_{i}=-{k}_{v}\mathrm{\Delta}V,$$

(14)

where Δ*V* is the volume change of the cell and Δ*P _{i}* is the internal pressure against this change. Δ

1. Hoffman J. Blood Cells Mol. Dis. 2001;27:57. [PubMed]

2. Lim GHW, Wortis M, Mukhopadhyay R. Proc. Natl. Acad. Sci. 2002;99:16767.

3. Khairy K, Foo J, Howard J. Celluar and Molecular Bioengr. 2008;1:173.

4. Li J, Dao M, Lim CT, Suresh S. Biophys. J. 2005;88:3707. 2005. [PubMed]

5. Zarda PR, Chien S, Skalak R. J. Biomechanics. 1977;10:211. [PubMed]

6. Hwang WC, Waugh RE. Biophys. J. 1997;72:2669. [PubMed]

7. Hochmuth R, Marcus WD. Biophys. J. 2002;82:2964. [PubMed]

8. Zhu Q, R. Asaro R. Biophys. J. 2008;94:2529. [PubMed]

9. Li J, Lykotrafitis G, Dao M, Suresh S. Proc. Natl. Acad. Sci. 2007;104(12):4937. [PubMed]

10. Mohandas N, Gallagher PG. Blood. 2008;112:3939. [PubMed]

11. Deutsch S, Tarbell JM, Manning KB, Rosenberg G, Fontaine AA. Annu. Rev. Fluid Mech. 2006;38:656.

12. Discher DE, Boal DH, Boey SK. Biophys. J. 1998;75:1584. [PubMed]

13. Feng F, Klug WS. J .Comput. Phys. 2006;220:394.

14. Ma L, Klug WS. J .Comput. Phys. 2008;227:5816.

15. Berk DA, Hochmuth RM. Biophys. J. 1992;61:9. [PubMed]

16. Knowles DW, Tilley L, Mohandas N, Chasis JA. Proc. Natl. Acad. Sci. USA. 1997;94:12969. [PubMed]

17. Waugh RE, R.G. Bauserman RG. Annals Biomed. Engr. 1995;23:308. [PubMed]

18. Butler J, Mohandas N, Waugh RE. Biophy. J. 2008;95:1826. [PubMed]

19. Mohandas N, Evans EA. Annu. Rev. Biophys. Biomol. Struct. 1994;23:787. [PubMed]

20. DeSilva TM, Peng K-C, Speicher KD, Speicher DW. Biochemistry. 1992;31:10872. [PubMed]

21. An X, Lecomte MC, Chasis JA, Mohandas N, Gratzer W. J. Bio. Chem. 2002;277(35):31796. [PubMed]

22. Anong WA, Weis TL, Low PS. J. Biol. Chem. 2006;281:22360. [PubMed]

23. An X, Takakuwa Y, Nunomura W, Manno S, Mohandas N. J. Biol. Chem. 1996;271:33187. [PubMed]

24. Rief M, Pascual J, Saraste M, Gaub HE. J. Mol. Biol. 1999;286:553. [PubMed]

25. Lee JCM, Discher DE. Biophys. J. 2001;81:3178. [PubMed]

26. Law R, Carl P, Harper S, Dalhaimer P, Speicher DW, Discher DE. Biophys. J. 2003;84:533. [PubMed]

27. Paramore S, Voth GA. Biophys. J. 2006;91:3436. [PubMed]

28. Rief M, Gautel M, Oesterhelt F, Fernandez JM, Gaub HE. Science. 1997;276:1109. [PubMed]

29. Sung LA, Vera C. Ann. Biomed. Eng. 2003;31:1314. [PubMed]

30. Vera C, Skelton R, Bossens F, Sung LA. Ann. Biomed. Eng. 2005;33:1387. [PubMed]

31. Zhu Q, Vera C, Asaro R, Sche P, Sung LA. Biophys. J. 2007;93:386. [PubMed]

32. Evans EA, Skalak R. Mechanics and thermodynamics of biomembranes. Boca Raton, FL: CRC Press; 1980.

33. Dao M, Li J, Suresh S. Mater. Sci. Engr. 2006;C26:1232.

34. Svoboda K, Schmidt CF, Branton D, Block SM. Biophys. J. 1992;63:784. [PubMed]

35. Seifert U, Berndl K, Lipowsky R. Physical Review. A. 1991;44(2):1182. [PubMed]

36. Peterson MA, Strey H, Sackmann E. J. Phys. II (France) 1992;2:1273.

37. Boal DH, Seifert U, Zilker A. Phys. Rev. Lett. 1992;69:3405. [PubMed]

38. Waugh RE, Evans EA. Biophys. J. 1979;26:11532.

39. Evans EA, Fung YC. Microvasc. Res. 1972;4:335. [PubMed]

40. Discher DE, Mohandas N. Biophys. J. 1996;71:1680. [PubMed]

41. Discher DE, Mohandas N, Evans EA. Science. 1994;266:1032. [PubMed]

42. Kozlov MM, Markin VS. Contemporary Problems of Biomechanics. CRC Press; 1990. p. 11.

43. Waugh RE, Hochmuth RM. Biophys. J. 1987;52:391. [PubMed]

44. Mukhopadhyay R, Lim GHW, Wortis M. Biophys. J. 2002;82:1756. [PubMed]

45. Svetina S, Kuzman D, Waugh RE, Ziherl P, Zêkŝ B. Bioelectrochemistry. 2004;62:107113. [PubMed]

46. Hochmuth RM, Mohandas N, Blackshear JR PL. Biophys. J. 1973;13:747. [PubMed]

47. Borghi N, Brochard-Wyart F. Biophys. J. 2007;93:1369. [PubMed]

48. Heinrich V, Ritchie K, Mohandas N, Evans E. Biophys. J. 2001;81(3):1452. [PubMed]

49. Cirak F, Ortiz M, Schröder P. Int. J. Num. Meths. Engr. 2000;47:2039.

50. Hughes TJR, Liu WK. Comput. Meths. Appl. Mech. Engr. 1981;26:331.

51. Hughes TJR, Liu WK. Comput. Meths. Appl. Mech. Engr. 1981;27:167.

52. Belytschko T, Liu W, Moran B. Nonlinear Finite Elements for Continua and Structures. New York: John Wiley; 2000.

53. Asaro RJ, Lubarda V. Mechanics of Solids and Materials. Cambridge Univ. Press; 2006.

54. Lin LC, Brown FL. Phys. Rev. E. 2005;72:011910. [PubMed]

55. Marcelli G, Parker KH, Winlove CP. Biophys. J. 2005;89:2473. [PubMed]

56. Bo L, Waugh RE. Biophys. J. 1989;55:509. [PubMed]

57. Helfrich W. Z. Naturforsch. 1973;C28:693. [PubMed]

58. Malone JG, Johnson NL. Int. J. Num. Meths. Engr. 1994;37(4):559.

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |