Ordering and arrangement of deformed red blood cells in flow through microcapillaries

The shapes and alignment of elastic vesicles similar to red blood cells (RBCs) in cylindrical capillary flow are investigated by mesoscopic hydrodynamic simulations. We study the collective flow behavior of many RBCs, where the capillary diameter is comparable to the diameter of the RBCs. Two essential control parameters are the RBC volume fraction (the tube hematocrit, HT), and the suspension flow velocity. Depending on HT, flow velocity and capillary radius, the RBC suspension exhibits a disordered phase and two distinct ordered phases, consisting of a single file of parachute-shaped cells and a zigzag arrangement of slipper-shaped cells, respectively. We argue that thermal fluctuations, included in the simulation method, coupled to hydrodynamic flows are important contributors to the RBC morphology. We examine the changes to the phase structures when the capillary diameter and the material properties (bending rigidity κ and stretching modulus μ) of the model RBCs are varied, constructing phase diagrams for each case. We focus on capillary diameters, which range from about 1.0 to about 1.4 times the RBC long diameter. For the smallest capillary diameter, the single-file arrangement dominates; for the largest diameter, the ordered zigzag arrangement begins to loose its stability and alternates with an asymmetric structure with two lanes of differently oriented cells. In simulations with long capillaries, the coexistence of different phases can be observed.


Introduction
Human blood is composed of red blood cells (RBCs) at about 45% volume fraction, white blood cells at about 1% volume fraction and blood plasma. The blood plasma is rich in biological components such as lipids, enzymes and platelets. In the absence of significant shear flow gradients, RBCs form aggregates called rouleaux where the RBCs loosely adhere to each other like a stack of coins [1,2]. This formation depends on the presence of fibrinogen in the cell membrane and globulin in the plasma. Subjected to sufficiently large shear gradients, the RBC rouleaux in a flow break-up and disperse into their constituent RBCs [3].
The goal of our numerical investigation is to elucidate the structures and arrangements of RBCs in flow through narrow cylindrical capillaries that arise from the coupling between the RBC deformability, hydrodynamics under flow and thermal fluctuations that induce RBC membrane undulations. Since we want to focus on the interplay of RBC deformation and hydrodynamics, we neglect any direct attractive interactions; this corresponds to sufficiently large flow rates, or to washed RBCs in Ringer's solution, where globulin is absent [4].
In blood flows with no significant shear gradients, human RBCs have a biconcave-disc shape with a maximum diameter and thickness of about 7.6 and 2 µm, respectively [1]. The RBC membrane consists of a lipid bilayer supported by an attached spectrin network that acts as a cytoskeleton. The bilayer resistance to bending is controlled by the curvature elasticity with a bending rigidity κ; the spectrin's resistance to a shear strain is characterized by the in-plane shear elasticity with shear modulus µ. At thermal equilibrium, the discocyte shape of an RBC can be predicted theoretically by minimizing the membrane's bending and stretching energy subject to the constraints of fixed surface area and volume. Under flowing conditions, however, the shapes and arrangements adopted by many RBCs is determined by the competition between these mechanical properties and the external hydrodynamic flow forces arising from the blood 4 has been quantified. With increasing flow velocity, the transition from a discocyte to a parachute shape was found to occur abruptly at a critical flow velocity, not gradually as in earlier studies based on the boundary integral method [13,30]. This difference in results can be attributed to the symmetry breaking for the discocyte state in [33], unlike in the earlier studies [13,30] where a symmetric, coaxial RBC orientation was imposed. Recently, we have generalized the MPC approach to treat many model RBCs; we investigated the hydrodynamically mediated clustering of parachute-shaped RBCs at high flow velocities and low hematocrits [36,37]. The results corroborate optical microscopy experiments [22,50] which showed that RBCs, moving in single-file in capillary flow, can form clusters which are distinct from the rouleaux described above. Unlike a rouleaux, there are no direct contacts between any two RBCs in these clusters, and the RBCs are not discocytes.
In this paper, we consider in more detail the flow-induced morphology of dense RBC suspensions with high H T 0.2, in narrow capillaries with diameters in the range 8.1-10.4 µm, i.e. with diameters slightly larger than the RBC diameter. In our previous study of capillary blood flow [36], we have focused on six RBCs in a capillary with R cap /R 0 = 1.4, where R 0 = √ S/4π is the effective radius of RBCs defined by its surface area S. The aim of the present investigation is to obtain a more detailed understanding of the RBC dynamics, in particular its dependence on capillary radius and membrane elasticity.
Recently, there has been an impressive body of work examining the structure and flow behavior of many RBCs (or deformable particles) under flow in channels whose dimensions are considerably larger than the maximum RBC diameter [31, 32, 34, 35, 38-40, 51, 52]. These simulation studies involve many hundred RBCs, sometimes in structured channels with choking points and bifurcations [31,32]. In contrast, we are interested in suspension with several RBCs per simulation cylinder (with periodic boundary conditions in the flow direction). We map out the structure and phase behavior with variations in the flow velocity, hematocrit H T , capillary diameter and membrane material properties, κ and µ.
This paper is organized as follows. In section 2, we outline our simulation method and membrane model. Results are given in section 3. First, we explain the dynamics of a single RBC in a cylinder with periodic boundary condition in section 3.1. Next, we construct a phase diagram for one choice of capillary diameter in section 3.2, then we observe changes to this initial diagram and associated phase behavior given first changes in the capillary diameter in sections 3.2 and 3.3, then changes in the membrane material parameters in section 3.4. We conclude with a summary and discussion in section 4.

