Crossover from tumbling to tank-treading-like motion in dense simulated suspensions of red blood cells

Via computer simulations, we provide evidence that the shear rate induced red blood cell tumbling-to-tank-treading transition also occurs at quite high volume fractions, where collective effects are important. The transition takes place as the ratio of suspension stress to the characteristic cell membrane stress exceeds a certain value, independent of volume fraction and cell deformability. This value coincides with that for a transition from an orientationally less ordered to a highly ordered phase. The average cell deformation does not show any signature of the transition, but rather follows a simple scaling law independent of volume fraction.


I. INTRODUCTION
Despite the large interest in a better understanding of the circulatory system and related diseases, there are still many open issues regarding the microscopic mechanisms which ultimately determine the rheology of suspensions of red blood cells (RBCs), the major particulate constituent of blood.Focusing on purely mechanical aspects, a RBC can be considered as a thin and flexible incompressible biconcave membrane which encloses an internal fluid (haemoglobin solution) and resists shearing and bending.
In the present work, we focus on a dynamic phenomenon in dense RBC suspensions.Via computer simulations, we provide evidence for the transition from a tumbling to a tank-treading-like state (henceforth called TB-TT transition) upon increasing the capillary number Ca.The capillary number is defined as the ratio of fluid stress to the characteristic membrane stress, where the latter is related to the shear elasticity of the cell.Transcending previous studies, we find that the TB-TT transition occurs not only in the well-known case of a dilute suspension, but up to haematocrit values as high as Ht = 65% (haematocrit, or volume fraction, is the ratio of the volume occupied by RBCs to the total suspension volume).It is shown that-for all studied shear rates, cell deformabilities, and volume fractions-this transition is characterised by an effective capillary number Ca * (ratio between effective suspension stress and the characteris-tic membrane stress) rather than by the bare capillary number Ca.A detailed analysis of the average RBC inclination angle θ (a measure of average cell orientation) and the corresponding order parameter Q > is also provided.When plotted versus Ca * , all values for Q > and θ tend towards a master curve as Ca * exceeds a certain threshold value Ca * cr which, remarkably, coincides with that for the onset of the TB-TT transition for a single cell.
Our results provide evidence that the change from tumbling to tank-treading-like dynamics is accompanied by a phase transition from weak to strong orientational ordering of the cells.Interestingly, the average deformation of a cell, quantified by the Taylor deformation parameter D a changes continuously in the investigated parameter range, without any signature of the observed transition.
The article is organised as follows.The setup of the simulations is discussed in section II.Section III contains the results and discussions of the TB-TT transition (section III A), the orientational order (section III B), and the RBC deformation (section III C).We conclude our work and provide an outlook in section IV.

II. SIMULATION SETUP A. Numerical model
A hybrid immersed-boundary-lattice-Boltzmannfinite-element model has been used [16,17].The liquid phase is fully considered and solved by the lattice-Boltzmann method [18,19], while the fluid-particle interaction is realised by the immersed-boundary method [20].The red blood cells are represented by a finite-element mesh consisting of 1620 triangular facets.Physically, the cell membrane is made of a lipid bilayer (being essentially an incompressible 2D fluid), giving rise to a finite resistance against bending and changes of surface area, and a cytoskeleton, leading to a finite shear In the central region (grey, from Lz/4 to 3Lz/4), the volume fraction is nearly constant, the velocity profile is linear (red solid line), and the stress is constant.
resistance [21][22][23].The ensueing effects on the collective and rheological properties of the RBC suspension can be captured by assuming a total membrane energy of the form where describes the energy penalty against shear and area dilation, is energy change due to bending, while and penalize changes of membrane surface area A and cell volume V over their undeformed counterparts, A (0) and V (0) , respectively.Eq. ( 2) represents Skalak's energy density [24] with κ S and κ α being the shear and area resistance, respectively.I 1 and I 2 denote the in-plane strain invariants, which are related to the local membrane deformation tensor (see Krüger et al. [16] for more

