INTRODUCTION

Top Abstract Introduction Conclusions References

Mechanical responses of cells originate in disparate physics over length scales ranging from intraprotein distances up through and beyond the characteristic units of organized assemblies. Cytoskeletal proteins may, for instance, unfold (Rief et al., 1997) or even dissociate (Evans and Ritchie, 1997) when a cell is extended. Network entanglements among thermally fluctuating filaments, in certain cases, may also contribute significantly to the elasticity of complex cells (MacKintosh et al., 1995). Such phenomena are among the many that reflect a novel hierarchy of scales in cytoskeletal mechanics. An illustration of some of the unique features that can arise in cell deformation after smaller scale details are integrated out is provided in the present paper, together with a companion work (Boey et al., 1998; referred to hereafter as paper I). Both papers focus on large deformation elasticity of the red cell membrane cytoskeleton. In paper I, Monte Carlo simulations have been presented for a submicron patch of several bead-and-tether idealizations of cross-linked spectrin chains. The present paper coarse-grains these quasi-ordered models, allowing a very general consideration of the submicron to cellular scales. Nonlinearities and associated anisotropies in large deformations of triangulated networks are explicitly revealed in stable nonhomogeneous states. Direct comparisons of ensemble-averaged computer "experiments" are thus made with published micromechanical tests on the red cell membrane's cytoskeleton.

Micropipette aspiration techniques have been applied to a range of cells, red cells in particular, for many years (e.g., Rand and Burton, 1964; subsequent work reviewed in Evans and Skalak, 1980). In standard analyses of such experiments, zero-temperature continuum notions and axisymmetry have been invoked to estimate elastic moduli and other constitutive responses. Such physical quantities no doubt have a basis in micro- and mesostructure. As a primary example of early success in identifying a material basis, the magnitude of resistance to aspiration of a flaccid red cell initially appeared to correlate moderately well with what was considered to be a molecularly thin structureless layer, of polymer-like chains, each with a size approximating that of spectrin (Evans and Skalak, 1980). In more recent fluorescence imaging measurements of red cells, detailed maps of membrane cytoskeleton deformation also correlated, to a degree, with some of the properties expected of a spectrin network (Discher et al., 1994; Discher and Mohandas, 1996). It was specifically shown that micropipette aspiration leads to a nonhomogeneous network deformation with stretching of the network as great as 250% and as small as 40%; spectrin, based on its contour length, is certainly expected to be capable of sustaining such a large range of deformation. Raw measurements such as these on red cells establish definitive benchmarks for more detailed models of the cytoskeleton. With the development of such models as a general aim, both this work and paper I focus on the erythrocyte's quasi-ordered triangulated meshwork of spectrin, a structure suggested by at least some electron microscopy studies (e.g., Byers and Branton, 1985; Liu et al., 1987).

In paper I, three variations of a polymer chain model are motivated and investigated as candidate descriptions of the erythrocyte cytoskeleton. In all of these models, temperature plays a central role, because the elasticity of the network arises from the configurational entropy, i.e., thermal fluctuations, of multisegmented and interconnected polymer chains. Furthermore, in all of these models, the chains are joined together at sixfold coordinated junction vertices, and the chains are attached by their midpoints to a flat plane representing the lipid bilayer. Where these microscopic models differ is in 1) the number of segments n_seg of the chain, 2) whether the cytoskeleton is under stress in its resting or reference state in situ, and 3) whether there are attractive interactions between chain elements. The computer simulations of paper I deal explicitly with the various complex steric interactions in the polymer chain networks and yield both elastic moduli and geometrical properties (such as the mean area per junction vertex). With either 12 or 26 segments per chain, the equivalent physical size of the plasma membrane that is directly simulated with ~10³ monomers is only several hundred nanometers on a side, more than an order of magnitude smaller than the human erythrocyte. In principle, there is nothing to prevent one from simulating larger systems, except that the number of segments in a model representation of a red cell would be 10⁶, which is beyond the reach of most researchers' computing resources. Two alternative strategies are available for developing simulations of whole cells based upon the models of paper I. The dimensionality of the problem could be reduced from three to two, or even a single curvilinear coordinate, by assuming that the network deformation has perfectly cylindrical symmetry. Such an assumption limits, by definition, both the symmetry of the deformation that can be examined, as well as the true nature of the network's response. A second approach, and that taken here, is to coarse-grain away some of the details of the polymer chain networks so as to decrease significantly the number of degrees of freedom.

In essence, our approach relies upon the fact that the motion of the junction vertices is similar to that of the nodes of a triangulated network at a low but nonzero temperature. Thus, the 3-D polymer chains can be replaced by two- and three-body 2-D interactions that effectively represent the many-segmented chains. This permits us to reduce the number of degrees of freedom by a factor of 3n_seg, a sufficient reduction that the simulation of whole cells is possible with modest computing resources. Fig. 1 illustrates an application of the effective representation technique: micropipette aspiration of an erythrocyte. In the figure, each bond segment represents a convoluted polymer chain which, in paper I, would have 12 or 26 segments. Such a whole-cell simulation can be used to investigate not only the global response of the cytoskeleton to imposed stresses, as exemplified by the length of the aspirated section of network in the pipette, but also the detailed responses, such as the average nodal density and fluctuating, anisotropic shapes of the triangular plaquettes near the entrance to the micropipette. The full statistical mechanical approach that we take is motivated by 1) actual experiments which show that micropipette-imposed network deformations are generally, as mentioned before, nonhomogeneous with very large strains; 2) theoretical complexity of the large deformation responses of triangulated structures emulating either polymer chain networks or, still more complicated, true Hookean spring networks in which a nonzero, force-free spring length exists; and 3) the nontrivial entropic contribution of thermal fluctuations at the biologically relevant temperature of ~300 K. Altogether, the general sort of mechanics problem that is considered within the present methodology and exemplified by Fig. 1 is rather complex: nonhomogeneous deformation of a nonlinear, anisotropic, and thermally fluctuating sparse material that may, in places, undergo hysteretic and finite-size dependent phase transitions.