Red blood cell (RBC) model and mesoscale hydrodynamics technique
The membrane of each RBC is modeled as a collection of N mb = 500 vertices of mass m mb , interconnected by two triangular networks of bonds [53,54]: a fixed network whose bonds are harmonic springs, and a dynamically triangulated network whose bonds undergo 'flips' [55,56]. The elastic, fixed network models the spectrin cytoskeleton of a real RBC, while the fluid, dynamic network models the viscous lipid bilayer. The spring constant, k el of a harmonic bond generates a shear modulus, µ = √ 3k el . In the fluid network, the fluctuations and shape changes of the membrane are controlled by the Helfrich curvature energy [57] with a bending rigidity, κ, set in most simulations to κ = κ 0 = 20k B T ; here k B and T are the Boltzmann constant and temperature, respectively. The membrane vertices have an excluded volume interaction, where neighboring vertices, whether bonded or not, experience a repulsion such that the minimum distance between any two vertices is l min = 0.67a on the same RBC, with the MPC lattice constant a introduced below. The maximum bond length of the fluid network is l max = 1.33a. In addition, global volume and local area constraint potentials are added that ultimately keep the volume and global area of the membrane fixed at V = 450a 3 and S = 407a 2 . The reduced volume of a given RBC is 7a. For this volume, the stable shape of a fluid vesicle (without an elastic network) at thermal equilibrium is a biconcave discocyte. In terms of this effective RBC radius, the shear modulus µ of the elastic vesicle (RBC) is set in most simulations to µR 2 0 /k B T = 110, but two variations are considered in section 3.4. RBCs have a surface area of about 140 µm 2 , which gives an effective RBC radius R 0 = 3.3 µm. A list of these and other symbols and abbreviations is provided in table 1.
The fluid is modeled by N s point particles of mass m s where m s = m mb /10. The motion of these particles is governed by MPC dynamics [44][45][46][47], which was designed to describe hydrodynamic flows on mesoscopic length scales. This dynamics consists of a series of two alternating steps: a free streaming step with time step h = t CD where the fluid particles do not interact with one another, followed by a collision step where the fluid particles exchange momentum. In the collision step, the fluid particles are first sorted into the cells of a cubic lattice with lattice constant a; then, for particles in a given collision box, the particle velocities relative to the center-of-mass velocity of the collision box, are rotated by an angle π/2 about a randomly chosen axis. During a collision, the momentum and kinetic energy of all particles in a collision box are conserved; consequently, MPC dynamics describes the hydrodynamic flows in a Newtonian fluid. Thermostats are added to extract heat generated by shear gradients in the flow.

