Effect of cytosol viscosity on the flow behavior of red blood cell suspensions in microvessels

The flow behavior of blood in microvessels is directly associated with tissue perfusion and oxygen delivery. Current efforts on modeling blood flow have primarily focused on the flow properties of blood with red blood cells (RBCs) having a viscosity ratio $C$ of unity between the cytosol and suspending medium, while under physiological conditions the cytosol viscosity is about five times larger than the plasma viscosity (i.e., $C\approx 5$). The importance of $C$ for the behavior of single RBCs in fluid flow has already been demonstrated, while the effect of $C$ on blood flow has only been sparsely studied. We employ mesoscopic hydrodynamic simulations to perform a systematic investigation of flow properties of RBC suspensions with different cytosol viscosities for various flow conditions in cylindrical microchannels. Our main aim is to link macroscopic flow properties such as flow resistance to single cell deformation and dynamics as a function of $C$. Starting from a dispersed cell configuration, we find that the flow convergence and the development of a RBC-free layer (RBC-FL) depend only weakly on $C$, and require a convergence length in the range of $25D-50D$, where $D$ is the channel diameter. The flow resistance for $C=5$ is nearly the same as that for $C=1$, which is facilitated by a slightly larger RBC-FL thickness for $C=5$. This effect is due to the suppression of membrane motion and dynamic shape deformations by a more viscous cytosol for $C=5$, resulting in a more compact cellular core of the flow in comparison to $C=1$. The weak effect of cytosol viscosity on the flow resistance and RBC-FL explains why cells can have a high concentration of hemoglobin for efficient oxygen delivery, without a pronounced increase in the flow resistance.

The flow behavior of blood in microvessels is directly associated with tissue perfusion and oxygen delivery.Current efforts on modeling blood flow have primarily focused on the flow properties of blood with red blood cells (RBCs) having a viscosity ratio C of unity between the cytosol and suspending medium, while under physiological conditions the cytosol viscosity is about five times larger than the plasma viscosity (i.e., C ≈ 5).The importance of C for the behavior of single RBCs in fluid flow has already been demonstrated, while the effect of C on blood flow has only been sparsely studied.We employ mesoscopic hydrodynamic simulations to perform a systematic investigation of flow properties of RBC suspensions with different cytosol viscosities for various flow conditions in cylindrical microchannels.Our main aim is to link macroscopic flow properties such as flow resistance to single cell deformation and dynamics as a function of C. Starting from a dispersed cell configuration, we find that the flow convergence and the development of a RBC-free layer (RBC-FL) depend only weakly on C, and require a convergence length in the range of 25D − 50D, where D is the channel diameter.The flow resistance for C = 5 is nearly the same as that for C = 1, which is facilitated by a slightly larger RBC-FL thickness for C = 5.This effect is due to the suppression of membrane motion and dynamic shape deformations by a more viscous cytosol for C = 5, resulting in a more compact cellular core of the flow in comparison to C = 1.The weak effect of cytosol viscosity on the flow resistance and RBC-FL explains why cells can have a high concentration of hemoglobin for efficient oxygen delivery, without a pronounced increase in the flow resistance.