View larger version (81K): [in this window] [in a new window]   FIGURE 1   Simulation of a small erythrocyte under aspiration. The micropipette, indicated by the solid gray shading, has an inside diameter of 12sR = 0.9 µm. The surface of the cell is triangulated with 6110 vertex nodes that represent the spectrin-actin junction complexes of the erythrocyte cytoskeleton. The volume of the cell is 0.6 times the fully inflated volume, and the simulation is drawn from the stress-free model in the free shape ensemble, as described in the text.

The contents of this paper are organized as follows. In the next section, effective representations are developed for each of the models in paper I. The representations are used in two different ensembles, referred to as the free shape and fixed shape ensembles, and we present computational details for each. As the principal applications of our approach, two types of aspiration experiments are simulated. The free shape ensemble is used to obtain the pressure-dependent aspiration of flaccid cells, i.e., a process in which the pressure inside the cell always approximates that outside the cell. The fixed shape ensemble is then used to map out the cytoskeleton's deformation in cells aspirated beyond the flaccid regime into a regime of whole cell pressurization. The main results are summarized in a concluding section.

	EFFECTIVE NETWORKS

In paper I, we determined the geometrical and elastic properties of three polymer chain networksall seemingly reasonable candidate models of the human erythrocyte cytoskeleton. Each polymer chain is represented by a series of hard beads with diameter a, linked together by tethers with a maximum extension of a. A total of n_seg tethers or segments make up each chain. This particular choice of tether length, combined with the hard core repulsion between all beads, enforces self-avoidance of the chains and gives rise to an effective bending resistance for each chain, because a next-nearest neighbor of a given bead in a chain is forbidden from passing between the given bead and its immediate neighbor. This bending resistance, at the scale of roughly two monomers, is employed later in a nonlinear elastic model for the chains.

In small deformation, although not in large deformation, the polymer chain models presented in paper I behave like low-temperature networks of linear Hooke's law springs. That is, at small stress only, the ratio of the area compression modulus to the shear modulus, and the stress dependence of the network area, are roughly those of a triangulated network of Hooke's law springs in two dimensions. To be clear, however, the complete properties of the polymer networks are not those of spring networks: the polymer networks display neither the collapse transition under compression (Discher et al., 1997) nor the unbounded expansion under tension (Boal et al., 1993) seen in spring networks. However, the averages of and the fluctuations in the positions of the junction vertices of the polymer networks are close to those of the junction vertices of a spring network whose spring constant k_sp is in the range k_sps_o² ~20-40, where  is the inverse temperature, (k_BT)¹, and s_o is the equilibrium spring length.

Because the fluctuations in the positions of the junctions correspond to a low-temperature system, it is tempting to search for an effective energy representation of the models in which the junctions are the primary degrees of freedom, and the effects of the chains are subsumed into two- and three-body interaction potentials between adjacent junctions. Obviously, a network composed solely of Hookean springs is not an appropriate description of the chain networks. A viable effective network must include the effects of both steric interactions among the chain elements (preventing network collapse) and the maximum chain length imposed by the tether constraints (preventing unbounded expansion). As a starting point though, it should certainly be well appreciated that a polymer chain exhibits elastic behavior because of its entropy, and an ideal chain of n_seg segments in three dimensions behaves like a polymer spring, with zero force-free length and an effective spring constant of 3/(n_sega²), where a is the segment length. However, the presence of steric interactions at short distances and the limitations imposed by the maximum length of the chain cause the stress versus strain relationship of real chains to be different from that of ideal chains. Although many simple but unmotivated functional forms can be constructed to represent the polymer chain models, the approach taken here is to faithfully model the results of paper I by first extending to networks a description of single chain elasticity recently proposed by Marko and Siggia (1995) and then, immediately thereafter, combining this with a Flory-like, mean field approach to chain elasticity balanced against sterics.

Marko and Siggia (1995) have proposed a simple interpolating model between ideal chain behavior at very small chain extension and divergent-tension behavior at large extension, i.e., as the end-to-end length approaches the full contour length. Their worm-like chain model appears to be particularly well suited to lambda DNA, which has a contour to persistence length ratio on the order of ~16 µm:0.1 µm, as well as the protein titin, which is composed of a sequence of separately folded immunoglobulin domains (Rief et al., 1997). Through a simple integration of the Marko and Siggia force-extension relationship, the end-to-end distance s of a chain with maximum extension s_max can be described by an effective potential,

(1)

where b is the persistence length of the chain segments and

(2)

This attractive potential has a minimum at s = 0, diverges at s = s_max, and has an effective spring constant near the minimum equal to (3/[2bs_max]), like that of an ideal spring. It should be remarked that the potential above truly represents a free energy, and that an expansion for small x differs most notably from the more classical, nonlinear freely jointed chain result (e.g., Discher et al., 1997) through the appearance of odd powers in x and a simple factor of 2 in the microscopic length scale.

The first element of the effective network, then, will be Eq. 1 as the two-body attractive potential acting along the triangulating "bonds" of a 2-D network. A minimum in Eq. 1 at s = 0, however, corresponds to a network collapsed to a point, and such a collapse is certainly not seen in the networks of paper I. Hence a second term added to the effective potential reflects those steric interactions between chain elements that prevent network collapse. For this repulsive energy arising from steric interactions, we choose a simple functional form, C/A, where A is twice the area of a single triangular plaquette. Thus the total energy of the effective network that we use to represent the polymer chain models of paper I is

(3)

where C > 0 prevents network collapse. The mean-field balance of sterics against chain elasticity should be recognized as Flory's classical approach to real chains. It is worth noting that Eq. 3 can be Taylor expanded for a triangulated network in the zero temperature limit, and this expansion yields an expression for E_net (see Eq. 19 of Discher et al., 1997) that is similar in form, at lowest order, to that employed in prior continuum analyses of red cell cytoskeleton deformation (specifically, Discher et al., 1994). It is also noteworthy that the worm-like chain part of the model employed here is algebraically simpler than the freely jointed chain model; in C₆-symmetrical networks, however, the two polymer models have similar linear regimes and can be shown to be anisotropic at higher order (Discher et al., 1997).