6
In MPC simulations, length and time are typically scaled according tox = x/a and t = t k B T /(m s a 2 ) so that the mean free path, l = h √ k B T /m s , is l/a =ĥ in reduced units. Given a MPC time step ofĥ = 0.025 and a fluid mass density of ρ m = 10m s /a 3 (on average 10 fluid particles per collision box), the shear viscosity of the fluid is η 0 = 20.1 √ m s k B T /a 2 ; we set this viscosity to be the same inside and outside of the RBC. The MPC time step is chosen to produce a Schmidt number, S c (the ratio between momentum and mass transport), consistent with liquids like water where S c ∼ 10 2 -10 3 . The membrane dynamics is integrated in time using a molecular dynamics (MD) algorithm; we employ a shorter time step than t CD , specifically Simulations are performed with a single (n RBC = 1) and many (n RBC = 2-6) RBCs in a cylindrical capillary with various simulation box (i.e. cylinder) lengths, L z ; periodic boundary conditions are employed along the flow (z) direction. In reporting our results, we scale all lengths by the RBC effective radius R 0 , energy by 8πκ 0 and time by τ = η 0 R 3 0 /κ 0 , the bending energy and the relaxation time of a nearly spherical fluid vesicle with a bending rigidity κ 0 = 20k B T , respectively. To mimic the pressure-driven flow typically employed in experiments on capillary flow with a pressure difference between an inlet and an outlet reservoir, the flow in our simulations is driven by a uniform 'gravitational' force in the z direction [49], ∇ z P = −ρ m g, which in our reduced units is g * = ρ m g R 4 0 /κ 0 . The advantage of such a driving force is the minimization of density gradients in the MPC fluid. The mean flow velocity is v m = ρ m g R 2 cap /8 0 = (5k B T R 2 cap /2η 0 R 4 0 )g * in the absence of RBCs. The MPC fluid particles interact with the capillary walls via a 'bounce-back' rule-where fluid particle velocities are inverted upon collision with the walls-that produces no-slip boundary conditions [48]. Similarly, the fluid particles interact with the membranes via a bounce-back rule that scatters fluid particles off the membrane triangles (fluid network), thereby also preventing the fluid from the exterior of the membrane penetrating into the interior and vice versa. In addition, the membrane vertices interact with the fluid in the MPC collision step, where the vertices in a given collision box are treated just like the fluid particles (except for their different mass). The membrane vertices of different RBCs interact with one another via an excluded volume potential, where the minimum distance between vertices on two different vesicles is l min = 0.77a. This minimum distance is larger than that between two vertices on the same membrane (l min = 0.67a) to prevent interpenetration of two neighboring RBCs.
More details of the model and simulation technique can be found in [33,41,42,45].

Results
The structure and spatial correlations of an RBC suspension depend on the flow velocity, the hematocrit H T , the capillary radius R cap and the membrane elasticity. We will examine the RBC structure and arrangements in a capillary with a diameter of R * cap = R cap /R 0 = 1.23 or R * cap = 1.4 in sections 3.1 and 3.2, where the number of RBCs per simulation cylinder, n RBC , is set to n RBC = 1 and 6, respectively. The hematocrit is varied by changing the tube length (per RBC) L * RBC = L * z /n RBC (where L * z = L z /R 0 ). The spatial correlations for these two capillary radii are found to be not substantially different. The RBC arrangement becomes different in wider capillaries, as will be discussed in section 3.3 for R * cap = 1.58. For real RBCs, the three radii R * cap = 1.23, 1.4 and 1.58 correspond to capillary diameters of 8.1, 9.2 and 10.4 µm, respectively. In sections 3.1-3.3, the bending rigidity and shear modulus are set to κ = κ 0 = 20k B T and µR 2 0 /k B T = 110, respectively. Results for the changes to this structure due to a variation of the membrane material parameters κ and µ will be presented in section 3.4.
3.1. Periodic RBC arrays in narrow capillaries with R * cap = 1.4 First, we briefly review simulation results for the dynamics of a regular array of RBCs in a microcapillary of radius R * cap = 1.4, which have been obtained in [36,37] from simulations of a single RBC (n RBC = 1) in a cylinder with periodic boundary conditions. At low flow velocities, corresponding to g * 8, RBCs assume a discocyte shape. The interactions with neighboring cells (periodic images) imposed by the periodic boundary conditions force the RBC to adopt a tilted orientation with respect to the direction of flow [37]. At higher flow velocities, the RBC adopts a parachute shape when the hematocrit is sufficiently low; however, as H T increases (and L * RBC decreases), the RBC gradually transforms from a parachute shape into a shallow 'bowl' shape at a critical value of H T 0.18, corresponding to L * z 2.2, which is slightly less than the cylinder diameter [37]. The average maximum radial distance, r * max , of the RBC membrane from the capillary axis increases, and the average mean curvature at the front of the parachute/bowl decreases with increasing H T . The critical distance L * z 2.2 indicates approximately where the shape of an RBC ceases to be hydrodynamically independent from its neighbors with decreasing cell-to-cell distance. This result is consistent with earlier solutions of the Stokes equation for a line of equally spaced and axisymmetrically located rigid spheroids, cylinders and discocytes in a cylindrical tube [13][14][15][16][17].

