INTRODUCTION

TOP ABSTRACT INTRODUCTION THE MODEL THEORY RESULTS DISCUSSION APPENDIX A REFERENCES

The normal resting shape of the human red blood cell (RBC or erythrocyte) is a flattened biconcave disc (discocyte) ~8 µm in diameter. Treatment of erythrocytes in vitro with a variety of amphipathic agents is known to transform this shape systematically and reversibly into various other shapes such as echinocytes (crenated shapes) and stomatocytes (cup-like shapes), which are further subdivided into classes labeled I, II, and III (Brecher and Bessis, 1972; Bessis, 1973; Chailley et al., 1973; Mohandas and Feo, 1975). In particular, the echinocyte III is a more-or-less spherical body covered evenly with 25 to 50 rounded protuberances, which we shall call spicules.

These shape transformationsfrom the discocyte to the stomatocyte on one side and from the discocyte to the echinocyte on the otherhave been studied experimentally for more than 50 years. Understanding these shapes and shape transformations is a classic problem in cell biology; over the past three decades it has also attracted the attention of physicists. What makes this problem so intriguing is that the structure of the red cell is remarkably simple. To a good approximation it is simply a bag of fluid, a concentrated solution of hemoglobin surrounded by a thin macroscopically homogeneous membrane. Thus, the diverse resting shapes that it exhibits can be thought of as (locally) stable mechanical equilibrium shapes of the membrane. This membrane is composite (Steck, 1989; Alberts et al., 1994). On the outside is the plasma membrane, a self-assembled fluid bilayer with a thickness of ~4 nm, composed of a complex mix of phospholipids, cholesterol, and dissolved proteins. Lipids such as phosphatidylcholine, sphingomyelin, and glycophospholipids, which are neutral at physiological pH, are common in the outer leaflet, whereas phosphatidylserine and phosphatidylethanolamine predominate in the inner leaflet (Gennis, 1989; Alberts et al., 1994). Due to the negative charge of phosphatidylserine, located in the inner leaflet, there is a significant difference in charge between the two leaflets of the bilayer. Inside the plasma membrane but linked to it by protein anchors is a thin, two-dimensionally cross-linked protein cytoskeleton, the membrane skeleton (Steck, 1989; Bennett, 1990). The skeleton is an hexagonally linked net, each unit of which is a filamentous spectrin tetramer with an extended length of ~200 nm; the end-to-end distance is reduced appreciably in the relaxed state by thermal fluctuations (Boal, 1994). The spectrin tetramer is negatively charged at physiological pH. In the red cell membrane, the distance between the vertices of the net is ~76 nm, whereas the offset between the plasma membrane and the spectrin network is 30 to 50 nm (Byers and Branton, 1985; Liu et al., 1987). Functionally, the plasma membrane serves as an osmotic barrier, passing water with relative ease but controlling, via a system of pumps and channels, the flow of ions and larger solute molecules. The membrane skeleton has the role of supporting and toughening the plasma membrane, which would otherwise disintegrate in circulatory shear flow.