Equation 3 contains three parameters, C, b, and s_max, which are determined by fitting the predictions of Eq. 3 for network area as a function of stress against the results found from the different polymer chain models. To simplify the fitting procedure, we derive a simple relationship between C and the bond length s at essentially zero temperature and stress. This is a state in which all bonds have the same length s_o, and the area per vertex is A_o = s_o²/2. The expression for C is obtained by demanding that E_net be a minimum at s = s_o, or

(4)

Solving Eq. 4 gives C in terms of the other two parameters, b and s_max:

(5)

(6)

By expressing the areas from the polymer chain models in terms of the equilibrium interjunction distance s_o, the fitting procedure involves only two parameters, b and s_max. In principle, the fit should reflect the effective temperature of the network through  appearing in Eq. 5; in practice, only the combination C is relevant in the Monte Carlo simulations of the effective networks.

To perform the fit, we construct a mean-field version of the effective network. The procedure is reasonably accurate, because we know from previous work (Boal et al., 1993) that the mean-field approximation provides a good description of triangulated spring networks at low temperature. Furthermore, because the mean-field approach is analytic, the fitting procedure is computationally trivial to perform. Once b and s_max have been fitted, the results are double-checked by performing a full simulation of the effective network and comparing again with the full polymer chain networks. The stress dependence of the polymer chain network is obtained from the mean field version of Eq. 3 at constant pressure, yielding the effective free energy per junction vertex _j,

(7)

where the in-plane pressure  is positive for a network under compression and s = s(). That is, for a given , a value of s can be found that minimizes _j for any choice of b and s_max; hence, the area per vertex of the effective network can be predicted as a function of the applied pressure  for any parameter set. By comparing the predicted areas from the effective network with those found in the model cytoskeleton simulations, a ² statistic can be calculated for each b and s_max. The fitting procedure searches for the values of b and s_max that minimize ².

To review, the three polymer chain models from paper I are characterized by the quantities n_seg, the number of segments in each polymer chain, and s_o, the average distance between junction vertices at zero stress. The average distance between junction vertices in the red blood cell, s_R, has the physical value of 75 nm, so that if the model network for the cytoskeleton is not under stress, then s_o is equal to s_R (see paper I). However, if the model network represents a cytoskeleton under stress in the erythrocyte, due, for example, to assembly onto the encapsulating bilayer, then s_o may be different from s_R. The area per vertex of the effective network at zero stress, A_o, is s_o²/2; however, the reference state value of the area per vertex in the cytoskeleton is A_R = s_R²/2. The chains are assumed to have a nominal contour length l_c = 1.2n_sega = 200 nm, with the factor 1.2 arising from the expectation value for tether length fluctuations between a and a. Corresponding to the contour length is a nominal contour area, A_c, of l_c²/2, and each of the models is constructed such that A_c/A_R  7, as estimated experimentally (Byers and Branton, 1985; Liu et al., 1987). The values of n_seg and s_o in each polymer chain model are as follows.

Stress-free model. This model fixes n_seg at 26, and assumes that the equilibrium state of the cytoskeleton in vivo is at zero stress; hence s_o = s_R.

Condensed model. Here, n_seg is fixed at 12, and the equilibrium state of the cytoskeleton in the cell is at zero stress (again, s_o = s_R). However, there is an attractive interaction between the vertices of the polymer chain that reduces the value of the interjunction separation below that for a network without interactions (see paper I for further details). The strength of the attractive interaction is set so that A_c/A_R  7.

Prestress model. The number of segments n_seg equals 12 in this model, as in the condensed model, but there are no attractive interactions between chain vertices. Because a network with n_seg = 12 has A_c/A_o  3.78, the network must be placed under stress to force A_c/A_R  7 in a simulation of an erythrocyte. That is, the junctions are forced to be closer together in the cell (s_R) than they would be if the cytoskeleton were extracted and viewed flat in isolation at zero temperature (s_o) (again, see paper I for further details). This can be accomplished by setting up the geometry of the cell such that s_o = 1.36s_R.

Fig. 2 shows a comparison of the polymer chain simulations from paper I with their effective representations in the mean field limit. Over 12-fold changes in network density or area, the effective potential tracks the full polymer chain simulation rather well. Given that the effective potential was fitted in the mean field limit, each parameter set was checked by constructing a full two-dimensional network in the isobaric-isothermal ensemble, using the effective potential with the appropriate parameters from the mean-field fit. The ensemble averages for geometrical quantities found in the full network agreed to within a few percent of their mean-field values. Specifically, the fits shown in Fig. 2 yield

(8a)

(8b)

(8c)

The quantities b and s_max have physical meaning for single chains as the persistence length and maximum extension of the chain. However, chains in a network have somewhat different characteristics compared to isolated chains because of the interchain interactions present in a network (see Boal, 1994). Thus we must treat b and s_max as quantities that have values specific to each polymer network, although we expect that their magnitudes are not greatly removed from those of single chains.

View larger version (23K): [in this window] [in a new window]   FIGURE 2   Area per vertex Aj as a function of two-dimensional pressure. Both quantities are quoted in the dimensionless combinations Aj/Ao and Ao. Comparison is made between the full cytoskeleton simulations from paper I against their mean field representations in two dimensions. The data points are taken from paper I. The solid lines indicate the effective representations with values for the parameters smax/so and b/so as given in the text. Ao is the average area per vertex at  = 0. In both the stress-free and condensed models, Ao is equal to the in situ area, or reference area AR, of the cytoskeleton. In the prestress model, Ao is greater than AR.

Consider b/s_o. For a single chain, the persistence length b can be extracted from the mean square end-to-end displacement of the chain r_ee² via the approximate expression

(9)

Replacing r_ee² with s_o² and substituting the contour length l_c = 1.2n_sega, Eq. 9 gives an approximate value for b in the polymer chain models. As can be seen from Table 1, which summarizes the results for s_o found in the polymer chain models, the estimated values for b/s_o are in the correct range, but are about a factor of 2 higher than the values from the fits displayed in Eq. 8. Note that, owing to their mutual interaction and interaction with the flat bilayer, chains in a network have a greater r_ee² than do chains in isolation. Significantly improved estimates for b/s_o are obtained by noting that worm-like chains are expected to have b = a/2.

