Deformation and dynamics of red blood cells in ﬂ ow through cylindrical microchannels †

The motion of red blood cells (RBCs) in microcirculation plays an important role in blood ﬂ ow resistance and in the cell partitioning within a microvascular network. Di ﬀ erent shapes and dynamics of RBCs in microvessels have been previously observed experimentally including the parachute and slipper shapes. We employ mesoscale hydrodynamic simulations to predict the phase diagram of shapes and dynamics of RBCs in cylindrical microchannels, which serve as idealized microvessels, for a wide range of channel con ﬁ nements and ﬂ ow rates. A rich dynamical behavior is found, with snaking and tumbling discocytes, slippers performing a swinging motion, and stationary parachutes. We discuss the e ﬀ ects of di ﬀ erent RBC states on the ﬂ ow resistance, and the in ﬂ uence of RBC properties, characterized by the F¨oppl – von K´arm´an number, on the shape diagram. The simulations are performed using the same viscosity for both external and internal ﬂ uids surrounding a RBC; however, we discuss how the viscosity contrast would a ﬀ ect the shape diagram.


Introduction
The behavior of so mesoscopic particles (e.g., polymers, vesicles, capsules, and cells) in ow has recently received enormous attention due to the wide range of applications of such suspensions and their rich physical properties. 1 To better control and/or manipulate the suspension's properties, for example in lab-on-chip applications, 2,3 a deeper understanding of the interplay among ow forces, elastic response and dynamics of so objects is required. An important example is the motion of red blood cells (RBCs) in microcirculation, which inuences many vital processes in microvasculature; 4 however, similar mechanisms and phenomena encompass a much wider class of capsule suspensions. 5,6 RBCs are extremely exible and experience strong deformation in microcirculation due to the ow and/or geometrical constraints. RBC deformation is important for the reduction of blood ow resistance 7-11 and for ATP release and oxygen delivery. 12,13 RBCs in microcirculation may attain various shapes including parachutes and slippers. 7,[14][15][16][17][18][19][20][21] Parachutes are characterized by a rather symmetric shape resembling a semispherical cap and are located at a position near the tube center. Slippers correspond to asymmetric RBC shapes, and therefore their membranes are typically in motion (e.g., tank-treading). Both non-centered slipper 20 and centered slipper shapes 17 have been observed experimentally, where the latter may only differ slightly from parachute shapes. The stable slipper shapes are well established at higher hematocrits due to hydrodynamic cell-cell interactions. 14,15,19 However, it is still not fully clear whether slippers are stable or transient states for single cells in ow. 7,9,17,20,21 The most convincing evidence so far comes from simulations in two dimensions (2D). 9,10,21 Simulations of 2D vesicles in slit channels have shown the existence of stable parachutes, slippers, and a snaking dynamics of discocytes 9,10 -where snaking refers to an oscillating RBC dynamics near the tube center. A phase diagram of various shapes was predicted, depending on relative connement and ow rate. Simulations of single RBCs in three dimensions (3D) 7,8,11,18 have been restricted so far to a limited number of studies, which only reveal (except for uctuations and transient states) stationary parachutes and discocytes. It is important to note that these shapes (averaged over thermal uctuations) are characterized by different symmetry classes, ranging from cylindrical symmetry (parachutes) to a single mirror plane containing the capillary axis and the RBC center (slippers). This raises several important questions: are slipper shapes also stable in 3D capillary ow? Is snaking dynamics around the center line possible in cylindrical microchannels in 3D? Do thermal uctuations destroy the regular snaking oscillations? What is the role of the membrane shear modulus (absent in 2D) in the phase diagram?
In this paper, we present a systematic study of single RBCs owing in microchannels. We construct diagrams of RBC shapes for different ow conditions and analyze RBC deformation. Changes in RBC properties (e.g., shear elastic modulus, bending rigidity) are also considered, since they are of importance in various blood diseases and disorders. 22 The presented 3D shape diagrams describe RBC deformation in microchannels, which mimic small vessels in microcirculation, and show that the parachute shape occurs mainly in small channels, while in large channels the slipper shape may occur. All simulations assume the same viscosity for both cytosol and suspending uid; however, we discuss how the viscosity contrast would affect the shape diagram. We also compare the 3D results with the 2D shape diagrams in ref. 9 and 10 and emphasize essential differences. For instance, due to membrane shear elasticity, there exists a region of RBC tumbling in 3D, which is absent in 2D. Finally, the simulation results are compared to available experimental data. 17,20,23 The article is organized as follows. In Section 2 we briey describe the mesoscopic method employed for uid ow, a RBC model, and simulation setup and parameters. In Section 3 we rst present the shape diagram for cell parameters typical of a healthy RBC. Then, RBC membrane properties are varied in order to elucidate their effects on the RBC shape and dynamics in microchannel ow. We also discuss the effect of viscosity contrast between external and internal uids on the shape diagram and the effects of different shapes on the ow resistance. We conclude briey in Section 4.

