New Research In
Physical Sciences
Social Sciences
Featured Portals
Articles by Topic
Biological Sciences
Featured Portals
Articles by Topic
 Agricultural Sciences
 Anthropology
 Applied Biological Sciences
 Biochemistry
 Biophysics and Computational Biology
 Cell Biology
 Developmental Biology
 Ecology
 Environmental Sciences
 Evolution
 Genetics
 Immunology and Inflammation
 Medical Sciences
 Microbiology
 Neuroscience
 Pharmacology
 Physiology
 Plant Biology
 Population Biology
 Psychological and Cognitive Sciences
 Sustainability Science
 Systems Biology
Shape transitions of fluid vesicles and red blood cells in capillary flows

Edited by Tom C. Lubensky, University of Pennsylvania, Philadelphia, PA (received for review May 23, 2005)
Abstract
The dynamics of fluid vesicles and red blood cells (RBCs) in cylindrical capillary flow is studied by using a threedimensional mesoscopic simulation approach. As flow velocity increases, a model RBC is found to transit from a nonaxisymmetric discocyteto an axisymmetric parachute shape (coaxial with the flow axis), while a fluid vesicle is found to transit from a discocyte to a prolate ellipsoid. Both shape transitions reduce the flow resistance. The critical velocities of the shape transitions are linearly dependent on the bending rigidity and on the shear modulus of the membrane. Slipperlike shapes of the RBC model are observed around the transition velocities. Our results are in good agreement with experiments on RBCs.
The rheological properties of red blood cells (RBCs) are among the key factors for the flow resistance of blood in microvessels, because the hematocrit H (the volume fraction of RBCs) of normal human blood is ≈45%. Human RBCs have a biconcavedisk shape with a diameter of 8 μm under physiological conditions. Because of the large surface area S compared with a sphere of equal volume V, RBCs are easily deformed by external forces. RBCs are well known to form parachute shapes under flow in microvessels or glass capillaries with diameters between 3 and 13 μm (1–5). Some diseases, such as diabetes mellitus and sickle cell anemia, change the mechanical properties of RBCs; a reduction of RBC deformability was found to be responsible for an enhanced flow resistance of blood (6). Therefore, it is very important to understand the relation of RBC elasticity and flow properties in capillaries. Another motivation to study deformable vesicles and RBCs in capillary flow is the recent design of microfluidic devices, which allow the manipulation of small amounts of suspensions of cells and particles, with many applications in biology and medicine (3, 7, 8).
The cytoplasm of RBCs behaves as a Newtonian fluid because RBCs do not have a nucleus and other intracellular organelles. The membrane of RBCs consists of a lipid bilayer with an attached spectrin network as cytoskeleton. The lipid bilayer is an areaincompressible fluid membrane; its shape is controlled by the curvature energy with a constant volume (9). The bending rigidity of RBC has been measured by micropipette aspiration (10) and atomic force microscopy (11) to be approximately κ ≃ 2 × 10^{19} J (corresponding to κ ≃ 50k _{B} T). The shear modulus of the composite membrane, which is induced by the spectrin network, has been determined by several techniques; it is found to be μ ≃ 2 × 10^{6} N/m from optical tweezers manipulation (12), while the value μ ≃ 6 × 10^{6} N/m is obtained from micropipette aspiration (10). The relative importance of bending and stretching is determined by the dimensionless Föppl–von Kármán number , where K _{0} = 4μK _{A}/(μ + K _{A}) is the 2D Young modulus (with area compression modulus K _{A}) of the composite membrane, and R _{0} is the mean cell radius (13). Because the fluid membrane is nearly incompressible, K _{A} >> μ and thus K _{0} = 4μ. Therefore, RBCs are characterized by a Föppl–von Kármán number of γ ≃ 400. Interestingly, this value is not far from the boundary γ_{b} ≃ 150, which separates the regimes where the shape of crystalline vesicles is controlled mainly by bending and stretching energies, respectively (13).
Theoretically, the shapes of vesicles (9) and RBCs (14, 15) in the absence of flow have been calculated very successfully on the basis of a mechanical model of membranes, which includes both curvature and shear elasticity. In particular, it was shown recently that the full stomatocyte–discocyte–echinocyte sequence of RBCs can be reproduced by this model (15). We employ a very similar membrane model in our studies of RBCs in flow.
Two macroscopic quantities of bloods, the flow resistance and the hematocrit ratio H _{T}/H _{D}, are well investigated experimentally (16). The hematocrit H _{T} in a narrow tube is lower than the “discharge” hematocrit H _{D} (the hematocrit of the suspension collected at the end of the tube), because RBCs concentrate in the center of the tube, so that their mean velocity v _{ves} is higher than the mean blood velocity v _{m} (Fåhraeus effect). The deformations of RBCs (1–3) and fluid vesicles (17) in flow have been measured by optical microscopy. However, the quantitative relation of cell deformation and membrane properties is still unclear.
The deformation of RBCs and fluid vesicles in capillary flows were studied theoretically by lubrication theories (4, 5, 18) and boundary integral methods (19–21). In most of these studies, axisymmetric shapes, which are coaxial with the center of the capillary were assumed and cylindrical coordinates were used. The effect of nonaxisymmetric shapes has been investigated (22) with the assumption that the vesicle has a fixed shape corresponding to a convex ellipsoidal front part and an offcenter ellipsoidal indentation at the rear, to mimic the “slipper”like shapes observed experimentally (1, 2). It was found that the nonaxisymmetric shape does not change the dynamics of the parachutes very much. Very recently, the boundaryintegral method was extended to study deformable elastic capsules, without a resistance to surface dilation and bending deformation, in capillary flow (21). Biconcave discs with initial coaxial and noncoaxial orientations were found to deform into parachute and slipperlike shapes, respectively. However, this method is only applicable to slightly deformed capsules and shorttime dynamics due to an incipient buckling instability.
In this work, we investigate the transitions from nonaxisymmetric to axisymmetric shapes using a 3D simulation approach without numerical instability. Recently, several mesoscopic simulation techniques for fluid flow have been developed. The RBC aggregation in capillary vessels was simulated in ref. 23; RBCs were modeled as solid elastic bodies, in which dissipative particles are connected harmonically in a 3D mesh, and embedded in a particlebased solvent (23). Such a description can only be reasonable for small deviation from the discocyte rest shape and would completely miss the tanktreading motion of a RBC in shear flow. In contrast, to study the deformation of RBCs in flow, we model a RBC as an elastic capsule with bending rigidity, in which the 2D membrane is impenetrable to the inside and outside fluids. We denote this model as an “elastic vesicle” below. We also simulate the flow of fluid vesicles in capillaries, both in its own right and to distinguish the effects of shear and bending elasticity of the membrane. Thermal fluctuations are taken into account consistently in our simulations.
Methods
We employ a mesoscopic approach, which combines a particlebased hydrodynamics model (24–29) for the solvent and a coarsegrained, dynamically triangulated surface model (30) for the membrane. The mesoscale hydrodynamics method has been called multiparticle collision dynamics (27–29) or stochastic rotation dynamics (25, 26). We briefly outline here the simulation technique. The multiparticle collision dynamics method is explained in detail in refs. 24–29; the combination with the membrane model is described in refs. 31 and 32. We have applied this method to fluid vesicles in simple shear flow and have found that it reproduces the vesicle dynamics very well (31, 32). Here, we also employ a membrane model with shear elasticity, which is based on the RBC membrane model by Discher et al. (14).
Solvent. The solvent is described by N _{s} pointlike particles of mass m _{s}. The algorithm consists of alternating streaming and collision steps. In the streaming step, the particles move ballistically and the position of each particle r _{i} is updated according to r _{i} (t + h) = r _{i} (t) + v _{i}(t)h, where v _{i} is the velocity of particle i and h is the time interval between collisions. In the collision step, the particles are sorted into cubic boxes of lattice constant a. The collision step consists of a stochastic rotation of the relative velocities of each particle in a box, , where v _{cm} is the velocity of the center of mass of all particles in the box. The matrix Ω rotates velocities by the angle π/2 around an axis, which is chosen randomly for each box.
Fluid Vesicle. The fluid membrane is described by N _{mb} vertices, which are connected by tethers to form a triangular network. The vertices have excluded volume and mass m _{mb}. Soft pairwise potentials are used for the tetherbond and excluded volume (32). The average bond length is chosen to equal the lattice constant a of the collision boxes. The shapes and fluctuations of the membrane are controlled by curvature elasticity with the energy H _{cv} = (κ/2)∫(C _{1} + C _{2})^{2} dS, where κ is the bending rigidity, and C _{1} and C _{2} are the principal curvatures at each point of the membrane (33). The curvature energy is discretized as described in refs. 30 and 34. To model the fluidity of the membrane, tethers can be flipped between the two possible diagonals of two adjacent triangles.
To keep the volume V and the surface area S constant in the presence of external forces, we used a global volumeconstraint potential U _{V} = (1/2)k _{V}(V  V _{0})^{2}, and a local areaconstraint potential, U_{Si} = (1/2)k_{Si} (S_{i}  S _{0})^{2}, where S_{i} is the area of the elementary triangles of the network. Here, we use the parameters k_{Si} = 4,000k _{B} T, k _{V} = k _{B} T, and S _{0} = 0.41 a ^{2}. With these potentials, the volume V = 450a ^{3} and surface area S = 407 a ^{2} of a vesicle are kept constant with a deviation of <1% for all simulated systems. The reduced volume is V ^{*} = V/(4πR ^{3} _{0}/3) = 0.59, whereR _{0} = = 5.7a is the effective vesicle radius. At this reduced volume, a biconcave discocyte is the equilibrium shape, and a prolate ellipsoid and stomatocyte are metastable in the absence of flow (9, 31).
The solvent particles interact with the membrane in two ways. First, the membrane vertices are included in the multiparticle collision dynamics collision procedure (26, 32). Second, the solvent particles are scattered with a bounceback rule from membrane triangles (32). These interactions together ensure that the fluid satisfies a noslip boundary condition on the membrane. The fluids in the interior and exterior of the vesicle are taken to be the same, in particular to have the same viscosity η_{0}. The relative membrane viscosity is η_{mb}/η_{0} R _{0} = 3 (31). We focus on steady states, which depend only very little on the viscosities of the internal fluid and the membrane.
Vesicle with Shear Elasticity. The elastic membrane is modeled as a composite network, which consists of a dynamically triangulated surface as in the case of fluid vesicles, coupled to an additional network of harmonic springs with fixed connectivity (no bondflip). The same number of the bonds is used for both the fluid and the tethered networks. The harmonic potential (1/2)k _{el}(r _{i}  r _{j} )^{2} is used for the tethers. This tether network generates the shear modulus μ = .‡
We present data for various shear and bending moduli in the range 0.5 < k _{el} a ^{2}/k _{B} T < 16 (corresponding to 28 < μR ^{2} _{0}/k _{B} T < 900) and 5 ≤ κ/k _{B} T ≤ 40. Thereby, we cover a range of Föppl–von Kármán numbers from γ = 11 to γ = 360. The upper limit corresponds to the case of RBCs, where γ ≃ 400.
Capillary Flow. We use cylindrical capillaries with radius R _{cap}. The solvent interacts with the wall by the procedure described in refs. 27 and 28, which implies noslip boundary conditions. Periodic boundary conditions are used in the flow (z) direction. A gravitation force is used to generate flow (28). In experimental conditions of RBCs, the Reynolds number Re = ρv _{ves} R _{0}/η_{0} is very small, typically Re ≃ 10^{2}, where v _{ves} is the mean velocity of the vesicle. Therefore, we chose parameters such that Re < 1 in all simulations (29). We use the capillary radius R _{cap} = 8a so that R _{0}/R _{cap} = 0.71, the system length L _{z} = 80a, the solvent density ρ = 10m _{s}/a ^{3}, the viscosity of the solvent η_{0} = 20.1/a ^{2}, the number of membrane vertices N _{mb} = 500, and the mass of membrane particles m _{mb} = 10m _{s}. This value corresponds to a volume fraction of vesicles of H _{T} = V/(πR ^{2} _{cap} L _{z}) = 0.028. The results are scaled with the intrinsic time τ = η_{0} R ^{3} _{cap}/k _{B} T. With a RBC radius of R _{0} = 3.3 μm and the viscosity η_{0} = 10^{3} Pa s of water, we obtain R _{cap} = 4.6 μm and a characteristic velocity R _{cap}/τ = 0.2 μm/s. The error bars are estimated from three or four independent runs.
Results and Discussion
Both fluid and elastic vesicles retain their discoidal shapes in slow capillary flows (see Fig. 1A ). The vesicles turn such as to align the longest axis of the momentofinertia tensor with the flow (z) direction, even if their initial conformations are discocytes, which are coaxial with the capillary. The discoidal shape is elongated in the flow direction, and its front–rear symmetry is broken, but the biconcave dimples and the mirror symmetry with respect to the plane determined by the two eigenvectors of the momentofinertia tensor with the largest eigenvalues are retained. This shape is stable because it obstructs solvent flows less than a discocyte in coaxial orientation. The shapes and orientation of vesicles thermally fluctuate around the stable average shape. As the flow velocity increases, the fluid vesicle transits into a prolate ellipsoidal shape (see Fig. 2). The freeenergy barrier between discocyte and prolate is overcome by the elongation force of the flow field. A related transition also occurs in simple shear flow (31). Conversely, the elastic vesicle transits into a parachute shape (compare Figs. 1B and 2; see also Movie 1, which is published as supporting information on the PNAS web site). The shear elasticity prevents the elongation of the vesicle into a prolate shape in this case, because a very large shear deformation of the membrane would be required. The elastic vesicle often forms a nonaxisymmetric slipperlike shape around the transition velocity, as shown in Fig. 1C . RBCs with similar slipper shapes have been observed in microvessels (1, 2). When the flow velocity is reduced, the parachute shape jumps back into a discocyte.
The fluid vesicle can remain in a metastable prolate state below the transition velocity. However, its fluctuations in orientation and overall curvature increase, and some of them return to discocytes by means of a Vshaped conformation, with the apex pointing into the flow direction (see Movie 2, which is published as supporting information on the PNAS web site). In fast flows, a fluid vesicle also can attain a parachute shape when the initial conformation is a coaxial discocyte. When the flow is reduced or stopped, the parachute shape of fluid vesicles changes into a stomatocyte, which exists as metastable state for V ^{*} = 0.59 in the absence of flow. Thus, the parachute shape of fluid vesicles can be understood as a deformed stomatocyte. Because a shear flow induces a lift force on a vesicle near a wall (36, 37), all vesicles stay around the center of the cylindrical capillary.
The mean fluid velocities v _{m} at both shape transitions are found to increase linearly with the bending rigidity κ of the membrane (see Fig. 2 A ). Here, the mean fluid velocity is given by v _{m} = (1  H _{T})v _{sol} + H _{T} v _{ves}, where v _{sol} is the mean velocity of the exterior solvent. Because the data extrapolate to a very small flow velocity in the limit κ/k _{B} T → 0, we conclude that κ is the dominant factor in determining the transition velocity in the regime μR ^{2} _{0}/k _{B} T ≤ 110, corresponding to γ ≤ 90. At fixed bending rigidity κ/k _{B} T = 10, the transition velocity is found to increase also roughly linearly with increasing shear elasticity, as shown in Fig. 2B . This finding implies that in the investigated regime of elastic parameters, the flow forces required to induce a shape transformation are linearly related to the elastic bending and stretching forces, v ^{c} _{m}τ/R _{cap} = 0.1μR ^{2} _{0}/k _{B} T + 4κ/k _{B} T. Our data can therefore be extrapolated to other combinations of elastic parameters. In particular, it suggests that parachute shapes of RBCs should appear for flow velocities larger than 800R _{cap}/τ ≃ 0.2 mm/s under physiological conditions. This result is consistent with the experimental results of ref. 2 and is in the range of microcirculation in the human body.
Fig. 2B also shows that there is a metastable region, where discocytes are seen for increasing flow velocity, but parachute shapes are seen for decreasing flow velocity. This hysteresis becomes more pronounced with increasing shear modulus (at κ/k _{B} T = 10). We believe that this effect is due to a suppression of thermal fluctuations with increasing μR ^{2} _{0}/k _{B} T. The reason for this metastability, and the pathway from the parachute to the discocyte shape, can be extracted from the simulations. After the flow velocity is sufficiently reduced, the parachute shape becomes unstable to a rotation of its symmetry axis. As the orientation changes, the flow forces pull the vesicle into the discocyte shape shown in Fig. 1 A . The tails of slipper and parachute shapes are found to sharpen as μ increases or κ decreases.
To investigate the dependence of the vesicle shapes on the elastic moduli, we consider several characteristic quantities. The curvature of the rear part of nonaxisymmetric discoidal vesicles in slow flows is smaller than of their front part (see Fig. 3 and snapshots in Figs. 1 A and 2). The difference between rear and front ends becomes more pronounced in faster flows, while the shear elasticity reduces the deformation.
Fig. 4 shows the κ dependence of parachute shapes. The length ratios are l _{max}/2r _{max} ≃ 1 and l _{dmp}/l _{max} ≳ 0.5. Thus, a sufficiently large deviation from the discocyte is required to stabilize the parachute shape. The relative fluctuations of the lengths l _{max} and r _{max} are small, typically ≈3–4%. As κ increases, the deformations of both elastic and fluid vesicles are reduced: the vesicles are less elongated in the z direction, and their front ends become more smoothly rounded. The curvature in the center of the dimples of fluid vesicles decreases, because larger κ suppresses large curvatures. Conversely, the curvature of the dimple of an elastic vesicle increases. When the shear modulus μ is varied at bending rigidities κ/k _{B} T = 10 and 20, we find that the length ratios of parachute shapes in Fig. 4 depend very weakly and monotonically on μ. Also, the front curvature remains essentially unchanged. The strongest effect is seen in the rear curvature, which decreases with increasing μ. In fact, data for the rear (dimple) curvature follow a single scaling curve, when they are plotted as a function μR ^{2} _{0}/κ (see also Fig. 7, which is published as supporting information on the PNAS web site).
Fig. 5 shows the profile of the solvent velocity v_{z} in the capillary. In the absence of the vesicle, the parabolic velocity profile of Poiseuille flow is obtained. A vesicle modifies the flow field. The solvent inside the vesicle moves with constant velocity v _{ves} because of the impermeability of the membrane. This movement reduces v_{z} near the center of the capillary, r < 0.7R _{0}, which implies that the pressure at the rear end of the vesicle becomes higher than at the front end. Then, the flow velocity in the annulus between the vesicle and the wall becomes faster, and a flow in the radial direction is generated near the rear and front parts of the vesicle.
Fig. 6 shows the hematocrit ratio H _{T}/H _{D} and the flow resistance ΔP _{ves}/v _{m} in the capillary, the two quantities which are mainly studied experimentally. Because the vesicle stays in the center of the capillary, the vesicle velocity v _{ves} is larger than the mean fluid velocity v _{m}. This value generates the hematocrit reduction according to the relation H _{T}/H _{D} = v _{m}/v _{ves} (16). The pressure drop ΔP _{ves} per vesicle is calculated by , where v _{0} is the mean fluid velocity in absence of vesicle. This expression is obtained from the wellknown result of the pressure drop in simple Poiseuille flow, , where for a parabolic velocity profile the maximum velocity v _{max} at the center of capillary is related to the mean velocity v _{0} by v _{max} = 2v _{0}. In the presence of the vesicle, there is an extra pressure drop, so that . This equation can be written in the same form as for Poiseuille flow when an apparent viscosity η_{app} is introduced, so that , which implies η_{app} v _{m} = η_{0} v _{0}.
We find that the hematocrit ratio and the flow resistance decrease with increasing v _{m}, with jumps at the shape transitions.§ The reduction of both quantities can be understood by the reduction of the effective vesicle radius. We compared our results with the “axialtrain” or “stackedcoins” model (4, 38). This model assumes a coaxial cylindrical core with radius r _{eff}, which represents a line of RBC, to move as a rigid body. It predicts H _{T}/H _{D} = (1 + λ^{2})/2 and apparent viscosity η_{app} = η_{0}/(1  λ^{4}), where λ = r _{eff}/R _{cap} v (38). This equation implies immediately ΔP _{ves}/v _{m} ∼ λ^{4}. Although this model applies to high H _{T}, whereas our simulation are carried out at very low H _{T}, we find that these equations reproduce the dependences of H _{T}/H _{D} and ΔP _{ves}/v _{m} on the mean fluid velocity v _{m} surprisingly well (compare Fig. 6 Inset). Here, we choose r _{eff} = r _{max}, where r _{max} is maximum radial distance from the center of the capillary (see Fig. 4).
Ellipsoidal fluid vesicles with larger reduced volume V ^{*} > 0.8 have also been studied experimentally (17). Therefore, we have simulated prolate fluid vesicles at V ^{*} = 0.78 and 0.9. Our results qualitatively reproduce the observed deformations (compare figures 4 and 5 of ref. 17).
In summary, we have found steady states of nonaxisymmetric discoidal vesicles in slow capillary flows. We have explored the dependence of the transitions from nonaxisymmetric to axisymmetric parachute shapes, as well as the geometry of the parachutes, on the elastic properties of the composite membrane. Under physiological conditions, RBC aggregates are observed in slow flows. Our simulation method can easily be applied to multiple vesicles in a capillary, as well as to more complex flow geometries. It will be interesting to study the coupling of shape transitions and aggregation of vesicles in the future.
Acknowledgments
We thank D. M. Kroll, K. Mussawisade, M. Ripoll, G. Vliegenthart, and R. G. Winkler for helpful discussions. H.N.'s stay at Forschungszentrum Jülich was supported by the Japan Society for the Promotion of Science. G.G. was supported in part by the Deutsche Forschungsgemeinschaft through the priority program “Nano and Microfluidics.”
Footnotes