Consider s_max/s_o. For the bead and tether chains in the cytoskeleton models, the nominal contour length l_c is less than the maximum distance to which the chains can be stretched under infinite tension. Numerically,

(10)

whereas l_c = 1.20n_sega. Because the stress dependence of the network area was determined for large stresses in paper I, the effective maximum chain length is closer to s_max than to l_c, although the two differ by only 15%. One can see from Table 1 that the expected values of s_max/s_o are within 20% of the values quoted in Eq. 8 from fitting the stress dependence of the area.

From the fits to the polymer chain cytoskeleton models, we draw a number of conclusions:

1.	The effective potential given by Eqs. 3 and 5, based on short-range steric repulsion and nonlinear entropic elasticity of polymer chains, provides a very good representation of the simulation results in paper I.
2.	The parameters of the effective potential are in the approximate range expected from the elementary structure of single chains.
3.	Because the basic network elements or bonds have a definitive maximum length, triangular networks of these nonlinear elements will be highly anisotropic in their stress versus strain relationship, as elaborated elsewhere (Discher et al., 1997).

	ASPIRATION SIMULATIONS

The nonlinearity of the effective models, the network anisotropy, the finite temperature nature of the physics, the possible prestress on the cytoskeleton, and other complicating aspects of structure all motivate the use of the above effective potential in direct simulations of whole cell deformation. Micropipette aspiration of red blood cells is the focus here, but the models can certainly be applied to other problems, such as the motion of blood cells through capillaries. An effective potential reduces the number of degrees of freedom by replacing multisegmented chains by few-body interactions, so that the deformation of a single cell with many thousands of junction vertices can be simulated using a conventional workstation. Two types of aspiration experiments are simulated in this paper, and we treat each type of experiment with a different type of ensemble. In one set of experiments, a flaccid erythrocyte is aspirated under moderate pressure and is simulated with a free shape algorithm. In a second set of experiments, a swollen cell is aspirated under large pressure, and a fixed shape algorithm is employed. Both codes use the Metropolis Monte Carlo technique to determine the ensemble-averaged detailed deformation of the cell under aspiration. The codes are described in some depth in this section before the results from the simulations are presented.

Free shape

A snapshot from this simulation is shown in Fig. 1. The system consists of a simple, closed surface decorated with 6110 vertices, all but 12 of which are joined to six neighboring vertices. The initial configuration has the shape of two parallel sheets in the form of a hexagon, with each vertex on the perimeter of one sheet connected to two vertices on the perimeter of the other sheet. Thus, there are six vertices on each sheet, or 12 in total, that are at the "corners" of the hexagons and therefore have only fivefold coordination. Note that the 12 fivefold defects in the model are the minimum number required by topology for the triangulation of a spherical surface. The connectivity is fixed, in that each vertex has a fixed set of neighbors, even though all vertices may move in space subject to a collection of energetic restraints. The total energy of the system has several components:

1.

The two- and three-body potential Veff of Eq. 3, which represents the in-plane properties of the polymer chain model of the cytoskeleton.

2.

A bending energy Ebend, which represents the bending resistance of the lipid bilayer of the plasma membrane, and which is constructed from a set of unit vectors n normal to each triangular element of the cell's surface. The bending energy then involves a sum of 3N  12 terms, each corresponding to one of the 3N  12 bonds between the N vertices that define the surface. Explicitly, (11) where i and j are labels for neighboring triangular elements of the surface. Experimentally, the continuum bending resistance kc for an erythrocyte bilayer (including cholesterol) has been measured to be ~20kBT, which corresponds to a value for the descretized bending resistance kbend of 69kBT (because kbend = 2kc; see Boal and Rao, 1992). 3.  A term Esurf, which enforces approximate surface area conservation. Because the lipid bilayer is relatively incompressible, we add a term to the energy to suppress fluctuations in surface area. The reference surface area Acell is defined to be (12) where sR is the reference bond length whose physical value is 75 nm. The form of the surface energy term is chosen to be (13) where ksurf is a parameter. Choosing the value ksurfsR2 = 600 constrains the root dispersion of the surface area to lie within a few percent of Acell for surfaces at zero stress. 4.  A term Evol, which enforces approximate volume conservation. Because the enclosed volume of a human erythrocyte does not change appreciably during deformation, either in vivo or in the aspiration experiments, a term (14) is introduced to constrain the volume, where kvol is a parameter. Choosing the value kvolsR3 = 600 limits the root dispersion of the volume to lie within a few percent of the desired Vcell.

In the simulation, we found it efficient to attempt to move four vertices simultaneously, then accept or reject the move on all four according to the change in summed energy of Eqs. 3, 11, 13, and 14. The initial pancake-like configuration of the cell is allowed to relax for ~10⁶ trial moves per vertex to generate an equilibrated configuration for the aspiration simulation. During the relaxation process, the volume of the cell is driven toward a predetermined value, as a consequence of the presence of Eq. 14 in the Boltzmann weight. A normal human erythrocyte is not spherical, but has a volume of ~60% of the volume of a sphere V_sphere, which is given by

(15)

for a surface of area A. Because the free shape calculations are meant to simulate the aspiration of a flaccid cell, we fix V_cell = 0.6V_sphere.

The shape of the computational micropipette is shown in cross section by the gray area in Fig. 1. The pipette is a hollow cylinder of radius R_P over most of its length, but its mouth is a semicircle in cross section with a radius of R_P/2. Because the cell has a discretized surface in the simulations, a rounded entrance to the simulation pipette is needed to reduce computational friction as the cell is drawn up the pipette. Our coordinate convention is that the bottom edge of the pipette is a ring in the xy plane (with radius 3R_P/2), and the symmetry axis of the pipette is the positive z axis. We choose R_P = 6s_R, corresponding to a pipette with an inside diameter of 0.90 µm. A pressure P is applied to the cell through the interior of the pipette after equilibration of the cell configuration is finished. The complete Boltzmann factor exp() for the Monte Carlo moves thus involves

(16)

where V_pip is the volume of the upper surface of the cell contained within the pipette region with z > 0. The sign convention of Eq. 16 is such that P < 0 acts to increase V_pip and pull the cell into the pipette.