Ca η/η0
Ht = 35%, s Ht = 45%, s Ht = 55%, s Ht = 65%, s Ht = 35%, r Ht = 45%, r Ht = 55%, r Ht = 65%, r FIG. 2 Relative suspension viscosity η/η0 as function of bare capillary number Ca = η0 γr/κS.η0 is the reference viscosity of the suspending fluid without red blood cells.For all volume fractions, the suspension is shear thinning.Additionally, the viscosity increases with increasing volume fraction.One observes that the data for soft and hard RBCs overlap, when plotted versus Ca.Error bars are of the order of the symbol size.It is also noteworthy that the viscosities can be approximated by a Herschel-Bulkley law (dashed lines) of the form η/η0 = a + b × Ca c with the parameters (a, b, c) = (0, 2.0, 0.78) for 35%, (0.004, 2.4, 0.73) for 45%, (0.015, 3.4, 0.72) for 55%, and (0.043, 5.0, 0.71) for 65%.FIG. 3 Cross-section through an undeformed red blood cell (red) together with its equivalent inertia ellipsoid (black).r is the large RBC radius (dotted), a, b are the two large semiaxes (a = b = 1.1r) and c is the small semiaxis (c = 0.36r) of the inertia ellipsoid.The vector ô is perpendicular to the a-b-plane and is given by the normalized eigenvector of the inertia ellipsoid corresponding to the smallest eigenvector (i.e., the shortest semiaxis, c).The intrinsic membrane orientation vector î, on the other hand, is computed from the positions of the vertices of the RBC mesh and connects the centre of masses (green dots) of the top and bottom halfs of the undeformed membrane.The identity of a mesh vertex of being in the "top" or "bottom" half remains unaltered during the course of the simulation, making î an intrinsic measure of the RBC orientation that is not available from the inertia tensor.
details).Note that both the parameters κ α and κ A penalize changes in the membrane surface area (physically, κ α is related the cytoskeleton, while κ A arises from the incompressibility of the the liquid bilayer).The bending energy in Eq. ( 3) has the classical Helfrich form [25], of which a simplified version is employed in the simulations [26].Here, κ B is the bending modulus and the sum runs over all pairs of neighbouring facets of the tessellated RBC surface with mutual equilibrium normal-to-normal angles θ eq ij .We remark that, in the present case, the biconcave shape of the RBC is not a result of the minimization of the bending energy subject to the constraint of fixed volume [23], but rather specified as input to the simulation via the construction of the membrane mesh (and ensured by the presence of the θ eq ij in Eq. ( 6)).In our simulations, a fixed number (494, 635, 776, and 917) of neutrally buoyant RBCs with equilibrium shape of a biconcave disk of radius of r = 9 (all quantities in lattice units) and equal internal and external viscosities have been distributed throughout a fixed fluid volume (L x × L y × L z = 100 × 100 × 160) resulting in volume fractions Ht = 35, 45, 55, and 65%, respectively.To improve numerical stability and avoiding particle overlap, a repulsive force inspired by Feng and Michaelides [27] is introduced.It is zero for distances larger than one lattice constant and behaves like 1/r 2 for smaller distances.The effect of the membrane energy, Eq. ( 1), is realized in terms of local forces derived from the principle of virtual work [16].
The shear flow (periodic in x-and y-directions) has been realised by moving the bottom and top walls at z = 0 and z = 160 in opposite directions along the xaxis with velocities ±u w , giving rise to an average shear rate of γ = 2u w /L z .Inertia is neglected.One layer of RBCs has been stuck to each wall to avoid slip.A typical simulation snapshot illustrating the geometry is shown in Fig. 1(a).
The elasticity parameter values κ S and κ B in simulation units can be found by matching the relevant dimensionless parameters of the problem (in this case the capillary number η 0 γr/κ S and the reduced bending modulus κ B /(κ S r 2 )).Two different RBC rigidities have been considered.The softer RBCs, (κ [23].The more rigid RBCs obey (κ S , κ B , κ α , κ A , κ V ) = (0.06, 0.012, 1, 1, 1).The former are denoted 's', the latter 'r'.The choice of κ α = κ A = κ V = 1 ensures that local area and volume deviations are restricted to a few percent.
In order to improve the statistics, all simulations have been repeated with various random initial RBC configurations.Ten and five independent runs have been performed for the softer and more rigid RBCs, respectively.All reported observables are averaged over independent runs and time in the steady state.Due to the absence of thermal fluctuations in the present model, all observed effects are shear-induced.

B. Characterisation of the flow
Since we are interested in bulk properties rather than wall-induced effects, we first examine RBC concentration, suspension velocity, and shear stress profiles to see whether a bulk-like behavior may be expected.
Fig. 1(b) reveals that, within the inner 50% of the volume between the walls (henceforth called the central region), the volume fraction is nearly constant and the velocity is linear, defining a constant shear rate.The shear stress is confirmed to be constant throughout the entire volume as expected for simple shear flow [17].Thus, an effective bulk shear viscosity η( γ) ≡ σ xz ( γ)/ γ can be defined as the ratio of shear stress σ xz and shear rate γ using the data in the central region (Fig. 2).For the present purposes, the viscosity is only needed to define an effective capillary number (see eq. ( 10) below), which will turn out to be the central quantity governing the cell dynamics.A detailed investigation of the rheology will be presented in another publication.
The results analysed in the central region are expected to be characteristic of the bulk properties of the system.In the following, system properties will thus be determined in the central region where all studied observables have been found to be rather independent of the transverse position z.

A. Transition from tumbling to tank-treading-like dynamics
We investigate the average RBC tumbling (rigid bodylike) frequency ω by tracking the orientation of the cells' inertia tensor in time.We recall the average tumbling frequency of a single rigid ellipsoid in a simple viscous shear flow, where p = a/c (a ≥ b ≥ c are the semi-axes of the ellipsoid, a and c are aligned with the shearing plane, see Fig. 3) [28].For a rigid sphere (p = 1), the tumbling frequency is ω/ γ = 1 2 .One finds for the inertia ellipsoid of an undeformed RBC a = b = 1.1r and c = 0.36r, and Jeffery's solution [28] predicts ω/ γ = 0.30.On the contrary, for purely tank-treading cells, the average is ω = 0, as they do not show any tumbling activity.
Additionally, we define the instantaneous reduced angular velocity of a RBC as ω * γ := where ac is approximately the RBC cross-sectional area in the shearing plane and   is the surface-averaged rotational component of the membrane velocity (r and v are location and velocity of a RBC surface element relative to the current RBC centre, respectively, and A is the RBC surface area).Note that ω * is sensitive to any form of membrane rotation, even if the cell's inertia tensor is not rotating in space (i.e., if there is no tumbling).
A first evidence for the onset of the TB-TT transition is provided in Fig. 4 where the reduced angular membrane velocity ω * / γ and the reduced tumbling frequency ω/ γ cr , a type of rotational motion is dominant which allows for a non-rotating tensor of inertia.Tank-treading-like dynamics provides such an alternative.Interestingly, for an isolated RBC we see a TB-TT transitionin the same region (Fig. 5), which is in qualitative agreement with previous works [29][30][31].Furthermore, at σ ≈ 0.4 Pa, which corresponds to Ca * ≈ 0.2 for healthy RBCs, the transition has been observed experimentally for single RBCs [10].Note that all data for Ca * > Ca * cr collapse onto a single master curve, which is not the case when plotted as function of the bare capillary number Ca (not shown).The decay of the tumbling frequency in the region Ca * ≈ 0.2-0.3 can be captured by a simple power law, ω/ γ ∝ Ca * −3 .We have to emphasize, however, that this fit is entirely ad-hoc.The establishment of such a power-law requires a more careful analysis which goes beyond the scope of the present work.
The picture of a TB-TT transition characterized by a critical capillary number without explicit dependence on haematocrit is further supported by a survey of the dynamics at the level of the individual cells.Fig. 6 exemplifies the typical temporal behaviour of an individ-ual RBC in moderately dense suspension, from which one can clearly identify tumbling motion at low effective capillary numbers, whereas tank-treading is performed at large Ca * .Interestingly, even at this comparatively large value of haematocrit (45%), the typical rotational motion of an individual cell does not appear to be significantly different from that of a single cell, apart from occasional hydrodynamical collisions with its neighbours.While a more detailed analysis of these aspects will be presented in a separate work, we point out here that a comparison of the rate of flips [32] of the equivalent inertia ellipsoid on the one hand and of a suitable intrinsic membrane orientation vector [47] on the other hand provides an independent means to distinguish tumbling from tank-treading behaviour.This is demonstrated in Fig. 7, where we plot the ratio of the total number of flips (obtained by integrating the flip rate over several decades in inverse shear rate) of the two orientation vectors.Consistent with the expectations from the analysis of the angular velocities, we find that for low effective capillary numbers, both the inertia ellipsoid and intrinsic membrane orientation vector flip at the same rate, indicative of tumbling-like motion.In contrast, at higher Ca * , only the intrinsic membrane orientation vector performs flips while the orientation of the inertia ellipsoid remains virtually fixed, as would be expected for tank-treading-like motion.
It deserves special notice that, while the TB-TT transition is fairly well understood in terms of flow gradient effects and cell deformability (it is expected to occur when the fluid stress is strong enough to push the cells through the maximum of the elastic energy) [2], its occurrence in a dense suspension, where collective effects play a major role, is a more complex phenomenon.Our results suggest that, even at a haematocrit of Ht = 65%, RBCs can perform a tank-treading-like motion in an effective medium whose viscosity is no longer the viscosity of the ambient fluid but the significantly higher effective suspension viscosity.It has also to be noted that, due to RBC collisions and related complicated dynamics, an unperturbed swinging or tank-treading motion for Ca Ca * cr is not expected.Instead, the cells' inclination and deformation fluctuate about their time averages, which can also be inferred from Fig. 6.Therefore, the term tank-treadinglike is used to emphasise that the cells are not tumbling, but neither show a perfect, unperturbed tank-treading motion.

B. Cell alignment and orientational ordering
The onset of the TB-TT transition seems to be accompanied by a transition in the orientational order parameter Q > of the system which is defined as the largest eigenvalue of the order tensor [34]  versus Ca * .The power-law fit, a/b|max − 1 = 0.89 × Ca * 0.9 , is shown as dashed line.There is a one-to-one correspondence between a/b and Da.The grey area denotes the TB-TT transition.Note the excellent agreement with the experimental data obtained for young RBCs reported in Fig. 2 in Pfafferott et al. [33] (we converted the shear stress reported therein to Ca * by assuming a shear elasticity of κS = 5 µN/m [23]).
The RBC orientation vector ô is defined as the inertia tensor eigenvector corresponding to the shortest semiaxis c (see Fig. 3).The related eigenvector of Q is called the director n indicating the average orientation of the cells.
As illustrated in Fig. 8, both Q > and the average inclination angle θ (average angle between the director and the flow axis) change rapidly at Ca * ∼ Ca * cr .We have also investigated the time evolution and spatial dependence of n and Q > .Both quantities show only slight fluctuations about their averages.Note that RBC ordering at high shear rates has been observed experimentally [1].Similarly to the behavior of the tumbling frequency, all simulated data follow a master curve for Ca * ≥ Ca * cr and can be described by one single rather than by two independent parameters (Ca * instead of Ca and Ht).
In qualitative agreement with analytical predictions [2,3,35] and previous simulations [36,37], the average inclination angle approaches 90 deg in the limit of vanishing capillary number [28] and grows with increasing capillary number toward the onset of tank-treading (it is bounded from above by 135 deg).In the tank-treading regime, the average inclination angle is expected to decrease again [35,36,38], as the capsule becomes more elongated and aligns further with the flow.For the limited range of capillary numbers accessible to present simulations, however, we do not observe such a behavior.Rather, the inclination angle is found to remain roughly constant or even slightly increases (which is, however, in line with the predictions of Skotheim and Secomb [2]).
A typical distribution of inclination angles is shown in Fig. 8(c).For comparison, the probability of finding the rigid inertia ellipsoid of an RBC (aspect ratio p = a/c) with a given orientation angle in shear flow is This relation can be extracted from the known expression for the angular velocity of the ellipsoid [28].As γ increases, the centre of P (θ) is shifted towards larger angles and its shape becomes narrower.While the former results in an increase of θ, the latter leads to a higher orientational order Q > since deviations from a given cell orientation become restricted to a narrower range.It is striking that, for Ca * Ca * cr , the average inclination angle equals that of the isolated cell (Fig. 8(b)).

C. Deformation behavior
In Fig. 9(a), some typical deformation probability distributions P (D a ) are shown.Here, is the Taylor deformation parameter as a measure for the RBC asymmetry in the a-b-plane, where the index zero refers to the undeformed state and a 0 = b 0 holds for an RBC.Using this data, we determine D max a , the deformation parameter with maximum probability.It is found (Fig. 9(b)) that the data on D max a collapse onto a single master curve when plotted as a function of Ca * for all studied values of control parameters.
In marked contrast to the behavior of tumbling frequency and orientational order parameter, D max a shows no signature of the TB-TT transition.A similar observation for dilute RBC suspensions has been reported before [10].The entire observed range of capillary numbers can be parameterised by a simple power law, a/b| max − 1 ∝ Ca * 0.9 , as shown in Fig. 9(b).

IV. CONCLUSIONS AND OUTLOOK
The focus of the present study is on the tumblingto-tank-treading transition (TB-TT transition) in suspensions of aggregation-free red blood cells in the limit of high volume fractions where collective effects become dominant.We provide evidence that this transition occurs when the ratio of the effective suspension stress to the characteristic membrane stress (effective capillary number) exceeds a threshold value Ca * cr .Note that all dependence on volume fraction is implicit in Ca * through its dependence on the effective stress.For a single cell, Ca * cr corresponds to the stress where the cell is driven through the maximum of its elastic energy, thus allowing tank-treading-like dynamics.
The average tumbling frequency and cell inclination angle signal the onset of the TB-TT transition and a scaling collapse above Ca * cr when expressed in dependence of the effective capillary number.Remarkably, the most probable cell deformation changes continuously in the investigated parameter range following a scaling law, ∝ Ca * 0.9 , without any signature of the observed transition.The value of Ca * cr coincides with that for a transition from a less ordered cell orientation distribution to a highly ordered phase.
The above findings support the following conclusions: (i) The cell dynamics is dominated by Ca * and therefore the suspension stress.The cells in the tank-treadinglike state behave more like isolated particles experiencing their environment only via the suspension stress.(ii) The large degree of orientational ordering at Ca * ≥ Ca * cr is related to the onset of tank-treading-like dynamics.(iii) Cell deformation alone does not seem to be a relevant factor for the TB-TT transition.
Interestingly, a scaling of static and dynamic quantities in terms of an effective, concentration-corrected dimensionless number is not uncommon in soft matter systems and has been observed previously, for instance, for polymer solutions [39][40][41].An important difference between these and the present system is, however, the volume and area conservation of the RBCs, which implies that excluded volume effects should play a much more dominant role for RBC suspensions.This might be one of the reasons why many quantities lack scaling with Ca * in the tumbling regime, where, in constrast to the tanktreading case, mutual disturbance of the cells is expected to be important.
During tank-treading, particles explore a smaller volume as compared to tumbling, thus providing a less strong hindrance to the motion of other particles.This is expected to affect other quantities, such as stresses and diffusivity, which are sensitive to the interactions between particles.Indeed, it is commonly expected that the suspension viscosity is reduced as the number of tanktreading particles increases [1,10,14,42].
Furthermore, our work opens a route to equipping coarse-grained blood models [43,44] with proper constitutive laws for the tumbling/tank-treading probability or the average particle deformation given an ambient stress value.
Finally, while we focused on the limit of high Ca * , many interesting open questions remain regarding the opposite limit of small shear rates such as the possible existence of a yield stress or flow heterogeneity [45,46].
FIG. 1 (a) Snapshot of a simulation with average haematocrit (volume fraction) Ht = 55% and capillary number Ca = 0.011.The bottom (top) wall is located at z = 0 (z = Lz = 160) and moves to the left (right) with velocity ±uw (black arrows).(b) Layer-resolved haematocrit, Ht, shear stress, σxz, and flow velocity, ux, for the same control parameters as in panel (a).Each curve is obtained as average over steady state and independent runs.In the central region (grey, from Lz/4 to 3Lz/4), the volume fraction is nearly constant, the velocity profile is linear (red solid line), and the stress is constant.

FIG. 4
FIG. 4 Reduced (a) angular velocity ω * / γ and (b) tumbling frequency ω/ γ versus effective capillary number Ca * .The legend applies to both panels.The grey area denotes the TB-TT transition.As seen in panel (b), for Ca * > 0.1, the average tumbling frequency of a single isolated cell is zero, whereas ω/ γ approximately decays like ω/ γ ∝ Ca * −3 (inclined solid line) in the dense suspensions.Contrarily, the angular membrane velocity in (a) is rather independent of Ca * .Note that ω * / γ and ω/ γ are nearly identical for Ca * 0.1.The analytic values of ω/ γ for a rigid sphere and ellipsoid with aspect ratio p = a/c = 1.1/0.36 are also shown in (b).

10 FIG. 5
FIG.5This sequence shows the time evolution of the rotational behavior of a single RBC in an external shear flow with rate γ.The cross-section is parallel to the shearing plane, and the vorticity of the shear flow is clockwise.The evolution as function of dimensionless time γt is shown for three different capillary numbers, Ca = 0.1 (loosely dashed line), 0.2 (densely dashed line), and 0.5 (solid line).For the smallest value (Ca = 0.1), the RBC performs a tumbling motion since it is not sufficiently deformed to undergo a tank-treading motion.However, already for Ca = 0.2, the RBC can rotate without tumbling, but shape oscillations are still observed.The tank-treading motion is fully developed for the highest capillary number (Ca = 0.5).Note that in the dilute limit Ca * ≈ Ca holds.

FIG. 6
FIG. 6 Exemplary cross-sections of a rotating RBC in a 45%-suspension at (a) Ca * = 0.05 and (b) Ca * = 0.49.The cross-sections correspond to the shearing plane and run through the cell centre.The shear flow is in clockwise direction.In the left panel, a flipping event (rigid body-like rotation) is shown (temporal order: red, blue, green; time difference between snapshots is 0.14/ γ).The right panel illustrates the irregular tank-treading-like rotation (temporal order: red, blue, green; time difference between snapshots is 0.35/ γ).Shape fluctuations due to collisions with neighbouring particles are clearly visible.Arrows indicate the local surface velocity.

FIG. 7
FIG.7Ratio of the number of flips (i.e., half-turns) of the inertia ellipsoid, N flips,IE , and intrinsic cell orientation vector (see Fig.3), N flips,orient , of RBCs in a suspension at different haematocrit values.At low effective capillary numbers, both inertia ellipsoid and orientation vector flip at the same rate (the flip number ratio being close to unity), characteristic for tumbling-like motion.In contrast, at high Ca * only the intrinsic orientation vector performs flips, while the flip rate of the inertia ellipsoid is strongly reduced, indicating tank-treading-like motion.

FIG. 8
FIG. 8 (a) Order parameter and (b) average inclination angle versus effective capillary number Ca * .The curves obtained for one single cell are also shown.The grey area denotes the TB-TT transition.The legend applies to both panels.(c) Inclination angle probability distribution for the softer cells with Ht = 55% and selected values of Ca * .The analytic curve for a single rigid ellipsoid (aspect ratio p = a/c = 1.1/0.36) is also shown.The single cell in panel (b) approaches 90 • for decreasing Ca * , as expected from the analytical solution.