RBC structures and correlations in narrow capillaries with R *
cap = 1.4 and 1.23 Hydrodynamic interactions between cells do not only affect their shapes, but also their spatial arrangements. In order to investigate the latter effect, several cells are needed in the simulation box. Therefore, we now consider six RBCs (n RBC = 6) in a tube with periodic boundary conditions. The length of the simulation cylinder is varied from L * RBC = 0.8 to 1.8. Figure 1 shows phase diagrams with several different flow phases; simulation snapshots of the corresponding RBC shapes and arrangements are displayed in figure 2. We first present the structure and flow behavior for R * cap = 1.4, and discuss the similarities and differences for R * cap = 1.23 at the end of this section. At low densities and low flow rates, a disordered-discocyte phase, D, appears. In this phase, no long-range correlations in the orientation of discocytes exist, although transient clusters are observed, in which the orientations of the discocytes are correlated; there are no direct contacts between the membranes of neighboring RBCs.
As expected, in faster flows with g * above a critical value g * = 13 (for R * cap = 1.4), the RBCs exist as a line of equally spaced parachutes, the aligned-parachute phase, Pc. Interestingly, at higher hematocrits of H T 0.4, corresponding to relatively low L * RBC = 0.8-1.2, the RBCs arrange as slippers into two parallel interdigitated rows [36], denoted the S phase: the centers of mass of two consecutive RBCs, each one in a different row, are displaced from one another by L * RBC . If the suspension of RBCs at L * RBC < 1.2 is initially prepared in a Pc conformation (by replicating the final fluid and membrane configurations of a simulation with n RBC = 1 at the same H T and g * ), then the suspension at some later moment in time spontaneously reorders into an S conformation. This reordering is triggered by an off-axis deviation in the position of one RBC as a result of thermal fluctuations. When this RBC moves radially off axis, it slows down and simultaneously deforms into a slipper shape; the subsequent RBC behind flow speeds up and then is pushed in the opposite radial direction. This sequence of events then repeats itself with the other RBCs farther back in the flow.
The different RBC shapes can be distinguished by a shape parameter, the asphericity, which quantifies the deviation from a spherical shape. It is defined as [42,58] where λ 1 λ 2 λ 3 are the eigenvalues of the gyration tensor; this implies α = 0 and α 0.18 for a sphere and a discocyte, respectively. The average asphericity, α , is shown in figure 3 for L * RBC = 1.75 and n RBC = 6. A sudden transition in shape would be reflected by a sudden drop in α with increasing g * ; the gradual and monotonic decrease of α for n RBC = 6 signals the absence of such a transition, qualitatively consistent with the n RBC = 1 results. However, the fluctuations in RBC shape, measured by the standard deviation, ( ) = 2 − 2 , display a pronounced peak at g * 11 (see figure 3).
The orientational order of the RBCs can be characterized by the orientational order parameter, Here θ z is the angle between the eigenvector of the RBC gyration tensor with the smallest eigenvalue and the z-axis. If this eigenvector is parallel or perpendicular to the flow direction, then S = 1 or −1/2, respectively. Figure 3 shows that at L * RBC = 1.75, S is a monotonically and gradually increasing function of g * , with mainly perpendicular orientation of RBC in the disordered phase. The disordered phase at low g * 11 is characterized by large orientational fluctuations, which are quantified by the standard deviation, σ (S); the fluctuations decrease with increasing g * once g * 12, i.e. in the Pc region of the phase diagram.
The pressure drop is the difference between the pressures required to maintain the fluid with RBCs and the original Poiseuille fluid with no RBCs at a particular mean flow velocity v m . Across a section of capillary of length L z , the pressure drop per RBC is given by where v 0 is the mean flow velocity in a capillary with no RBCs that is generated by the same uniform pressure gradient, g * . Figure 4 shows for H T = 0.37 (L * RBC = 1.05), that the pressure drop per RBC in simulations with n RBC = 1-which generate axisymmetric bowl-shaped RBCs (similar to those in the Pc phase)-is less than P drp in simulations with n RBC = 6, which generate S conformations; in addition, the internal energy of a slipper-shaped RBC is larger than of a bowl-shaped RBC. Thus, the system favors an evolution towards the S phase at sufficiently high H T , although the S phase creates a larger resistance to flow than the Pc phase, and RBCs in the S phase have a higher internal energy than the RBCs in Pc phase.
We consider next the sensitivity of our results to finite-size effects. At L * RBC = 1.05, the S-phase occurs if the number of RBCs in the simulation box is even, n RBC = 2, 4 or 6. For n RBC = 5, the S-phase conformation can still be recognized, but now contains a defect, see figure 5. Two interdigitated rows each with two slipper-shaped RBCs are found; however, the fifth RBC acts as a defect with its center of mass near the capillary axis. For n RBC = 3, the  suspension adopts two possible conformations, a 'sliding' conformation where two RBCs are pressed up against the capillary walls while the remaining RBC, centered about the capillary axis, slides between them, and a 'jammed' conformation where the middle RBC cannot pass its two neighboring RBCs. Figure 6 shows that the pressure drop per RBC converges to a value independent of n RBC for n RBC 6; therefore, we conclude that to examine the structural correlations of the suspension, the considered number of RBCs is appropriate for this channel diameter.
Spatial pair correlations along the z direction between two RBCs are characterized by the axial pair distribution function, G(z * nb ), of RBC centers of mass, where z * i j = |z * i − z * j | (z * i being the z coordinate of the center of mass of RBC i) and density ρ B = (n RBC − 1)/L * z . We employ a discretized -function with (z nb ) = 1/ z if z nb is within the interval [− z/2, z/2], and 0 otherwise (in the numerical analysis, we use z * = 0.03). In equation (4), the time average is denoted by . . . t . This function is proportional to the conditional probability of finding the center of mass of an RBC a distance z * nb away from the center of mass of a given, tagged RBC. Results are shown in figure 7. In both S and Pc phases, simulations with n RBC = 6 show pronounced peaks in G(z * nb ), evenly spaced at a distance L * RBC apart; G(z * nb ) nearly vanishes at z nb = 0 and between two peaks. Corresponding data for n RBC = 3, 5 show smaller peaks in G(z nb ), indicating a less ordered structure, in agreement with the snapshots of figure 5. By comparison, for the D phase, the peaks in G(z * nb ) are not as pronounced as in the S and Pc phases, and G(z * nb ) does not vanish at any z nb . The Pc and S phases can also be distinguished by comparing distributions of the radial position, R, of an RBC center of mass, where ρ 2D b = n RBC /(π R * 2 cap ), and δ(R * ) = 1/ R if R * is within the interval [− R/2, R/2] and 0 otherwise. Figure 7 shows that in the S phase, P(R * ) has a peak away from R * = 0 and essentially vanishes at R * = 0, while the Pc and D phases display a strong peak at R * = 0. The phase structure for the more narrow capillary, R * cap = 1.23, closely resembles the structure for R * cap = 1.4 (see figure 1). For R * cap = 1.23, the suspension still exists in a D, Pc or S phases depending on L * RBC and g * . However, if the suspension is prepared initially in an alignedparachute (Pc) conformation by replicating the steady-state fluid and membrane configurations from a n RBC = 1 simulation at L * RBC 1.4, then the Pc conformation can remain intact for a substantial duration of time before eventually reordering into an S conformation, particularly at L * RBC = 1.4. In other words, the suspension appears more stable against thermal fluctuations, which could trigger a structural change, at R * cap = 1.23 than at 1.4. The RBC energies in the S and the Pc phases at R * cap = 1.23 are almost identical to these energies (shear elastic and bending) at R * cap = 1.4, but the pressure drops, reflecting the resistance to flow, are much higher at R * cap = 1.23, see figure 4.