The free shape simulation has been run for all three effective networks from paper I, at five aspiration pressures per model. The simulation times required for the cells to reach their equilibrium configuration under aspiration are rather long: although the cell boundary is drawn to 90% of its aspirated length within a moderately short time, the approach to the final configuration is typically about (1-2) × 10⁶ Monte Carlo attempts per vertex. Thus we allow at least this number of moves per vertex for the cell to reach equilibrium, depending upon the model, then collect 20 samples of the configuration separated by 10⁴ moves per vertex. For a total of 2 × 10⁶ moves per vertex, the deformation of a 6000-vertex cell obeying Eq. 16 can be simulated in a week's cpu time on a 200-MHz workstation.

Fixed shape

The fixed shape ensemble, in which a model network is allowed to relax over a specified shape, is related to previous continuum mechanics analyses of cell network deformation that assumed axial symmetry (Evans and Mohandas, 1994; Discher et al., 1994). This reduced the three-dimensional cell shape to a specified curve in two dimensions. The approach is extended here by dropping the assumption of axisymmetry and performing a full three-dimensional, finite-temperature simulation on a specified or fixed surface. The surface over which the network nodes move is broken down into several simple geometric elements: a spherical surface (representing the main body of the cell) is joined to the interior surface of the micropipette, whose geometry is the same as that employed in the free shape ensemble. The network follows the cylindrical interior of the pipette, before ending in a fixed hemispherical cap of radius R_P. The distance along the symmetry axis from the entrance of the micropipette to the tip of the cap is defined as L. Two sample configurations are later shown in Fig. 5. Within the accuracy of double precision, the network nodes are required to move on the fixed surface. The computational advantage of this approach is that only the E_net term in Eq. 16 need be evaluated; the resulting gain in the execution time of the simulation allows much larger systems to be investigated compared to the free shape algorithm. Furthermore, the relaxation time of the network is significantly shorter than the free shape approach. Although the cell shape itself is an input to the fixed shape ensemble, the free shape ensemble results will demonstrate the appropriateness of the assumed geometry.

The discrete network in the fixed shape ensemble consists of 18,434 nodes, of which 12 are symmetrically distributed fivefold centers and the remainder are sixfold centers. Because a human red cell, in comparison, has ~30,000 spectrin-actin nodes, the model's sphered diameter is within 22% of the average for sphered human red cells.

For a first ensemble, the nodes are all placed on the surface of a sphere, and ~2 × 10⁶ Monte Carlo sweeps are taken with all nodal motions attempted only along lines of principal curvature of the defined surface. Correct Boltzmann sampling with curvilinear motions requires a supplemental weight factor given by a simple ratio of the node's radial distance from the surface generator's axis, r_new/r_old. In addition, a nodal move is rejected outright if any of the node's six (or five) local normals makes an angle greater than /2 with respect to the surface normal at the node. This constraint takes the place of E_bend in Eq. 16 by establishing a signed-area sterics similar in form to that adopted previously in simulations of strictly planar nets (e.g., Boal et al., 1993; Discher et al., 1997). In a second ensemble, the spherical surface is progressively deformed in a simple sequence of reversible but nonequilibrated mappings from a sphere to a new equi-area shape having a projection in the z direction of L = 2R_P. The equi-area constraint reflects the relative incompressibility of the lipid bilayer and necessitates a decrease in the radius of the sphere outside the micropipette as the projection length is increased. As a consequence, no constraint is thus placed on the volume enclosed by the surface. In experiments, such cell volume adjustment is readily achieved osmotically, with the result being that the projection length is essentially set by osmotic rather than aspiration pressure (Discher et al., 1994). Once L = 2R_P is achieved, ~2 × 10⁶ Monte Carlo steps are taken to relax the network from the stretched state imposed by the initializing transformation. Longer projections were incrementally achieved by subsequent equi-area mappings followed by extended relaxation intervals. Each simulation reported here took ~1 week on an R10000 processor of an SGI-Cray Origin 2000.

	SIMULATION RESULTS

As a first step in simulation, we revisit the type of aspiration experiments represented in Waugh and Evans (1979), in which flaccid red blood cells were aspirated under modest pressure. By measuring the length of the aspirated projection L (in the simulation shown in Fig. 1, L is the distance from the bottom of the pipette to the top of the network), an apparent shear modulus of 6-9 × 10⁶ J/m for erythrocytes had been extracted. It is shown in paper I that the three polymer chain models have shear moduli in this range, and so a simulation of the full cell should approximately reproduce such L versus P results. The free shape ensemble is used to simulate the experiments; the results are shown in Figs. 3 and 4.

View larger version (23K): [in this window] [in a new window]   FIGURE 3   Aspiration length L/RP as a function of in-plane tip tension in the free shape ensemble. The representative experimental data of Waugh and Evans (1979) for flaccid human erythrocytes are indicated by the solid gray triangles. Between different cells, a standard deviation of 10-20% is typically seen in the slope of the experimental data. In all situations, the in-plane stress is obtained from  = RPP/2, where P < 0 is the aspiration pressure; the sign convention employed is that  < 0 corresponds to positive tension in the network.

View larger version (19K): [in this window] [in a new window]   FIGURE 4   Ratio of the cap area per vertex Acap to the reference area per vertex AR = sR2/2 in the free shape aspiration simulations of two effective models. The in-plane pressure  in the aspiration simulations is found from the law of Laplace:  = RPP/2. Shown for comparison by the solid curves are the area ratios expected for the same networks confined to two dimensions and subject to an in-plane stress (calculated in the mean field limit).

The predictions of the three models for L/R_P as a function of pressure are plotted in Fig. 3. For both the simulation and experiment, the pressure in the figure is actually the two-dimensional pressure, , at the tip of the aspirated section as calculated from the law of Laplace, namely,

(17)

The conversion to physical units in Fig. 3 uses s_R = 75 nm and ¹ = 4.0 × 10²¹ J. Given the uncertainties in the experiment and simulation, the agreement is acceptable for all three models. The stress-free model tends to be stretched further into the pipette than the other two models for a given pressure.

