Deformation and dynamics of red blood cells in flow through cylindrical microchannels

Abstract

The motion of red blood cells (RBCs) in microcirculation plays an important role in blood flow resistance and in the cell partitioning within a microvascular network. Different 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 confinements and flow rates. A rich dynamical behavior is found, with snaking and tumbling discocytes, slippers performing a swinging motion, and stationary parachutes. We discuss the effects of different RBC states on the flow resistance, and the influence of RBC properties, characterized by the Föppl–von Kármán number, on the shape diagram. The simulations are performed using the same viscosity for both external and internal fluids surrounding a RBC; however, we discuss how the viscosity contrast would affect the shape diagram.

---

1 Introduction

The behavior of soft mesoscopic particles (e.g., polymers, vesicles, capsules, and cells) in flow has recently received enormous attention due to the wide range of applications of such suspensions and their rich physical properties.¹ 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 flow forces, elastic response and dynamics of soft objects is required. An important example is the motion of red blood cells (RBCs) in microcirculation, which influences many vital processes in microvasculature;⁴ however, similar mechanisms and phenomena encompass a much wider class of capsule suspensions.^5,6 RBCs are extremely flexible and experience strong deformation in microcirculation due to the flow and/or geometrical constraints. RBC deformation is important for the reduction of blood flow 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–21 Parachutes are characterized by a rather symmetric shape resembling a semi-spherical 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²⁰ and centered slipper shapes¹⁷ 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 flow.^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 confinement and flow 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 fluctuations and transient states) stationary parachutes and discocytes. It is important to note that these shapes (averaged over thermal fluctuations) 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 flow? Is snaking dynamics around the center line possible in cylindrical microchannels in 3D? Do thermal fluctuations 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 flowing in microchannels. We construct diagrams of RBC shapes for different flow 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.²² 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 fluid; 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 briefly describe the mesoscopic method employed for fluid flow, a RBC model, and simulation setup and parameters. In Section 3 we first 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 flow. We also discuss the effect of viscosity contrast between external and internal fluids on the shape diagram and the effects of different shapes on the flow resistance. We conclude briefly in Section 4.

2 Models & methods

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