Models & methods
To model tube ow, we employ the smoothed dissipative particle dynamics (SDPD) method 24 for the suspending uid. SDPD is a particle-based mesoscopic simulation technique, where each SDPD particle represents a small volume of uid rather than individual atoms or molecules. The RBC membrane is represented by a triangulated network model 7,8,25,26 and coupled to a uid through friction forces.

Red blood cell model
The RBC membrane is modeled by a triangulated network of springs, 7,8,25,26 which includes elastic, bending, and viscous properties. A RBC is represented by a collection of N v particles connected by N s ¼ 3(N v À 2) springs with the potential where l j is the length of the j-th spring, l m is the maximum spring extension, x j ¼ l j /l m , x is the persistence length, and k p is the spring constant. Note that this spring denition allows us to dene a nonzero equilibrium spring length l 0 . Then, we employ a stress-free model for the membrane 8 so that each spring has its own l 0j set to the spring length of an initially triangulated cell surface. In addition to the elastic contribution in eqn (1), each spring may also have dissipative and random force terms 8 in order to incorporate membrane viscosity. The bending rigidity of a membrane is modeled by the bending energy where k b is the bending constant, q j is the instantaneous angle between two adjacent triangles having the common edge j, and q 0 is the spontaneous angle, which is set to zero in all simulations. Finally, to maintain a constant cell area and volume which mimic area-incompressibility of the lipid bilayer and incompressibility of the inner cytosol, we introduce two potentials where k a , k d , and k v are the global area, local area, and volume constraint coefficients, respectively. A and V are the instantaneous cell area and volume, while A j is the instantaneous area of an individual face within a triangulated network. A r , A 0 j , and V r are the desired total RBC area, area of the j-th face (set according to the initial triangulation), and total RBC volume, respectively.
To relate the model parameters in the spring potential (1) (e.g., x, k p ) and the bending potential (2) to the macroscopic membrane properties (e.g., Young's modulus Y r , bending rigidity k r ), we use analytic relationships derived for a regular hexagonal network. 8,27 The ratio l m /l 0 is set to 2.2 for all springs. 26 To relate simulation parameters to the physical properties of RBCs, we need a basic length and energy scale. Therefore, we dene an effective RBC diameter D r ¼ ffiffiffiffiffiffiffiffiffiffi A r =p p with A r being the RBC membrane area. From A r we can also calculate the average bond length l 0 for a given number of is the total number of triangular elements on a membrane. Experimental results 28 for the RBC area imply that D r ¼ 6.5 mm. Table 1 summarizes the parameters for the RBC model in units of D r and the thermal energy k B T, and the corresponding average values for a healthy RBC in physical units. The global area (k a ) and volume (k v ) constraint coefficients are chosen large enough to approximate closely the area-incompressibility of the lipid bilayer and incompressibility of the inner cytosol. Finally, a relationship for time scale is based on the characteristic RBC relaxation time, which is dened further below in the text.

Smoothed dissipative particle dynamics
SDPD 24 is a mesoscopic particle method, which combines two frequently used uid-dynamics approaches: the smoothed particle hydrodynamics 29,30 and dissipative particle dynamics 31,32 methods. The SDPD system consists of N point particles of mass m i , position r i and velocity v i . SDPD particles interact through three pairwise forces: conservative (C), dissipative (D), and random (R), so that the force on particle i is given by where e ij ¼ r ij /|r ij | and v ij ¼ v i À v j . p i and p j are particle pressures assumed to follow the equation of state p ¼ p 0 (r/r 0 ) a À b, where p 0 , r 0 , a, and b are selected parameters. The particle density is calculated locally and dened as r i ¼ 1 À r r c 2 and w ij ¼ w(r ij ). The (distancedependent) coefficients g ij and s ij dene the strength of dissipative and random forces and are equal to g ij ¼ 5h 0 3 w ij r i r j and where h 0 is the desired uid's dynamic viscosity and k B T is the thermal energy unit. The notation tr[dW ij ] corresponds to the trace of a random matrix of independent Wiener increments dW ij , and dW -S ij is the traceless symmetric part.
The time evolution of velocities and positions of particles is determined by Newton's second law of motion The above equations of motion are integrated using the velocity-Verlet algorithm. The SDPD uid parameters are given in Table 2. A natural length scale in the uid is the cut-off radius r c ; however, since we investigate the dependence of uid properties on r c , we use the membrane bond length l 0 instead, which is very similar in magnitude to r c . In addition, the exponent a in the equation of state is chosen to be a ¼ 7, and r 0 ¼ n, where n is the uid's number density (in particles per l 0 3 ). A relatively large value of a provides a good approximation of uid incompressibility, since even small changes in local density may lead to strong local pressure changes. Furthermore, the speed of sound, c, for the selected equation of state can be given as c 2 ¼ p 0 a/r 0 . The corresponding Mach numbers have been kept below 0.1 in all simulations providing a good approximation for an incompressible uid ow.
To span a wide range of ow rates, we employed different values of uid viscosities h in simulations with an input parameter h 0˛ [ 15; 120] ffiffiffiffiffiffiffiffiffiffiffiffi mk B T p =l 0 2 (m is the uid's particle mass), since the uid viscosity modies linearly the RBC relaxation time scale dened further below. Large values of viscosity were used to model high ow rates of the physical system in order to keep the Reynolds number in simulations, based on characteristic RBC size, sufficiently low (see also an argument at the end of Section 2.3). Even though we can directly input the desired uid viscosity h 0 in SDPD, the assumption that h 0 equals the actual uid viscosity h is reliable only when each SDPD particle has large enough number of neighboring particles, which may require a large enough cutoff radius and/or a density of uid particles. Therefore, it is always advisable to calculate the uid viscosity directly (e.g., in shear-ow setup) to check validity of the approximation of the simulated uid viscosity by h 0 . Note that for the uid viscosity we have always used the precalculated values of h rather than input values of h 0 . The tube wall has been modeled by frozen particles which assume the same structure as the uid, while the wall thickness is equal to r c . Thus, the interactions of uid particles with wall particles are the same as the interactions between uid particles, and the interactions of a RBC with the wall are identical to those with a suspending uid. The wall particles also provide a contribution to locally calculated density of uid particles near a wall, while the local density of wall particles is set to n. To prevent wall penetration, uid particles as well as vertices of a RBC are subject to reection at the uid-solid interface. We employed bounce-back reections, because they provide a better approximation for the no-slip boundary conditions in comparison to specular reection of particles. To ensure that no-slip boundary conditions are strictly satised, we also add a tangential adaptive shear force 33 which acts on the uid particles in a near-wall layer of a thickness r c . Table 1 RBC parameters in units of the effective RBC diameter D r and the thermal energy k B T, and the corresponding average values for a healthy RBC in physical units. N v is the number of membrane vertices, A r is the RBC membrane area, l 0 is the average bond length, V r is the RBC volume, T is the temperature, Y r is the membrane Young's modulus, k r is the membrane bending rigidity, and k d , k a , and k v are the local area, global area, and volume constraint coefficients, respectively. In all simulations, we have chosen A r ¼ 133.5 and k B T ¼ 0.4, which implies that D r ¼ 6.5 and l 0 ¼ 0.4

RBC parameters
Scaled units Physical units Table 2 SDPD fluid parameters in simulation and physical units. Mass and length for the SDPD fluid are measured in units of the fluid particle mass m and the membrane bond length l 0 . p 0 and b are parameters for the pressure equation, and h is the fluid's dynamic viscosity. In all simulations, we have set m ¼ 1, l 0 ¼ 0.4, and the thermal energy Fluid parameters Scaled units Physical units Coupling between the uid ow and RBC deformation is achieved through viscous friction between RBC nodes and surrounding uid particles, which is implemented via dissipative particle dynamics interactions. 8 Each membrane vertex interacts with uid particles within a spherical volume with a radius r 0 c using dissipative and random forces similar to those in SDPD. The strength of dissipative (friction) coupling depends on the uid viscosity and particle density as well as on the choice of r 0 c . The RBC membrane also separates inner and outer uids, which is implemented through bounce-back reections of uid particles on a membrane surface. 8 Finally, the local density of uid particles near the membrane includes contributions of both inner and outer uid particles.

Simulation setup
The simulation setup consists of a single periodic tube-like channel characterized by a diameter D and the length L ¼ 10D r , lled with a uid and a single suspended RBC. For simplicity, the uid viscosity inside a RBC is set to be the same as that of blood plasma. The ow is driven by a constant force f applied to each uid particle, which is equivalent to a constant pressure gradient DP/L ¼ fn, where DP is the pressure drop. To characterize the ow we dene a non-dimensional shear rate given by is the average shear rate (or pseudoshear rate), v ¼ Q/A is the average ow velocity with crosssectional area A ¼ pD 2 /4 and volumetric ow rate Q ¼ pD 4 fn/ (128h), and s ¼ hD r 3 /k r is a characteristic relaxation time of a RBC. Note that we dene c g based on the Poiseuille ow solution for a Newtonian uid without a RBC, since a single RBC does not signicantly affect the total ow rate. This assumption for Q also results in _ g* to be proportional to the pressure drop, which would be a convenient parameter to control in experiments. Furthermore, the ow behavior is determined by the cell connement c ¼ D r /D and the Föppl-von Kármán number G ¼ Y r D r 2 /k r ¼ 2662 (average value for a healthy RBC), which characterizes relative importance of cell elasticity to bending rigidity.
To interpret the non-dimensional shear rate with respect to experimental measurements, we can compute the characteristic RBC relaxation time from eqn (6) to be s ¼ 1.1 s. Thus, the pseudo-shear rate c g used in experiments is roughly equivalent in magnitude to _ g* in inverse seconds. It is important to note that since we employ distinct viscosity values in different simulations, the RBC relaxation time in simulation units also changes. Therefore, the same shear rate in simulations with different viscosities corresponds to different shear rates in physical units. This approach allows us to keep the Reynolds number low in the simulations, while a large range of shear rates in physical units can be spanned.

Sensitivity of simulation results to the discretization of uid and RBC membrane
A too coarse discretization of uid and RBC membrane may affect simulation results. To check whether our RBC discretization is ne enough, a number of simulations have been performed using signicantly different numbers of membrane vertices, N v ¼ 1000 and N v ¼ 3000. The comparison reveals that N v ¼ 1000 is sufficient to obtain accurate results for the investigated range of ow rates. Much larger ow rates may require ner RBC discretization due to strong membrane deformation.
Another potential source of error arises from the discretization of uid ow. There are two main parameters here, which are related to each other: the particle density n and the cutoff radius r c within the SDPD uid. The value of r c cannot be arbitrarily small, since the SDPD method properly functions only if each particle has a large enough number of neighboring particles. Thus, the choice of r c is directly associated with the particle density and can be selected smaller in magnitude for higher number densities. To study the sensitivity of the results to uid discretization, the uid density has been varied between n ¼ 0.2l 0 À3 and n ¼ 0.8l 0 À3 , while the corresponding r c values were between 3.8l 0 and 2.3l 0 . Simulation results show that values of r c ( 3l 0 and n T 0.4l 0 À3 are small and large enough, respectively, to properly reproduce the ow around a RBC for the studied ow rates. Note that the cutoff radius does not directly reect strong local correlations, since local interactions are scaled by the weights w ij , which decay to zero at distance r c . Finally, coupling between the RBC and uid ow is also performed over a smoothing length r 0 c . Even though generally there are no restrictions on the choice of r 0 c , it has to be small enough to impose properly the coupling between RBC vertices and local uid ow. To test the sensitivity of our simulation results to the choice of this parameter, we varied the coupling radius between 1.2l 0 and 2.4l 0 . A comparison of simulation results indicated that r 0 c ( 1.9l 0 appears to be sufficient to obtain results independent of r 0 c for _ g* ( 100. All results in the paper are obtained using the discretization parameters which comply with the estimations made above.

Results and discussion
3.1 Shapes and dynamics of a healthy RBC Fig. 1 shows several RBC shapes for c ¼ 0.58 (corresponding to a channel diameter D ¼ 1.72D r ), which are typically encountered in microcirculatory blood ow; see also Movies S1-S4. † For slow ows, see Fig. 1(a), the RBC shape is similar to the biconcave discocyte shape in equilibrium. For higher ow rates, see Fig. 1(b), a slipper shape may be observed, which is characterized here by an off-center position within the tube so that the membrane displays a tank-treading motion due to local shear gradients resembling the tank-treading in shear ow. 34,35 At the highest ow rate, see Fig. 1(c), a parachute shape is obtained, where the cell ows at the channel center and the membrane is practically not moving in the lab frame. Fig. 2 presents our main result, the shape diagram for different ow rates and connements, where the cell parameters are similar to those of a healthy RBC. The parachute RBC shape is predominantly observed in the region of strong connements and high enough ow rates, where large ow forces are able to strongly deform a RBC. Here, it should be noticed that parachutes are also stable for weak or no connement when the curvature of the parabolic ow in the center exceeds a critical value. 9,36 At weak connements, we nd offcenter slippers with tank-treading motion for higher ow rates, and discocytes with tumbling motion for lower ow rates. Both regions arise from the transition from strongly deformed parachute to more relaxed (discocyte and slipper) shapes, similar to the transition seen in the diagram for 2D vesicles. 9,10 However, the boundary between slippers and discocytes is governed by the critical shear rate _ g * ttt of the tumbling-to-tanktreading transition of a RBC; 37,38 tumbling occurs off the tube center, when the local shear rate drops below _ g * ttt . In the case of small viscosity contrast between inner and outer uids (equal to unity here), the origin of the tumbling-to-tank-treading transition is the anisotropic shape of the spectrin network, which requires stretching deformation in the tank-treading state, 37,38 and therefore cannot be captured by simulations of 2D vesicles. In addition, near the tumbling-slipper boundary, tumbling motion of a RBC exhibits a noticeable orbital dri so that the tumbling axis is not xed and oscillates in the vorticity direction (see Movie S2 †). This effect is qualitatively similar to a rolling motion (also called kayaking) found in experiments 39 and in simulations 40 of a RBC in shear ow. Orbital oscillations of a tumbling RBC are attributed to local membrane stretching deformation due to small membrane displacements whose effect becomes reduced if a RBC transits to a rolling motion. 39 At small shear rates _ g*, there also exists a so-called snaking region, rst observed for 2D vesicles in ref. 9 and 10, where a RBC performs a periodic oscillatory motion near the center line. In contrast to snaking in 2D, the snaking motion in 3D is fully three dimensional and exhibits an orbital dri (see Movie S1 †), which is similar to that for a RBC rolling motion in shear ow occurring in a range of shear rates between RBC tumbling and tank-treading. 39,40 The origin of orbital oscillations in the snaking regime might be similar to that for a rolling RBC; however, this issue requires a more detailed investigation. Note that at very low _ g* ( k B T/k r , the rotational diffusion of RBCs becomes important, and RBC dynamics is characterized by random cell orientation. Another striking difference between the phase diagrams in Fig. 2 and in ref. 9 and 10 is that at high connements the "conned slipper" found in the 2D vesicle simulations is suppressed in 3D. The conned slipper in 2D found for c T 0.6 is qualitatively similar to a slipper at low connements, which is also called "unconned slipper" in ref. 9 and 10, since this vesicle state exists in unbound parabolic ow. Note that the regions of conned and unconned slippers in 2D have no common boundary. The absence of slippers at high connements in 3D is due to the cylindrical shape of a channel, which would cause the conned slipper to conform to the wall curvature, which is energetically unfavorable. To better understand the differences between various RBC states, we now analyze the cell orientational angle, displacement from the channel center, and asphericity. The RBC orientational angle is dened as an angle between the eigenvector of the gyration tensor corresponding to the smallest eigenvalue (RBC thickness) and the tube axis. The RBC displacement r is computed as a distance between the RBC center of mass and the tube center. The RBC asphericity characterizes the deviation of a cell from a spherical shape and is dened as [(l 1 À l 2 ) 2 + (l 2 À l 3 ) 2 + (l 3 À l 1 ) 2 ]/(2R g 4 ), where l 1 # l 2 # l 3 are the eigenvalues of the gyration tensor and R g 2 ¼ l 1 + Fig. 1 Simulation snapshots of a RBC in flow (from left to right) for c ¼ 0.58. (a) A biconcave RBC shape at _ g* ¼ 5; (b) an off-center slipper cell shape at _ g* ¼ 24.8; and (c) a parachute shape at _ g* ¼ 59. 6. See also Movies S1-S4. † Fig. 2 A phase diagram for G ¼ 2662 (Y r ¼ 18.9 Â 10 À6 N m À1 , k r ¼ 3 Â 10 À19 J), which mimics average membrane properties of a healthy RBC. The plot shows various RBC dynamics states depending on the flow strength characterized by _ g* and the confinement c. The symbols depict performed simulations, with the RBC states: parachute (green circles), slipper (brown squares), tumbling (red diamonds) and snaking (blue stars) discocytes. The phase-boundary lines are drawn schematically to guide the eye. l 2 + l 3 . The asphericity for a single RBC in equilibrium is equal to 0.15. Fig. 3 presents the temporal dependence of these properties for different RBC states, including snaking and tumbling discocyte, slipper, and parachute. For the snaking dynamics, the orientational angle oscillates between 40 and 90 degrees ( Fig. 3(a)), the cell remains close to the channel center ( Fig. 3(b)), and it shows only slight deformation compared to the equilibrium shape (Fig. 3(c)). The parachute shape is characterized by a small orientational angle (aligned with the tube axis), a cell position right in the tube center, and a small asphericity which indicates that the RBC shape attains a more spherical shape. Both tumbling discocytes and tank-treading slippers are displaced further from the channel center ( Fig. 3(b)) than snaking discocytes and parachutes, and show an oscillating orientational angle. However, tumbling discocytes clearly show cell rotations, while slippers display a swinging motion characterized by small orientational oscillations around the tank-treading axis. 37 Moreover, a tumbling RBC does not experience strong deformation (Fig. 3(c)), while a slipper shows large oscillations in cell asphericity. Note that to determine the RBC shape under given conditions, we used both visual assessment of the corresponding RBC shapes in ow and the analysis of the characteristics discussed above.

Comparison with experiments
There exist several experimental studies of a RBC in microchannel ow. 17,20,23,41 In the experiments of ref. 23, the rotation of single RBCs and of their rouleaux structures in tube ow with radii ranging from 30 mm to 100 mm has been investigated. Even though the tube diameters in the experiments were larger than those used in our simulations, these experiments provide direct evidence of the existence of RBC off-center tumbling dynamics for low connements and low ow rates with c g ( 50 s À1 , in agreement with the simulation results in Fig. 2.
In the experiments of ref. 17, the imposed ow velocities were very large, ranging from 1 cm s À1 to 30 cm s À1 in a capillary with the diameter of 9 mm. This is much faster than the typical ow velocities in microcirculation, where for venules and arterioles with a similar diameter ow velocities are in the range of 0.2-7 mm s À1 . 4,42 The range of ow velocities we span in our simulations is about 0.2-1.0 mm s À1 , and is therefore comparable with that in microcirculation. The results in ref. 17 show the existence of parachute and slipper shapes, where a weak connement favors non-centered slipper shapes, which is in qualitative agreement with the simulations in Fig. 2. Furthermore, a good agreement between experiments and simulations is found for low connements, where parachute shapes are observed for ow rates with the velocities lower than 4-7 cm s À1 . 17 At ow velocities larger than approximately 7 cm s À1 and at low connements, centered slippers are observed which resemble parachutes, but become slightly asymmetric. 17 Recent 2D simulations 10 have also found that at high enough ow rates centered slippers and parachutes may coexist. Currently, we are not able to reach such high ow rates in 3D simulations due to numerical limitations. In addition, this region might be of limited interest, since the corresponding ow rates are far beyond the physiologically relevant values.
Experimental data in ref. 20 were obtained for narrow capillaries with diameters ranging from 4.7 mm to 10 mm; however, the ow velocities are considerably smaller than those in ref. 17, in the range of 1-40 mm s À1 . For strong connements c T 1, centered bullet-like shapes are observed which resemble an elongated cylindrical shape with a semi-spherical cap at the front end and a semi-spherical dip at the rear end. A comparison with our simulation results for the weaker connement of c ¼ 0.65 shows a good agreement since centered parachute shapes are found in both experiments and simulations. The transition to the parachute shape for c ¼ 0.65 occurred at the RBC velocity of about 0.5 mm s À1 in the experiments, 20 while our simulations predict the transition velocity of about 0.45 mm s À1 .

The effects of membrane properties on the shape diagram
In order to investigate the effects of membrane elastic parameters on the RBC phase diagram in capillary ow, we calculated Fig. 3 Characteristics of different RBC states: blue trianglesnaking discocyte ( _ g* ¼ 9.9, c ¼ 0.72), red diamondtumbling discocyte ( _ g* ¼ 14.9, c ¼ 0.44), brown squareslipper ( _ g* ¼ 49.7, c ¼ 0.44), and green circleparachute ( _ g* ¼ 64.6, c ¼ 0.65). (a) Cell orientational angle between the eigenvector of the gyration tensor corresponding to the smallest eigenvalue (RBC thickness) and the tube axis. (b) Distance between the RBC center of mass and the tube center normalized by D r . (c) RBC asphericity, which characterizes the deviation from a spherical shape. The asphericity is defined as [(l 1 À l 2 ) 2 + (l 2 À l 3 ) 2 + (l 3 À l 1 ) 2 ]/(2R g 4 ), where l 1 # l 2 # l 3 are the eigenvalues of the gyration tensor and R g 2 ¼ l 1 + l 2 + l 3 . The asphericity for a single RBC in equilibrium is equal to 0.15. See also Movies S1-S4. † phase diagrams for both reduced Young's modulus and increased bending rigidity. The state diagram for membrane bending rigidity increased by a factor of ve, which leads to k r ¼ 1.5 Â 10 À18 J and G ¼ 532, is shown in Fig. 4. We expect a dependence of the transition lines on the membrane parameters and geometry of the channel to be of the form _ g * c ¼Ũ(G, c, k r /k B T), whereŨ is a universal function for each transition. In addition, a more general form of _ g * c should also include the strength of thermal uctuations characterized by ambient temperature as an independent variable. As an example, different temperatures would affect the rotational diffusion of RBCs, which may become important at very low _ g* and lead to random cell orientation modifying potentially a snaking region. Also, recent experiments on vesicles 43 and a corresponding theory 44 suggest that the transition lines might be affected by thermal noise due to the sensitivity of a nonlinear dynamics to small perturbations near transition lines. A comparison of the results of Fig. 4 with those of Fig. 2 shows that for stiffer cells, the parachute and slipper/tumbling regions shi to lower values of _ g*. This behavior is consistent with the roughly linear dependence of the shear rate c g at the parachute-to-discocyte transition on RBC bending rigidity and shear modulus reported in ref. 7, which is equivalent to _ g* ¼ c 1 (c) + c 2 (c)G. Similarly, the snaking region shrinks towards lower _ g* values. Fig. 5 shows the state diagram for a RBC in tube ow with the membrane Young's modulus reduced by the factor of three (Y r ¼ 6.3 mN m À1 ) in comparison to that of a healthy RBC. This diagram should be compared with Fig. 2 and 4, which are for healthy RBCs and cells with an increased bending rigidity, respectively. In Fig. 5, the transition from snaking discocytes and swinging slippers to the parachute shape, as well as the transition from tumbling discocytes to swinging slippers occur at lower ow rates than those for a healthy RBC (Fig. 2). As a consequence, the snaking region shrinks substantially and is observed primarily at very low _ g*. Another feature is that the tumbling region is signicantly reduced in comparison to that in the diagram for healthy cells (Fig. 2). These results are consistent with the fact that the transition from discocyte to parachute shapes in Poiseuille ow depends roughly linearly on the RBC elastic properties and bending rigidity, 7 as well as that the tumbling-to-tank-treading transition of a RBC in shear ow depends nearly linearly on the RBC elastic properties such as Young's modulus. 37,38 This implies that the transition lines for a xed ambient temperature are linear functions of the Föpplvon Kármán number G, as discussed above. A comparison with the results of ref. 7 for RBCs with considerably smaller bending and shear rigidities indicates that strong thermal uctuations destroy the regular snaking oscillations.
Correspondence of different systems with a xed Föppl-von Kármán number G is also supported by the argument that the RBC relaxation time s is a linear function of RBC membrane elastic parameters (k r or equivalently Y r for a xed G). Thus, a simulation with the parameters {D r , Y r , sk r , _ g*} should lead to identical results as those obtained from a simulation with {D r , Y r /s, k r , _ g*/s}, where s is a scaling constant. Even though this argument is quite general, other characteristics of a system also need to be properly considered including thermal uctuations and Reynolds number of the ow. Finally, the assumption of linear dependence of s on the membrane properties may become invalid for strong enough ow rates, which may lead to non-linear membrane deformation. A semi-quantitative comparison of cell shape regions can be done by looking at the state diagram of Fig. 5 for the reduced membrane Young's modulus and the state diagram of Fig. 4 for an increased Fig. 4 A phase diagram for G ¼ 532 (Y r ¼ 18.9 Â 10 À6 N m À1 , k r ¼ 1.5 Â 10 À18 J) such that the bending rigidity of a cell membrane is five times larger than k r ¼ 3 Â 10 À19 J of a healthy RBC. The diagram depicts RBC behavior in tube flow with respect to the dimensionless shear rate _ g* and confinement c. The symbols correspond to performed simulations and RBC states are the parachute (green circles), slipper (brown squares), tumbling (red diamonds) and snaking (blue stars) discocytes. The phase-boundary lines are drawn schematically to guide the eye.
where the Young's modulus of a cell membrane is three times lower than that assumed for a healthy RBC. The plot shows different RBC dynamical states in tube flow with respect to the non-dimensional shear rate _ g* and the confinement c. The symbols depict performed simulations and RBC states include the parachute (green circles), slipper (brown squares), tumbling (red diamonds) and snaking (blue stars) discocytes. The phase-boundary lines are drawn schematically to guide the eye. membrane bending rigidity. The corresponding Föppl-von Kármán numbers are not the same, but similar in magnitude and both are considerably smaller than those for healthy RBCs. Thus, we expect that the two state diagrams for G ¼ 532 (Fig. 4) and G ¼ 887 (Fig. 5) should be similar, and show the same trends in comparison with the diagram for healthy cells (Fig. 2). This comparison further supports our argument about the existence of universal functions _ g * c to describe the transition lines.

The effect of cytosol viscosity on the shape diagram
We have focused in our simulations on the investigation of the effect of membrane elasticity on RBC dynamics in microcapillary ow, and have therefore employed same viscosity for the suspending uid and RBC cytosol. For a healthy RBC, the viscosity contrast l, dened as the ratio of cytosol over blood plasma viscosity, is approximately l ¼ 5. Therefore, we want to discuss briey the possible effect of the viscosity contrast on the dynamical states. Clearly, cytosol viscosity only plays a role when the internal uid (or equivalently membrane) is in motion, i.e., only when the cell tank-treads in the slipper state. Thus, we mainly need to discuss how the slipper region will be modied. Recent experiments 45 and simulations 46 indicate that the viscosity contrast indeed strongly affects the tumbling-totank-treading transition characterized by _ g * ttt . An increase in l is known to shi _ g * ttt of the tumbling-to-tanktreading transition to higher shear rates. 46 Therefore, we expect that the boundary between the tumbling and slipper regions in Fig. 2 would shi to higher values of _ g* for a real RBC leading to the expansion of the tumbling region. A very large viscosity of either RBC membrane or cytosol, which may occur in some bloodrelated diseases or disorders, 22 may also lead to a complete disappearance of the slipper region, which would be replaced by the RBC tumbling state. Recent simulations of 2D vesicles with l ¼ 1 in ref. 9 and l ¼ 5 in ref. 10 have shown an expansion of the snaking region toward the slipper region, since the RBC tanktreading becomes less favorable. Note that since RBC tumbling due to membrane elastic anisotropy is not possible in 2D, the snaking region in ref. 9 and 10 has a large common boundary with the slipper region, while in 3D the snaking region has practically no boundary with the dynamical slipper region and is mainly connected to the tumbling region (Fig. 2). Thus, no signicant changes in the snaking region due to the viscosity contrast is expected in 3D. In analogy with the expansion of the snaking region in 2D, an expansion of the tumbling region in 3D can be expected.
An effect of the viscosity contrast on the boundary between the slipper and parachute regions is not obvious. 2D simulations for different viscosity contrasts 9,10 indicate that the boundary is slightly altered that makes the slipper region get mildly expanded. Thus, it is plausible to expect a similarly weak effect in 3D; however, a denite statement on this issue is only possible aer a systematic numerical investigation of the effect of viscosity contrast has been performed in 3D.

Snaking, tumbling, and swinging frequencies
The simulations also provide interesting information about the dependence of the snaking, tumbling (both in the discocyte state) and swinging (in the slipper state) frequencies u on shear rate, as shown in Fig. 6 for three different connements extracted from the state diagram of Fig. 4. In all regimes, the frequencies increase linearly with shear rate. However, the prefactors of this linear dependence are very different. The strongest dependence is found in the swinging regime, somewhat larger than that in the snaking regime, while the frequency is almost independent of _ g* in the tumbling regime. Increasing connement signicantly reduces the snaking, tumbling and swinging frequencies. Note that the swinging frequencies are likely to be overpredicted in Fig. 6 in comparison with those of real RBCs due to the assumption of having the same viscosity of the suspending medium and RBC cytosol. Fig. 7 presents the volumetric ow rate Q RBC with a RBC measured in simulations and normalized by the ow rate Q without a RBC. Data are shown for different _ g* for the case of G ¼ 2662, corresponding to the shape diagram in Fig. 2. The volumetric ow rate at low connements remains essentially unaffected by the presence of a RBC. As the connement is increased, the ow rate decreases and the effect of a RBC on the ow rate appears to be stronger at low _ g* values. Note that Fig. 7 also presents the change in ow resistance, since the apparent viscosity is inversely proportional to the volumetric ow rate. Thus, the apparent viscosity increases with increasing connement. A reduction in the ow resistance with the increase of _ g* (represented by different curves in Fig. 7) is due to the transition from discocyte to parachute and slipper shapes, Fig. 6 Normalized snaking, tumbling, and swinging (slipper state) frequencies u of RBCs as a function of the dimensionless shear rate _ g*. The elastic membrane parameters are the same as in Fig. 4 (Y r ¼ 18.9 Â 10 À6 N m À1 , k r ¼ 1.5 Â 10 À18 J), i.e. the bending rigidity is five times larger than that used for a healthy RBC, and G ¼ 532. Three confinements are shown, as indicated. The symbols depict simulation results for various confinements. The regions of _ g* corresponding to different RBC states can be distinguished in the plot by colors and line typesswinging (brown, dashed line), tumbling (red, dotted line), and snaking (blue, solid line). since for large enough _ g* no signicant further changes in ow resistance occur.

Flow resistance
Even though Fig. 7 shows the effect of connement on Q RBC , the ow resistance for different c may be also affected by the change in hematocrit. In all simulations, the length of the channel was kept the same, while the tube diameter was varied, which means that the tube hematocrit (H t ), calculated as the ratio of the RBC volume to total tube volume, is inversely proportional to D 2 . For the strongest connement c ¼ 0.79 (D ¼ 8.2 mm) H t is equal to 0.027, while for the lowest connement c ¼ 0.37 (D ¼ 17.8 mm) H t ¼ 0.0057.

Conclusions
Even though the phase diagrams of RBC shapes show some qualitative similarities with the corresponding diagram of 2D vesicles, 9,10 there are several qualitative and quantitative differences. In 3D, slippers are essentially absent at high connements (c T 0.7) and low ow rates due to the cylindrical channel geometry, which requires the RBC slipper to comply with the channel wall curvature and cannot be captured by a 2D model. Therefore, parachute shapes are preferred at high cin agreement with the 3D results of ref. 7. At even higher connements than those in this study, RBCs are expected to attain bullet shapes. 20 Another evident difference between RBCs in 3D and vesicles in 2D 9,10 is the existence of a RBC tumbling region, which appears in 3D due to anisotropic elastic properties of RBCs. 37,38 For vesicles in 2D, the only possibility to include tumbling is to introduce a high enough viscosity contrast between the uids inside and outside a RBC.
The calculated state diagrams in 3D provide a better description of RBC shapes and dynamics in microvascular ow than the previous 2D results. The ow resistance is affected only weakly by the transition from discocyte to parachute and slipper shapes, most signicantly at high connements (c T 0.6). The geometrical complexity of microvasculature induces non-trivial partitioning of RBCs, which oen leads to very low cell volume fractions within various vessel structures, so that our simulation results provide an important step towards an understanding of blood ow and RBC behavior in microcirculation. Finally, similar physical mechanisms are expected for capsule suspensions. Fig. 7 The volumetric flow rate Q RBC with a RBC normalized by the flow rate Q without a RBC as a function of the confinement c for different _ g*. The Föppl-von Kármán number G ¼ 2662 is representative of healthy RBCs.