Strictly speaking, Eq. 17 applies only if the cap is hemispherical in shape. We therefore use a self-consistency check to verify its applicability. Fig. 2 illustrates how a given network model's area per junction changes as a function of applied pressure. If the law of Laplace is applicable to the aspiration simulations, then the area per junction at the cap of the aspirated segment as a function of  from Eq. 17 should be the same as in Fig. 3. Suitably reduced units facilitate the comparison shown for the stress-free and condensed models in Fig. 4. The agreement is excellent, confirming that the law of Laplace is appropriate to the cap region. Good agreement is found for the prestress model as well.

More detailed studies of erythrocyte deformation have been made possible by advances in imaging techniques. In recent red cell aspiration experiments (Discher et al., 1994), the density profile of the cytoskeleton was viewed on micropipette-pressurized cells containing fluorescently labeled proteins. Osmotic adjustment of cell volume was used to control the length of the aspirated projection subject to the constraint of nearly constant cell area as imposed by the lipid bilayer. These experiments demonstrated that the surface density of cytoskeleton at the entrance to the micropipette is higher than the mean density of the undeformed cell, and that the density decreases quasilinearly along the length of the aspirated projection toward the cap. Our simulation of these experiments is based upon the fixed shape ensemble, and Fig. 5 displays two configurations characterized by projection lengths of L/R_P = 2 and 8. The specific model parameters used for the simulations of Fig. 5 are those of the Prestressed model.

View larger version (143K): [in this window] [in a new window]   FIGURE 5   Simulations of an effective cytoskeletal network in two fixed shape ensembles of different aspiration lengths (A and B); also shown (C) are the corresponding relative density profiles. The network consists of 18,434 nodes, predominantly sixfold coordinated, and the nodes are confined in their motions to four smoothly connected and well-defined surfaces: a large spherical body, a right circular cylinder of radius RP, a quarter-sector of a torus (minor radius 0.6RP) that connects the cylinder and sphere, and a hemispherical cap. The latter three surfaces define the aspirated projection, of length L, from z = 0. The projection is formed incrementally through a sequence of equi-area transformations from an initial sphere of radius Rs. In physical units, if sR = 75 nm is assumed, then RP = 668 nm and the initial Rs = 2.67 µm. The density profiles, as projected into bins along the z axis, are normalized in C by the homogeneous density of a network on a sphere. The particular model shown is the prestress model.

To begin the discussion of the fixed shape results, we describe several generic features of the aspiration experiments represented in Fig. 5. After the initial establishment of an aspirated configuration, a steady state is reached in which nodal fluxes between the four geometric surfaces composing the aspirated shape, along with other monitored quantities, appear to fluctuate about stable averages. One such monitored quantity is the relative density distribution of nodes. With _R denoting the mean density of network nodes over the surface of the undeformed cell, the relative density,  = /_R, is obtained by averaging and suitably area-normalizing the number of nodes to obtain the average nodal density in fixed increments of width z = R_P/4 along z. In physical dimensions this interval closely corresponds to the maximum resolution of an optical microscope. For the two projection lengths shown in Fig. 5, A and B, the ensemble-averaged density profiles along the z direction in Fig. 5 C indicate that the relative density of nodes over much of the sphere outside of the micropipette is uniform (within 5%) and essentially unchanged from the nonaspirated spherical state. Toward the entrance of the micropipette, the relative surface density of nodes generally increases only to decrease, quasilinearly, to a minimum value near the tip of each projection's cap, i.e., z/R_P = 2 or 8 in Fig. 5. These features, namely a uniform relative density of network over much of the sphere outside the micropipette and a gradient in relative density on the aspirated projection, are the key signatures of red blood cell cytoskeletal deformation, as revealed in recent micropipette experiments employing fluorescence imaging (Discher et al., 1994). Quantitative comparisons of relative density profile characteristics allow a discrimination between network models that fit experiments and those that do not.

As a consistency check between the fixed shape and free shape ensembles, the work required to incrementally lengthen the fixed shape projection is computed. The total energy of an aspirated configuration is calculated by simply summing the effective E_net over the entire network and subtracting the same quantity for the nonaspirated sphere. The ensemble-average energy of the prestressed worm-like chain network is shown as an increasing function of projection length L/R_P in Fig. 6 A. As concluded in paper I, the characteristic energy scales of the various network models are all consistent with a low-temperature state. Consequently, it is assumed that free energy changes with projection length are dominated by the sum total of the effective energy, at least for small projection lengths, so that changes in the entropy associated with, say, restrictions on nodal motions contribute negligibly. The work done by the aspiration pressure in displacing the cytoskeletal projection can then be equated with the additional energy stored in deformation of the cytoskeleton, ignoring for this analysis the constraints imposed by cell area and volume. Thus, to extend the cytoskeletal projection from L/R_P = 2 to L/R_P = 3, corresponding to an incremental length change of R_P and volume change R_P³, the necessary increase in pressure is simply given in the energy plot (Fig. 6 A) by the [initial slope]/(R_P³) = 8100k_BT/(668 nm)³ = 350 dyn/cm². Previous discussions in this paper of the aspiration of flaccid red cells (e.g., Waugh and Evans, 1979) suggest a value somewhat closer to 300 dyn/cm², a value not significantly different, given the uncertainties in experiment. Exact agreement is not to be expected, in any case, because the present ensemble is based not on flaccid cells as in the first part of this paper, but on cells pressurized by a strong aspiration leading to a distinct overall geometry of deformation. Note in Fig. 6 A that the energy scale is in units of ~10⁴ k_BT, which, together with the thousands of molecular degrees of freedom involved in deformation, indicates a work per molecule on the order of several k_BT. Such an energy scale is obviously well within the regime where Boltzmann sampling and thermal fluctuations are both extremely relevant, as further elucidated below. One last feature of note in Fig. 6 A is that the strain energy in the network grows more rapidly than the projection length, reflecting, in part, the intrinsic nonlinearity of the deformation.