2.1 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 potentialwhere l_j is the length of the j-th spring, l_m is the maximum spring extension, x_j = l_j/l_m, ξ is the persistence length, and k_p is the spring constant. Note that this spring definition allows us to define a nonzero equilibrium spring length l₀. Then, we employ a stress-free model for the membrane⁸ 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⁸ in order to incorporate membrane viscosity. The bending rigidity of a membrane is modeled by the bending energywhere k_b is the bending constant, θ_j is the instantaneous angle between two adjacent triangles having the common edge j, and θ₀ 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 potentialswhere 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⁰_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., ξ, k_p) and the bending potential (2) to the macroscopic membrane properties (e.g., Young's modulus Y_r, bending rigidity κ_r), we use analytic relationships derived for a regular hexagonal network.^8,27 The ratio l_m/l₀ is set to 2.2 for all springs.²⁶ To relate simulation parameters to the physical properties of RBCs, we need a basic length and energy scale. Therefore, we define an effective RBC diameter with A_r being the RBC membrane area. From A_r we can also calculate the average bond length l₀ for a given number of membrane vertices where N_t = 2N_v − 4 is the total number of triangular elements on a membrane. Experimental results²⁸ for the RBC area imply that D_r = 6.5 μm. Table 1 summarizes the parameters for the RBC model in units of D_r and the thermal energy k_BT, 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 defined further below in the text.

Table 1 RBC parameters in units of the effective RBC diameter D_r and the thermal energy k_BT, 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₀ is the average bond length, V_r is the RBC volume, T is the temperature, Y_r is the membrane Young's modulus, κ_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_BT = 0.4, which implies that D_r = 6.5 and l₀ = 0.4

RBC parameters	Scaled units	Physical units
N v	1000
A r		133.5 × 10^−12 m^2
D r		6.5 × 10^−6 m
l 0	0.061Dr	3.93 × 10^−7 m
V r	0.34Dr^3	93 × 10^−18 m^3
T		310 K
Y r		18.9 × 10^−6 N m^−1
κ r	70kBT	3 × 10^−19 J
k d		4.3 × 10^−6 N m^−1
k a		2.1 × 10^−4 N m^−1
k v		220 N m^−2

---

2.2 Smoothed dissipative particle dynamics

SDPD²⁴ is a mesoscopic particle method, which combines two frequently used fluid-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 bywhere 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₀(ρ/ρ₀)^α − b, where p₀, ρ₀, α, and b are selected parameters. The particle density is calculated locally and defined as with weight function being the Lucy function, where r_c is the cutoff radius. Note that ∇W_L(r) = −rw(r) so that and w_ij = w(r_ij). The (distance-dependent) coefficients γ_ij and σ_ij define the strength of dissipative and random forces and are equal to and , where η₀ is the desired fluid's dynamic viscosity and k_BT 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 fluid parameters are given in Table 2. A natural length scale in the fluid is the cut-off radius r_c; however, since we investigate the dependence of fluid properties on r_c, we use the membrane bond length l₀ instead, which is very similar in magnitude to r_c. In addition, the exponent α in the equation of state is chosen to be α = 7, and ρ₀ = n, where n is the fluid's number density (in particles per l₀³). A relatively large value of α provides a good approximation of fluid 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² = p₀α/ρ₀. The corresponding Mach numbers have been kept below 0.1 in all simulations providing a good approximation for an incompressible fluid flow.

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₀. p₀ and b are parameters for the pressure equation, and η is the fluid's dynamic viscosity. In all simulations, we have set m = 1, l₀ = 0.4, and the thermal energy k_BT = 0.4

Fluid parameters	Scaled units	Physical units
p 0	16	1.07 Pa
b	12.8	0.86 Pa
η		1.2 × 10^−3 Pa s
k B T		4.282 × 10^−21 J

---

To span a wide range of flow rates, we employed different values of fluid viscosities η in simulations with an input parameter η₀ ∈ [15; 120] (m is the fluid's particle mass), since the fluid viscosity modifies linearly the RBC relaxation time scale defined further below. Large values of viscosity were used to model high flow 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 fluid viscosity η₀ in SDPD, the assumption that η₀ equals the actual fluid viscosity η 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 fluid particles. Therefore, it is always advisable to calculate the fluid viscosity directly (e.g., in shear-flow setup) to check validity of the approximation of the simulated fluid viscosity by η₀. Note that for the fluid viscosity we have always used the precalculated values of η rather than input values of η₀.

The tube wall has been modeled by frozen particles which assume the same structure as the fluid, while the wall thickness is equal to r_c. Thus, the interactions of fluid particles with wall particles are the same as the interactions between fluid particles, and the interactions of a RBC with the wall are identical to those with a suspending fluid. The wall particles also provide a contribution to locally calculated density of fluid particles near a wall, while the local density of wall particles is set to n. To prevent wall penetration, fluid particles as well as vertices of a RBC are subject to reflection at the fluid–solid interface. We employed bounce-back reflections, because they provide a better approximation for the no-slip boundary conditions in comparison to specular reflection of particles. To ensure that no-slip boundary conditions are strictly satisfied, we also add a tangential adaptive shear force³³ which acts on the fluid particles in a near-wall layer of a thickness r_c.

Coupling between the fluid flow and RBC deformation is achieved through viscous friction between RBC nodes and surrounding fluid particles, which is implemented via dissipative particle dynamics interactions.⁸ Each membrane vertex interacts with fluid particles within a spherical volume with a radius r^′_c using dissipative and random forces similar to those in SDPD. The strength of dissipative (friction) coupling depends on the fluid viscosity and particle density as well as on the choice of r^′_c. The RBC membrane also separates inner and outer fluids, which is implemented through bounce-back reflections of fluid particles on a membrane surface.⁸ Finally, the local density of fluid particles near the membrane includes contributions of both inner and outer fluid particles.

2.3 Simulation setup

The simulation setup consists of a single periodic tube-like channel characterized by a diameter D and the length L = 10D_r, filled with a fluid and a single suspended RBC. For simplicity, the fluid viscosity inside a RBC is set to be the same as that of blood plasma. The flow is driven by a constant force f applied to each fluid particle, which is equivalent to a constant pressure gradient ΔP/L = fn, where ΔP is the pressure drop. To characterize the flow we define a non-dimensional shear rate given bywhere = v̄/D = Dfn/(32η) is the average shear rate (or pseudo-shear rate), v̄ = Q/A is the average flow velocity with cross-sectional area A = πD²/4 and volumetric flow rate Q = πD⁴fn/(128η), and τ = ηD_r³/κ_r is a characteristic relaxation time of a RBC. Note that we define based on the Poiseuille flow solution for a Newtonian fluid without a RBC, since a single RBC does not significantly affect the total flow rate. This assumption for Q also results in [small gamma, Greek, dot above]* to be proportional to the pressure drop, which would be a convenient parameter to control in experiments. Furthermore, the flow behavior is determined by the cell confinement χ = D_r/D and the Föppl–von Kármán number Γ = Y_rD_r²/κ_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 τ = 1.1 s. Thus, the pseudo-shear rate used in experiments is roughly equivalent in magnitude to [small gamma, Greek, dot above]* 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.

2.4 Sensitivity of simulation results to the discretization of fluid and RBC membrane

A too coarse discretization of fluid and RBC membrane may affect simulation results. To check whether our RBC discretization is fine enough, a number of simulations have been performed using significantly 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 flow rates. Much larger flow rates may require finer RBC discretization due to strong membrane deformation.

Another potential source of error arises from the discretization of fluid flow. 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 fluid. 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 fluid discretization, the fluid density has been varied between n = 0.2l₀⁻³ and n = 0.8l₀⁻³, while the corresponding r_c values were between 3.8l₀ and 2.3l₀. Simulation results show that values of r_c ≲ 3l₀ and n ≳ 0.4l₀⁻³ are small and large enough, respectively, to properly reproduce the flow around a RBC for the studied flow rates. Note that the cutoff radius does not directly reflect 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 fluid flow is also performed over a smoothing length r^′_c. Even though generally there are no restrictions on the choice of r^′_c, it has to be small enough to impose properly the coupling between RBC vertices and local fluid flow. To test the sensitivity of our simulation results to the choice of this parameter, we varied the coupling radius between 1.2l₀ and 2.4l₀. A comparison of simulation results indicated that r^′_c ≲ 1.9l₀ appears to be sufficient to obtain results independent of r^′_c for [small gamma, Greek, dot above]* ≲ 100. All results in the paper are obtained using the discretization parameters which comply with the estimations made above.

3 Results and discussion

3.1 Shapes and dynamics of a healthy RBC

Fig. 1 shows several RBC shapes for χ = 0.58 (corresponding to a channel diameter D = 1.72D_r), which are typically encountered in microcirculatory blood flow; see also Movies S1–S4.† For slow flows, see Fig. 1(a), the RBC shape is similar to the biconcave discocyte shape in equilibrium. For higher flow 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 flow.^34,35 At the highest flow rate, see Fig. 1(c), a parachute shape is obtained, where the cell flows at the channel center and the membrane is practically not moving in the lab frame.

Fig. 1 Simulation snapshots of a RBC in flow (from left to right) for χ = 0.58. (a) A biconcave RBC shape at [small gamma, Greek, dot above]* = 5; (b) an off-center slipper cell shape at [small gamma, Greek, dot above]* = 24.8; and (c) a parachute shape at [small gamma, Greek, dot above]* = 59.6. See also Movies S1–S4.†

Fig. 2 presents our main result, the shape diagram for different flow rates and confinements, where the cell parameters are similar to those of a healthy RBC. The parachute RBC shape is predominantly observed in the region of strong confinements and high enough flow rates, where large flow forces are able to strongly deform a RBC. Here, it should be noticed that parachutes are also stable for weak or no confinement when the curvature of the parabolic flow in the center exceeds a critical value.^9,36 At weak confinements, we find off-center slippers with tank-treading motion for higher flow rates, and discocytes with tumbling motion for lower flow 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 [small gamma, Greek, dot above]^*_ttt of the tumbling-to-tank-treading transition of a RBC;^37,38 tumbling occurs off the tube center, when the local shear rate drops below [small gamma, Greek, dot above]^*_ttt. In the case of small viscosity contrast between inner and outer fluids (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 drift so that the tumbling axis is not fixed and oscillates in the vorticity direction (see Movie S2†). This effect is qualitatively similar to a rolling motion (also called kayaking) found in experiments³⁹ and in simulations⁴⁰ of a RBC in shear flow. 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.³⁹

Fig. 2 A phase diagram for Γ = 2662 (Y_r = 18.9 × 10⁻⁶ N m⁻¹, κ_r = 3 × 10⁻¹⁹ J), which mimics average membrane properties of a healthy RBC. The plot shows various RBC dynamics states depending on the flow strength characterized by [small gamma, Greek, dot above]* and the confinement χ. 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.

At small shear rates [small gamma, Greek, dot above]*, there also exists a so-called snaking region, first 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 drift (see Movie S1†), which is similar to that for a RBC rolling motion in shear flow 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 [small gamma, Greek, dot above]* ≲ k_BT/κ_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 confinements the “confined slipper” found in the 2D vesicle simulations is suppressed in 3D. The confined slipper in 2D found for χ ≳ 0.6 is qualitatively similar to a slipper at low confinements, which is also called “unconfined slipper” in ref. 9 and 10, since this vesicle state exists in unbound parabolic flow. Note that the regions of confined and unconfined slippers in 2D have no common boundary. The absence of slippers at high confinements in 3D is due to the cylindrical shape of a channel, which would cause the confined 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 defined 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 defined as [(λ₁ − λ₂)² + (λ₂ − λ₃)² + (λ₃ − λ₁)²]/(2R_g⁴), where λ₁ ≤ λ₂ ≤ λ₃ are the eigenvalues of the gyration tensor and R_g² = λ₁ + λ₂ + λ₃. 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.³⁷ 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 flow and the analysis of the characteristics discussed above.

Fig. 3 Characteristics of different RBC states: blue triangle – snaking discocyte ([small gamma, Greek, dot above]* = 9.9, χ = 0.72), red diamond – tumbling discocyte ([small gamma, Greek, dot above]* = 14.9, χ = 0.44), brown square – slipper ([small gamma, Greek, dot above]* = 49.7, χ = 0.44), and green circle – parachute ([small gamma, Greek, dot above]* = 64.6, χ = 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 [(λ₁ − λ₂)² + (λ₂ − λ₃)² + (λ₃ − λ₁)²]/(2R_g⁴), where λ₁ ≤ λ₂ ≤ λ₃ are the eigenvalues of the gyration tensor and R_g² = λ₁ + λ₂ + λ₃. The asphericity for a single RBC in equilibrium is equal to 0.15. See also Movies S1–S4.†

3.2 Comparison with experiments

There exist several experimental studies of a RBC in microchannel flow.^17,20,23,41 In the experiments of ref. 23, the rotation of single RBCs and of their rouleaux structures in tube flow with radii ranging from 30 μm to 100 μm 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 confinements and low flow rates with ≲ 50 s⁻¹, in agreement with the simulation results in Fig. 2.

In the experiments of ref. 17, the imposed flow velocities were very large, ranging from 1 cm s⁻¹ to 30 cm s⁻¹ in a capillary with the diameter of 9 μm. This is much faster than the typical flow velocities in microcirculation, where for venules and arterioles with a similar diameter flow velocities are in the range of 0.2–7 mm s⁻¹.^4,42 The range of flow velocities we span in our simulations is about 0.2–1.0 mm s⁻¹, and is therefore comparable with that in microcirculation. The results in ref. 17 show the existence of parachute and slipper shapes, where a weak confinement 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 confinements, where parachute shapes are observed for flow rates with the velocities lower than 4–7 cm s⁻¹.¹⁷ At flow velocities larger than approximately 7 cm s⁻¹ and at low confinements, centered slippers are observed which resemble parachutes, but become slightly asymmetric.¹⁷ Recent 2D simulations¹⁰ have also found that at high enough flow rates centered slippers and parachutes may coexist. Currently, we are not able to reach such high flow rates in 3D simulations due to numerical limitations. In addition, this region might be of limited interest, since the corresponding flow rates are far beyond the physiologically relevant values.

Experimental data in ref. 20 were obtained for narrow capillaries with diameters ranging from 4.7 μm to 10 μm; however, the flow velocities are considerably smaller than those in ref. 17, in the range of 1–40 mm s⁻¹. For strong confinements χ ≳ 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 confinement of χ = 0.65 shows a good agreement since centered parachute shapes are found in both experiments and simulations. The transition to the parachute shape for χ = 0.65 occurred at the RBC velocity of about 0.5 mm s⁻¹ in the experiments,²⁰ while our simulations predict the transition velocity of about 0.45 mm s⁻¹.

3.3 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 flow, we calculated phase diagrams for both reduced Young's modulus and increased bending rigidity. The state diagram for membrane bending rigidity increased by a factor of five, which leads to κ_r = 1.5 × 10⁻¹⁸ J and Γ = 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 [small gamma, Greek, dot above]^*_c = [capital Omega, Greek, tilde](Γ, χ, κ_r/k_BT), where [capital Omega, Greek, tilde] is a universal function for each transition. In addition, a more general form of [small gamma, Greek, dot above]^*_c should also include the strength of thermal fluctuations 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 [small gamma, Greek, dot above]* and lead to random cell orientation modifying potentially a snaking region. Also, recent experiments on vesicles⁴³ and a corresponding theory⁴⁴ 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 shift to lower values of [small gamma, Greek, dot above]*. This behavior is consistent with the roughly linear dependence of the shear rate at the parachute-to-discocyte transition on RBC bending rigidity and shear modulus reported in ref. 7, which is equivalent to [small gamma, Greek, dot above]* = c₁(χ) + c₂(χ)Γ. Similarly, the snaking region shrinks towards lower [small gamma, Greek, dot above]* values.

Fig. 4 A phase diagram for Γ = 532 (Y_r = 18.9 × 10⁻⁶ N m⁻¹, κ_r = 1.5 × 10⁻¹⁸ J) such that the bending rigidity of a cell membrane is five times larger than κ_r = 3 × 10⁻¹⁹ J of a healthy RBC. The diagram depicts RBC behavior in tube flow with respect to the dimensionless shear rate [small gamma, Greek, dot above]* and confinement χ. 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.

Fig. 5 shows the state diagram for a RBC in tube flow with the membrane Young's modulus reduced by the factor of three (Y_r = 6.3 μN m⁻¹) 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 flow rates than those for a healthy RBC (Fig. 2). As a consequence, the snaking region shrinks substantially and is observed primarily at very low [small gamma, Greek, dot above]*. Another feature is that the tumbling region is significantly 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 flow depends roughly linearly on the RBC elastic properties and bending rigidity,⁷ as well as that the tumbling-to-tank-treading transition of a RBC in shear flow depends nearly linearly on the RBC elastic properties such as Young's modulus.^37,38 This implies that the transition lines for a fixed ambient temperature are linear functions of the Föppl–von Kármán number Γ, as discussed above. A comparison with the results of ref. 7 for RBCs with considerably smaller bending and shear rigidities indicates that strong thermal fluctuations destroy the regular snaking oscillations.

Fig. 5 A phase diagram for Γ = 887 (Y_r = 6.3 × 10⁻⁶ N m⁻¹, κ_r = 3 × 10⁻¹⁹ J) 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 [small gamma, Greek, dot above]* and the confinement χ. 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.

Correspondence of different systems with a fixed Föppl–von Kármán number Γ is also supported by the argument that the RBC relaxation time τ is a linear function of RBC membrane elastic parameters (κ_r or equivalently Y_r for a fixed Γ). Thus, a simulation with the parameters {D_r, Y_r, sκ_r, [small gamma, Greek, dot above]*} should lead to identical results as those obtained from a simulation with {D_r, Y_r/s, κ_r, [small gamma, Greek, dot above]*/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 fluctuations and Reynolds number of the flow. Finally, the assumption of linear dependence of τ on the membrane properties may become invalid for strong enough flow 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 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 Γ = 532 (Fig. 4) and Γ = 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 [small gamma, Greek, dot above]^*_c to describe the transition lines.

3.4 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 flow, and have therefore employed same viscosity for the suspending fluid and RBC cytosol. For a healthy RBC, the viscosity contrast λ, defined as the ratio of cytosol over blood plasma viscosity, is approximately λ = 5. Therefore, we want to discuss briefly the possible effect of the viscosity contrast on the dynamical states. Clearly, cytosol viscosity only plays a role when the internal fluid (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 modified. Recent experiments⁴⁵ and simulations⁴⁶ indicate that the viscosity contrast indeed strongly affects the tumbling-to-tank-treading transition characterized by [small gamma, Greek, dot above]^*_ttt.

An increase in λ is known to shift [small gamma, Greek, dot above]^*_ttt of the tumbling-to-tank-treading transition to higher shear rates.⁴⁶ Therefore, we expect that the boundary between the tumbling and slipper regions in Fig. 2 would shift to higher values of [small gamma, Greek, dot above]* 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 blood-related diseases or disorders,²² 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 λ = 1 in ref. 9 and λ = 5 in ref. 10 have shown an expansion of the snaking region toward the slipper region, since the RBC tank-treading 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 significant 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 definite statement on this issue is only possible after a systematic numerical investigation of the effect of viscosity contrast has been performed in 3D.

3.5 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 ω on shear rate, as shown in Fig. 6 for three different confinements 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 [small gamma, Greek, dot above]* in the tumbling regime. Increasing confinement significantly 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. 6 Normalized snaking, tumbling, and swinging (slipper state) frequencies ω of RBCs as a function of the dimensionless shear rate [small gamma, Greek, dot above]*. The elastic membrane parameters are the same as in Fig. 4 (Y_r = 18.9 × 10⁻⁶ N m⁻¹, κ_r = 1.5 × 10⁻¹⁸ J), i.e. the bending rigidity is five times larger than that used for a healthy RBC, and Γ = 532. Three confinements are shown, as indicated. The symbols depict simulation results for various confinements. The regions of [small gamma, Greek, dot above]* corresponding to different RBC states can be distinguished in the plot by colors and line types – swinging (brown, dashed line), tumbling (red, dotted line), and snaking (blue, solid line).

3.6 Flow resistance

Fig. 7 presents the volumetric flow rate Q_RBC with a RBC measured in simulations and normalized by the flow rate Q without a RBC. Data are shown for different [small gamma, Greek, dot above]* for the case of Γ = 2662, corresponding to the shape diagram in Fig. 2. The volumetric flow rate at low confinements remains essentially unaffected by the presence of a RBC. As the confinement is increased, the flow rate decreases and the effect of a RBC on the flow rate appears to be stronger at low [small gamma, Greek, dot above]* values. Note that Fig. 7 also presents the change in flow resistance, since the apparent viscosity is inversely proportional to the volumetric flow rate. Thus, the apparent viscosity increases with increasing confinement. A reduction in the flow resistance with the increase of [small gamma, Greek, dot above]* (represented by different curves in Fig. 7) is due to the transition from discocyte to parachute and slipper shapes, since for large enough [small gamma, Greek, dot above]* no significant further changes in flow resistance occur.

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 χ for different [small gamma, Greek, dot above]*. The Föppl–von Kármán number Γ = 2662 is representative of healthy RBCs.

Even though Fig. 7 shows the effect of confinement on Q_RBC, the flow resistance for different χ 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². For the strongest confinement χ = 0.79 (D = 8.2 μm) H_t is equal to 0.027, while for the lowest confinement χ = 0.37 (D = 17.8 μm) H_t = 0.0057.

4 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 confinements (χ ≳ 0.7) and low flow 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 χ – in agreement with the 3D results of ref. 7. At even higher confinements than those in this study, RBCs are expected to attain bullet shapes.²⁰ 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 fluids inside and outside a RBC.

The calculated state diagrams in 3D provide a better description of RBC shapes and dynamics in microvascular flow than the previous 2D results. The flow resistance is affected only weakly by the transition from discocyte to parachute and slipper shapes, most significantly at high confinements (χ ≳ 0.6). The geometrical complexity of microvasculature induces non-trivial partitioning of RBCs, which often 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 flow and RBC behavior in microcirculation. Finally, similar physical mechanisms are expected for capsule suspensions.

Acknowledgements

A CPU time grant by the Jülich Supercomputing center is gratefully acknowledged. Dmitry A. Fedosov acknowledges funding by the Alexander von Humboldt Foundation. We also thank the DFG Research Unit FOR 1543 “SHENC – Shear Flow Regulation in Hemostasis” for support.

References

1. R. G. Larson, The structure and rheology of complex fluids, Oxford University Press, Oxford, NY, 1999.
2. J. A. Davis, D. W. Inglis, K. M. Morton, D. A. Lawrence, L. R. Huang, S. Y. Chou, J. C. Sturm and R. H. Austin, Proc. Natl. Acad. Sci. U. S. A., 2006, 103, 14779–14784.
3. H. Craighead, Nature, 2006, 442, 387–393.
4. A. S. Popel and P. C. Johnson, Annu. Rev. Fluid Mech., 2005, 37, 43–69.
5. T. Omori, Y. Imai, T. Yamaguchi and T. Ishikawa, Phys. Rev. Lett., 2012, 108, 138102.
6. A. Kumar and M. D. Graham, Phys. Rev. Lett., 2012, 109, 108102.
7. H. Noguchi and G. Gompper, Proc. Natl. Acad. Sci. U. S. A., 2005, 102, 14159–14164.
8. D. A. Fedosov, B. Caswell and G. E. Karniadakis, Biophys. J., 2010, 98, 2215–2225.
9. B. Kaoui, N. Tahiri, T. Biben, H. Ez-Zahraouy, A. Benyoussef, G. Biros and C. Misbah, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2011, 84, 041906.
10. N. Tahiri, T. Biben, H. Ez-Zahraouy, A. Benyoussef and C. Misbah, Microvasc. Res., 2013, 85, 40–45.
11. D. A. Fedosov, H. Noguchi and G. Gompper, Biomech. Model. Mechanobiol., 2014, 13, 239–258.
12. R. S. Sprague, M. L. Ellsworth, A. H. Stephenson, M. E. Kleinhenz and A. J. Lonigro, Am. J. Physiol., 1998, 275, H1726–H1732.
13. A. M. Forsyth, J. Wan, P. D. Owrutsky, M. Abkarian and H. A. Stone, Proc. Natl. Acad. Sci. U. S. A., 2011, 108, 10986–10991.
14. R. Skalak, Science, 1969, 164, 717–719.
15. P. Gaehtgens, C. Dührssen and K. H. Albrecht, Blood Cells, 1980, 6, 799–812.
16. U. Bagge, P.-I. Branemark, R. Karlsson and R. Skalak, Blood Cells, 1980, 6, 231–237.
17. M. Abkarian, M. Faivre, R. Horton, K. Smistrup, C. A. Best-Popescu and H. A. Stone, Biomed. Mater., 2008, 3, 034011.
18. I. V. Pivkin and G. E. Karniadakis, Phys. Rev. Lett., 2008, 101, 118105.
19. J. L. McWhirter, H. Noguchi and G. Gompper, Proc. Natl. Acad. Sci. U. S. A., 2009, 106, 6039–6043.
20. G. Tomaiuolo, M. Simeone, V. Martinelli, B. Rotoli and S. Guido, Soft Matter, 2009, 5, 3736–3740.
21. B. Kaoui, G. Biros and C. Misbah, Phys. Rev. Lett., 2009, 103, 188101.
22. M. Diez-Silva, M. Dao, J. Han, C.-T. Lim and S. Suresh, MRS Bull., 2010, 35, 382–388.
23. H. L. Goldsmith and J. Marlow, Proc. R. Soc. London, Ser. B, 1972, 182, 351–384.
24. P. Espanol and M. Revenga, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2003, 67, 026705.
25. D. E. Discher, D. H. Boal and S. K. Boey, Biophys. J., 1998, 75, 1584–1597.
26. D. A. Fedosov, B. Caswell and G. E. Karniadakis, Comput. Methods Appl. Mech. Eng., 2010, 199, 1937–1948.
27. H. S. Seung and D. R. Nelson, Phys. Rev. A, 1988, 38, 1005–1018.
28. E. A. Evans and R. Skalak, Mechanics and thermodynamics of biomembranes, CRC Press, Inc., Boca Raton, Florida, 1980.
29. L. B. Lucy, Astron. J., 1977, 82, 1013–1024.
30. J. J. Monaghan, Annu. Rev. Astron. Astrophys., 1992, 30, 543–574.
31. P. J. Hoogerbrugge and J. M. V. A. Koelman, Europhys. Lett., 1992, 19, 155–160.
32. P. Espanol and P. Warren, Europhys. Lett., 1995, 30, 191–196.
33. D. A. Fedosov and G. E. Karniadakis, J. Comput. Phys., 2009, 228, 1157–1171.
34. R. Tran-Son-Tay, S. P. Sutera and P. R. Rao, Biophys. J., 1984, 46, 65–72.
35. T. M. Fischer, Biophys. J., 2007, 93, 2553–2561.
36. J. L. McWhirter, H. Noguchi and G. Gompper, Soft Matter, 2011, 7, 10967–10977.
37. M. Abkarian, M. Faivre and A. Viallat, Phys. Rev. Lett., 2007, 98, 188302.
38. J. M. Skotheim and T. W. Secomb, Phys. Rev. Lett., 2007, 98, 078301.
39. J. Dupire, M. Socol and A. Viallat, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 20808–20813.
40. D. Cordasco and P. Bagchi, Phys. Fluids, 2013, 25, 091902.
41. G. Tomaiuolo, V. Preziosi, M. Simeone, S. Guido, R. Ciancia, V. Martinelli, C. Rinaldi and B. Rotoli, Ann. Ist. Super. Sanita, 2007, 43, 186–192.
42. A. R. Pries, T. W. Secomb and P. Gaehtgens, Am. J. Physiol., 1995, 269, H1713–H1722.
43. M. Levant and V. Steinberg, Phys. Rev. Lett., 2012, 109, 268103.
44. D. Abreu and U. Seifert, Phys. Rev. Lett., 2013, 110, 238103.
45. V. Vitkova, M.-A. Mader, B. Polack, C. Misbah and T. Podgorski, Biophys. J., 2008, 95, L33–L35.
46. A. Z. K. Yazdani and P. Bagchi, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2011, 84, 026314.

---

Footnote

† Electronic supplementary information (ESI) available: Four movies of different RBC shapes in capillary flow: snaking discocytes, tumbling discocytes, swinging slippers, and parachutes. Blue particles on the RBC are firmly attached to the membrane and serve as tracers in order to visualize the membrane dynamics; these particles have no physical meaning. See DOI: 10.1039/c4sm00248b

---