RBC structures and correlations in wider capillaries with R
The phase diagram as a function of flow velocity g * and capillary length, L * RBC , obtained from simulations with n RBC = 6, is shown in figure 8. The comparison of the phase diagrams for different R * cap , figures 1 and 8, indicates that the transition from the D to the Pc phase at low hematocrit shifts to larger g * and thus larger flow velocities with increasing R * cap at a given L * RBC (i.e. fixed line density of RBCs), as expected from simulations at low hematocrit in [37]. The arrangement of RBCs for R * cap = 1.58 at high hematocrit (low L * RBC ) is not as well defined as in the more narrow capillaries (R * cap = 1.4 and 1.23). In the high-hematocrit region of the phase diagram, where the S phase appears in the narrower capillaries, the suspension alternates in time between the interlocked state of the zigzag-slipper (S) phase and a disordered state with no immediately discernible structure, i.e. no significant long-range spatial and orientational correlations. Note that this region extends to lower H T -values for larger R * cap , compare figures 1 and 8. At the highest hematocrit examined (H T = 0.37), the S conformation  does not occur, but rather an asymmetric-lane (Al) conformation consisting of two parallel rows: one with four slipper-shaped RBCs, and the other with two discocyte-shaped RBCs (see figure 9); the suspension alternates between this ordered Al state and a disordered state. Figure 10 shows examples of the disordered state. The dynamics of the transformations between different phases in the high-Hematocrit region of the phase diagram of figure 8 is shown in movies 1 and 2 (available from stacks.iop.org/NJP/14/085026/mmedia). Here, movie S1 shows the transformation from the Pc to the S state and vice versa, movie S2 the transformation from the S state to the Al state and vice versa. To characterize the 'mixed state' with its alternating assembly and disassembly of ordered rows of RBCs at low L * RBC , we calculate spatial correlation functions and radial distributions. Figure 11 shows G(z * nb ) and P(R * ) at five points in the phase diagram, in the range L * RBC = 1.05-1.4, and g * = 13.2-18.5. In all cases, G(z * nb ) has peaks located at integer multiples of L * RBC . These peaks are not as pronounced as those observed for the S phase in figure 7 for R * cap = 1.4 (with n RBC = 6), and there is no z * nb where G(z * nb ) nearly vanishes. The reason is that the peaks are broadened by the intermittent presence of the disordered state. In the radial distribution function P(R * ), see figure 11, configurations in the disordered state mainly contribute to the peak located at R * = 0, while configurations in the interlocked S state contribute to the offcenter peak at R * 0.4-0.6.
At L * RBC = 1.05, the peaks in G(z * nb ) and the off-center peak in P(R * ) all decrease in magnitude with increasing g * (see figure 11), which implies that the suspension assumes an ordered, off-axis state less frequently as g * increases. This result is confirmed by an estimation of the percentage of the total simulation time in which all RBCs are found with their centers of mass at a radial distance greater than * R = 0.35 (a value corresponding to the minimum of P(R * ) in figure 11) from the capillary axis, which yields about 44, 24 and 18% at g * = 13.2, 15.9 and 18.5 respectively.
When R * > * R for all RBCs, we define the times for which the suspension exists in either the S or Al conformation as follows. We let N nb be the maximum number of neighbors j to a given RBC i that are a distance r x y = (x * i j ) 2 + (y * i j ) 2 > * x y away from i (where [x i , y i , z j ] are the coordinates of the center-of-mass position vector of an RBC). If N nb = 3 or 4, then the suspension exists in an S or an Al conformation, respectively. We set * x y = 0.88. At L * RBC = 1.05 and g * = 13.2, the suspension spends 93 and 6% of the time in the S and the Al states, respectively; however, it spends 83 and 15% at g * = 15.9, and 72 and 27% at g * = 18.5. Thus, the suspension is found much more often in the S state than in the Al state for L * RBC = 1.05, but the probability of the Al state increases with increasing g * . The suspension also becomes more disordered as the hematocrit decreases and L * RBC increases for fixed g * = 13.2; at L * RBC = 1.23 and 1.4, the suspension spends 38 and 26% of the time, respectively, in the S state and never occurs in the Al state.
In the mixed-state region, the displacement of RBCs from the center of the capillary increases with decreasing L * RBC at fixed g * ; figure 11 shows this trend at g * = 13.2. As the suspension is squeezed, there is an increased chance of finding three RBCs in a 'stack' whose axis is on average perpendicular to the flow direction; the existence of such stacks implies G(z * nb 0) > 0 (see figure 10(c)). Figure 12 shows orientational distributions and correlation functions at the same five points in the mixed-state region of the phase diagram, figure 8, discussed above. Let u i be the unit Here, u i and u j are orientation vectors of two neighboring RBCs. For the orientational distribution, the angle θ is defined by θ = arccos(|u i z |) where u i z is the z-component of u i . The angle distribution P(θ ), displays two trends, see figure 12. At fixed g * , as the average RBC distance L * RBC decreases, the RBC orientation vectors rotate towards the capillary axis, becoming more parallel with this axis; when the suspension exists in the S state, the deformation of the slipper-shaped RBCs increases as L * RBC decreases. At fixed L * RBC , the suspension becomes more disordered with increasing g * ; as a result, the distribution P(θ) at L * RBC = 1.05 broadens. For the suspension in the S state, the orientational correlation function, G uu (z * nb ), exhibits oscillation with a period of 2L * RBC starting at z * nb = L * RBC as shown in figure 12. The first minimum (negative value) and maximum (positive value) appear at z * nb 1.2 and 2.4, respectively. As with G(z * nb ) and P(θ), at fixed L * RBC the maxima in G uu (z * nb ) become less pronounced with increasing g * , reflecting a decrease of order in the interlocked zigzag conformations with increasing g * , compare figure 12.
By the analysis of these correlation and distribution functions, the various arrangements of RBCs in our microchannels can be clearly distinguished. For the cells with a smaller bending rigidity and γ = 8γ 0 (top panel of figure 13), the phase structure in the region, where the S phase occurs in figure 1, is less well defined or more loose compared to the S phase for (κ 0 , µ 0 ). For 1.0 L * RBC 1.2, the suspension alternates in time between a disordered state and an ordered S state, consisting of an interlocked zigzag of slipper-shaped RBCs. The transitions between these two states are seen during a simulation run. The deformation of the RBCs away from a discocyte shape becomes weaker with decreasing g * , so that at the lowest g * considered the zigzag-slipper conformations are replaced by zigzagdiscocyte conformations. At the lowest L * RBC and sufficiently high g * , the suspension also spends time in an Al state (see figure 9), alternating in time between the disordered, S, and Al states. Stacked conformations also appear, like for the wider capillary with R * cap = 1.58 discussed in section 3.3.
In an earlier study [33], the critical velocity, where an 'isolated' discocyte undergoes a shape transition from a discocyte to a parachute with increasing g * , was found to increase with increasing bending and shear elastic rigidity of the RBC membrane. The bottom panel of figure 13, where both moduli, κ and µ, are larger compared to those in figure 1, displays the same qualitative behavior; the transition from the disordered-discocyte phase to the aligned parachute phase is shifted to higher g * compared to this critical g * in figure 1. In addition, only the ordered zigzag state is found for g * 22 and L * RBC 1.15 below the transition line. However, at g * 22, the suspension is observed to alternate between an ordered, zigzag-slipper state and a disordered state. The RBCs in the interlocked zigzag conformations are discocytes at g * 22, and become slippers at higher g * .
As discussed above, the change of the membrane properties not only shifts the phase boundaries but also yields qualitative changes on the phase behavior. Some changes such as the existence of the Al phase are similar to the results of the variation of R cap . However, others such as discocytes with interlocked zigzag arrangement are not observed in the variation of R cap . The changes of the membrane properties and the capillary radius R cap modify the RBC deformation and hydrodynamic confinement effects, since the stability of each deformed shape depends on the membrane properties.
The pressure drop of suspensions of RBCs with different membrane parameters are displayed and compared in figure 14. For the suspension with a larger bending and shear rigidity, and γ = 2γ 0 , the resistance to flow, measured by the pressure drop, is higher than this drop for the suspensions with γ 0 and 8γ 0 . In addition, the pressure drop for the suspension with bending rigidity κ 0 and Föppl-von Kármán number γ 0 is almost the same as for the suspension with the reduced bending rigidity κ 0 /2 and 8γ 0 ; increasing the shear modulus by a factor of 4 is compensated by decreasing the bending rigidity by a factor of 2. Such a compensation is not unexpected qualitatively, because for Föppl-von Kármán numbers of about 100, bending and stretching energies are comparable, as shown for icosahedral shells in [59].