Although the RBC membrane is certainly heterogeneous at the length scale of individual lipid molecules, lipid patches or the basic spectrin tetramer, on scales longer than 100 nm it is reasonably homogeneous in its properties, and it may make sense to treat it as a mechanical continuum. In this spirit and following early work by Canham (1970), Helfrich (Helfrich, 1973; Deuling and Helfrich, 1976), and others, we imagine replacing the (fluid-phase) plasma membrane by an ideal uniform two-dimensional surface characterized by a bending rigidity _b, taken for the RBC to be (Waugh and Bauserman, 1995; Strey et al., 1995) 2.0 × 10¹⁹ J (roughly 50 k_BT_room, where k_B is Boltzmann's constant), a spontaneous curvature C₀, which recognizes the outside-inside asymmetry of the leaflet composition, and an area-difference-elasticity (ADE) term (Helfrich, 1973; Evans, 1974; Svetina et al., 1985; Boi et al., 1992; Wiese et al., 1992; Miao et al., 1994), which reflects the fact that a difference in relaxed area, A₀, between the outer and inner leaflets produces a "bilayer couple" (Sheetz and Singer, 1974) tending to make the membrane become convex with the larger-area side on the outside of the convexity (like a bimetallic strip). We assign to the modulus associated with the ADE term a value of (Waugh and Bauserman, 1995) 1.27 × 10¹⁹ J. These terms taken together constitute what is called the ADE model of the plasma membrane. In addition, we represent the RBC membrane skeleton as a two-dimensional elastic continuum characterized by a bulk stretching modulus K and a shear modulus µ with µ  2.5 × 10⁶ J/m² and K  3 µ (see Discussion). These two mechanical subsystems, constrained to the same mathematical surface, will constitute our model of the RBC shape at length scales above 100 nm (for full details, see The Model).

The basic hypothesis of the mechanical approach is that the observed RBC shapes must be shapes that minimize the energy subject to appropriate constraints on volume and area, i.e., that all the observed RBC shapes must emerge as local energy minima of this model and that the observed shape transformations must come about as a response of the minimum-energy shapes to changes in the "control parameters," which characterize the model. There are relatively few such control parameters. In addition to the known area and volume of the RBC, we include in the list the three mechanical moduli noted above, the spontaneous curvature C₀, the relaxed area difference A₀ between the bilayer leaflets, and parameters describing the effective relaxed size and shape of the membrane skeleton. If all these control parameters were known or easily measurable, we could predict the stable RBC shapes using continuum mechanics. The fact that they are not restricts us, as we shall see, to somewhat more generic predictions.

This program has already been carried out, albeit in somewhat restricted form, with respect to the discocyte and stomatocyte shapes. The shapes and shape transitions of laboratory-prepared fluid-phase phospholipid vesicles (lacking any skeleton) have been successfully described by the ADE model (Miao et al., 1994; Döbereiner et al., 1997). In particular, plausible discocytic and stomatocytic shapes do occur as energy-minimizing shapes of the ADE model. Furthermore, the bilayer-couple mechanism (Sheetz and Singer, 1974) accounts qualitatively for many of the chemically induced discocyte-stomatocyte transitions observed in the laboratory, in that stomatocytogenic agents tend to segregate to the inner leaf of the bilayer, thus making A₀ negative and favoring the invaginated stomatocyte shape (see more below). There has been discussion in the literature of the way in which introducing the elastic energy of the membrane skeleton modifies or improves the detailed prediction for the shape of the normal resting discocyte (Zarda et al., 1977; Evans and Skalak, 1980). However, most treatments in the literature tend to regard the cytoskeletal elastic energy as providing only perturbative modification of the basic shapes, which emerge from the pure-bending-energy model. We shall argue that this point of view fails completely in describing echinocyte shapes.

Until recently, it was a major difficulty for a fully mechanical picture of the RBC shapes that echinocytic shapes have not been found in the catalog of minimum-energy shapes of the ADE model; however, several authors have now pointed out what we believe is the correct resolution of this problem (Waugh, 1996; Igli, 1997; Igli et al., 1998a,b; Wortis, 1998). A pure-lipid vesicle with a sufficiently positive spontaneous curvature (or, equivalently, sufficiently large positive A₀) adopts a vesiculated or budded shape, in which one or more quasispherical buds appears on the outside of the main body of the vesicle, attached to it via a narrow neck or necks (Fig. 1). Such necks are allowed as low-energy structures for pure-lipid vesicles, because the bending energy depends on the sum of the two principal curvatures (see The Model), so that large curvatures of opposite sign (characteristic of necks) can still lead to a small mean curvature and, thus, to low energy (Fourcade et al., 1994). But (and this is the crucial point) small necks involve large elastic strains in the membrane skeleton and are, thus, inhibited by elastic energy cost. It is not surprising that they are not seen for undamaged RBCs. What our calculations show is that the formation of echinocytes is driven by positive C₀ (equivalently, positive A₀) and that the spiculated shapes seen in experiments can arise from a competition between the bending energy of the plasma membrane, which by itself would promote budding, and the stretching and shear elastic energies of the membrane skeleton, which prevent it. To make this more precise, we note that the ratio of the bending modulus to the shear modulus (or the stretching modulus) defines an elastic length scale,

(1)

which yields _el  0.28 µm. We shall see below that this number sets the scale of spicule size and, thus, fixes the spicule density of the fully developed echinocyte (echinocyte III).

View larger version (9K): [in this window] [in a new window]   FIGURE 1   (a) Budded or vesiculated shape, which would be a low-energy configuration of a lipid membrane without cytoskeleton at sufficiently positive spontaneous curvature. (b) Spiculated shape into which (a) is transformed by the elastic-energy cost of cytoskeletal deformations in the high-shear neck regions.

The success of this approach strengthens the hypothesis that RBC shapes can be understood on the basis of continuum mechanics; however, in themselves, such calculations cannot answer the question, "Why does a particular shape occur under given experimental conditions?" Such a question has two parts. Why do the control parameters have the values they do? How does this set of control parameters lead to the observed shape? The second part is the proper domain of continuum mechanics. The first part is primarily biological or biochemical in content and revolves around the mechanisms that control the lipid composition of the leaflets of the plasma membrane, local electrostatic effects, the composition and environment of the spectrin network, and its coupling to the plasma membrane, and any other determinants of C₀, A₀, the elastic moduli, and other control parameters. Hopefully, this separation of the question is useful.

The plan of this paper is as follows. In the remainder of this section, we briefly review some history of echinocyte shapes and shape calculations. The Model section presents the details of the continuum-mechanics model. Theory section discusses the theory required to solve the model numerically for echinocyte shapes, including a full account of the scaling argument suggested above. The Results section presents our results for the predicted echinocyte shapes and spicule density. The Discussion section contains discussion of our results and further speculations.

The discocyte-to-echinocyte transformation was precisely identified by Ponder (1948, 1955) and has since been studied extensively by many other authors, as reviewed, for example, by Bessis (1973), Lange et al. (1982), and Steck (1989). Although this transformation can be driven in many ways, the final spiculated shapes produced are apparently the same (Smith et al., 1982; Mohandas and Feo, 1975). This quasi-universal behavior is evidence for the kind of mechanism suggested above, in which quite diverse biochemical processes may set (a few) important control parameters to the same or similar values. In the same spirit, we note that RBC ghosts (formed by hemolysis) behave quite similarly to intact red cells (Lange et al., 1982).

The simplest picture of this typeconsistent with some (but not necessarily all) experimental findingsis the so-called bilayer-couple hypothesis of Sheetz and Singer (1974), which focuses attention on the control parameter A₀. Thus, adding exogenous phospholipids to the exterior solution is observed to promote echinocyte formation (Ferrell et al., 1985; Christiansson et al., 1985). This is explained in the bilayer-couple picture by arguing that such addition is expected to produce an area increase in the outer leaflet only, as phospholipids do not readily flip-flop from one leaflet to the other. This increases the area difference A₀ and is expected to promote echinocyte formation. Similarly, cholesterol does equilibrate across the bilayer but apparently prefers the outer leaflet, presumably because of the particular lipids there present. Thus, adding cholesterol will tend to increase A₀ and promote echinocytosis, whereas depleting cholesterol will tend to decrease A₀ and to promote stomatocytosis (Lange and Slayton, 1982). More generally, anionic amphipaths all tend to produce echinocytes, whereas cationic amphipaths tend to produce stomatocytes (Deuticke, 1968; Weed and Chailley, 1973; Mohandas and Feo, 1975; Smith et al., 1982). The bilayer-couple explanation of these observations is that, when incorporated into the plasma membrane, the cationic compounds segregate preferentially to the inner leaflet and the anionic compounds, to the outer leaflet because of the predominantly negative charge of the inner-leaflet lipids. Outer-leaflet segregation increases A₀, thus promoting echinocytosis; conversely, inner-leaflet segregation promotes stomatocytosis. Time-dependent effects have even been seen, where a molecule initially intercalates into the outer leaflet, causing echinocytosis, but then slowly migrates to the inner leaflet, following the electrostatic gradient, causing return to the discocyte and subsequent stomatocytosis (Isomaa et al., 1987).

A similar hypothesis is made to explain the observed tendency of ghosts to become echinocytes at high salt concentrations but stomatocytes at low salt concentrations (Lange et al., 1982). The argument relies on electrostatics. Cationic species neutralize the negative lipid charges on the inner leaflet and salt decreases the Debye length, more effectively screening such charges. Both effects lead to a contraction of the inner leaflet and consequent increase in A₀, promoting crenation.

Other observations can be linked to the bilayer-couple effect but somewhat less directly. Thus, it is found that hemolyzed echinocytes produce discocytic ghosts (Lange et al., 1982) and that electroporation suppresses shape changes (Schwarz et al., 1999a,b). To explain these effects, it is argued that hemolysis and electroporation both involve perforation of the membrane and resulting in contact via the pore surface between the lipids in the inner and outer leaflets. This provides a mechanism for lipid transport that may reduce the magnitude of A₀ by partial symmetrization of the lipid composition of the leaflets, thus (equivalently) reducing C₀. Similarly, it is known that ATP depletion drives discocytes to crenate (Nakao et al., 1960, 1961, 1962; Backman, 1986). The hypothesis here is that ATP is required in some way to maintain the lipid asymmetry of the leaflets (related to C₀), perhaps for the operation of ATP-driven translocases (Steck, 1989).

Some effects are not readily explained by the bilayer-couple mechanism. Weed and Chailley (1972) and Gedde et al. (1995, 1997a,b, 1999) report that RBC shape can be controlled experimentally by varying the external pH with high pH promoting echinocyte formation and low pH, stomatocyte formation (the effect of proximity to a glass surface in promoting echinocytosis (Furchgott and Ponder, 1940) is probably related to this effect). The mechanism for this effect does not seem to be well established. Some authors argue for a mechanism involving membrane proteins. Band 3 has been specifically implicated (Gimsa and Ried, 1995; Jay, 1996; Gimsa, 1998; Hägerstrand et al., 2000). For example, it is known that high pH induces dissociation of ankyrin and band 3, which plays a role in the linkage of the membrane skeleton to the bilayer (Low et al., 1991). Indeed, it was shown that in vitro the membrane skeleton expands at high pH and contracts at low pH (Elgsaeter et al., 1986; Stokke et al., 1986). It has recently been proposed that pH change may induce conformational transformation in band 3 protein, leading to a change in A₀. Similarly, Wong (1994, 1999) has proposed that the folding of the spectrin in the cytoskeletal net is controlled by a conformational change of the band 3 protein. Such mechanisms might shift the emphasis of shape determination from the dominantly lipid-related control parameters, C₀ and A₀, to the elastic parameters and the relaxed shape of the membrane skeleton. We shall comment briefly on these issues in the Discussion section.

Numerous shape calculations have previously been done using variants of the model we describe in The Model. Early calculations (Fung and Tong, 1968; Zarda et al., 1977; Evans and Skalak, 1980) omitted the ADE term ( = 0), set C₀ = 0, and concentrated entirely on the shape of the resting discocyte and the modifications caused by osmotic swelling. They solved numerically the full mechanical shape equations. Landman (1984) treated a viscoelastic shell surrounding a viscous droplet and modeled the formation of protrusions as sudden local addition of shell material. Waugh (1996) studied a full ADE model with local shear elasticity (but treating the compressibility modulus as infinite) and studied the instability of a flat membrane to formation of a local "bump" (axisymmetric spicule) for positive values of C₀ or A₀ or both. He assumed a specific spicule-shape parametrization, focussed on a single spicule, calculated its energy, and in this way was able to discuss generically the conditions under which spicule formation might become energetically favorable. Most recently, Igli and collaborators (Igli, 1997; Igli et al., 1998a,b) have used the same model as Waugh (1996) but parametrize the spicule more simply, using a hemispherical cap on a cylindrical body and joining this to a spherical surface via a base shaped like a toroidal section. They correctly identify the emergence of echinocyte shapes as a competition between the bending energy and the skeletal elasticity. They model the echinocyte as a collection of n_s spicules attached to a spherical body and then determine the radius of the sphere, the parameters of the spicules, and the preferred number of spicules by minimizing the mechanical energy subject to constraints on cell area and volume. In this way, they are able to estimate the expected number of spicules and their shape as a function of the relaxed area difference A₀.

The present paper builds on the work of Waugh (1996) and Igli (1997, 1998a,b). We extend the model to include the area modulus K of the membrane skeleton. We calculate spicule shape for the first time from the Euler-Lagrange equations (these calculations are closely related to the earlier work of Zarda et al. (1977)). Although we continue with Igli et al. (1998a,b) to treat the echinocyte as a sphere decorated with spicules, we deal seriously (although still approximately) with the mechanical boundary conditions that must be met where the spicule joins the sphere.

	THE MODEL

TOP ABSTRACT INTRODUCTION THE MODEL THEORY RESULTS DISCUSSION APPENDIX A REFERENCES

The central assumption of our work is that the red blood cell assumes a shape that minimizes its overall membrane energy subject to the appropriate constraints. The RBC membrane is a composite of the plasma membrane and the cytoskeletal network; correspondingly, we take the membrane energy to be a sum of two terms,

(2)

the bending energy,

(3)

associated with the bilayer and an elastic energy of stretching and shear,

(4)

associated with the skeleton. In doing this, we treat the plasma membrane as incompressible. An estimate of its stretching modulus is K ~ 10¹ J/m², comparable with that of a pure-phospholipid bilayer and some four orders of magnitude larger than that of the skeleton. Correspondingly, we ignore any bending rigidity associated with the isolated membrane skeleton (the dimensional estimate  ~ K × (thickness)² leaves it two orders of magnitude smaller than that of the bilayer).

Eq. 3 for the plasma membrane's contribution to the overall membrane energy is the now-standard ADE Hamiltonian (Svetina et al., 1985; Boi et al., 1992; Wiese et al., 1992; Miao et al., 1994). The first term was originally proposed by Helfrich (Helfrich, 1973; Deuling and Helfrich, 1976). C₁ and C₂ are the two local principal curvatures, C₀ is the spontaneous curvature, and the integral is over the membrane surface. The two leaflets of a closed bilayer of fixed interleaflet separation D are required by geometry to differ in area by an amount,

(5)

In calculations, we shall take D  3 nm, corresponding to the distance between the neutral surfaces of the leaflets. The second term in Eq. 3 is the so-called area-difference-elasticity energy and represents the cost in stretching (or compressional) energy of the individual leaflets necessary to force the change from the relaxed area difference A₀ so as to conform to this geometric requirement. This effect occurs because a strong hydrophobic barrier prevents lipids in one leaflet from "flip-flopping" passively to the other on the time scales of mechanical shape changes. A is the area of the plasma membrane. The modulus associated with the ADE term is generically of the same scale as _b. Note, finally, that we can substitute Eq. 5 to rewrite Eq. 3,

(6)

in which

(7)

 = /_b, and the constant term is shape independent. This shows that, in determining minimal shapes at given values of the control parameters, C₀ and A₀ do not enter independently but only in the form of an effective spontaneous curvature ₀ or, equivalently, an effective relaxed area difference,

(8)

which we shall often quote in the dimensionless, reduced form

(9)

Note for future reference that increasing a₀ by 0.01 is equivalent to increasing C₀ by 6.7 µm¹.

Eq. 4 measures the elastic-energy cost of the spectrin network. It depends both on the relaxed shape of the membrane skeleton and on the way it is actually distributed over the membrane surface. (The notion of a relaxed shape is, of course, somewhat nominal, because removing the skeleton from the plasma membranewhich can certainly be done (Sheetz, 1979)would radically change its local biochemical environment, in a way which would modify its shape and elastic constants). In this redistribution, each element of the network will be stretched or compressed. The quantities ₁ and ₂ are the local principal extension ratios of each membrane element (Evans and Skalak, 1980). K and µ are the (two-dimensional) moduli for stretch and shear, respectively, as introduced in the Introduction, and the integrals are over the undeformed shape. In writing Eq. 4, we are assuming, as appears to be typical, that the membrane skeleton does not plastically deform in the course of the experimental shape transformation being described. Note that Eq. 4 makes a particular choice of terms in the elasticity beyond those quadratic in the weak-deformation parameters, _i = _i  1. For the large elastic strains, which are present at narrow necks and for small spicules, these nonlinearities do play a role. As far as we are aware, RBC elasticity has not been measured well enough to produce a clear preference for a particular form of the nonlinearities. Various authors have used Eq. 4 for the elastic energy in other contexts with apparent success; however, alternative forms of the elasticity have also been proposed for dealing with problems involving large local deformations (Evans and Skalak, 1980; Discher et al., 1994). We shall make some additional comments in the Discussion.

One attractive feature of Eq. 4 is that the effect of changing the size of the relaxed membrane skeleton by a pure, uniform dilation is particularly simple. Suppose that the skeleton is decreased in linear scale by a factor b. This means that undeformed areas are reduced by a factor 1/b², whereas the extension ratios _1,2 are each increased by a factor b. Note that the integrand of the shear term in Eq. 4 is invariant under this change but the factor ₁₂ in the stretching integrand increases by b². Only the quadratic part of the stretching term influences the membrane shape because dA₀₁₂ is just the membrane area A, which is fixed by the incompressibility of the plasma membrane. It follows that the effect of decreasing the size of the skeleton is completely equivalent to keeping the size fixed but, instead, changing the elastic moduli according to

(10)

(11)

i.e., this "prestress" in the membrane skeleton is completely equivalent to no prestress, a harder stretching modulus, and a softer shear modulus. It is not known with certainty what "prestress" exists in the typical RBC skeleton. Computer simulation has suggested that the relaxed skeleton may be as much as 10 to 20% smaller in area than the plasma membrane (Boal, 1994). On the other hand, the experiment by Svoboda et al. (1992) shows that isolated skeletons are expanded. In our calculations, we shall assume zero net prestress (relaxed skeletal area equal to RBC area, i.e.,  A/A₀ = 1). Any deviation from this would result in modified, effective elastic constants according to Eqs. 10 and 11.

In our calculations, we will identify shapes that are minima of the (free) energy defined by Eqs. 2 through 4 subject to constraints of fixed membrane area A (bilayer incompressibility) and volume V (as the RBC volume is typically set by osmotic balance). Constrained minimization is achieved variationally by introducing the functional,

(12)

in which  and P are Lagrange multipliers used to enforce the surface-area and volume constraints. P is the pressure difference across the bilayer, whereas  has the dimension of a surface tension. Making  stationary with respect to variations of membrane shape and cytoskeletal distribution leads to a set of coupled Euler-Lagrange equations. These equations can then be solved to give the shapes of mechanical equilibrium. In principle, such shapes are expected to correspond to observed equilibrium shapes at temperature T = 0 only; however, because the energy scale _b is large on the scale of k_BT, thermal fluctuations are generally negligible and play an important role only for "soft" modes, especially near instabilities (Wortis et al., 1997). The Euler-Lagrange equations are in general nonlinear and can have multiple solutions. The lowest-energy solution for given A and V is automatically stable to small fluctuations; higher-energy solutions must be tested for stability. All local-minimum solutions are candidates for stable observable shapes, except in exceptional cases (near instability) where energy barriers become comparable to k_BT.

The stable shapes produced by this process depend on the control parameters: the geometrical parameters, A and V; the determinants of curvature, C₀ and A₀, which only enter in the combination found in Eq. 7 or equivalently, Eqs. 8 and 9, the moduli, , , K, and µ; plus, finally, any parameters required to characterize the relaxed shape of the membrane skeleton. The area and volume of a typical RBC we take as (Bessis, 1973) A = 140 µm² and V = 90 µm³. The moduli, given in the Introduction, are less well determined by experiment (we shall have more to say on this point in the Discussion Section). This leaves as unknown control parameters the curvature variables and the cytoskeletal parameters.

The Euler-Lagrange variational equations derived from (12) are not numerically tractable except in the case of axisymmetry, which clearly does not apply to echinocytic shapes. In this paper, we aim to treat only fully developed echinocytic shapes (echinocyte III). For this case, individual spicules are identical and axisymmetric in shape to a good approximation and, in addition, the central body, which they decorate is approximately spherical with radius R₀. The observed distribution of spicules is rather regular, and we shall approximate the local spicule packing as a triangular array (except for special numbers of spicules, there must, of course, be some defects), which will look increasingly like an hexagonal close-packed structure, as the number of spicules, n_s, becomes large.

The base of each individual spicule, where it meets the sphere tangentially, is a circular contour L of radius r_L. If _L is the angle subtended by L at the center of the sphere, then

(13)

Because of the close-packed structure, the circles L meet tangentially. We explain in the Theory section how to derive appropriate boundary conditions where the spicules meet the sphere. Solving the axisymmetric Euler-Lagrange equations then determines a spicule shape, including an individual spicule volume V_s and area A_s. In terms of these variables, the overall area and volume of the RBC are taken to be

(14)

and

(15)

with the spicule number

(16)

In writing these relations, we have assumed that 10% of the spherical surface is not covered by the circular spicule bases (close packing on a flat surface would give 9.3%; curvature effects and packing defects both increase this number). The spicule volume V_s is calculated with respect to a plane through L, and Eq. 15 neglects curvature corrections. These approximations are good when the number of spicules is large, as it will turn out to be for the echinocyte shapes.

	THEORY

TOP ABSTRACT INTRODUCTION THE MODEL THEORY RESULTS DISCUSSION APPENDIX A REFERENCES

Membrane mechanics

In our treatment, spicules are assumed to be axisymmetric, and individual spicules are joined to the central body along the contour L, as illustrated in Fig. 2. Thus, our calculation involves finding a family of energy-minimizing axisymmetric spicule shapes and selecting from that family those shapes consistent with appropriate mechanical boundary conditions applied along L.

View larger version (18K): [in this window] [in a new window]   FIGURE 2   Parametrization of a single-spicule shape, showing the boundary L where the spicule joins the spherical central body in constructing the full echinocyte shape. The bases of two adjacent spicules are sketched in to illustrate the points A of spicule tangency and the point B, which is symmetrically situated between three neighboring spicules.

Parametrization of the axisymmetric spicule shape is illustrated in Fig. 2. The variables z and r measure distances along and perpendicular to the symmetry axis, respectively, whereas s measures the arclength from the pole. The function z(r) determines the shape.  is the angle between the local normal and the symmetry axis; C_m(r), C_p(r) are the principal curvatures,

(17)

In calculating the spicule shape, we shall assume that the relaxed cytoskeleton is locally flat, a good approximation as long as the number of spicules is large, so that the spicule size is small compared with R₀. Thus, in forming each spicule, a flat circular patch of relaxed cytoskeleton must be elastically deformed to fit the spicule contour. We assume that this deformation is axisymmetric, so that the center of the patch remains at the apex of the spicule, and each point of the patch at relaxed radius s₀ maps to a point at arclength s and radius r on the spicule contour. The principal extension ratios can then be written,

(18)

When these expressions for ₁ and ₂ are inserted into Eq. 4 we obtain an explicit form for the elastic energy in terms of the functions s₀(s) and r(s),

(19)

Using Eqs. 2, 3, 17, and 19, we can implement the stationarity condition for the free-energy functional (12) with respect to variations of membrane shape and cytoskeletal strain. This leads to a set of five coupled first-order differential equations, which are written down explicitly in the Appendix. Note that the ADE term in Eq. 3 enters only through the appearance of an effective spontaneous curvature

(20)

which must be determined self-consistently via Eq. 5. Although these Euler-Lagrange equations look complicated, they can be related to simple mechanical force-balance conditions in a way that makes their content entirely transparent (Mukhopadhyay and Wortis, in preparation). We have used these equations to calculate (approximate) spicule shapes.

We now turn to the boundary condition where the spicules meet the sphere (see Fig. 2). Along each element of L, the spicule membrane exerts a tension in the plane of the membrane and perpendicular to L, and a tension normal to the plane of the membrane. There is no tension in the third ("shear") direction because of the membrane fluidity. In addition, the membrane skeleton exerts an independent tension that must be in-plane, because the skeleton lacks bending rigidity. In general, these tensions would be different at different points of L; however, in the axisymmetric approximation, the tensions are uniform around L and the skeletal tension is directed radially. To maintain mechanical equilibrium, these tensions must balance where the spicules meet at point A. Note that A is a point of bilateral symmetry between the two adjacent spicules, so that in-plane tensions, which act in opposite directions, always balance by symmetry. On the other hand, the normal tensions must vanish individually, as they act in the same direction. The normal tension is related to the isotropic bending moment per unit length,

(21)

by (Evans and Skalak, 1980; Mukhopadhyay and Wortis, in preparation)  = dM/ds. Thus, an appropriate boundary condition along L is

(22)

This is the boundary condition that we shall apply in our calculations. Of course, this boundary condition is approximate, because the real spicule is only approximately axisymmetric. Indeed, with equal logic, we could argue by symmetry that should vanish at point B, where three adjacent spicules meet at a radius from the spicule axis some 15% larger than r_L. We shall use this observation in the Results to get a rough measure of the error introduced by the approximation (22).

We now outline the algorithm used to solve numerically the Euler-Lagrange equations for stationary spiculated shapes. A "shooting method" was used to integrate the Euler-Lagrange equations from the apex to the edge of the spicule, the contour L along which the normal tension vanishes. The set of solutions can be characterized by five parameters. Three of them are the global parameters: the pressure difference P, the Lagrange multiplier , and the effective spontaneous curvature C, Eq. 20. The remaining two parameters correspond to initial values of variables in the integration procedure, i.e., the curvature (C_m = C_p) and the stretching ratios (₁ = ₂) of the cytoskeleton at the pole. Integrating the Euler-Lagrange equations to a point where Eq. 22 is satisfied determines r_L and _L and, therefore, R₀ and n_s from Eqs. 13 and 16. The initial conditions may be adjusted iteratively to fit four parameters, A (Eq. 14), V (Eq. 15), ₀ (Eq. 8), and the cytoskeletal stretching ratio, , which measures the ratio of the area of the calculated spicule to that of the corresponding relaxed cytoskeletal patch (roughly the same as the ratio of the area of the full echinocyte to that of the full relaxed membrane skeleton, as the corrections in Eq. 14 for the spherical segments are small). We are left, finally, with a one-parameter family of mechanical-equilibrium solutions, which we take to be labeled by the spicule number n_s. The predicted equilibrium configuration is the solution n, which minimizes the overall energy. Of course, for the exact solution, n_s is a discrete variable. Probably, there exists, for given initial conditions, a set of distinct branches of solutions characterized by different values of n_s and different arrangements of spicules on the surface, as discussed more fully in the Discussion section.

Scaling analysis

Before proceeding to solutions, we wish to identify the important energy and length scales of the problem. Because we are working at a temperature that is effectively low, one overall energy scale drops out of the mechanical problem. We take this scale to be set by the bending modulus _b (_b and are similar in magnitude). There are three length scales. The first is the overall scale of the RBC, which we denote R. The second, we may take to be 1/₀ from Eq. 7 (which is equivalent to using a length based on ₀, Eq. 8); alternatively, we can take the second length to be 1/C, where C is the effective spontaneous curvature from Eq. 20 (we assume that ₀ and C are comparable in magnitude). The third length is the elastic length scale _el, Eq. 1. In general, there may be other length scales associated with the relaxed cytoskeletal shape; however, we shall assume that any such lengths are comparable with R. For the RBC, R _el. The length scale 1/₀ is a control parameter whose magnitude can be tuned across the range defined by R and _el, resulting in the observed RBC shape classes.

The size of the RBC sets the total area and volume in the problem. Note that the RBC has a relatively smaller volume for a given surface area than a sphere. If a sphere with the same surface area has a volume V_sphere, then the ratio of RBC volume to that of the equivalent sphere defines a reduced volume, v  V_RBC/V_sphere. For the RBC, v  0.6, significantly less than unity, which allows the RBC to assume its sequence of distinct nonspherical shapes. In its usual discocytic phase, the red blood cell is expected to have a magnitude of 1/₀ that is of the order of the red blood cell dimensions or larger. Tuning ₀ to sufficiently large positive values generates the sequence of evaginated echinocytic shapes; tuning it to sufficiently large negative values generates the sequence of invaginated stomatocytic shapes.

The significance of the elastic length scale _el is that it determines how strong the effect of elasticity is for structures with some characteristic length, say L. To see this, we imagine rescaling the problem with the energy scale _b and the length scale L in such a way as to make all the parameters appearing in the Euler-Lagrange equations dimensionless. In this way, we arrive at the rescaled parameters,

(23)

Equivalent rescaling of the variables, C'_m = C_mL, etc., leaves the form of the Euler-Lagrange equations unchanged, although the parameter values (Eq. 23) vary with L. For example, suppose L is the echinocyte spicule radius. If L > _el, then µ' is greater than unity, and the elastic terms from Eq. 2 dominate the bending-energy terms; conversely, if L < _el, then µ' is smaller than one, and the bending-energy terms dominate. Thus, generically, elasticity is strong at length scales larger than _el and weak at length scales smaller than _el. This statement may seem surprising to those who think of the discocyte as effectively shaped by bending energy. The key point here is the plasticity of the membrane skeleton, albeit on time scales much longer than the induced shape changes we are discussing. If the relaxed membrane skeleton closely matches the shape preferred by bending energy alone, then (of course) elastic effects are small. This is apparently the case for the resting discocyte (but definitely not for the echinocyte).

A similar scaling argument gives insight into the characteristic spicule size in the echinocyte region. Suppose we let L = _el in Eq. 23. At this length scale, the bending and elastic terms ('_b, K', and µ') are comparable and the only remaining parameters are P', ', and C (or, equivalently, '₀). If, in addition, P' and ' are small, then the Euler-Lagrange equations are equivalent to those for the unconstrained minimization of the combined bending-plus-elastic problem (2). But, this is a problem that we understand. When C is small in magnitude, then the membrane remains flat (at the length scale of _el); on the other hand, when C is large in magnitude, then bending energy (which dominates the shape at scales below _el) will favor budded shapes on the scale of 1/C, either invaginated or evaginated depending on the sign of C. Between these limits, when _el and 1/C are comparable, we expect spicules (for positive C) scaled by 1/C, tending towards smoother shapes when _el is somewhat smaller than 1/C and towards sharper, more columnar shapes when _el is somewhat larger than 1/C. Note that, as long as we rescale to the length _el, the expected spicule shape is predicted to be independent of cell size R and to depend only on the scaled value of C and the ratio of the two elastic constants, µ and K.

The above argument depends on the condition that the scaled values of P and  be small at the elastic length scale. For smooth and flaccid shapes like the discocyte (and when the cytoskeleton is not significantly stretched or compressed) the pressure difference P is generated by the bending energy, so P' is of order unity on the scale of the cell size R, and similarly for '. It then follows from Eq. 23 that the condition is well satisfied for the discocyte and other nearby shapes (as we shall see, the initial echinocyte shapes fall into this class). Of course, if the system is pushed too hard, then this condition breaks down and other length scales enter the spicule-shape problem.

In summary, the shape of a spicule arises as a result of the interplay between the two length scales, _el and 1/C. When 1/C is positive and much larger than _el, the effect of elasticity is strong and we can expect spicules that are at most gentle bumps. It is only when 1/C becomes of the order _el or smaller that we may expect sharp spicules. Thus, if we look at the series of shape transformations, from discocytes through discocytes with gentle bumps to sharply spiculated structures, the sharp spicules will first arise with a radius comparable with _el. We shall verify this scenario in the next section by explicit solution for the stable shapes.

	RESULTS

TOP ABSTRACT INTRODUCTION THE MODEL THEORY RESULTS DISCUSSION APPENDIX A REFERENCES

In this section we report echinocyte shapes calculated according to the program outlined in the Theory section and using the elastic moduli (_b, , µ, and K) given in Introduction and the RBC parameters (A and V) given in the Model section. The cytoskeletal stretching ratio  is not known to good accuracy for the RBC. In this section, we take  = 1; however, we have tested values in the range 0.7 to 1.2, and it is one of our results that variation within this range produces no appreciable changes in the number of spicules or their shape. The spontaneous curvature and relaxed area difference are not known and presumably vary somewhat over a typical RBC population; however, they enter the mechanics problem only in the combination ₀, Eq. 9, which we take as our principal control parameter.

Fig. 3 displays the equilibrium spicule shapes that we obtain for a sequence of values of ₀. Figs. 4 and 5 plot the corresponding calculated number of spicules n_s, which appear on the fully developed echinocyte III at equilibrium and the corresponding values of C as functions of ₀.

View larger version (22K): [in this window] [in a new window]   FIGURE 3   Sequence of spicule shapes for increasing values of the reduced effective relaxed area difference a0, Eq. 9, which measures the combined effect of spontaneous curvature and additional area in the outer leaflet in driving the membrane to bend outward. (a-d) Correspond to 0 = 0.014, 0.018, 0.020, and 0.022, respectively. Note the bar shows the elastic length scale el. Note how the spicule sharpens as 0 increases and how the neck begins to form when the spicule dimension (set by 1/C) falls below el so that the bending energy begins to dominate.

View larger version (12K): [in this window] [in a new window]   FIGURE 4   Calculated equilibrium number of spicules ns plotted as a function of the reduced effective relaxed area difference 0.

View larger version (12K): [in this window] [in a new window]   FIGURE 5   Calculated effective spontaneous curvature C, Eq. 20, plotted as a function of 0.

How well do the calculated shapes and spicule dimensions agree with experimental observations (Bessis, 1973; Brecher and Bessis, 1972; Chailley et al., 1973). As shown in Fig. 3, the calculated shapes fall into three distinct categories. For n_s less than ~20, corresponding to ₀ appreciably smaller than 0.016, the spicule shapes are broad and rather hat shaped, rather different from those seen in experiment. As ₀ and C increase, the number of spicules increases and the individual spicules become smaller, just as predicted by the scaling argument at the end of the previous section. By the time ₀ is at or above 0.016, for n_s in the range 30-60, the spicule shapes are columnar with rounded tops, in good general agreement with the kinds of shapes seen in experiment. Beyond n_s = 70, the spicules become narrow-necked, a shape not typically seen in experiment but illustrating the dominance of the bending energy, as discussed in the Theory section.

In the "good," central region, the spicule dimensions and numbers are very comparable with those seen in experiment. A typical example (₀ = 0.018) is shown in Fig. 6 a. For this case, we find n_s = 41, a radius of the central, spherical body of R₀ = 2.57, and an aspect ratio, height over width at base (which also corresponds approximately to the distance between adjacent spicules), of ~0.8, all of which are in reasonable agreement with observation. The only significant discrepancy is that the spicule height is close to 1.0 µm, somewhat larger than the range 0.6 to 0.8 µm seen in experiment, as is the aspect ratio. We shall comment on possible reasons for this discrepancy in the Discussion section. We emphasize that these shapes are a prediction of the model, not an assumption as they are in other recent work (Igli, 1997; Igli et al., 1998a,b). The only other spicule-shape predictions are those given by Waugh (1996), whose calculation suggests cross-sections that are less steep at the sides and more pointed at the apex than what is observed. Whether this is a consequence of the variational form assumed or of the neglect of stretching elasticity in the model is not clear. The calculations of Igli et al., based on a postulated spicule shape and assuming skeletal incompressibility (K  ), produce somewhat higher values of the aspect ratio. At a given reduced area difference ₀, they find a spicule number, which is larger than ours by almost a factor of two. We will return to some of these issues in the Discussion section.

View larger version (29K): [in this window] [in a new window]