↵ † To whom correspondence should be addressed. Email: hi.noguchi{at}fzjuelich.de.

Author contributions: H.N. and G.G. designed research, performed research, and wrote the paper.

This paper was submitted directly (Track II) to the PNAS office.

↵ ‡ This value differs by a factor of 4 from the result of ref. 35 because we employ a harmonic potential with a minimum at vanishing distance.

↵ § The pressure difference between front and rear end of a vesicle generates a solventdensity gradient, because the multiparticle collision dynamics fluid obeys a idealgas equation of state. The density exhibits a maximum at the dimple of the parachute shape in our simulations and is ≈50% higher than the average density at v _{m}τ/R _{cap}=218. To check whether this method causes any artifacts, we have simulated parachuteshaped elastic and fluid vesicles with a twice as large solvent density, ρ = 20m _{s}/a ^{3}, where the relative density increase is only 25%. The peak of v _{z} in the annulus between the vesicle and the wall (see Fig. 5) is slightly lower, and ΔP _{ves} is ≈10% lower than for ρ = 10 m _{s}/a ^{3}. Thus, we overestimate ΔP _{ves} by ≈10... 20%. Because the vesicle shapes are not changed, the effect on the vesicle dynamics should be very weak, however.
 Copyright © 2005, The National Academy of Sciences
References

↵
Skalak, R. (1969) Science 164 , 717719. pmid:5778020
 ↵
 ↵
 ↵
 ↵
 ↵

↵
Shelby, J. P., White, J., Ganesan, K., Rathod, P. K. & Chiu, D. T. (2003) Proc. Natl. Acad. Sci. USA 100, 1461814622. pmid:14638939
 ↵

↵
Seifert, U. (1997) Adv. Phys. 46 , 13137.
 ↵
 ↵
 ↵
 ↵
 ↵

↵
Lim, H. W. G., Wortis, M. & Mukhopadhyay, R. (2002) Proc. Natl. Acad. Sci. USA 99 , 1676616769. pmid:12471152
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵

↵
Gompper, G. & Kroll, D. M. (2004) in Mechanics of Membranes and Surfaces, eds. Nelson, D. R., Piran, T. & Weinberg, S. (World Scientific, Singapore), 2nd Ed., pp. 359426.
 ↵
 ↵

↵
Helfrich, W. (1973) Z. Naturforsch. 28c, 693703.
 ↵
 ↵
 ↵
 ↵

↵
Whitmore, R. L. (1967) J. Appl. Physiol. 22 , 767771. pmid:6023191