View larger version (38K): [in this window] [in a new window]   FIGURE 6   (A) Cumulated strain energy in the prestress model as an increasing function of projection length L/RP. (B) Relative density of aspirated networks at the entrance of the micropipette and at the projection's cap: experiments (gray; Discher et al., 1994), simulations of the prestress model (black), and simulations of the stress-free model (plus signs). Results for the condensed model are within a few percent of the stress-free model. Horizontal error bars on the simulation results all have a length equal to the height of the quarter-sector of the torus at the micropipette entrance, as motivated by the observation that the comparable dimension on actual micropipettes used in experiments is vanishingly small.

Density profiles of aspirated model networks are obtained, as described above, by z axis projection for direct comparison to experiment. Fig. 6 B shows very clearly that simulations of both the stress-free and condensed model underestimate the relative density at the cap of the projection, even with the significant uncertainty in experiments. The prestressed network model, however, provides a better fit of relative density at both the cap and the entrance. The close agreement with experiment provided by this latter model thus suggests that at least one set of the microscopic simulations reported in paper I correlates well with some of the available micromechanical measurements. In addition to the three network models focused upon here, we simulated the aspiration of both unappended Hookean spring networks, i.e., linear springs having a nonzero force-free length (e.g., Hansen et al., 1996; note that polymers have a force-free length of zero), and lipid bilayer models in which individual plaquettes strongly resist area changes. The latter model of the bilayer gave, consistent with experiment (Discher et al., 1994), a homogeneous density distribution over the aspirated projection as well as the rest of the cell. The Hookean spring model, in contrast, exhibited cap densities that were closest to the (poorly fitting) stress-free model, despite the complicating appearance of predicted broken symmetry states (Discher et al., 1997) in the compressed region at the micropipette entrance. Importantly, only by the sort of finite temperature simulation approach taken here can one deal rigorously with the nonhomogeneous deformation in micropipette aspiration of networks susceptible to hysteretic phase transitions and finite size effectstwo demonstrated characteristics of Hookean spring networks.

Reasonable agreement between the one set of simulations and experiment may possibly be due to the successful model's area modulus in the reference state, K_A^Ref, being roughly four- to eightfold greater than the reference state shear modulus, µ^Ref. By reference state, we again mean the state of the cytoskeleton in the undeformed model cell; the reference state is unstressed in two of the models and prestressed in the third. This explanation in terms of moduli was first suggested by previous continuum mechanics analyses (see footnote on p. 809 of Mohandas and Evans, 1994), with moduli suitably transformed into those lowest order moduli commonly employed and defined in paper I. Combined with the thorough study of microscopic models in paper I, it is clear that 1) the zero surface pressure ( = 0) value for K_A/µ of the microscopic bead-and-tether networks and Hookean spring networks is invariably close to 2, and 2) only in a reference state that is compressed, i.e., prestressed, is this characteristic ratio larger. Beyond the reference state values for moduli, it is important to emphasize that the present network approach rigorously includes, through the worm-like chain model, nonlinear chain mechanics that have heretofore been omitted from all cell deformation analyses, continuum mechanics or otherwise. Partly because of such omissions and partly because of an absence in prior analyses of complicating structure as described below, it is remarkable that significant guidance in simulation has been provided by the sort of lowest order moduli identified in previous analyses.

Beyond comparisons to well-documented experimental results, additional descriptions of structure are readily garnered from the present discrete simulation approaches. What follows are results that would be exceedingly difficult to obtain in detail by most other methods of analysis. The focus of this final results section is the prestress network in a moderate length projection of L/R_P = 4. First, spectrin tethers are not stretched the same along the length of the projection, as is apparent in Fig. 5 B. Toward the micropipette entrance, spectrin is very clearly compressed in the circumferential direction and extended in the axial direction; near the projection tip, in contrast, spectrin tethers are more isotropically stretched in being part of equilateral triangles. This stretching is quantitatively elaborated in Fig. 7 A, where ensemble-averaged distributions of spectrin stretch between network nodes are illustrated for three intervals along the projection. Note that the maximum spectrin stretch in this model had been fitted to be s_max/s_R = 3.12, so that a significant fraction of the chains at the pipette entrance are strongly stretched to approximately two-thirds of their maximum. Within the worm-like chain model (Eq. 1), such extensions are associated with forces in excess of the linear or Gaussian regime forces by at least a factor of 2, i.e., the nonlinearities of the model become increasingly important. In addition, as noted above, a significant fraction of chains are also strongly compressed in the vicinity of the entrancehence the bimodal distribution for spectrin stretching. Near the very tip of the same projection, in contrast, the stretch distribution is single-peaked albeit broad enough to include a small fraction of compressed chains.

View larger version (34K): [in this window] [in a new window]   FIGURE 7   Mesostructure in network deformation for a projection length L/RP = 4. (A) Average distributions of spectrin stretched between network nodes, s/sR, at three intervals along the projection. Peaks of the distributions are approximated by two geometric stretch ratios denoted in the text and in Fig. 8 by m and , and calculated essentially from the relative density profile alone. (B) Anisotropic thermal motions of two representative nodes of the network near the entrance of the pipette and near the very tip of the projection. (C) Projected profiles for spectrin orientation components, Nx2, Ny2, Nz2, and the relative density for an entire aspirated cell network. The orientation components here are simply the squared projections of spectrin chain unit vectors, e.g., cos2x, ensemble-averaged over all of the chains in each z bin. Note that the aspirated projection is to the right of z = 0.

The peaks of all of the distributions in Fig. 7 A can be located fairly accurately by two "mean" stretch values, one parallel and one orthogonal to the surface generator. These mean stretch ratios, denoted by _m and , are illustrated in Fig. 8 in a rectangular distortion of an initially square piece of elastic surface on a flaccid cell. Estimates of such mean stretches are obtainable directly from an integral over the projection's relative density profile and with respect to z' = L  z; i.e.,

(18a)

(18b)