Summary and discussion
Using a mesoscopic simulation technique, which includes hydrodynamics and thermal fluctuations, we have explored the alignment and shape changes of a dense suspension of RBC-like vesicles-of discocyte shape with bending and stretching energies-in microcapillary flows. Three phases were observed, depending on the hematocrit H T , the cylinder diameter and the flow velocity. At relatively low H T , a disordered-discocyte (D) phase is found at low flow velocities; at velocities above a critical velocity an aligned-parachute (Pc) phase emerges. The transition between these phases is characterized not only by the change in RBC shape, but also by the spatial and orientational alignment of successive RBCs.
The zigzag-slipper (S) phase, predicted by our many-RBC simulations at high H T [36], was unexpected given the results of single-RBC simulations, which show a phase of coherently tilted discocytes at low flow velocities which deform gradually with increasing flow velocity into an aligned arrangement of axisymmetric RBCs shaped as shallow bowls or parachutes depending on the H T . However, earlier experiments on the flow of suspensions of human blood cells through narrow glass capillaries (with diameters in the range of 5-10 µm) did exhibit this phase (or 'zipper flow') for a particular range of tube hematocrits [60]. When we prepare the RBC suspension initially by replicating the steady-state results from the corresponding single RBC simulations, the suspension, initially in the Pc state, remains in this state for a length of time before quickly reordering into the S state. We believe that thermal fluctuations can trigger this structural change which is subsequently amplified by the hydrodynamics. Simple two-dimensional model studies of pressure-driven flow of a row of discs in a capillary have shown that the single-file arrangement is indeed unstable to initial periodic perturbations in the disc positions depending on the initial lateral position of the single-file row relative to the capillary axis [13,61,62]. For unstable motion, the initial, single-file row reordered into two interlocked rows of rotating discs (the 'zipper flow'). For RBCs, there is in addition the effect of a hydrodynamic lift force due to the tank-treading motion in the slipper shape, which pushes RBCs back to the center. The origin of a critical hematocrit separating the Pc and S phases is that if an RBC in the Pc phase fluctuates off-axis and slows down, then provided the L * RBC is sufficiently large this RBC can recover from this fluctuation, returning to being on-axis without disturbing other RBCs farther back in the flow.
Real RBCs have a larger bending rigidity and shear modulus than our model RBCs (approximate values are κ/k B T 50 and R 2 0 /k B T 5000, corresponding to a Föppl-von Kármán number 400). Therefore, a larger flow velocity is needed for a real RBC to induce the transformation from a discocyte to a parachute; for a single cell in [33], only a (linear) extrapolation of the simulation results to the measured values of bending rigidity and shear modulus produced good agreement with the experimental results. It may therefore not be too surprising that for typical values of R cap /R 0 , the critical value of L * RBC for the S phase determined experimentally for real RBCs in [60] is about three times larger than the one determined for our suspensions of model RBCs. In Stokes flow and in the absence of thermal fluctuations, the critical flow velocity is a linear function of the elastic moduli for a fixed Föppl-von Kármán number-but characteristic distances and shapes should be unaffected. Therefore, two factors could be responsible for this difference in the critical values of L * RBC . First, the smaller flow velocities in the simulations yield a smaller Péclét number, i.e. diffusion is more important; this implies that the arrangement of the cells is disordered, and higher hematocrits are necessary for cells to interact sufficiently strongly. Second, the smaller Föppl-von Kármán numbers in the simulations imply that the parachute and slipper shapes are more smoothly rounded; this probably makes them less favorable in flow, and again favors the parachute phase. Furthermore, since the shear-flow gradients are larger, for an initial Pc conformation to remain intact a much larger distance is needed between a parachute-shaped real RBC undergoing an off-axis fluctuation in position and the next real RBC farther back in the flow.
The Pc state is more stable at lower R * cap . At higher R * cap , the structure of the suspension in this region is more 'loose'; the suspension alternates in time between different states, making sudden transitions between ordered and disordered arrangements; these disordered arrangements occur more frequently as the flow velocity increased. This 'mixed state' can also be achieved by decreasing the bending rigidity while keeping the capillary diameter fixed.
The experiments of [60] also displayed a more disordered arrangement of RBCs at larger capillary diameters.
Our simulation results have been limited to a relatively small number of cells, usually containing six RBCs-although some exploratory simulations were performed with less and more RBCs. The existence of a mixed state in our simulations implies that in much longer simulation cylinders, containing orders of magnitude more RBCs at the same H T , there can be ordered and disordered arrangements of RBCs coexisting in different spatial regions at the same time. A simulation with n RBC = 16 RBCs in the simulation cylinder indeed demonstrates such a behavior, as can be seen from the configurations displayed in figure 15, as well as movie 3 (available from stacks.iop.org/NJP/14/085026/mmedia). Here, a region of Pc phase nucleates in S phase, the two phases coexist for some time, then the Pc phase disappears again. Furthermore, for state parameters where our simulations yielded the S phase (i.e. where there were no transitions to other states), it is possible that in much longer simulation cylinders there can be defects in the interlocked arrangement of slipper shaped RBCs, such as the one illustrated by five RBCs in the simulation cylinder in figure 5. The existence of the Al state suggests the possibility of non-interdigitated, non-axisymmetric ordered structures along some spatial regions in larger-scale suspension simulations. The present work should provide insight and inspiration for such future studies.