I. INTRODUCTION
Blood is a multi-component suspension which consists of plasma (≈ 55%) and cells (red blood cells ≈ 45%, white blood cells and platelets < 1%).Flow properties of blood are mainly governed by red blood cells (RBCs) [1][2][3][4], which play an important role in many physiological processes.For instance, RBCs are responsible for oxygen delivery and mediate margination (or migration) of platelets [5][6][7] and leukocytes [8][9][10] toward vessel walls, thus affecting the hemostatic process and immune response.Blood-flow properties are also crucial in many applications, such as the enrichment or separation of rare circulating tumor cells from blood [11,12] and the effectiveness of drug carriers for delivery to the targeted sites [13][14][15].
Suspensions of blood cells reveal complex flow properties [4,16] and rheology [17,18].Two representative examples of the flow behavior of RBC suspension in microvessels or glass capillaries are the Fahraeus [19] and Fahraeus-Lindqvist [20,21] effects.The former effect concerns RBC volumetric flux (or the so-called discharge hematocrit H d ), which appears to be larger than the tube (bulk) hematocrit H t in vessels with a diameter D in the range of 7 − 200 µm.The latter effect describes a minimum of blood-flow resistance in a tube with a diameter ≈ 8 µm, such that the resistance to flow increases for both smaller and larger vessel diameters.The main mechanism governing these phenomena is the formation of a RBC depleted region next to the wall, called RBC-free layer (RBC-FL), as the suspension flows [22,23].The thickness of RBC-FL is directly associated with blood flow resistance [24][25][26][27] and plays a crucial role in the adhesion of leukocytes, platelets, and drug-delivery carriers to vessel walls [6,15].
The thickness of RBC-FL is governed by a competition between the hydrodynamic lift force acting on RBCs in the direction away from the wall and cell-cell interactions which disperse them and drive them toward the wall [4,28,29].Particle migration in the Stokes flow regime (i.e., no inertia) relies on breaking the time-reversal symmetry, which is achieved through RBC dynamics and/or deformations.For instance, a tumbling rigid spheroidal particle does not experience any migration on average in the presence of a wall [30], but a stable tank-treading motion of a spheroidally shaped membrane leads to the migration away from the wall [31,32].As shown by theoretical analysis [31,33], experiments [32,34] and simulations [35], the migration velocity v l is proportional to f (C) γR 3 /d α w , where f (C) is a function of the ratio C = η in /η ex of internal η in and external η ex fluid viscosities, γ is the local shear rate, R is the characteristic particle size, d w is the distance away from the wall, and α is an exponent whose value is often reported to be close to two.Time reversibility of RBC motion is generally broken by cell orientation and deformation in flow [29,35,36].
The lift force is counterbalanced by RBC dispersion due to cell-cell hydrodynamic interactions in flow [4,28,37].It is intuitive that fluid flow can significantly enhance such interactions or collisions between particles, which are often referred to as shear-induced dispersion forces and depend on local shear rate, particle size, deformation, and dynamics [28,32,38].The response of RBCs to fluid stresses is known to be sensitive to the viscosity ratio C between internal and external fluids [18,[39][40][41].For example, at low viscosity contrasts C 3, RBCs tumble at low shear rates and tank-tread at high shear rates [42][43][44], while for C 3, the tank-treading motion is suppressed and replaced by dynamic multi-lobed shapes [18,39].These differences in RBC deformation and dynamics as a function of C are expected to affect the lift force, cell-cell interactions in blood flow, and local structure of RBC suspensions, which influence blood-flow resistance.
The main focus of our investigation is the effect of C on the behavior of RBC suspensions in microvessels and the dependence of RBC-FL thickness and flow resistance on the viscosity contrast.Even though recognized, the importance of C > 1 for blood flow is not well studied so far.Katanov et al. [28] have investigated the formation of RBC-FL and flow convergence to steady state for C = 1, starting from an initially dispersed configuration of RBCs, and found that the full flow convergence requires a length of about 25D (D is the tube diameter), which is nearly independent of flow rate, RBC hematocrit, and channel size for 10 µm < D < 100 µm.A recent numerical investigation of blood flow in a tube with a diameter of D = 70 µm [45] has predicted a comparable RBC-FL thickness and flow resistance for suspensions with C = 1 and C = 5, where the differences become more pronounced at high flow rates.Another simulation study of blood flow in a slit [46] has suggested a domination of the lift force on RBCs over cell-cell interactions for C = 1, such that a slightly smaller RBC-FL was found for C = 5 suspension in comparison to C = 1.To clarify the importance of C for the behavior of RBC suspensions, we have performed a systematic investigation using mesoscopic hydrodynamic simulations, which include suspensions with C ∈ [1,20] as well as rigidified RBCs, several different flow rates, hematocrits, and tube diameters.In particular, we investigate the development of RBC-FL for the various conditions, and connect it to the flow resistance, deformation and dynamics of single RBCs.Our results show that the RBC-FL develops faster for C = 1 in comparison with C = 5 due to a larger lift force on RBCs with C = 1.The flow-convergence length becomes larger for elevated C values, but remains within approximately 50D.The RBC-FL thickness for C = 5 is slightly larger than that for C = 1, resulting in a nearly negligible effect of C on the flow resistance.This property is due to a smaller dispersion of cells for C = 5 in comparison to C = 1, since a larger internal viscosity dampens shape changes and membrane dynamics of RBCs.Suspensions with stiffened RBCs, which approximate the case of C → ∞, exhibit the smallest RBC-FL and the largest flow resistance.The robustness of flow resistance with respect to C ∈ [1,20] permits RBCs to contain a cytosol with a high concentration of hemoglobin, which maximizes oxygen delivery and does not strongly affect flow resistance.
The paper is organized as follows.Simulation methods, models, and setup are introduced in Section II.Section III presents simulation results, where the behavior of RBC suspensions with C ∈ [1,20] is investigated.The analysis of single cell characteristics is performed in Sections III B and III C, in order to explain differences in RBC-FL for various C. Finally, the dependence of RBC-FL thickness on several parameters, including flow rate, tube hematocrit and diameter, is investigated in Section III D. Our main results are discussed and summarized in Section IV.

II. METHODS AND MODELS
Fluid flow is modeled by the smoothed dissipative particle dynamics (SDPD) method [47] with angular momentum conservation [48], which is a mesoscopic particle-based hydrodynamics approach.The conservation of angular momentum is crucial for the proper representation of cellular motion when distinct fluid viscosities inside and outside the cell are employed [48].RBCs are represented by a spring-network model [49][50][51], and coupled to fluid flow through dissipative forces.Below, we briefly review several model ingredients with an emphasis on the implementation of the viscosity contrast between internal and external fluids separated by the membrane.More details about the methods and models can be found in Refs.[28,48,52].

A. RBC membrane model
The RBC membrane is represented by a spring-network model [49][50][51] with N v vertices distributed at a biconcave cell shape.Potential energy of the membrane, consists of several contributions.U sp corresponds to the spring's energy, which mimics elasticity of the spectrin network attached to the back side of the lipid membrane.U bend is the bending energy, representing bending resistance of the lipid bilayer.U area and U vol impose area and volume conservation constraints, which mimic area incompressibility of the lipid bilayer and incompressibility of a cytosol, respectively.The biconcave shape of a RBC at rest is imposed by setting the reduced volume V * = 6V r /(πD 3 r ) = 0.64, where V r is the RBC volume and D r = A r /π, with A r being the area of a RBC.The RBC membrane is characterized by the shear modulus µ and bending rigidity κ, which are implemented through in-plane elastic forces from the modeled springs and out-of-plane bending forces acting on each pair of adjacent triangles, respectively.The membrane parameters are set to mimic average properties of a healthy RBC with µ = 4.8 µN/m and κ = 70k B T .The effective size of a RBC is D r = 6.51 µm, the surface area is A r = 132.9µm 2 , and the total volume is V r = 92.45µm 3 .The stress-free shape of a RBC elastic network is assumed to be an oblate spheroid with a reduced volume of 0.96.

B. Fluid-membrane interactions
To model the viscosity ratio C = η in /η ex = 1, internal and external fluids have to be separated by the membrane.An impenetrable membrane is implemented through bounce-back boundary conditions (BCs) for both internal and external fluid particles at every triangular face of the membrane.Thus, internal fluid particles are subject to bounced-back BCs from inside the cell, while external fluid particles are bounced back from the outer membrane surface.Different fluid viscosities are implemented through different friction coefficients of dissipative forces in the SDPD method.Dissipative interactions between internal and external fluids assume the average of the two friction coefficients.
The frictional (dissipative) coupling between fluid and membrane particles is implemented through a dissipative force [51], similar to that in the dissipated particle dynamics method [53,54], where γ is the friction coefficient, α = 0.2 is an exponent of the weight function, r m is the cutoff radius, v ij = v i − v j is the velocity difference, and The value of γ is computed as [51], where η is the fluid viscosity, ρ is the fluid density, and V h represents a half sphere of radius r m in the positive z direction.This estimation of the friction coefficient assumes a linear flow-velocity profile within a distance r m near the membrane surface, so that the local shear rate cancels out.The cutoff for fluid-membrane coupling is set to r m = 0.75 µm.
To prevent an overlap of two membranes, a short-ranged repulsive Weeks-Chandler-Anderson (WCA) potential is applied between pairs of membrane particles belonging to different cells within the cutoff distance r W CA .Note that the thickness of RBC-FL is slightly sensitive to the choice of r W CA , which depends on membrane resolution.For N v = 1000, the average spring length at rest is approximately l ave = 0.4 µm and r W CA is set to 0.3 µm to prevent overlap.Doubling the membrane resolution allows a slight decrease of r W CA , which may result in a slight change of the RBC-FL thickness.Nevertheless, N v = 1000 is large enough to properly represent membrane deformation and nearly eliminate the effect of the cutoff distance r W CA on blood flow properties.

C. Simulation setup
The computational domain is a periodic cylindrical tube with a diameter D and a length L = 60 µm (≈ 9.2D r ), which is long enough to avoid finite-size effects.Initially, RBCs are introduced into the computational domain with an ordered structure.Their number is determined by hematocrit H t (the volume fraction of RBCs in a tube), which is assumed to be H t = 30% in most simulations.Then, fluid particles are randomly placed with a uniform distribution into the computational domain.The number density of fluid particles is n = 9 µm −3 , the smoothing length for SDPD pair interactions is r c = 1.04 µm, resulting in about 30 particles within the interaction range r c .The SDPD fluid particles inside the cells are set to represent the internal fluid with viscosity η in , while particles outside the cells correspond to the external fluid with viscosity η ex .To relax the initially ordered structure of RBCs, the cellular suspension is mixed by applying a flow within the tube, which is driven by a force exerted on all fluid particles.When the cells are well mixed, the flow is stopped, and the RBCs are let to diffuse and fill up the whole tube, resulting in a dispersed RBC configuration illustrated in Fig. 1(a).After the dispersed cell configuration is reached, a constant force f applied on all fluid particles in the x direction is turned on again to drive the fluid flow.Note that this driving force represents a uniform pressure gradient ∆P L = n • f , where ∆P is the pressure drop over the length L. The generated flow leads to RBC migration away from the wall and the formation of a denser cellular region near the tube center, as shown in Figs.1(b) and 1(c) for two different C values.
The flow strength is characterized by the dimensionless capillary number which represents the ratio between fluid stresses and elastic membrane stresses.Here, γ = v/D is the average shear rate with the average velocity v for a Newtonian fluid with viscosity η ex driven by the same pressure gradient.Thus, the capillary number directly characterizes the applied driving force or the pressure gradient.Note that τ = η ex D r /µ represents a characteristic relaxation time of a RBC, which is equal to approximately 1.6 × 10 −3 s for µ = 4.8 µN/m and blood-plasma viscosity η ex = 1.2 mP a • s.Simulations are performed at low enough Reynolds numbers, such that the largest investigated flow rate corresponds to Re = ρvD r /η ex = 2.34 (ρ is the fluid density).Furthermore, the Mach number M a = v/c s with the speed of sound c s is less than 0.1 in all simulations, so that the fluid flow can be considered incompressible.Solid wall BCs are modeled by a layer of immobile SDPD particles with a width r c .The structure and density of the wall layer are identical to those of equilibrated SDPD fluid in a periodic box.The pair forces between fluid and wall particles are the same as those for fluid-fluid interactions.Furthermore, an adaptive shear force is added to fluid particles within a near-wall layer of thickness r c to fully ensure no slip conditions [55].To prevent wall penetration, both fluid and membrane particles are reflected back inside the tube domain.

A. Development of the flow and RBC-FL
To examine flow development for various viscosity contrasts C ∈ [1,20], the viscosity of internal fluid η in is varied.This range of C covers RBC physiological conditions as well as some diseased states.After the pressure gradient is applied, the flow develops and reaches a terminal velocity profile with a position-averaged velocity v T at long times.Note that the velocity of individual cells depends on their location within the tube and can fluctuate in time even after the average flow velocity has reached steady state.v T is inversely proportional to the flow resistance, and can be used to obtain the relative suspension viscosity η rel = v/v T , which compares the volumetric flow rate of RBC suspension with that of external fluid (without RBCs) for the same pressure gradient.Thus, η rel quantifies an increase in the flow resistance due to the presence of RBCs.The dependence of η rel on C for H t = 30% is shown in Table I.Interestingly, the relative viscosity increases only by about 8% when C is increased 20 times.Furthermore, Table I also shows η rel = 1.7 for a suspension of stiffened cells (SC), whose shear modulus µ is 100 times larger than that of healthy RBCs.Stiffened RBCs do not exhibit significant deformation in fluid flow and represent a limit of very large viscosity contrast.
As the flow develops, RBCs migrate toward the tube center, resulting in the formation of RBC-FL near the wall whose thickness is directly associated with the flow resistance.To measure the thickness of RBC-FL, RBC suspension at a fixed time is projected onto the y-z plane [26], which is similar to taking a snapshot from experimental movie [56].Then, distances between the tube wall and projected RBC-core edge are extracted at several positions along the tube.Averaging these distances yields an average thickness of the RBC-FL δ.To improve RBC-FL statistics, we also employ flow axisymmetry, such that the thickness is sampled for 10 different angles by rotating a simulated configuration before the projection onto the y-z plane is performed.Furthermore, the RBC-FL data are accumulated over a certain time window, which is chosen long enough to avoid large deviations in RBC-FL thickness measurements and short enough to resolve the dynamics of RBC-FL development.This time window corresponds to about 8 ms, within which RBCs move on average 3 µm.The RBC-FL thickness for the initial configuration without flow (Fig. 1(a)) is δ ≈ 2.2 µm.This non-zero RBC-FL thickness is due to finite H t = 30%, the biconcave RBC geometry, which affects cell close-packing, and entropic repulsion from the wall that originates from the rotational diffusion of RBCs. Figure 2 shows the development of RBC-FL as a function of the average flow-convergence length L e = vt.At short times (L e < 5D = 100 µm), δ increases faster for the suspension with C = 1 than for those with C > 1, see Fig. 2(a).However, the RBC-FL thickness for C = 1 suspension saturates at a smaller value than δ for C > 1 with a difference of about 200 − 300 nm, as shown in Fig. 2(b).As C increases from 5 to 20, the converged RBC-FL thickness shows a similar plateau value of δ ≈ 3.1 µm.Converged or final RBC-FL thicknesses δ f for different C values are given in Table I.For the SC suspension with stiffened RBCs, δ fluctuates within the range of 2.0 − 2.3 µm, which is close to the RBC-FL thickness without flow.The differences in RBC-FL thicknesses for various C values are due to cell deformation and dynamics in flow, which will be discussed later.The critical convergence length L c e required for the development of RBC-FL is close to 25D (i.e., L c e ≈ 500 µm here) for both C = 1 and C = 5 suspensions, which is consistent with the previous investigation for C = 1 [28].Nevertheless, L c e becomes longer with increasing C, as L c e is approximately 50D for C = 20.The decrease of δ from its initial thickness for the SC suspension makes it difficult to define L c e due to a weak migration strength in this case.The development and dependence of RBC-FL thickness on C is determined by the migration and interactions between RBCs as cellular core forms.

B. Structure and mobility of the flowing RBC suspension
To better understand the dependence of final RBC-FL thickness δ f on the viscosity contrast, it is instructive to take a look at the structural properties of RBC suspensions after the flow has fully developed.The distribution of RBCs within the tube can be characterized by local hematocrit H t (r) obtained from simulations through spatial averaging of the density of fluid particles inside RBC membranes.Figure 3(a) shows the normalized local hematocrit H t (r)/H t within the tube for various viscosity contrasts.All H t (r) distributions contain a depletion zone near the wall, whose thickness is directly associated with δ f .For instance, the difference in δ f can clearly be seen for C = 1, C = 5 and SC suspensions.Inside the RBC-rich region, H t (r) is nearly uniform with a small peak near the center, which is typical for small tube diameters [27].For C ≥ 5, the H t (r) distributions are nearly independent of the viscosity contrast.
Radial cell density can also be characterized by the distribution of RBC centers of mass denoted as COM (r) and shown in Fig. 3(b).COM distributions have a peak near the RBC-FL.As the tube center is approached, COM (r) first decreases and then slightly increases for both C = 1 and C = 5 suspensions.As C is increased, spatial inhomogeneity in COM distribution becomes stronger.This is related to the RBC dynamics in flow and a decrease in radial migration of RBCs with increasing C, which will be discussed later.The SC suspension shows the strongest variations in COM (r), which is consistent with the H t (r) distribution in Fig. 3 Figure 4(a) presents flow velocity profiles v(r) normalized by the average velocity v for different RBC suspensions and blood plasma (i.e., for a Newtonian fluid).All velocity profiles for RBC suspensions are flattened at the tube center due to the presence of the cellular core.The velocity profiles do not exhibit large differences for various C values and overlap with the Newtonian-fluid case only in the RBC-FL.An increase of C results in slight flattening of v(r) near the tube center.Figure 4(b) shows radial profiles of local shear rates γ(r) normalized by the wall shear rate γw = 8 γ, which are computed from the velocity profiles v(r) in Fig. 4(a).Inside the cellular core, γ(r) for all RBC suspensions is close to zero (i.e., plug flow) and much less than that for the Newtonian case.Within the RBC-FL, γ(r) quickly increases from nearly zero to the wall shear rate.The SC suspension exhibits a steeper increase in γ(r) within the RBC-FL in comparison to soft-RBC suspensions.Even though fully developed velocity profiles are stable, RBCs within the cellular core are mobile, as can be seen from the trajectories of selected individual cells in Fig. 5(a).Thus, cells can migrate between different fluid layers laterally, even after the RBC-FL has fully developed (L e > 50D).To characterize the lateral mobility of RBCs, we define a dimensionless lateral mobility coefficient R T = |∆r com |/∆x com , where the absolute value of the cell velocity in radial direction, ∆r com /∆t, is normalized by its translational velocity along the x direction, ∆x com /∆t (here ∆t ≈ 8 ms).R T (r) is shown in Fig. 5(b) and can be interpreted as a measure of cell-cell collisions.R T (r) remains nearly constant within the cellular core and decreases slightly near the RBC-FL.Interestingly, R T (r) or fluctuations in lateral cell motion become smaller as C increases.Nevertheless, a large enough increase in C should eventually lead to an increase in R T (r), as for the SC suspension, where R T (r) values are slightly larger than those for C = 1.
For comparison, we have also performed simulations of a single RBC migrating from the wall to the tube center for C = 1 and C = 5 (the dashed lines in Fig. 5(b)).For a single cell, lateral migration due to the hydrodynamic lift force near the wall is much faster than that in a suspension, because cell-cell collisions are not present.However, in the tube center, cell-cell interactions within RBC suspensions enhance lateral migration of cells in comparison to the case of a single RBC.Note that R T (r) of a single RBC for C = 1 in Fig. 5(b) is larger than for C = 5, indicating a stronger lateral migration.Figure 6 presents the lift velocity v l of a single migrating RBC normalized by local shear rate γ(r) for C ∈ [1,10].Clearly, the RBC with C = 1 migrates faster than that with C = 5.This result is consistent with previous simulation studies [35,57], where the lift velocity on single RBCs in pure shear flow has been found to decrease with increasing C. The ratio v l / γ(r) is expected to be proportional to 1/d α w , where d w is the distance from the wall to the cell's COM.We find that α 2 for C =≤ 2, and 1 < α < 2 for C = 5, in agreement with experimental measurements [32,34,58].For the RBC with C = 10, α is nearly zero.Note that these simulations are performed in a tube with diameter D = 20 µm (D r /D = 0.33), representing a rather strong confinement with varying local shear rates.Furthermore, initial cell migration might be affected by the flow development, as we start from no-flow conditions.
The faster migration velocity of RBCs for C = 1 in comparison with C = 5 is due to differences in cell dynamics for different C. At high enough shear rates (e.g., near the wall), RBCs generally exhibit a tank-treading-like motion with a preferred alignment in flow for C = 1, while a tumbling-like dynamics with multi-lobed shapes is found for C = 5 [18,39].The tank-treading dynamics of RBC membrane near a wall leads to a larger lift force than for the tumbling dynamics [30,31].Note that the slower migration velocity of RBCs for C = 5 in comparison to C = 1 is consistent with a slower development of the RBC-FL for C = 5 in Fig. 2(a).

C. Shape and dynamics of single cells inside the RBC core
RBC dynamics inside the cellular core involves frequent changes in the shape, orientation, and membrane tanktreading motion, which often cannot be decoupled completely.Figure 7 shows several representative snapshots of RBCs near the RBC-FL and in the tube center.To better understand the behavior of single RBCs within the cellular core, we look at several cell characteristics.Shape changes of RBCs are quantified by the asphericity λ 1 ≤λ 2 ≤λ 3 are the eigenvalues of the gyration tensor.For a biconcave RBC shape at rest, A = 0.15, while A = 0 corresponds to a sphere and A = 1 to a long thin rod.For example, the largest eigenvalue λ 3 can be interpreted as the extension of a RBC in the flow direction (λ 3 = 4.77µm for a RBC at rest). Figure 8(a) shows that RBCs are significantly stretched in most locations within the cellular core region, except near the tube center.Close to the center, different shapes with A< 0.15 are generally observed (Fig. 8(b)), whose representative snapshots are shown in Fig. 7.Both λ 3 and A in Fig. 8 increase nearly linearly from the tube center to the RBC-FL.At the edge of the cellular core, elongated slipper shapes prevail, as shown in Fig. 7. SC cells show a crumpled configuration, which is due to not only fluid-flow stresses, but also residual elastic stresses within the membrane, since the RBC stress-free shape corresponds to an oblate spheroid with a reduced volume of 0.96 and the SC cells have a high shear modulus.Therefore, stiffened cells even without flow show some degree of membrane roughness in comparison with the smooth surface of soft RBCs.Nevertheless, the overall rest shape of SC cells remains biconcave.The shape characteristics of SC cells in Fig. 8 are nearly independent of the radial position within the tube, indicating that cell deformation can nearly be neglected.As C increases from 1 to 5, λ 3 decreases and its slope reduces as well, indicating that RBCs at C = 5 are stretched less than those at C = 1.Suspensions with C > 5 show similar shapes as those for C = 5.In comparison to RBCs in a suspension, single cells in tube flow have a larger (smaller) elongation than for C = 1 (C = 5).FIG.9: Definition of the inclination angle θ and the material angle θ M .Eigenvector v 1 of the gyration tensor (red arrow), the vector r c from the tube central axis to the cell's center of mass (blue arrow), and the material vector r M from the cell's center of mass to a specific material point (pink arrow) are shown.
In addition to deformation in flow, cell orientation changes and membrane exhibits a tank-treading motion.To quantify RBC orientation with a varying shape, we include in the analysis only the shapes with an asphericity larger than that of a biconcave shape at rest (A > 0.15), where a well-defined axis from the eigenvector v 1 corresponding to λ 1 can always be obtained.v 1 often does not lie within the flow-velocity gradient plane, due to cell-cell interactions and complex cell relaxation under fluid stresses.To simplify the quantification of cell dynamics, we define an inclination angle θ (marked in red in Fig. 9) as the angle between the vector v 1 and the flow direction x within the flow-velocity gradient plane (i.e.within x − r c plane, where r c is the vector from the tube central axis to the cell's center of mass).Figure 10(a) shows θ as a function of cell position in the flow direction for one selected cell from each suspension.θ may frequently exhibit discontinuous jumps of 180 degrees due to the symmetric disk-like rest shape, so that even small deformations can cause a reversal of the v 1 vector.Despite these jumps, RBC tumbling in the SC suspension can clearly be identified from the green trajectories in Fig. 10(a).Tumbling motion is infrequent for soft-RBC suspensions, where jumps between +90 o and −90 o are found instead and the departure of θ from these two values is generally within 30 degrees.Rare tumbling of RBCs for C = 20 suspension can be seen by the brown line in Fig. 10(a), but the tumbling period is extremely long (> 120D).Thus, cell-cell interactions in a crowded environment strongly hinder solid-like tumbling motion and RBCs are forced to relax fluid stresses via shape deformation instead.
The average orientation angle of RBCs as a function of radial position can be tracked through the absolute value of |θ|, as shown in Fig. 10(b).Near the RBC-FL, the average angle θ approaches 90 o , so that RBCs are aligned along the flow direction.θ decreases toward the tube center, and near the center the orientation angle is often not well defined as RBCs attain shapes with A < 0.15.The average angle θ for C ≥ 5 is slightly larger than for C = 1, indicating that RBCs are more aligned with the flow for large C values, which can also be seen in Fig. 7.
Strong dynamic changes in RBC shape significantly complicate cell orientation analysis and do not always allow decoupling between tumbling motion with shape rotation and tank-treading motion with membrane circulation.To analyze relative membrane motion, we monitor the material angle θ M defined as the angle between the material vector r M and the flow direction, see Fig. 9.Note that we cannot fully distinguish between tank-treading and tumbling motion by using θ M .Changes in θ M for a selected RBC are shown in Fig. 11(a) for various viscosity contrasts.θ M exhibits a nearly monotonic increase with increasing C. For C = 1, a smooth dependence of θ M is observed, indicating continuous tank-treading motion of the membrane.For C ≥ 5, an increase in θ M is slower than for C = 1, and membrane rotation shows short frequent pauses followed by periods of rapid θ M increase.Membrane motion can be quantified further by an average change in θ M as RBCs flow, which is defined as the ratio ∆θ M /∆x com between the change in θ M and the corresponding change in x com calculated for a fixed time interval ∆t ≈ 8 ms.where η rel increases as Ca decreases.Furthermore, η rel for C = 5 suspension is only slightly larger than η rel for C = 1 suspension and the difference is most prominent at intermediate flow rates.The converged RBC-FL thickness δ f decreases as the flow slows down.The C = 1 suspension has a smaller RBC-FL thickness than the C = 5 suspension with a difference of about 350 nm at Ca = 0.61 and about 100 nm at Ca = 0.06.
Table III shows that η rel increases as H t increases.The difference in η rel caused by C is most prominent for dense suspensions.The development of RBC-FL for different H t values is shown in Fig. 12, where δ f decreases for increasing H t , consistently with the increase of η rel .δ f is nearly the same for C = 1 and C = 5 suspensions.The difference between C = 1 and C = 5 for η rel and δ f remains similar for hematocrits H t = 15%, 30%, and 45%.

IV. SUMMARY AND DISCUSSION
The main focus of our study is the effect of the viscosity ratio C on the behavior of RBC suspensions.Interestingly, C ∈ [1,20] has only a very weak effect on the resistance of converged flow, quantified by the relative viscosity η rel .We have systematically analyzed the RBC-FL thickness, which is closely related to the flow resistance.When C increases from 1 to 5, the change in η rel is within 2%, while δ f increases by about 10 − 15% (or by approximately 300 nm).A further increase of C from 5 to 20 leads to an increase of η rel by about 6%, while δ f stays nearly unaffected.Thus, a larger viscosity of the RBC core for C ≥ 5 in comparison to C = 1 is complemented by a larger RBC-FL thickness, such that the overall flow resistance remains nearly independent of C.
The development of RBC-FL is governed by the two main mechanisms [28]: (i) RBC migration away from the wall due to the hydrodynamic lift force and (ii) the dispersion of RBCs within the cellular core due to cell-cell interactions in flow.The migration of RBCs away from the wall can be attributed to their shape deformation and membrane dynamics.For instance, RBCs with C = 1 migrate faster than those with C = 5, which is consistent with a slower development of the RBC-FL for C = 5 in Fig. 2(a).This results from the fact that both shape deformation and membrane dynamics are suppressed as C increases from 1 to 5, since a larger internal viscosity dampens RBC shape changes and dynamics.The importance of cell deformability for migration and the formation of RBC-FL is further illustrated by the results for SC suspension.Stiffened RBCs exhibit tumbling dynamics, which leads to a weak lift force, small RBC-FL thickness, and a large flow resistance.
The dispersion of RBCs within the cellular core counterbalances the lift force and is governed by local shear rate, shape deformations, and membrane dynamics.Both shape deformations and membrane dynamics are attenuated for C = 5 in comparison to C = 1.For example, the shape asphericity A for C = 5 is smaller than for C = 1, indicating that RBCs are less stretched.Furthermore, membrane dynamics for C = 5 is slower than for C = 1, because membrane tank-treading is more pronounced for the low viscosity contrast.This leads to a weaker hydrodynamic repulsion between RBCs in the C = 5 suspension, whose origin is similar to the lift force near a wall.To further confirm that the dispersion of RBCs is larger for C = 1 than for C = 5, we have performed several simulations of the collision of two cells in tube flow.After subtracting the migration effect, we find that the collision of two RBCs leads to a larger change in lateral cell displacement for C = 1 than for C = 5, indicating a stronger dispersion effect for C = 1.Even though cell collisions in dense suspensions involve more complex multicellular interactions, the lateral mobility coefficient R T in Fig. 5(b) is smaller for C = 5 than for C = 1, which is consistent with the simulations of binary collisions.Therefore, the cellular core for C = 5 remains more compactly packed than for C = 1, which is documented by a slightly thicker RBC-FL.The importance of RBC dynamics in flow is further emphasized by the behavior of stiffened RBCs.Hardened cells exhibit tumbling dynamics even in the cellular core, which results in their significant dispersion within the tube and a thin RBC-FL.
The dependence of the converged RBC-FL thickness δ f on the flow rate is consistent with the discussion of cell deformation and dynamics above, such that a decrease in driving pressure gradient weakens cellular dynamics and the exerted lift force, resulting in a reduction of δ f and an increase of flow resistance.Thus, the flow resistance is larger in the venular part of microvasculature, where blood-flow velocities are significantly smaller than those in the arteriolar part of the microvasculature.As expected, an increase in hematocrit leads to a decrease of δ f due to increased cell crowding, accompanied also by a slightly attenuated RBC dynamics.Note that the difference in η rel between C = 1 and C = 5 is more pronounced for H t = 45% than for smaller H t values.However, this difference can likely be neglected, because microvascular hematocrit values are generally smaller than 35%.A decrease of tube diameter D causes more pronounced cell-cell and cell-wall interactions, which result in a reduction of flow resistance, consistently with the Fahraeus-Lindqvist effect [20,21].Note that the viscosity ratio C has a more significant effect on flow resistance for large tube diameters (D > 40 µm), where η rel for C = 5 can be larger by more than 10% than for C = 1.
An interesting observation from our investigation is that the final RBC-FL thickness first increases with increasing C and then decreases for stiffened RBCs (this can be approximated by C → ∞).As discussed above, the initial increase in δ f for increasing C is the main reason that the flow resistance is nearly unaffected by C ∈ [1,20].In fact, this property of flowing blood allows the maximization of hemoglobin content in RBCs, and therefore oxygen delivery, without negative effects on the flow resistance.Furthermore, the effect of C in this range on blood-flow convergence is rather weak, such that RBC-FL convergence is obtained after a distance of about 25D − 50D.The flow-convergences distance can easily be compared with an average length of vessels (0.5 − 1 mm) within a branching network-like microvasculature [1,4], e.g.L e = 50D = 1 mm for D = 20 µm.This means that a converged flow within the microvasculature can only be expected in small vessels such as capillaries, while in microvessels with a diameter larger than about 10 − 20 µm, blood flow would likely always correspond to a transient (non-converged) flow.Even though the dependence of η rel on C ∈ [1,20] is nearly negligible, the structure and dynamics of RBC suspension are different for various C values.These differences in flow behavior for different viscosity contrasts are likely to be important for the margination of particles (e.g., platelets, drug-delivery carriers) in blood flow as well as for partitioning of RBCs within a complex microvascular network.In particular, the slow flow convergence behind branching points implies that many flow properties, such as resistance, particle margination, and oxygen delivery, can be highly inhomogeneous in complex vessel networks, depending on vessel diameters, branching lengths, and distance from the previous branching point.These aspects of blood flow need still to be addressed in future research.

FIG. 1 :
FIG. 1: Simulation snapshots for D = 20 µm and H t = 30%.(a) Configuration of dispersed RBCs, before the flow is applied.(b) Converged flow of a RBC suspension for C = 1 and Ca = 0.61 ( γ = 377.4s −1 ).(c) Converged flow for C = 5 and the same Ca as in (b).The flow is from left to right.

FIG. 2 :FIG. 3 :
FIG. 2: Evolution of the RBC-FL thickness in flow as a function of average convergence length L e = vt for C = 1, 5, 10, 20 and a suspension of stiffened RBCs.(a) Transient behavior at the beginning of RBC-FL development.(b) Dynamics of the RBC-FL thickness over the total simulation time.H t = 30%, D = 20 µm, and Ca = 0.61 ( γ = 377.4s −1 ).

FIG. 5 :
FIG. 5: Characteristics of RBC mobility perpendicular to the flow direction.(a) Trajectories of selected individual cells whose initial positions are closer to the tube center (top) and closer to the wall (bottom).(b) Distributions of the dimensionless lateral mobility coefficient R T (r) as a function of radial position within the tube.H t = 30%, D = 20 µm, and Ca = 0.61 ( γ = 377.4s −1 ).

FIG. 6 :
FIG. 6: Lift velocity v l normalized by local shear rate γ(r) for a single RBC migrating away from the wall as a function of the distance d w from the wall to the cell's COM in tube flow with D = 20 µm.Dashed lines indicate the power-law functions of d −1 w and d −2 w in the log-log plot.

FIG. 8 :
FIG. 7: Representative snapshots of RBC shapes near the edge of the cellular core and in the tube center for different RBC suspensions.Each column denotes a specific C value.H t = 30%, D = 20 µm, and Ca = 0.61 ( γ = 377.4s −1 ).

Figure 11
(b) presents ∆θ M /∆x com , where the rotation tendency is most prominent near the RBC-FL and decreases toward the tube center.A slight enhancement of ∆θ M /∆x com can be observed near the tube center.Large C values suppress membrane rotation, since a large difference in ∆θ M /∆x com is observed between C = 1 and C = 5 suspensions.The dependence of ∆θ M /∆x com is similar for all C ≥ 5 suspensions.In comparison to a single cell dynamics in tube flow, the membrane rotation for cell suspensions is suppressed by the surrounding cells.D. Dependence of the RBC-FL thickness on other flow parametersTo study the dependence of RBC-FL thickness on other flow parameters, we compare suspensions with C = 1 and C = 5 for different flow conditions Ca ∈ [0.06, 0.61] ( γ ∈ [37, 378] s −1 ), hematocrits H t ∈ [15, 45] %, and tube diameters D ∈[10,40] µm.

TABLE I :
Relative suspension viscosity η rel and final RBC-FL thickness δ f as a function of C. Ht = 30%, D = 20 µm, and Ca = 0.61 ( γ = 377.4s −1 )."SC" denotes stiffened cells, whose shear modulus is increased 100 times in comparison to that of a healthy RBC.
Table II presents relative viscosity η rel and final RBC-FL thickness δ f for various Ca,

TABLE II :
Relative viscosity η rel and final RBC-FL thickness δ f for various flow conditions characterized by Ca.D = 20 µm and Ht = 30%.