A boundary condition at the very tip or pole of the cap takes the form _m =  = 1/, and the radius of the spherically swollen, simulated cell before aspiration is R_s. These simple expressions, with r(z') being the radial distance from the surface at z' to the symmetry axis of the projection, have been derived previously in the context of analyzing experimentally determined density profiles (Discher and Mohandas, 1996). In the present theoretical work, we see that these same quantities, which presuppose axisymmetry, correspond very closely to local means of the roughly bimodal distributions of chain stretch. The distributions of stretch reflect an anisotropic structure in a nonaxisymmetrical deformation. It is also clear that the maximum forces in chains of the deformed network are achieved in those chains counted in the high stretch tails of the distributions. In the development of field theories of network failure, such considerations may prove crucial, because even a small underestimation of stretching of nonlinear polymer chains can lead to a gross miscalculation of the sustained force. With the present aspirated projection and nonlinear worm-like chain model, for example, a chain force of ~2 pN is estimated from the maximum mean field stretch given by _m = 2, whereas a number of chains in the network actually sustain much greater relative stretching of 2.5, corresponding to a significantly higher force of ~6 pN. Last, these numerical values for force help support a central, underlying premise here and in paper I, which is that red cell deformation, even very large deformation, is generally within the elastic regime; this is because protein dissociation and unfolding under force, despite possible loading-rate dependencies (Evans and Ritchie, 1997), is not expected to occur up to moderate experimental time scales, except at forces that are perhaps an order of magnitude higher, as found for titin (Rief et al., 1997). Direct measurements of spectrin dissociation under force are, however, certainly needed.

View larger version (58K): [in this window] [in a new window]   FIGURE 8   Mean stretch ratios, m and , as surface distortion metrics in the torsion-free deformation of an initially square surface element.

On top of the strong, nonlinear stretching of chains between network nodes, significant thermally driven displacements of nodes are observed. Such stochastic motions are explicitly illustrated in Fig. 7 B for two representative nodes of the network: one node near the entrance of the pipette, and one node near the very tip of the projection. At the pipette entrance, the ensemble of fluctuations about some stationary average position is seen to be different in magnitude between two orthogonal directions on the surface; the typical anisotropy ratio is between 1:2 and 1:3. This certainly reflects anistropy in the deformed local structure and is undoubtedly associated with the strong circumferential compression evident in the leftmost distributions of Fig. 7 A. Accordingly, a node at the tip of the projection should and indeed does exhibit more isotropic fluctuations. Moreover, because of the softness of the network, the amplitude of the fluctuations is often a very significant fraction of s_R = 75 nm. With modern approaches to nanometer-scale particle tracking, experiments on such local scales seem feasible and should provide insights into local structure in deformed cells. Obviously, the thermal fluctuations of network nodes is a natural and very physical feature of the systema distinctive feature excluded from finite-element and other sorts of continuum mechanics approaches.

Another signature of structural heterogeneity in the deformed network is provided by the spectrin orientation fields, N_x², N_y², and N_z², shown in Fig. 7 C as ensemble-averaged, z-binned profiles. The orientation components here are simply the squared projections, e.g., cos²_x, of unit spectrin chain vectors pointing along the "spectrin" bonds between nodes. Polarization microscopy experiments may prove particularly useful in correlations with these structural simulations. The simulated projections, for example, exhibit strong chain alignment in the pipette axis direction, as indicated by the fact that the largest orientation component is N_z². This is a feature of structure related to the dominating high extension peak in the two leftmost distributions of Fig. 7 A. For reference, the average squared projection of a unit vector constrained to a surface but otherwise isotropically oriented is simply 1/2; in contrast, the squared projection of a unit vector isotropically oriented in three dimensions is 1/3, as is well known in the field of liquid crystals. Furthermore, a unit vector that is, respectively, parallel or orthogonal to a fixed coordinate axis has a squared projection of either 1 or 0. These three simple physical limits, i.e., parallel: 1; orthogonal: 0; surface isotropic: 1/2, are, respectively, the relevant limits for N_x², N_y², and N_z². On the pole of the sphere antipodal to the pipette, for instance, N_x² and N_y² = 0.5, and N_z² = 0. Toward the center of the sphere, N_z² increases to near 0.5, consistent with thermal averaging of the lowest order isotropy of a sixfold structure. Importantly, over most of the projection, N_z² is very nearly 1, whereas N_x² and N_y² almost vanish. Such a set of simulation results unambiguously establishes a quantitative basis for experimental assessment of chain alignment induced by deformation.

	CONCLUSIONS

Top Abstract Introduction Conclusions References

In paper I, the geometrical and elastic properties of polymer chain networks at the intended level of spectrin persistence length were determined. Although the large number of chain segments makes such models computationally unwieldy when considered for the simulation of whole cells, the present paper shows that the segmented chains can be faithfully replaced by effective potentials that substantially reduce the number of degrees of freedom in the models and permit the simulation of whole cells on conventional workstations. The effective potentials include a worm-like chain two-body term representing the individual spectrin molecules, and a three-body term representing the steric interaction between different chains. These terms provide a better description of the polymer chain models than does a network of Hooke's law springs, because the latter spring networks possess two instabilities (in collapse and expansion) that are absent from the full chain networks of paper I.

The three parameters of the effective potential can be reduced to two by demanding that the network energy be a minimum at a particular value of the intervertex separation. The remaining two parametersthe persistence length and maximum extension of the chainsare then fit by the condition that the stress dependence of the area in the effective network should reproduce that of the polymer chain network. The fitting procedure uses a mean-field approach for computational efficiency, although it is subsequently verified that the parameters from the fit are not sensitive to the mean-field approximation used in their determination. There is a unique parameter set for each polymer chain model, and the values of the parameters are found to be in the range expected, given their interpretation as persistence lengths and maximum extensions. Thus the microscopic ingredients of the effective potentialsentropic elasticity and steric repulsion of the chainsare supported by the rough agreement of the parameters with expectations.

As a first application, the effective potentials are employed to simulate the micropipette aspiration of erythrocytes. Two separate simulation codes were developed, mirroring the two principal aspiration techniques and necessitated by the different relaxation modes of the computational networks. In the free shape ensemble, networks based upon effective representations of all three polymer chain models do a credible job of reproducing the aspiration of flaccid erythrocytes, as measured by Waugh and Evans (1979). This result is not surprising, given that the three polymer chain models have shear moduli in the range estimated from flaccid cell experiments. Furthermore, the shape of the aspirated section of the cell within the micropipette was found to be cylindrical and capped with a hemisphere as confirmed