## Spatio-temporal dynamics of dilute red blood cell suspensions in a microchannel flow at low Reynolds number

Qi Zhou a,1 , Joana Fidalgo b,1 , Lavinia Calvi , Miguel O. Bernabeu , Peter R. Hoskins , Mónica a c d S. N. Oliveira b, GLYPH&lt;3&gt; , Timm Krüger a, GLYPH&lt;3&gt;

a School of Engineering, Institute for Multiscale Thermofluids, University of Edinburgh, Edinburgh EH9 3FB, UK b James Weir Fluids Laboratory, Department of Mechanical and Aerospace Engineering, University of Strathclyde, Glasgow G1 1XJ, UK

c Centre for Medical Informatics, Usher Institute, University of Edinburgh, Edinburgh EH16 4UX, UK d Centre for Cardiovascular Science, University of Edinburgh, Edinburgh EH16 4SB, UK

## Abstract

Microfluidic technologies are commonly used for the manipulation of red blood cell (RBC) suspensions and analyses of flow-mediated biomechanics. To enhance the performance of microfluidic devices, understanding the dynamics of the suspensions processed within is crucial. We report novel aspects of the spatio-temporal dynamics of RBC suspensions flowing through a typical microchannel at low Reynolds number. Through experiments with dilute RBC suspensions, we find an o GLYPH&lt;11&gt; -centre two-peak (OCTP) profile of cells contrary to the centralised distribution commonly reported for low-inertia flows. This is reminiscent of the well-known 'tubular pinch e GLYPH&lt;11&gt; ect' which arises from inertial e GLYPH&lt;11&gt; ects. However, given the conditions of negligible inertia in our experiments, an alternative explanation is needed for this OCTP profile. Our massively-parallel simulations of RBC flow in real-size microfluidic dimensions using the immersed-boundary-lattice-Boltzmann method (IB-LBM) confirm the experimental findings and elucidate the underlying mechanism for the counterintuitive RBC pattern. By analysing the RBC migration and cell-free layer (CFL) development within a high-aspect-ratio channel, we show that such a distribution is co-determined by the spatial decay of hydrodynamic lift and the global deficiency of cell dispersion in dilute suspensions. We find a CFL development length greater than 46 and 28 hydraulic diameters in the experiment and simulation, respectively, exceeding typical lengths of microfluidic designs. Our work highlights the key role of transient cell distribution in dilute suspensions, which may negatively a GLYPH&lt;11&gt; ect the reliability of experimental results if not taken into account.

## 1. Introduction

Red blood cells (RBCs) occupy around 45% of the blood volume and play a crucial role in maintaining the normal function of tissues and organs in humans and animals. They transport oxygen and carbon dioxide, regulate blood viscosity, and a GLYPH&lt;11&gt; ect the circulation of nutrients and immune cells. Pioneered by Chien and Skalak in the 1960s-1980s [1, 2, 3], there has been

GLYPH&lt;3&gt; Corresponding author e-mail: monica.oliveira@strath.ac.uk or timm.krueger@ed.ac.uk

1 These authors contribute equally to this work.

May 8, 2020

extensive research carried out theoretically, experimentally and numerically to elucidate the behaviour of single RBC and characterise the dynamics of RBC suspensions (see recent reviews [4, 5, 6, 7, 8]). However, despite major progress, a rigorous and quantitative connection between microscale RBC dynamics and macroscale haemorheology is still lacking [9].

Up to date, various factors a GLYPH&lt;11&gt; ecting blood microrheology have been determined by theoretical numerical models of RBCs or biomimetic vesicles, including (but not limited to): (i) viscosity / contrast (between inner and outer fluids) [10], deformability [11], orientation [12] and initial position [13] of the cell; (ii) shear component (distinct in linear quadratic flows) [14], suspending / viscosity [15], flowline curvature [16] and wall confinement [17] of the flow; (iii) two-body or multi-body hydrodynamic interactions between cells [18]. Most studies focus on the individual dynamics pairwise interaction of cells or the characteristic signatures of the overall suspension / system, e.g., e GLYPH&lt;11&gt; ective viscosity and normal stress di GLYPH&lt;11&gt; erences.

Much less attention has been paid to the collective behaviour of cells, e.g., their spatiotemporal organisation or local microstructures, at a realistic length scale in three-dimensional space. Common reasons for avoiding such studies include their complexity of analytical modelling or prohibitiveness of computational cost. In fact, the RBC behaviour at sub-suspension level plays a key role in bridging the RBC dynamics and haemorheology. From an experimental perspective, the maturing of microfluidic technologies over the last two decades has provided a powerful tool to explore the otherwise inaccessible behaviour of complex biofluids at the microscale. It is now widely employed for blood flow experiments aimed at various applications [19, 20, 21]. For practical reasons in microfluidics, dilute (concentration GLYPH&lt;30&gt; &lt; 10%) or semidilute RBC suspensions are commonly used, where the clogging of samples is minimal and the measurements of individual cells are quantitatively accurate [22, 23, 24]. Through observation of dilute semi-dilute RBCs in polydimethylsiloxane (PDMS) microfluidic channels, many / intriguing properties of RBC flow at the microscale have been demonstrated, such as uneven cellplasma partitioning [25], long-range hydrodynamic interaction [26], wall-induced lift of cells in tumbling motion [27], cytoplasmic viscosity-promoted polylobes on cell membrane [28] and di GLYPH&lt;11&gt; usion-governed collective spreading of cell suspensions [29].

From experiments with a semi-dilute RBC suspension ( GLYPH&lt;30&gt; = 16%), Sherwood et al. [30] revealed strongly heterogeneous distribution of cells, which presented locally skewed rather than axisymmetric haematocrit ( Ht ) profiles. With much lower RBC concentration ( GLYPH&lt;30&gt; GLYPH&lt;20&gt; 5%), Shen et al. [31] found an even more heterogeneous RBC distribution eventually leading to inversion of the classic Zweifach-Fung e GLYPH&lt;11&gt; ect (governing cell partitioning) downstream of bifurcations. More recently, Iss et al. [32] and Abay et al. [33] discovered emergent self-organisation or focusing behaviour of cells in dilute RBC suspensions arising from geometric features of the microfluidic channel, such as confinement and constriction. These studies imply that the spatial arrangement of RBCs in microchannels can be highly dynamical and varies from one location to another. Such local characteristics, if ignored, may give rise to inaccurate estimation of the haemodynamics in the microfluidic channels, and a GLYPH&lt;11&gt; ect the physiological relevance of the in vitro experiments designed to imitate in vivo systems.

The present work aims to elucidate the spatio-temporal distribution of RBCs in a dilute suspension flowing through a typical microfluidic channel of rectangular cross-section. For a dilute suspension ( Ht GLYPH&lt;20&gt; 1%) of horse RBCs in the microchannel, cell density profiles are obtained by a direct cell-counting routine developed in-house . Resultantly, we observe a counterintuitive cell distribution with o GLYPH&lt;11&gt; -centre two-peak ordering rather than a single peak near the centreline, even in the negligible inertia regime. Immersed-boundary-lattice-Boltzmann simulations performed for an equivalent system confirm the experimental finding and provide insight into the underly-

Figure 1: (a) Schematic of the experimental setup. (b-c) Inlet configuration of the microfluidic device, showing optical images of the flow entry region where the inlet port is (b) aligned and (c) misaligned relative to the channel centreline. In (b), the coordinate axes are shown. (b-c) share the same scale bar of 500 µ m. (d) Schematic showing designated regions of interest (ROIs) for the measurement of probability density distribution (PDD) of RBCs in the suspension. (e) An individual frame of RBCs flowing in the microchannel and (f) a composite image combining 300 frames of the same region, obtained using the Z-projection method based on minimum intensity. In (e-f), Q = 0 2 : µ L min, and they share / the same scale bar of 100 µ m.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00250v2-figure-48e7fb4b-202f-4450-aadd-a77eeb076231.jpg)

ing mechanisms. The paper is organised as follows: the experimental and numerical methods are introduced in Sec. Materials and Methods . In Sec. Counterintuitive RBC distribution , cell density profiles from experiments and simulations are presented. Sec. Flow properties and RBC dynamics provides analysis of the observed cross-streamline migration and explores the mechanism behind the density profiles. Sec. Discussion compares the present study with existing works and states the key implications of our findings. Sec. Conclusion concludes the paper.

## 2. Materials and Methods

The materials and methods used for the experiments and the simulations involved in this study are outlined in Sec. Microfluidic experiments and Sec. Computational modelling , respectively. More details are provided in the Supplemental Information (SI).

## 2.1. Microfluidic experiments

RBC suspensions in Dextran40 were prepared using horse blood (see Fig. S1, SI for Photomicrography of the RBC suspension) provided by TCS Biosciences, in anticoagulant EDTA (1.5 mg mL). The procedures for sample processing are described in Fig. S2 and Sec. / Preparation of RBC suspensions . The steady shear rheology of Dextran 40 was determined by a rotational rheometer (DHR2, TA Instruments, with minimum resolvable torque 0.1 nNm) using a cone-plate geometry of 60 mm diameter and angle GLYPH&lt;18&gt; = 1 GLYPH&lt;14&gt; (CP60 1), 30 / µ m truncation gap. The density and shear viscosity values used to calculate dimensionless numbers (e.g., Reynolds number Re and Capillary number Ca ) for experiments are 1.07 g mL and 5.0 mPa s, respectively. /

The microfluidic device was manufactured by soft-lithography using an SU8 mould produced by photo-lithography (microLIQUID). The microfluidic channel was made in PDMS (polydimethylsiloxane, Sylgard 184, Dow Corning) and bonded to a glass slide substrate with oxygen plasma (Zepto, Diener). The channel has a rectangular cross-section, with the width and

depth dimensions ( W = 87 µ m, H = 31 µ m) determined by optical microscopy. Fig. 1a shows a schematic of the experimental setup, where the sample is infused into the inlet of the microdevice using the withdraw function of the syringe pump located at the outlet. With this proposed system, the sample-containing Eppendorf can be gently stirred at regular time intervals to minimise cell sedimentation, which typically limits the duration of experiments. To test the e GLYPH&lt;11&gt; ect of inlet conditions, inlet ports are manually punched onto the PDMS, so that the tubing connection can be either aligned (Fig. 1b) or misaligned (Fig. 1c) with the centreline of the microdevice. The results presented in this study are for the aligned case, unless otherwise stated.

The experimental setup used an inverted microscope (Olympus IX71) for flow visualisation and a syringe pump (Nemesys, Cetoni) for fluid control. The microfluidic device was illuminated using a halogen lamp. Videos were acquired by a sensitive monochrome CCD camera (Olympus, XM10) using a 10 GLYPH&lt;2&gt; magnification objective with numerical aperture of 0.25 to guarantee that the depth of field and image depth are large enough to observe cell trajectories across the entire channel height (Fig. 1e and Sec. Data acquisition and image analysis ). An algorithm was developed in-house based on ImageJ software to detect and count the cells automatically (Fig. S3 and Fig. S4 in Sec. Data acquisition and image analysis ). The RBC probability density distribution (PDD) along the width direction of a given region of interest (ROI, see Fig. 1d) is then calculated using a minimum of 300 frames, containing more than 8,000 cells (Fig. S5 in Sec. Probability density distribution (PDD) of RBCs ). The cell-free layer (CFL) is determined as half the di GLYPH&lt;11&gt; erence between the overall channel width and the RBC-core width, using composite images (see Fig. 1f) generated via the Z-projection method. The detailed procedure is described Fig. S8 and Sec. Cell-free layer (CFL) .

## 2.2. Computational modelling

The immersed-boundary-lattice-Boltzmann method (IB-LBM [34]) is employed to model blood flow as a suspension of deformable RBCs. The fluid flow governed by the Navier-Stokes equations is solved by the lattice-Boltzmann method with standard D3Q19 lattice [35], BGK collision operator [36], Guo's forcing scheme [37] and the Bouzidi-Firdaouss-Lallemand (BFL) implementation of no-slip boundary condition at the walls [38]. To precisely control the volume flow rate within the channel, inflow outflow open boundaries are implemented with the Ladd / velocity boundary condition [39]. The RBCs are modelled as Lagrangian membranes using a finite element approach. The fluid flow and RBC membrane are coupled through the immersedboundary method [40], which tackles the interpolation of velocities and the spreading of forces. The model is implemented in the parallel lattice-Boltzmann flow simulator HemeLB [41] (open source at http://ccs.chem.ucl.ac.uk/hemelb ) by incorporating a new module for discrete deformable cells. Details about the model configuration (Fig. S10), simulation setup (Fig. S11) and numerical analysis (Fig. S12-S14) can be found in Sec. Model configuration , Sec. Simulation setup and Sec. Numerical analysis . See Table 1 for the key simulation parameters.

Each RBC is modelled as a closed finite-element membrane with 720 triangular facets (meshing resolution matching the lattice size GLYPH&lt;1&gt; x for numerical stability and accuracy [42]), with its unstressed shape assumed as a biconcave discoid. The RBC membrane is elastic, with its mechanical properties controlled by five moduli which govern di GLYPH&lt;11&gt; erent energy contributions. They are strain modulus GLYPH&lt;20&gt; s , bending modulus GLYPH&lt;20&gt; b , dilation modulus GLYPH&lt;20&gt; GLYPH&lt;11&gt; , surface modulus GLYPH&lt;20&gt; S and volume modulus GLYPH&lt;20&gt; V . For more details of the RBC model, please refer to [34]. To tackle close contact among cells or between the cell and the wall (albeit rare in the dilute suspensions we consider), a repulsion potential decaying with inverse distance between neighbouring RBC surfaces is numerically implemented with interaction intensities (see GLYPH&lt;20&gt; cc and GLYPH&lt;20&gt; cw in Table 1) comparable to the

Table 1: Key simulation parameters. The symbol ' GLYPH&lt;24&gt; ' represents simulation units (dimensionless).

| Parameter                | Symbol                                               | Value                         |
|--------------------------|------------------------------------------------------|-------------------------------|
| channel geometry         |                                                      |                               |
| width                    | W                                                    | 96 µ m                        |
| depth                    | H                                                    | 30 µ m                        |
| length                   | L                                                    | 1320 µ m                      |
| hydraulic diameter       | D h                                                  | 45.7 µ m                      |
| RBC model                |                                                      |                               |
| cell radius              | r rbc                                                | 4 µ m                         |
| viscosity contrast       | GLYPH<21>                                            | 1.0                           |
| feeding haematocrit      | H F                                                  | 1%                            |
| local haematocrit        | Ht                                                   | variable                      |
| simulation setup         |                                                      |                               |
| particle Reynolds number | Re p                                                 | 0.03                          |
| capillary number         | Ca                                                   | 0.6                           |
| LB relaxation parameter  | e GLYPH<28> BGK                                      | 1.0                           |
| voxel size               | GLYPH<1> x                                           | 0.6667 µ m                    |
| time step                | GLYPH<1> t                                           | 7.41 GLYPH<2> 10 GLYPH<0> 8 s |
| shear modulus            | e GLYPH<20> s                                        | 5.00 GLYPH<2> 10 GLYPH<0> 4   |
| bending modulus          | e GLYPH<20> b                                        | 4.50 GLYPH<2> 10 GLYPH<0> 5   |
| reduced bending modulus  | e GLYPH<20> b GLYPH<1> x 2 = ( e GLYPH<20> s r 2 rbc | 1 / 400                       |
| cell-cell interaction    | e GLYPH<20> cc                                       | e GLYPH<20> b                 |
| cell-wall interaction    | e GLYPH<20> cw                                       | 10 e GLYPH<20> b              |

bending elasticity of the RBC membrane. The behaviour of RBCs in flow is determined by two non-dimensional numbers: the particle Reynolds number Re p and the capillary number Ca Re . p represents the ratio of particle inertia to fluid viscous force and is defined based on the channel Reynolds number Re c :

$$R e _ { \mathrm p } = R e _ { \mathrm c } \left ( 2 r _ { \mathrm t b c } / D _ { \mathrm h } \right ) ^ { 2 }, \quad R e _ { \mathrm c } = \rho _ { \mathrm e x } \bar { u } D _ { \mathrm h } / \eta _ { \mathrm e x }$$

where r rbc is the RBC radius and D h is the hydraulic diameter of the rectangular channel; ¯ is the u cross-sectional mean velocity of the unperturbed flow; GLYPH&lt;26&gt; ex and GLYPH&lt;17&gt; ex are the density and dynamic viscosity of the suspending fluid, respectively. Ca represents shear-induced deformation of the RBC membrane and is defined using its shear elasticity (governed by the strain modulus GLYPH&lt;20&gt; s ):

$$C a = \tau _ { r b c } \dot { \gamma } _ { c }, \quad \tau _ { r b c } = \eta _ { e x } r _ { r b c } / \kappa _ { s }, \quad \dot { \gamma } _ { c } = 8 \bar { u } / D _ { h }$$

where GLYPH&lt;28&gt; rbc is a characteristic time scale for the RBC membrane to relax to its equilibrium shape from a transient state and ˙ GLYPH&lt;13&gt; c is the characteristic wall shear rate in equivalent Poiseuille flow. Both Ca and the viscosity contrast GLYPH&lt;21&gt; (between RBC cytosol and the suspending fluid) are kept the same for experiments and simulations to ensure comparability (see values in Table 1).

## 3. Results

The experimental and simulation results on the RBC distribution in the channel are shown in Sec. Counterintuitive RBC distribution . Additional data of RBC dynamics and the mechanisms

Figure 2: RBC distribution across the channel obtained by projecting cells either horizontally or vertically. (a) Experimental PDD profiles (lines) at x = 65, 565, 1220, 2065 µ mversus simulation haematocrit Ht histograms (bars) at x = 150, 450, 750, 1050 µ m, showing the widthwise distribution of RBCs at H F = 1%. Note that the experimental PDDs have been scaled by a constant to be superimposed onto the simulation Ht histograms. (b) (Simulation) Widthwise Ht profiles at H F = 2%. (c) (Simulation) Widthwise Ht profiles at H F = 5%. (d-f) (Simulation) Corresponding depthwise Ht profiles at H F = 1%, 2%, 5%, respectively.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00250v2-figure-7d843f48-8d78-4933-a729-b289fb7cad83.jpg)

leading to the observed RBC distribution is provided in Sec. Flow properties and RBC dynamics .

## 3.1. Counterintuitive RBC distribution

The main results of the study are reported here, focusing on the dynamical distribution of RBCs in a dilute suspension. Experimentally, by continuously feeding RBCs at fixed haematocrit ( H F GLYPH&lt;20&gt; 1%) into the microfluidic device, the transverse motion of cells can be evaluated by monitoring the probability density distribution (PDD, see Sec. Probability density distribution (PDD) of RBCs ) of RBCs at sequential locations along the axial direction. This process is also replicated numerically to corroborate the experimental observations and provide further insights into the 3D motion of cells in a microchannel of rectangular cross-section.

O -centre two-peak (OCTP) profile. GLYPH&lt;11&gt; In experiments with negligible inertia ( Re p &lt; 2 GLYPH&lt;1&gt; 10 GLYPH&lt;0&gt; 4 ), we observe an evident increase of cell concentration close to each wall, manifested by the augmentation of two o GLYPH&lt;11&gt; -centre peaks in the PDD profile depicting RBC fractions across the channel width (Fig. 2a, Q = 0 2 : µ L min). / Locations of the two peaks are approximately symmetric about the channel centreline and keep moving inwards until a certain distance about 0.2 GLYPH&lt;24&gt; 0.3 W = 2 from the wall is reached (i.e., y = (0 5 : W ) = 0 7 : GLYPH&lt;24&gt; 0 8). Beyond this distance, the inward migration of cells : ceases (compare x = 1220 ; 2065 µ m in Fig. 2a); this eventually causes an o GLYPH&lt;11&gt; -centre two-peak (OCTP) ordering in the density profile. Contrarily, in the central region, the cell concentration remains nearly unchanged. By comparing the PDD profiles along the channel length, a cell population originally close to each wall is found to cross streamlines and progressively shift inwards (see x = 65 ; 565 ; 1220 µ m). For an initial cell distribution with significant left-right asymmetry induced by misaligned flow inlet (see Fig. 1d), the same inward motion of cells near the wall occurs (Fig. S6, SI). The observed quasi-steady OCTP distribution of RBCs clearly deviates from

the commonly reported cell distribution (or haematocrit profile) in microscale channels, which typically exhibits a core of cells near the channel centre accompanied by a cell-depleted layer near the channel wall [43, 44]. To rule out the possibility that such an OCTP distribution is exclusive to the volume flow rate Q = 0 2 : µ L min, / Q is also varied up to 4.0 µ L min for the aligned / configuration (Fig. S7, SI). Qualitatively similar evolution of the PDD profile is identified for all cases, featuring the development enhancement of density peaks at two symmetric locations near / the lateral channel walls.

To understand the peculiar cell distribution observed in experiments, numerical simulations for a H F = 1% RBC suspension flowing at low Reynolds number ( Re p = 0.03) are performed in a long straight microchannel of identical hydraulic diameter ( Dh = 45 7 : µ m) and similar aspect ratio ( AR = W H = = 3 2). : In simulations, we have access to the three-dimensional organisation of RBCs within the channel and can observe the RBC distribution from not only the width direction W but also the depth direction H , contrary to common experimental set-ups where only one dimension is observable. To quantify the RBC distribution, we analyse the time-averaged haematocrit ( Ht ) distribution across respective directions of the channel (Fig. 3a). In the W direction, the RBC distribution agrees well with experimental findings of the OCTP pattern and equivalent 'pseudo-equilibrium' position of the peaks are found, i.e., y eq = (0 5 : W ) = 0 7 (Fig. 2a). :

Adistinct initial state of the suspension is found in the H direction upon the entry of cells into the rectangular channel (at x = 150 µ m GLYPH&lt;25&gt; 3 3 : D h ). While the widthwise distribution presents a more or less centralised pattern (Fig. 2a), the depthwise distribution features two primary density peaks (Fig. 2d). Furthermore, a substantial inward march of cells towards the channel centreline can be observed, with marked shifts of both density peaks. Downstream the channel, a pronounced five-layered ordering develops, featuring a cell-depletion band next to each wall, two cell-enrichment bands neighbouring the centre, and one cell-suppression band right at the centre (see x = 1050 µ m GLYPH&lt;25&gt; 23 D h in Fig. 2d). An arguable equilibrium position for the primary density peaks exists at about 0.7 times the half channel depth away from the wall, i.e., z eq = (0 5 : H ) = 0 3. :

Figure 3: (Simulation) Cross-sectional visualization of RBC patterns. (a) Example of the haematocrit ( Ht ) analysis showing superimposed RBCs from multiple simulation snapshots (at designated time intervals) in sequential sampling boxes (blue) measuring 50 µ mlong each at target locations (ROIs) downstream of the channel entrance. Cells within each sampling box are projected and allocated into 10 bins across the width ( W ) and depth ( H ) directions for counting, respectively, with their position determined by centre of mass. (b) Cross-sectional slices showing the accumulative transverse pattern of RBCs at x = 0 D h ; 10 D h ; 23 D h away from the entrance, combining snapshots from 50 consecutive time steps. The 1st, 2nd and 3rd row represent snapshots of the H F = 1% RBC suspension in a straight rectangular channel with aspect ratios AR = 3 2 : ; 1 7 : ; 1 0 (varying : W = 96 ; 50 5 : ; 30 µ mwhile keeping H = 30 µ mfixed), respectively.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00250v2-figure-045bc9c6-6b55-4dd1-86de-797271ff8e71.jpg)

E ect of feeding haematocrit. GLYPH&lt;11&gt; The OCTP pattern can still be observed in simulated RBC suspensions of moderately higher concentration, where the feeding haematocrit H F is increased to 2% (Fig. 2b,e) and 5% (Fig. 2c,f), respectively. However, the magnitude of the near-wall density peaks relative to the central region are in general weakened. Notably, in the depthwise distribution of H F = 5% (Fig. 2f), the two peaks (for which the local haematocrits can reach as high as 8-9 % and presumably violating the dilute limit) are progressively smoothed out downstream the channel and the OCTP pattern vanishes at x = 1050 µ m. Meanwhile, the central band is growingly enriched with cells and the over Ht profile becomes flattened, suggesting considerable RBC dispersion taking place to even out regional density heterogeneities.

E ect of channel aspect ratio. GLYPH&lt;11&gt; The occurrence of OCTP patterns across both dimensions of the channel suggests the phenomenon of RBC focusing. Indeed, by investigating the transverse organisation of RBCs in sequential cross-sections ( x = 0 D h ; 10 D h ; 23 D h ), we observe the gradual formation of an RBC loop reminiscent of the classic 'tubular pinch e GLYPH&lt;11&gt; ect' arising from inertia (Fig. 3b), yet our experiments and simulations are all under negligible inertia as indicated by the low Reynolds number. Reducing the channel aspect ratio from the original AR = 3 2 to : AR = 1 7 : and AR = 1 0 does not remove the focusing behaviour, but brings the RBC loop closer to channel : centre. Examination of the Ht profiles in AR = 1 7 : ; 1 0 channels against the original : AR = 3 2 : channel reveals an inward shift of the OCTP equilibrium position in the width direction from y eq = (0 5 : W ) = 0 7 to : y eq = (0 5 : W ) = 0 5 : ; 0 3 (Fig. S12a,c in SI). In the case of : AR = 1 0 where the : two dimensions are identical, the equilibria from the W and H directions become equivalent, i.e., y eq = (0 5 : W ) = z eq = (0 5 : H ) = 0 3 (Fig. S12c,d in SI). :

## 3.2. Flow properties and RBC dynamics

To reveal the underlying physics behind the intriguing OCTP distribution of RBCs presented in Sec. Counterintuitive RBC distribution , we investigate several key aspects of the problem, including the origin of transverse cell migration, the e GLYPH&lt;11&gt; ect of channel geometry, the distribution of cell velocity, the formation of the cell-depletion layer and the variation of tube haematocrit. The two main questions we aim to address are: 1) How do the two density peaks in the PDD profile come into being? 2) Why do the peaks keep building up instead of being dispersed?

Mechanism for cell migration. Given its small Reynolds number, our RBC flow practically falls into the low-inertia regime, where viscous e GLYPH&lt;11&gt; ects dominate. For deformable objects in viscous flow, the phenomenon of their lateral drift away from the wall is well-documented, dating back to the 1830s when Poiseuille first observed a plasma layer (or cell-free layer, CFL) in frog arterioles and venules [45]. The origin of such motion of RBCs was later determined to be the non-inertial hydrodynamic lift arising from unbalanced pressure forces on the cell. This lift can drive cells in microscale vessels channels to cross streamlines and migrate towards the central area, ultimately / leading to the formation of a two-phase flow featuring a cell core and a fluid periphery [43, 44]. Following the discovery by Goldsmith and Mason [46], the lateral migration of deformed particles in non-inertial flows has been confirmed by other experimental studies [47, 48, 49, 50] and numerous modelling reports [51, 52, 53, 54, 55, 56, 57, 58, 59]. In principle, there are three mechanisms for the cross-streamline motion of RBCs in low-Reynolds-number shear flows: (1) wall-repulsion-induced lift force on RBCs, (2) shear-gradient-induced asymmetric deformation of RBCs, and (3) hydrodynamic di GLYPH&lt;11&gt; usion (or shear-induced di GLYPH&lt;11&gt; usion) between RBCs due to cell-cell interactions. Depending on the flow properties (e.g., velocity profile, shear distribution and wall confinement) and cell conditions (e.g., position, shape and rigidity), the direction and

Figure 4: (Simulation) Statistical* lift velocities of RBCs within x = (50 GLYPH&lt;6&gt; 8) µ m, characterising the RBC migration along the (a) y -direction (widthwise) and (b) z -direction (depthwise), respectively. For statistical analysis in both (a) and (b), the channel is fictitiously folded along the channel centreline and then divided into 10 bins, followed by the assignment of cells according to the position of their center of mass. A positive lift velocity here indicates an inward motion of RBCs towards the channel centreline, vice versa

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00250v2-figure-e11ab2ea-f3e1-401d-8752-a9657e7a770e.jpg)

*The short line and triangle in the boxplots represent the median and mean lift velocity for each bin, respectively. The horizontal caps show the statistics from Welch's t-tests between neighbouring groups, where *P &lt; 0.05 is considered as statistical significance; **P &lt; 0.01, ****P &lt; 0.0001.

strength of RBC motion under the combined e GLYPH&lt;11&gt; ect of the above mechanisms can be di GLYPH&lt;14&gt; cult to determine.

Decay of lift velocity. The presence of RBCs at H F = 1% (or below) does not significantly disturb the apparent velocity profiles from their unperturbed states (compare Fig. S15c,d and Fig. S16 of the SI), we can therefore analyse the transverse motion of RBCs with respect to the unperturbed flow. Subject to the complex flow pattern in the channel as elaborated by Fig. S15, the motion of RBCs within are likely to be a GLYPH&lt;11&gt; ected by both the solid wall and the fluid shear gradients. At low Reynolds number, this motion can be conveniently decomposed into two orthogonal directions, hereby the W and H directions, on assuming negligible inertia and invoking the Stokes equations. Fig. 4a-b show the statistical RBC lift velocities V l (relative to the unperturbed flow mean velocity ¯) in 10 subdivisions of a cross-section, divided either widthwise u ( W ) or depthwise ( H ). Compared with the depthwise lift, the widthwise lift is found to decay much faster as the cell-wall distance increases. While V z l ; vanishes at 0.75 times of H = 2 away from the wall (i.e., z = (0 5 : H ) = 0 25), : V y l ; already becomes negligible at 0.35 times of W = 2 away from the wall (i.e., y = (0 5 : W ) = 0 65). The contrast of lift decay between the two directions can : be explained by the shear-rate and shear-rate-gradient profiles in the cross-section (Fig. S15e-f, SI). As the cells migrate towards the channel centreline, the local shear rate diminishes to zero, with much faster drop in W than in H . Meanwhile, the local shear gradient also quickly vanishes in W direction, whereas that in H direction decreases in a relatively mild way. Consequently, contributions of both the wall-induced lift and shear-gradient-induced lift in the W direction are weakened more quickly.

The distribution of lift velocities in both the width and depth directions coincides with the overall evolution of the Ht profiles presented in Fig. 2a,d: whereas a 'pseudo-equilibrium' position close to the wall exists for the widthwise migration (beyond which the lateral motion of cells ceases), the depthwise cell migration is more substantial and results in a 'pseudo-equilibrium' position much closer to the channel centre. Therefore, a correlation between the decaying lift of cells and their spatial distribution in the cross-sectional direction can be inferred.

Growth of cell-free layer. We now move on to investigate how the transverse motion of RBCs contributes to the spatial development of CFL along the channel axis. In the lower-flow-rate experiment Q = 0 2 : µ L min ( / H F GLYPH&lt;20&gt; 1%), the growth of CFL thickness (recording the averaged CFL from the two lateral walls) over 46 D h manifests an evident power-law behaviour with a best-fit exponent of 0.39 ( R 2 = 0 93, Fig. 5). : For the other experiment with higher flow rate

Figure 5: (Empty circles) Experimental CFL measured in the width direction of the channel (averaged between opposite walls) under Q = 0 2 : µ L min and / Q = 0 8 : µ L min. (Solid symbols) Numerical CFLs calculated along each / wall of the channel (Right Bottom Left Top-walls). Both the experimental and numerical CFL values / / / GLYPH&lt;14&gt; cfl have been normalised by the maximum CFL thickness GLYPH&lt;14&gt; max detected within the whole investigated range, and the development lengths L are normalised by the channel hydraulic diameter D h. The main frame and the inset of the figure contain the same experimental data, but di GLYPH&lt;11&gt; erent numerical data from two simulations: the former with normal human RBCs where Ca = 0 6 and the latter with ageing human RBCs where : Ca = 0 1. :

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00250v2-figure-5cff5dc6-d1ee-4be1-b8ea-17381199f7f3.jpg)

Q = 0 8 : µ L min (hence larger / Re ), nearly identical power-law behaviour is found for the CFL development. Similarly, the numerical CFLs over 28 D h (recording separate CFLs for individual lateral vertical walls in simulation) also present power-law trends, but with certain discrepancy in / the fitting exponents when compared to the experimental value of 0.39: for the widthwise CFLs at the right and left walls, the values are 0.26 ( R 2 = 0 91) and 0.29 ( : R 2 = 0 94), respectively; for : the depthwise CFLs at the bottom and top walls, they are 0.29 ( R 2 = 0 91) and 0.25 ( : R 2 = 0 95), : respectively (Fig. 5). Part of the discrepancy in CFL growth rates between the experiment and simulation is attributed to the usage of di GLYPH&lt;11&gt; erent RBC types, namely horse RBC and human RBC ( Ca = 0 6 here), the latter of which has higher deformability in nature when exposed to : shear. By tuning down the shear elasticity of the modelled cells in simulation ( Ca = 0 1), which : imitates ageing human RBCs of increased sti GLYPH&lt;11&gt; ness, the discrepancy between the experimental and numerical data is significantly reduced with numerical exponents of 0.29-0.34 achieved (inset of Fig. 5).

Despite the di GLYPH&lt;11&gt; erence between the growth rates of experimental and simulation CFL, a neutral finding is that no critical saturation of the CFL increase can be identified throughout the investigated range (see raw data of the CFLs in Fig. S9b and Fig. S13a of the SI). This means that the development of CFL in a typical microfluidic channel (rectangular cross-section) for a dilute suspension can take much longer than earlier estimates provided for microvessels or cylindrical tubes [60, 61, 62, 63]. In particular, Katanov et al. [61] reported a universal length L c GLYPH&lt;20&gt; 25 D (independent of flow rate and haematocrit) for full CFL development in vessels between D = 10 GLYPH&lt;24&gt; 100 µ m, whereas we find L c &gt; 28 D h in simulation (only serving as a lower-bound estimate because we have limited our flow domain to such a length for tractable computational

Table 2: (Simulation) Variation of tube haematocrit along the channel axis.

| Inlet type   | H F   | Ht (150 µ m)   | Ht (450 µ m)   | Ht (750 µ m)   | Ht (1050 µ m)   | Ht (1150 µ m)   |
|--------------|-------|----------------|----------------|----------------|-----------------|-----------------|
| CTRAC-L      | 1%    | 0.810%         | 0.749%         | 0.703%         | 0.696%          | 0.685%          |
| CTRAC-S      | 1%    | 0.705%         | 0.661%         | 0.659%         | 0.622%          | 0.619%          |

cost) and L c &gt; 46 D h in experiment (hydraulic diameter of the channel D h GLYPH&lt;25&gt; 45 7 : µ m). Based on our numerical findings of the 'pseudo-equilibrium' positions for the widthwise cell migration (until 0.3 times of half channel width from wall) and depthwise cell migration (until 0.7 times of half channel depth from wall), the CFL thickness in the simulation may continue to grow until reaching 10 GLYPH&lt;24&gt; 14 µ m, which are considerably larger than the values of 6 GLYPH&lt;24&gt; 8 µ mwehave monitored within L c GLYPH&lt;20&gt; 28 D h (Fig. S13a, SI).

Reduction of tube haematocrit. The long-range development of CFL will inevitably a GLYPH&lt;11&gt; ect the local rheological properties of the suspension. As shown in Table 2, the tube haematocrit Ht measured from either of our simulated geometries (CTRAC-L: large-contraction inlet, all simulation data and plots presented elsewhere in the paper are from this geometry; CTRAC-S: smallcontraction inlet. See more details of CTRAC-L and CTRAC-S in Sec. Model configuration ) drops continuously along the channel axis compared to the feeding haematocrit H F, even after 1050 µ m ( GLYPH&lt;25&gt; 23 Dh ) away from the entrance. This reduced volume fraction of blood flow in microchannels versus the discharge haematocrit ( H D GLYPH&lt;25&gt; H F) is known as the Fåhræus e GLYPH&lt;11&gt; ect, where Ht H = D = ¯ U B = ¯ V x &lt; 1 ( ¯ U B and ¯ V x are average velocities of the overall blood flow and the cells, respectively) [64, 65]. It stems from the cross-streamline motion of RBCs towards the highvelocity region around channel centreline, which empowers the cells within the channel to travel faster on average than the ambient fluid in the streamwise direction. Owing to the fairly long CFL development length ( &gt; 28 D h in simulation), the cross-sectional distribution of RBC fluxes and consequently the average ratio of cell-to-fluid velocities (Fig. S17 and Fig. S18 in Sec. RBC lateral migration and axial acceleration ) continues to change till the very end of the channel, leaving consistent measurement of tube haematocrits within the microchannel di GLYPH&lt;14&gt; cult.

## 4. Discussion

In the following we will discuss and compare our findings in Sec. Results with existing studies and interpretations, and elaborate on their implications for experimental studies of dilute blood flow.

## 4.1. Comparison with existing studies

Hydrodynamic lift. The mechanism of hydrodynamic lift for tank-treading particles in shear flows has been analysed theoretically by Olla and Seifert in the 1990s [66, 67]. A simplified scaling law V l GLYPH&lt;24&gt; ˙ GLYPH&lt;13&gt;= h 2 was derived for the lift velocity V l under constant shear, as a function of the shear rate ˙ GLYPH&lt;13&gt; and cell-wall distance h [68]. This quadratic scaling has been widely employed in the literature as an approximation to estimate the lateral migration of particles in various flow conditions. However, recent studies show that in pressure-driven channel flows where non-constant shear (or non-zero shear gradient) exists, e.g., Poiseuille flow, the curvature of the velocity profile modifies the scaling law in a non-trivial way and would significantly slow down the spatial decay of the lift velocity [54, 69, 59, 70]. Indeed, through systematic measurement of

Figure 6: Theoretical prediction of RBC lift velocities (line curves and contour plots) given by V l GLYPH&lt;24&gt; ˙ GLYPH&lt;13&gt;= h across (a) the width direction and (b) the depth direction of the channel, respectively. z and y values in the legends represent positions relative to the mid-plane in each direction (see Fig. S15c-d in the SI), with z = 0 and y = 0 denoting the widthwise and depthwise mid-planes, respectively. The black dots represent rescaled values of the numerical lift velocities extracted from Fig. 4a-b. The black vertical lines in (a-b) indicate the center of the outmost layer of cells at x = 50 µ mdetected in the RBC simulation.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00250v2-figure-2b257bac-6ebc-4c6b-92f3-5fd5cab10b35.jpg)

the lateral migration of various vesicles in bounded planar Poiseuille flow under di GLYPH&lt;11&gt; erent confinement degrees, Coupier et al. revealed that the scaling V l GLYPH&lt;24&gt; ˙ GLYPH&lt;13&gt;= h 2 would fail to predict the trajectory of the vesicle, while a markedly di GLYPH&lt;11&gt; erent empirical law V l GLYPH&lt;24&gt; ˙ GLYPH&lt;13&gt;= h held [50]. Losserand et al. confirmed the scaling V l GLYPH&lt;24&gt; ˙ GLYPH&lt;13&gt;= h using RBC experiments and further refined the empirical law [71].

To examine whether the lateral migration of RBCs in our simulation obeys the lift scaling that Coupier et al. [50] and Losserand et al. [71] proposed for Poiseuille flow, we adopt theoretical shear rates from the asymptotic solution of pressure-driven flow in rectangular channel [72] and calculate the lift velocities applying V l GLYPH&lt;24&gt; ˙ GLYPH&lt;13&gt;= h (abandoning all prefactors for generality). Further with weighted normalisation, we obtain the dimensionless theoretical lifts ( ˆ V y l ; and ˆ V z l ; ) for both W and H directions (Fig. 6a-b). For comparison with the theoretical prediction, numerical medians of our numerical lift velocities are extracted from Fig. 4a-b. After rescaling and overlapping them onto the theoretical curves, we observe satisfactory data agreement across the H direction but only partial agreement in the W direction for the near-wall region only. Approaching the central region in the W direction, the theoretical prediction breaks down. This is actually expected since the velocity profile there is abruptly blunted and significantly deviates from a parabolic shape. We have also examined against our simulation data the theoretical prediction given by V l GLYPH&lt;24&gt; ˙ GLYPH&lt;13&gt;= h 2 , which would overestimate the spatial decay of the lift in general and failed to yield an agreement with the simulation data (results not shown here).

Cell-free layer. CFL characterisation is an active research field owing to its prominent role in regulating blood viscosity and local haemodynamics. Limited by the complexity of in vivo imaging and measurement [73], glass capillaries and PDMS microchannels have been widely used to mimic simplified micro-circulatory environment for in vitro observation of CFLs under di GLYPH&lt;11&gt; erent rheological conditions [74]. Meanwhile, versatile computational tools are developed and plenty of in silico studies on CFL emerge [44]. However, most in vitro in silico studies assume / a steady CFL for a given channel segment and hence fail to capture the highly dynamical nature of realistic CFLs as in the microcirculation [75, 76]; albeit some authors made improvements

Figure 7: Regression analysis (log-log scale) of the (a) experimental and (b) numerical CFLs as in Fig. 5. The CFL values GLYPH&lt;14&gt; cfl are normalised by the maximum CFL thickness GLYPH&lt;14&gt; max detected within the whole investigated range, and the development lengths L are normalised by the channel hydraulic diameter D h. The inset of (b) shows a zoom-in view of the numerical CFL trends within L = 10 GLYPH&lt;24&gt; 28 Dh .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00250v2-figure-c25ff67c-5d72-48fc-a99e-d9cd6bf998f2.jpg)

by proposing a nominal CFL development length by multiplying the mean flow velocity and characteristic flow time [15, 61]. It is only recently that exact analyses of the spatio-temporal heterogeneity of CFL were presented [77, 62, 78], detailing the spatial development recovery of / the CFL alongside microchannels in the presence of upstream disturbances.

In the present study, having confirmed the spatial variation of lift velocity V l GLYPH&lt;24&gt; ˙ GLYPH&lt;13&gt;= h (rather than V l GLYPH&lt;24&gt; ˙ GLYPH&lt;13&gt;= h 2 ) against the cell-wall distance h for individual RBCs, estimation can be made for the longitudinal development of GLYPH&lt;14&gt; cfl in a dilute suspension using an analytical approach (see Sec. Scaling analysis for further details). Our derivation shows that considering a = 1 (corresponding to V l GLYPH&lt;24&gt; ˙ GLYPH&lt;13&gt;= h ) in Eq. S1 yields a power-law scaling for the CFL as expressed in Eq. S11: GLYPH&lt;14&gt; cfl GLYPH&lt;24&gt; ' GLYPH&lt;12&gt; ; GLYPH&lt;12&gt; = 1 3, where = ' is the axial distance along the channel. The exponent of 1 3 here is / surprising as it coincides with the power-law scaling for vesicle RBC motion in simple shear flow / with constant shear rate (see Eq. 3 in reference [49]). Note that because GLYPH&lt;14&gt; cfl GLYPH&lt;24&gt; ' 1 3 = is derived on the basis of a parabolic velocity profile under the assumption of planar Poiseuille flow (constant shear-rate gradient), whereas the flow patterns in a real microchannel are far more complex (see Fig. S15, SI), its accuracy for the prediction of spatial CFL development remains to be verified by our experimental measurements and numerical results.

The power-law growth GLYPH&lt;14&gt; cfl GLYPH&lt;24&gt; ' 1 3 = predicted by the scaling analysis assuming ideal Poiseuille flow turns out to agree well with both our experimental and numerical CFLs at H F = 1%, for which the best-fitted power-law exponents were earlier shown to be 0.39 and 0.25 GLYPH&lt;24&gt; 0.29, respectively (Fig. 5). In further regression analyses of the experimental and numerical CFLs, both data sets are found to collapse on GLYPH&lt;14&gt; cfl GLYPH&lt;24&gt; ' 1 3 = (Fig. 7a-b). Furthermore, the instantaneous slopes for the numerical data are nearly independent of ' at the early stage of CFL development, presenting a roughly steady increase of 1 3 (see Fig. 7b). / This implies the predominant role of a single mechanism for the transverse motion of cells, which is sensibly the hydrodynamic lift as our scaling analysis suggests. Moving downstream of the channel, a continual increase of the CFL is observed upon approaching ' = 28 D h , but with slower growth indicated by decreasing slopes noticeably smaller than the original 1 3 (see inset of Fig. 7b). / This deceleration in CFL growth is expected down to the fact that the lift forces are weakened as the cells migrate further away from the wall and enter the low-shear region with blunted velocity profiles, where the lateral motion of cells slows down. In simulations of higher feeding haematocrit H F = 2% and H F = 5%, the scaling of GLYPH&lt;14&gt; cfl GLYPH&lt;24&gt; ' 1 3 = still holds for early-stage development of the CFL until reaching far downstream of the channel (see Fig. S14e and Fig. S14f). Notably for H F = 5%, the slope approximates zero between ' = 20 GLYPH&lt;0&gt; 25 D h and suggests a saturated CFL. This is in line with what was found in the evolution of Ht profiles (Fig. 2c,f), where exceedingly amplified local haematocrits (about 8-9%) break the dilute regime and lead to considerable RBC dispersion that resists further drift of cells away from the wall.

From previous simulations of much denser RBC suspensions ( Ht = 15% GLYPH&lt;24&gt; 45%) in cylindrical tubes [61], a similar power-law behaviour was reported for the steep increase of CFL thickness at an initial development stage, before a second time scale induced by relaxation of the RBC core came into play. However, the power-law exponents under several shear rates were considerably smaller than 1 3 even at the very beginning. / The deviation from 1 3 was caused by the cell-/ cell interactions existing in their dense suspensions, which obstructed the initial migration of RBCs and led to slower CFL growth than in a dilute suspension like ours. Arguably, the authors reasoned for the ideal growth factor of 1 3 by invoking the lift force scaling / F l GLYPH&lt;24&gt; V l GLYPH&lt;24&gt; 1 = h 2 (where ˙ is a constant and does not a GLYPH&lt;13&gt; GLYPH&lt;11&gt; ect the scaling) derived from simple shear flows [52, 56]. Albeit such reasoning does explain the qualitative trend of their CFL growth, the simple shear flow approximation itself for quadratic flows in a microchannel (where ˙ varies) is questionable, GLYPH&lt;13&gt; which may fail to capture important features of the RBC distribution in the presence of shear rate gradients as recent theories of particle migration at low Reynolds number demonstrate [79, 80].

## 4.2. Implications and limitations

What we have covered so far not only elucidates the 'aberrant' RBC distribution observed in our experiment and simulation contradicting earlier findings, but also has wider implications for experimentalists working on dilute suspensions, especially in high-aspect-ratio microfluidic devices. In brief, microfluidic designs need to be longer if their purpose is to investigate the microscopic behaviour of a dilute suspension after completed lateral migration. Experimentally, it is notoriously di GLYPH&lt;14&gt; cult to accurately control the inlet distribution of a cell particle suspension / under various entrance e GLYPH&lt;11&gt; ects . Therefore, local observations are frequently reported based on the assumption that the suspension organisation is free of transient e GLYPH&lt;11&gt; ects at the chosen region of interest (ROI) as long as the ROI is adequately far away from the entrance. However, the required length as revealed by our study turns out more demanding than commonly believed [61, 60, 62, 63]. Some key messages are:

(I) For experiments of dilute RBC suspensions in high-aspect-ratio channels, e.g., AR = W H = &gt; 3, the hydrodynamic lift of RBCs are substantially suppressed due to the existence of low-shear-rate and low-shear-gradient zones in the centre. Consequently, the full development of a CFL in these channels may require a much longer length ( L c &gt; 28 D h in simulations and L c &gt; 46 D h in experiments) than what is normally expected for cylindrical tubes ( L c GLYPH&lt;20&gt; 25 D according to reference [61]). To achieve full CFL development within L c = 20 GLYPH&lt;0&gt; 25 D h in rectangular channels, RBC suspensions with H F GLYPH&lt;21&gt; 5% should be adopted. Furthermore, the relevance of in vitro measurements of cellular behaviour using such high-aspect-ratio microchannels to the realistic microcirculatory blood flow needs reappraisal. After all, the vessels in the microcirculation are more likely to have circular cross-sections and considerably di GLYPH&lt;11&gt; erent velocity profiles, which the radial distribution of RBCs relies on.

(II) When determining the organisation of a dilute suspension (e.g., haematocrit profile) in a typical microchannel of rectangular cross-section, local measurements may be misleading if the channel is short or only moderately long. This is because heterogeneities emerging in the suspension (e.g., density peaks) will persist and hinder the formation of an equilibrium cell distribution owing to the lack of RBC dispersion under weak cell-cell interactions . In particular, if certain upstream disturbance exists (e.g., geometric constriction or expansion), history e GLYPH&lt;11&gt; ects need to be considered when recording local suspension properties such as the tube haematocrit. Comparing the cases CTRAC-S and CTRAC-L in Table 2, with channel geometry, unperturbed flow rate and feeding haematocrit all being the same, yet the tube haematocrits measured at sequential loca-

tions constantly di GLYPH&lt;11&gt; er by more than 6% (relative change). The reason is simply their di GLYPH&lt;11&gt; erently configured flow inlets, causing unequal degrees of disturbance to the initial distribution of cells.

(III) Phenomena violating empirical laws may occur downstream of bifurcations if the upstream RBC suspension is at an intermediate stage of development and has special cell distribution, e.g., the OCTP profile observed in our study. Recently, counterintuitive inversion of haematocrit partition against the classic Zweifach-Fung e GLYPH&lt;11&gt; ect in a bifurcating channel was reported by Shen et al. [31], arising from the formation of consecutive layers of high and low RBC haematocrits in the parent branch. The juxtaposed configuration of RBCs discovered by that study is in line with the five-layered ordering of cells we find in the depth direction of our RBC flow (Fig. 2d-f). Therefore, the phenomenon Shen et al. captured might also have been a consequence of inadequate channel length for the suspension to be fully developed in the parent branch.

Some limitations of the present study are: (1) While the numerical model simulates human RBCs, horse RBCs are used in experiments due to laboratory restrictions. This causes a certain degree of discrepancy between the simulation and experiment, but its influence on the suspensions behaviour we studied, such as cell focusing and CFL development, is shown to be quantitative instead of qualitative. (2) The e GLYPH&lt;11&gt; ect of inlet configuration on suspension organisation has not been characterised in a quantitative manner, which can potentially provide more specific guidance for experimentalists to improve their microfluidic designs by reducing transient e GLYPH&lt;11&gt; ects. We aim to address this issue in future studies, involving geometric factors such as constrictions, expansions and bifurcations. (3) The experimental data is insu GLYPH&lt;14&gt; cient for analysis of instantaneous growth rates of the CFL thickness as has been performed for the simulation. More measurements using a denser arrangement of ROIs, though quite challenging, should be made in future experiments in order to enable precise comparison with simulations. (4) To make the real-scale (channel length &gt; 1 mm) simulations tractable, we have largely reduced the parameter space and refrained from an exhaustive investigation combining a range of flow scenarios and cell conditions. A parametric study integrating all these e GLYPH&lt;11&gt; ects will lead to enhanced understanding of the intricate behaviour of dilute suspensions processed in microfluidic devices.

## 5. Conclusion

We performed both microfluidic experiments and numerical simulations of dilute blood suspensions in a low-Reynolds-number channel flow ( Re p GLYPH&lt;28&gt; 1). Surprisingly, an o GLYPH&lt;11&gt; -centre twopeak (OCTP) ordering of RBCs was found, which is reminiscent of the 'tubular pinch e GLYPH&lt;11&gt; ect' typical of the radial distribution of particles in inertial microfluidics [81, 82]. However, the transverse motion of cells in our case has an entirely di GLYPH&lt;11&gt; erent origin: the non-inertial hydrodynamic lift of deformable particles [46].

Back to the two primary questions we set out to address (see the opening paragraph of Sec. Flow properties and RBC dynamics ): how do the two density peaks in the PDD profile come into being and why do the peaks keep building up instead of being dispersed? The reason behind the persistence of significant density peaks is the deficiency of hydrodynamic interactions among cells in the dilute limit, where shear-induced di GLYPH&lt;11&gt; usion remains weak and falls short to smooth out density heterogeneities brought by the lateral migration of cells. Consequently, the evolution of the density profiles is predominantly determined by the decay of hydrodynamic lift within the suspension, which exhibits distinct patterns in the width and depth directions of a high-aspect-ratio microchannel and contributes to a counterintuitive (OCTP) profile of cells in the larger dimension. Additionally, depending on the inflow configuration, the initial distribution

of cells upon their entry into the channel varies substantially, which brings extra complexity to the restoration of a converged density profile and the development of an equilibrium cell free layer.

Our findings highlight the importance of local and transient characteristics at the dilute limit for an RBC-laden flow in microchannels. Experimentalists should therefore be cautious when working with dilute suspensions in microfluidics and make judicious assumptions about suspension properties with the presence of upstream disturbance to the flow, especially when the microscopic behaviour of the suspension is the focus of research.

Author contributions. QZ and JF designed the research, analysed the data and wrote the article. QZ performed the simulations. JF performed the experiments. LC assisted in validation of the computational model. MOB and PRH discussed the results and revised the manuscript. MOB supervised the study. MSNO and TK conceived the research idea, supervised the study and edited the manuscript.

Acknowledgements. QZ thanks the University of Edinburgh for the award of a Principal's Career Development Scholarship and an Edinburgh Global Research Scholarship. TK's and MOB's contributions have been funded through two Chancellor's Fellowships at the University of Edinburgh. MOB is supported by grants from EPSRC (EP R029598 1, EP R021600 1), Fondation / / / / Leducq (17 CVD 03), and the European Union's Horizon 2020 research and innovation programme under grant agreement No 801423. Supercomputing time on the ARCHER UK National Supercomputing Service (http: // www.archer.ac.uk) was provided by the 'UK Consortium on Mesoscale Engineering Sciences (UKCOMES)' under the EPSRC Grant No. EP R029598 1. / / The authors declare that they have no competing financial interests.

## References

- [1] S. Chien, S. Usami, R. J. Dellenback, M. I. Gregersen, L. B. Nanninga, M. M. Guest, Blood Viscosity: Influence of Erythrocyte Aggregation, Science 157 (1967) 829-831.
- [2] S. Chien, R. G. King, R. Skalak, S. Usami, A. L. Copley, Viscoelastic properties of human blood and red cell suspensions, Biorheology 12 (1975) 341-346.
- [3] R. Skalak, S. Chien, Rheology of blood cells as soft tissues, Biorheology 19 (1982) 453-461.
- [4] P. M. Vlahovska, D. Barthes-Biesel, C. Misbah, Flow dynamics of red blood cells and their biomimetic counterparts, Comptes Rendus Physique 14 (2013) 451-458.
- [5] D. Abreu, M. Levant, V. Steinberg, U. Seifert, Fluid vesicles in flow, Advances in Colloid and Interface Science 208 (2014) 129-141.
- [6] J. Freund, Numerical Simulation of Flowing Blood Cells, Annual Review of Fluid Mechanics 46 (2014) 67-95.
- [7] T. M. Geislinger, T. Franke, Hydrodynamic lift of vesicles and red blood cells in flow - from Fåhræus &amp; Lindqvist to microfluidic cell sorting, Advances in Colloid and Interface Science 208 (2014) 161-176.
- [8] A. Yazdani, X. Li, G. E. Karniadakis, Dynamic and rheological properties of soft biological cell suspensions, Rheologica Acta 55 (2016) 433-449.
- [9] G. Gompper, D. A. Fedosov, Modeling microcirculatory blood flow: current state and future perspectives, Wiley Interdisciplinary Reviews: Systems Biology and Medicine (2015) n a-n a. / /

[10] G. Danker, C. Misbah, Rheology of a Dilute Suspension of Vesicles, Physical Review Letters 98 (2007).

- [11] T. Krüger, M. Gross, D. Raabe, F. Varnik, Crossover from tumbling to tank-treading-like motion in dense simulated suspensions of red blood cells, Soft Matter 9 (2013) 9008-9015.

[12] D. Cordasco, P. Bagchi, Orbital drift of capsules and red blood cells in shear flow, Physics of Fluids 25 (2013) 091902.

[13] A. Guckenberger, A. Kihm, T. John, C. Wagner, S. Gekle, Numerical-experimental observation of shape bistability of red blood cells flowing in a microchannel, Soft Matter (2018).

- [14] P. Olla, The behavior of closed inextensible membranes in linear and quadratic shear flows, Physica A: Statistical Mechanics and its Applications 278 (2000) 87-106.

| [15]      | J. Zhang, E GLYPH<11> ect of Suspending Viscosity on Red Blood Cell Dynamics and Blood Flows in Microvessels, Micro- circulation 18 (2011) 562-573.                                                                                                                                                                                                                                                            |
|-----------|----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
| [16]      | G. Ghigliotti, A. Rahimian, G. Biros, C. Misbah, Vesicle Migration and Spatial Organization Driven by Flow Line Curvature, Physical Review Letters 106 (2011).                                                                                                                                                                                                                                                 |
| [17]      | B. Kaoui, J. Harting, C. Misbah, Two-dimensional vesicle dynamics under shear flow: E GLYPH<11> ect of confinement, Physical Review E 83 (2011) 066319.                                                                                                                                                                                                                                                        |
| [18]      | O. Aouane, A. Farutin, M. Thiébaud, A. Benyoussef, C. Wagner, C. Misbah, Hydrodynamic pairing of soft particles in a confined flow, Physical Review Fluids 2 (2017).                                                                                                                                                                                                                                           |
| [19]      | M. Toner, D. Irimia, Blood-on-a-Chip, Annual Review of Biomedical Engineering 7 (2005) 77-103.                                                                                                                                                                                                                                                                                                                 |
| [20] [21] | G. Tomaiuolo, Biomechanical properties of red blood cells in health and disease towards microfluidics, Biomi- crofluidics 8 (2014) 051501. B. Sebastian, P. S. Dittrich, Microfluidics to Mimic Blood Flow in Health and Disease, Annual Review of Fluid                                                                                                                                                       |
| [22]      | Mechanics 50 (2018) 483-504. R. Lima, T. Ishikawa, Y. Imai, M. Takeda, S. Wada, T. Yamaguchi, Radial dispersion of red blood cells in blood flowing through glass capillaries: The role of hematocrit and geometry, Journal of Biomechanics 41 (2008) 2188- 2196.                                                                                                                                              |
| [23]      | T. Yaginuma, M. S. N. Oliveira, R. Lima, T. Ishikawa, T. Yamaguchi, Human red blood cell behavior under homogeneous extensional flow in a hyperbolic-shaped microchannel, Biomicrofluidics 7 (2013) 054110.                                                                                                                                                                                                    |
| [24]      | S. Roman, A. Merlo, P. Duru, F. Risso, S. Lorthois, Going beyond 20 um-sized channels for studying red blood cell phase separation in microfluidic bifurcations, Biomicrofluidics 10 (2016) 034103.                                                                                                                                                                                                            |
| [25]      | S. Yang, A. Ündar, J. D. Zahn, A microfluidic device for continuous, real time blood plasma separation, Lab on a Chip 6 (2006) 871-880.                                                                                                                                                                                                                                                                        |
| [26]      | V. Kantsler, E. Segre, V. Steinberg, Dynamics of interacting vesicles and rheology of vesicle suspension in shear flow, EPL (Europhysics Letters) 82 (2008) 58005.                                                                                                                                                                                                                                             |
| [27]      | X. Grandchamp, G. Coupier, A. Srivastav, C. Minetti, T. Podgorski, Lift and Down-Gradient Shear-Induced Di GLYPH<11> u- sion in Red Blood Cell Suspensions, Physical Review Letters 110 (2013).                                                                                                                                                                                                                |
| [28]      | L. Lanotte, J. Mauer, S. Mendez, D. A. Fedosov, J.-M. Fromental, V. Claveria, F. Nicoud, G. Gompper, M. Abkar- ian, Red cells' dynamic morphologies govern blood shear thinning under microcirculatory flow conditions, Pro- ceedings of the National Academy of Sciences (2016) 201608074.                                                                                                                    |
| [29]      | C.-H. Chuang, K. Kikuchi, H. Ueno, K. Numayama-Tsuruta, T. Yamaguchi, T. Ishikawa, Collective spreading of red blood cells flowing in a microchannel, Journal of Biomechanics 69 (2018) 64-69.                                                                                                                                                                                                                 |
| [30]      | J. M. Sherwood, D. Holmes, E. Kaliviotis, S. Balabani, Spatial Distributions of Red Blood Cells Significantly Alter Local Haemodynamics, PLOS ONE 9 (2014) e100473.                                                                                                                                                                                                                                            |
| [31]      | Z. Shen, G. Coupier, B. Kaoui, B. Polack, J. Harting, C. Misbah, T. Podgorski, Inversion of hematocrit partition at microfluidic bifurcations, Microvascular Research 105 (2016) 40-46.                                                                                                                                                                                                                        |
| [32] [33] | C. Iss, D. Midou, A. Moreau, D. Held, A. Charrier, S. Mendez, A. Viallat, E. Helfer, Self-organization of red blood cell suspensions under confined 2D flows, Soft Matter 15 (2019) 2971-2980. A. Abay, S. M. Recktenwald, T. John, L. Kaestner, C. Wagner, Cross-sectional focusing of red blood cells in a                                                                                                   |
| [34]      | constricted microfluidic channel, Soft Matter 16 (2020) 534-543. T. Krüger, Computer Simulation Study of Collective Phenomena in Dense Suspensions of Red Blood Cells under Shear, Vieweg + Teubner Verlag, 2012.                                                                                                                                                                                              |
| [35]      | Y. H. Qian, D. D'Humières, P. Lallemand, Lattice BGK Models for Navier-Stokes Equation, EPL (Europhysics Letters) 17 (1992) 479.                                                                                                                                                                                                                                                                               |
|           | P. L. Bhatnagar, E. P. Gross, M. Krook, A Model for Collision Processes in Gases. I. Small Amplitude Processes                                                                                                                                                                                                                                                                                                 |
| [36]      | in Charged and Neutral One-Component Systems, Physical Review 94 (1954) 511-525.                                                                                                                                                                                                                                                                                                                               |
| [38]      | Review E 65 (2002) 046308. M. Bouzidi, M. Firdaouss, P. Lallemand, Momentum transfer of a Boltzmann-lattice fluid with boundaries, Physics                                                                                                                                                                                                                                                                     |
| [39]      | of Fluids 13 (2001) 3452-3459. A. J. C. Ladd, Numerical simulations of particulate suspensions via a discretized Boltzmann equation. Part 1. Theoretical foundation, Journal of Fluid Mechanics 271 (1994) 285-309. C. S. Peskin, The immersed boundary method, Acta Numerica 11 (2002) 479-517. M. O. Bernabeu, M. L. Jones, J. H. Nielsen, T. Krüger, R. W. Nash, D. Groen, S. Schmieschek, J. Hetherington, |
| [40]      |                                                                                                                                                                                                                                                                                                                                                                                                                |
| [41]      | H. Gerhardt, C. A. Franco, P. V. Coveney, Computer simulations reveal complex distribution of haemodynamic forces in a mouse retina model of angiogenesis, Journal of The Royal Society Interface 11 (2014) 20140543.                                                                                                                                                                                          |
|           | T. Krüger, F. Varnik, D. Raabe, E GLYPH<14> cient and accurate simulations of deformable particles immersed in a using a combined immersed boundary lattice Boltzmann finite element method, Computers & Mathematics Applications 61 (2011) 3485-3505.                                                                                                                                                         |
| [42] [43] | fluid with G. R. Cokelet, H. L. Goldsmith, Decreased hydrodynamic resistance in the two-phase flow of blood through small                                                                                                                                                                                                                                                                                      |

vertical tubes at low flow rates, Circulation Research 68 (1991) 1-17.

- [44] D. A. Fedosov, B. Caswell, A. S. Popel, G. E. Karniadakis, Blood Flow and Cell-Free Layer in Microvessels, Microcirculation 17 (2010) 615-628.
- [45] S. P. Sutera, R. Skalak, The History of Poiseuille's Law, Annual Review of Fluid Mechanics 25 (1993) 1-20.
- [46] H. L. Goldsmith, S. G. Mason, Axial Migration of Particles in Poiseuille Flow, Nature 190 (1961) 1095-1096.
- [47] B. Lorz, R. Simson, J. Nardi, E. Sackmann, Weakly adhering vesicles in shear flow: Tanktreading and anomalous lift force, EPL (Europhysics Letters) 51 (2000) 468.
- [48] M. Abkarian, C. Lartigue, A. Viallat, Tank Treading and Unbinding of Deformable Vesicles in Shear Flow: Determination of the Lift Force, Physical Review Letters 88 (2002).
- [49] N. Callens, C. Minetti, G. Coupier, M.-A. Mader, F. Dubois, C. Misbah, T. Podgorski, Hydrodynamic lift of vesicles under shear flow in microgravity, EPL 83 (2008) 24002.
- [50] G. Coupier, B. Kaoui, T. Podgorski, C. Misbah, Noninertial lateral migration of vesicles in bounded Poiseuille flow, Physics of Fluids 20 (2008) 111702.
- [51] I. Cantat, C. Misbah, Lift Force and Dynamical Unbinding of Adhering Vesicles under Shear Flow, Physical Review Letters 83 (1999) 880-883.
- [52] S. Sukumaran, U. Seifert, Influence of shear flow on vesicles near a wall: A numerical study, Physical Review E 64 (2001) 011916.
- [53] T. W. Secomb, B. Styp-Rekowska, A. R. Pries, Two-Dimensional Simulation of Red Blood Cell Deformation and Lateral Migration in Microvessels, Annals of Biomedical Engineering 35 (2007) 755-765.
- [54] B. Kaoui, G. H. Ristow, I. Cantat, C. Misbah, W. Zimmermann, Lateral migration of a two-dimensional vesicle in unbounded Poiseuille flow, Physical Review E 77 (2008) 021903.
- [55] S. K. Doddi, P. Bagchi, Lateral migration of a capsule in a plane Poiseuille flow in a channel, International Journal of Multiphase Flow 34 (2008) 966-986.
- [56] S. Meßlinger, B. Schmidt, H. Noguchi, G. Gompper, Dynamical regimes and hydrodynamic lift of viscous vesicles under shear, Physical Review E 80 (2009).
- [57] L. Shi, T.-W. Pan, R. Glowinski, Lateral migration and equilibrium shape and position of a single red blood cell in bounded Poiseuille flows, Physical Review E 86 (2012).
- [58] D. S. Hariprasad, T. W. Secomb, Prediction of noninertial focusing of red blood cells in Poiseuille flow, Physical Review E 92 (2015).
- [59] A. Farutin, C. Misbah, Analytical and Numerical Study of Three Main Migration Laws for Vesicles Under Flow, Physical Review Letters 110 (2013) 108104.
- [60] A. R. Pries, K. Ley, M. Claassen, P. Gaehtgens, Red cell distribution at microvascular bifurcations, Microvascular Research 38 (1989) 81-101.
- [61] D. Katanov, G. Gompper, D. A. Fedosov, Microvascular blood flow resistance: Role of red blood cell migration and dispersion, Microvascular Research 99 (2015) 57-66.
- [62] S. S. Ye, M. Ju, S. Kim, Recovery of cell-free layer and wall shear stress profile symmetry downstream of an arteriolar bifurcation, Microvascular Research 106 (2016) 14-23.
- [63] Y. C. Ng, B. Namgung, S. L. Tien, H. L. Leo, S. Kim, Symmetry recovery of cell-free layer after bifurcations of small arterioles in reduced flow conditions: e GLYPH&lt;11&gt; ect of RBC aggregation, American Journal of Physiology - Heart and Circulatory Physiology 311 (2016) H487-H497.
- [64] R. Fåhraeus, The Suspension Stability of the Blood, Physiological Reviews 9 (1929) 241-274.
- [65] K. H. Albrecht, P. Gaehtgens, A. Pries, M. Heuser, The Fahraeus e GLYPH&lt;11&gt; ect in narrow capillaries (i.d. 3.3 to 11.0 GLYPH&lt;22&gt; m), Microvascular Research 18 (1979) 33-47.
- [66] P. Olla, Simplified Model for Red Cell Dynamics in Small Blood Vessels, Physical Review Letters 82 (1999) 453-456.
- [67] U. Seifert, Hydrodynamic Lift on Bound Vesicles, Physical Review Letters 83 (1999) 876-879.
- [68] P. Olla, The lift on a tank-treading ellipsoidal cell in a bounded shear flow, Journal de Physique II 7 (1997) 1533-1540. ArXiv: chao-dyn 9610017. /
- [69] G. Danker, P. M. Vlahovska, C. Misbah, Vesicles in Poiseuille Flow, Physical Review Letters 102 (2009) 148102.
- [70] S. Nix, Y. Imai, T. Ishikawa, Lateral migration of a capsule in a parabolic flow, Journal of Biomechanics 49 (2016) 2249-2254.
- [71] S. Losserand, G. Coupier, T. Podgorski, Migration velocity of red blood cells in microchannels, Microvascular Research 124 (2019) 30-36.
- [72] H. Bruus, Theoretical Microfluidics, Oxford Master Series in Physics, OUP Oxford, 2008.
- [73] S. Kim, P. K. Ong, O. Yalcin, M. Intaglietta, P. C. Johnson, The cell-free layer in microvascular blood flow, Biorheology 46 (2009) 181-189.
- [74] S. Tripathi, Y. V. B. V. Kumar, A. Prabhakar, S. S. Joshi, A. Agrawal, Passive blood plasma separation at the microscale: a review of design principles and microdevices, Journal of Micromechanics and Microengineering 25 (2015) 083001.

- [75] S. Kim, R. L. Kong, A. S. Popel, M. Intaglietta, P. C. Johnson, Temporal and spatial variations of cell-free layer width in arterioles, American Journal of Physiology - Heart and Circulatory Physiology 293 (2007) H1526-H1535.
- [76] P. K. Ong, S. Jain, S. Kim, Spatio-temporal variations in cell-free layer formation near bifurcations of small arterioles, Microvascular Research 83 (2012) 118-125.
- [77] O. Oulaid, J. Zhang, Cell-free layer development process in the entrance region of microvessels, Biomechanics and Modeling in Mechanobiology 14 (2014) 783-794.
- [78] P. Balogh, P. Bagchi, The cell-free layer in simulated microvascular networks, Journal of Fluid Mechanics 864 (2019) 768-806.
- [79] Q. M. Qi, E. S. G. Shaqfeh, Theory to predict particle migration and margination in the pressure-driven channel flow of blood, Physical Review Fluids 2 (2017).
- [80] R. G. Henríquez Rivera, X. Zhang, M. D. Graham, Mechanistic theory of margination and flow-induced segregation in confined multicomponent suspensions: Simple shear and Poiseuille flows, Physical Review Fluids 1 (2016).
- [81] G. Segré, A. Silberberg, Radial Particle Displacements in Poiseuille Flow of Suspensions, Nature 189 (1961) 209-210.
- [82] J.-S. Park, S.-H. Song, H.-I. Jung, Continuous focusing of microparticles using inertial lift force and vorticity via multi-orifice microfluidic channels, Lab on a Chip 9 (2009) 939-948.

## Spatio-temporal dynamics of dilute red blood cell suspensions in low-inertia microchannel flow

## Supplemental Information (SI)

Qi Zhou 1 ; + , Joana Fidalgo 2 ; + , Lavinia Calvi 1 , Miguel O. Bernabeu , Peter R. Hoskins , 3 4 Mónica S. N. Oliveira 2 ; GLYPH&lt;3&gt; , Timm Krüger 1 ; GLYPH&lt;3&gt;

- 1. School of Engineering, Institute for Multiscale Thermofluids, University of Edinburgh, Edinburgh EH9 3FB, UK
- 2. James Weir Fluids Laboratory, Department of Mechanical and Aerospace Engineering, University of Strathclyde, Glasgow G1 1XJ, UK
- 3. Centre for Medical Informatics, Usher Institute, University of Edinburgh, Edinburgh EH16 4UX, UK
- 4. Centre for Cardiovascular Science, University of Edinburgh, Edinburgh EH16 4SB, UK + These authors contribute equally to this work.
- * Correspondence: monica.oliveira@strath.ac.uk or timm.krueger@ed.ac.uk

## S1. Scaling analysis

RBC migration. Consider an RBC whose center of mass is initially located at a position ( x 0 ; y 0 ) in close proximity to the wall ( y = 0) upon entry into a microchannel of hydraulic radius R h . To make the asymptotic analysis tractable, we consider a unidirectional flow along the axial direction of the channel with velocity u x ( y ), which is a function of the distance to the channel wall. Driven by the flow, the RBC will gradually migrate towards the channel centreline ( y = R h ) under hydrodynamic lift while travelling along the streamwise direction (i.e., x -direction), for which the motion can be described by a two-component velocity v v ( x ; v y ). Note that v is deliberately used here to represent cell velocity while u is reserved for fluid velocity. For leadingorder estimation, we assume v x : = u x ( y ). Based on the empirical law by Losserand et. al (1), we propose an expression for the lift velocity of the RBC v y as

$$v _ { y } = \dot { y } = K \dot { \gamma } ( y ) y ^ { - a }, \quad \dot { \gamma } ( y ) = \frac { \partial u _ { x } ( y ) } { \partial y },$$

where ˙ is the shear rate of the unperturbed flow, and GLYPH&lt;13&gt; K and a are constants depending on flow conditions and cell properties. For a given period of time, e.g., t 2 [0 ; t 0 ], the RBC travels to a new position ( x 0 ; y 0 ), and the axial distance ' can be written as

$$\ell = x ^ { \prime } - x _ { 0 } = \int _ { 0 } ^ { \prime } v _ { x } \, \mathrm d t = \int _ { y _ { 0 } } ^ { \prime } u _ { x } ( y ) \, \frac { \mathrm d t } { \mathrm d y } \mathrm d y = \int _ { y _ { 0 } } ^ { \prime } \frac { u _ { x } ( y ) } { \dot { y } } \, \mathrm d y.$$

Combining Eqs. (S1) and (S2), one finds:

$$\ell = \frac { 1 } { K } \int _ { y _ { 0 } } ^ { y ^ { \prime } } y ^ { a } v _ { x } ( y ) \left ( \frac { \partial u _ { x } ( y ) } { \partial y } \right ) ^ { - 1 } \, d y.$$

On assuming Poiseuille flow, the form of u x is known and v x may be expressed as

$$v _ { x } = u _ { x } ( y ) = \hat { u } _ { x } \left ( 1 - \left ( \frac { R _ { h } - y } { R _ { h } } \right ) ^ { 2 } \right ), \quad \hat { u } _ { x } = 2 \bar { u },$$

where ˆ u x and ¯ are the maximum and mean flow velocity, respectively. By replacing u v x in Eq. S3 with the above, one gets:

$$\alpha \gtrless. & & \ell = \frac { 1 } { 2 K } \int _ { y _ { 0 } } ^ { y ^ { \prime } } y ^ { a + 1 } \left ( 2 + \frac { y } { R _ { h } - y } \right ) \, d y. & & ( S 5 )$$

Given that lateral migration of the RBC away from the wall happens very slowly, we may further assume that y GLYPH&lt;28&gt; Rh over the time scale of our investigation. This means y = ( Rh GLYPH&lt;0&gt; y ) is small, which allows us to arrive at a straightforward expression of the axial distance ' that can be integrated over the lateral displacement:

$$\ell \approx \frac { 1 } { K } \int _ { y _ { 0 } } ^ { y ^ { \prime } } y ^ { a + 1 } \, d y = \, \frac { y ^ { a + 2 } } { ( a + 2 ) K } \Big | _ { y _ { 0 } } ^ { y ^ { \prime } }.$$

After some rearrangement, we may write the equation as

$$y ^ { \prime } = \left ( \begin{matrix} ( a + 2 ) K \ell + y _ { 0 } ^ { a + 2 } \end{matrix} \right ) ^ { \frac { 1 } { a + 2 } }.$$

CFL development. Now let us consider an arbitrary RBC cloud in the microchannel flow setting o GLYPH&lt;11&gt; from the wall. It is actually the outermost cells of the cloud that determine the CFL thickness GLYPH&lt;14&gt; cfl , representing the distance between the RBC membrane and the wall. At high Ca numbers (e.g., Ca = 0 6), the RBCs are highly deformable and significantly elongated along the flow, and : GLYPH&lt;14&gt; cfl can be approximated by y , namely the distance of the RBC's centre of mass to the wall:

$$\delta _ { \text{eff} } ( 0 ) & \approx y _ { 0 }, \quad \delta _ { \text{eff} } ( t ^ { \prime } ) & \approx y ^ { \prime }.$$

Combining Eqs. (S7) and (S8), we arrive at an expression describing the growth of GLYPH&lt;14&gt; cfl along the channel axis at sequential positions of ' :

$$\delta _ { \text{eff} } ( t ^ { \prime } ) = \left ( ( a + 2 ) K \ell + \delta _ { \text{eff} } ( 0 ) ^ { a + 2 } \right ) ^ { \frac { 1 } { a + 2 } }.$$

Ideally, if the outermost cells are directly in contact with the wall, GLYPH&lt;14&gt; cfl (0) = 0. This simplification results in a power-law scaling between GLYPH&lt;14&gt; cfl and ' , with an exponent of GLYPH&lt;12&gt; :

$$\delta _ { \text{eff} } \sim \ell ^ { \beta }, \quad \beta = \frac { 1 } { a + 2 }.$$

Specifically, for a = 1, corresponding to V l GLYPH&lt;24&gt; ˙ GLYPH&lt;13&gt;= h , the exponent becomes 1 3 and the below / scaling holds:

$$\delta _ { \text{eff} } & \sim \ell ^ { 1 / 3 }.$$

## S2. Experiments

S2.1. Preparation of RBC suspensions

Figure S1: Photomicrography of the horse RBC suspension in PS (NaCl 0.9% w v) at / di GLYPH&lt;11&gt; erent magnifications. (a) 40x, scale bar 10 µ m; (b) 10x, scale bar 50 µ m.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00250v2-figure-5bab6f5a-c421-42a5-8707-a6bcb31cdeb5.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00250v2-figure-e7233e9f-b964-43f8-9156-084d90d3f839.jpg)

The RBC suspensions were prepared using horse blood (Fig. S1). The blood sample was provided by TCS Biosciences (UK), in anticoagulant ethylenediaminetetraacetic acid (EDTA) 1.5 mg mL, and stored at 4 / GLYPH&lt;14&gt; C until further use. Sample processing and disposal were performed following the Ethics Committee for Animal Experimentation at University of Strathclyde. Since the blood was provided from di GLYPH&lt;11&gt; erent animals in distinct periods, the packed cell volume (PCV) or haematocrit ranges from 30 to 50% according to the blood provider. This variability is overcome by producing cell suspensions from the cell sediment after washing and centrifugation. The haematocrit ( Ht ), defined by the cell volume fraction, is controlled to be Ht GLYPH&lt;20&gt; 1%.

Fig. S2 presents the steps involved in sample processing. First, an Eppendorf with 1.5 mL of whole blood (WB) was centrifuged at 6000 rpm for 1 min (miniFUGE), after which the supernatant containing the yellow plasma and bu GLYPH&lt;11&gt; y coat (mainly white blood cells and platelets) was discarded. Second, the RBCs were washed twice with physiological saline PS (NaCl 0.9%

w v), with the transparent supernatant discarded. / Third, the final sample was prepared by suspending the desired volume of RBC in Dextran40 solution (Dx40 Mw = 40000 g mol, 0.1 g mL / / in isotonic medium). This suspending medium will avoid the phenomenon of cell sedimentation during experiments, which is a common problem for suspensions prepared in PS or in original blood plasma.

## S2.2. Data acquisition and image analysis

An inverted microscope (Olympus, IX71) was used to observe the device with bright field, employing a small magnification objective (10x 0.25 NA) in order to inspect cells across the / entire channel depth and throughout a relatively large section of the channel length. We make use of a high precision syringe pump (neMESYS, Cetoni GmbH) of the low-pressure and 14:1 gear ratio type, which is designed to deliver accurate flow rates (on the order of nanoliter min) for / microfluidic applications. SGE gastight syringes of appropriate volume were used to ensure that the specified pulsation-free minimum dosing rate is always exceeded. For each flow condition tested, the system was allowed to stabilize before image acquisition by a sensitive monochrome CCD camera (Olympus, XM10) at a frame rate of 15 Hz. The exposure time was set to 100 µ s allowing well-defined images of the cell boundaries. Videos of at least 300 frames each, corresponding to approximately 20 seconds of experiment, were acquired for four regions of interest (ROI). All the images were obtained at the channel centreplane ( z = 0). All ROIs have a length of 130 µ m. The images acquired need post-processing for image enhancement and cell detection. For this purpose, we applied an automatic cell counting routine developed in-house on ImageJ software.

## S2.2.1. Background correction

The bright field images obtained from experiments usually present an unbalanced intensity distribution, in particular close to the channel inlets and outlets, where the fittings are connected to the device. To enhance the images, they were first corrected by subtracting the background. For this purpose, the Rolling ball radius filter was used, applying a radius of 10 pixels, which corresponds to an approximate diameter of a horse RBC ( DRBC GLYPH&lt;25&gt; 6 µ m). This procedure allowed obtaining a clear and balanced image without losing information of the RBC (Fig. S3), especially for those cells located far away from the focal plane.

## S2.2.2. Cell detection

The RBCs in images with corrected background are then detected and counted using an automatic routine developed on ImageJ open software. Fig. S4 describes the intermediate steps

Figure S2: Steps involved in the preparation of RBC suspensions. (a) Eppendorf with whole blood; (b) separation (by centrifugation) of the supernatant plasma and bu GLYPH&lt;11&gt; y coat from the RBCs sediment; (c) transparent supernatant after RBC washing and centrifugation with PS; (d) RBC suspension in Dx40 solution ( Ht GLYPH&lt;20&gt; 1%).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00250v2-figure-93c76fc8-f8dc-47c6-b022-a401b2080d68.jpg)

Figure S3: Background correction using Rolling Ball Radius function in ImageJ. (a) Image sections before and after background subtraction; (b) Corresponding intensity profiles to the image sections in (a).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00250v2-figure-01e4e965-3322-4b64-8a0e-e65872680e5b.jpg)

Figure S4: Procedures of detecting cells from the experimental images. (a) Image with corrected background; (b) The corresponding binary image; (c) function Dilate; (d) function Fill Holes; (e) function Watershed; (f) function Outline; (g) detected cells after function Analyse Particles on applying the Circularity filter.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00250v2-figure-62f56454-d066-4aa2-bc0b-873bdd595462.jpg)

for further image enhancement and cell detection, based on the image with corrected background. First, a binary image (Fig. S4b) is generated from the image with corrected background (Fig. S4a), followed by function Dilate (Fig. S4c) and function Fill holes (Fig. S4d). This procedure allows recovering the cell contour for brighter cells that are located away from the focal plane. Subsequently, function Watershed is applied to generate separation among slightly overlapped cells (Fig. S4e). Then, function Outline allows delimiting each foreground object in the binary image, by generating a one-pixel wide outline (Fig. S4f). Based on this image, function Analyse Particles provides information of each delimited object, following which the Circularity filter is applied to detect objects whose circularity is in the range of 0.5 GLYPH&lt;24&gt; 1.0, defined by 4 GLYPH&lt;25&gt; A P = 2 (where A is the cell area and P the cell perimeter). The output is the image presented in Fig. S4g, providing numbered cells together with a list of values including the center of mass ( XCoM , YCoM ) and perimeter of each cell. A second filter is applied to the perimeter of the automatic detected object, in order to delete the ones presenting even smaller sizes than an undersized cell ( DRBC GLYPH&lt;25&gt; 4 µ m).

## S2.3. Probability density distribution (PDD) of RBCs

## S2.3.1. Calculation of local RBC fractions

Local RBC fractions in experiments are determined using the center of mass information from detected cells within each ROI. We divide the whole channel width into a number of bins ( GLYPH&lt;21&gt; 26) and count the number of cells that fall into each bin. The RBC fraction for each bin is then calculated as corresponding cell percentage (%) relative to the total number of cells detected across the whole channel width. By plotting the RBC fractions on the cross-section of the

designated ROI, the Probability density distribution (PDD) of RBCs is obtained, which can be regarded as an approximate representation of the haematocrit profile.

## S2.3.2. Accuracy of the PDD profile

To verify the accuracy of PDD profiles depicting RBC distribution obtained from the automatic counting method, two validation tests were performed. First, we examined the minimum number of frames needed to ensure a converged distribution. For a given ROI, the RBC distribution was determined using four di GLYPH&lt;11&gt; erent image stacks, containing 50, 100, 200 and 300 image frames, respectively. The PDD profiles resulting from a relatively small stack, e.g., 50 and 100 frames, were found to present significant noises and di GLYPH&lt;11&gt; er much from each other (Fig. S5a). In contrast, the distributions obtained from a larger stack, namely 200 and 300 frames, were much smoother and appeared similar to each other (Fig. S5a). Therefore, we conclude that a stack with 300 frames is enough to generate a converged RBC distribution.

Second, we compared the RBC distribution obtained from our automatic counting using a stack of 300 frames with one based on manual counting. The result showed a high-level similarity between the two distributions (Fig. S5b), hence validating the accuracy of our proposed method for the analysis of RBC distribution in experiments using low haematocrit ( Ht GLYPH&lt;20&gt; 1%), where the boundaries of individual cells can be well defined.

## S2.3.3. E ect of cell overlapping GLYPH&lt;11&gt;

The e GLYPH&lt;11&gt; ect of cell overlapping is also investigated (Fig. S5c). In Fig. S4g, it is possible to observe false 'agglomerates' (given the RBC concentration and the type of Dextran used in our experiments, we expect no physical aggregation) caused by the overlapping of cells across the channel depth direction. To study the impact of these 'agglomerates' on the calculated RBC distribution, they were identified and excluded in a test case. The resulting RBC distribution was then compared to the original RBC distribution, where the two distributions obtained are found virtually identical (Fig. S5c). This suggests that the statistical e GLYPH&lt;11&gt; ect of cell overlapping is negligible on the final cell distribution obtained using our automatic counting method.

Figure S5: (a) E GLYPH&lt;11&gt; ect of the number of frames on the RBC distribution calculated by the automatic cell counting method. The PDD profiles are calculated using a stack of 50 (black line), 100 (green line), 200 (blue line) and 300 (red line) frames, respectively. (b) PDD profiles (using a stack of 300 images) obtained from automatic cell counting (red line) and manual cell counting (black line), respectively. (c) E GLYPH&lt;11&gt; ect of cell overlapping on the calculated RBC distribution, showing PDD profiles obtained from automatic counting with (red line) and without (black line) overlapped cells in the images, respectively.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00250v2-figure-d4952bf3-8f14-4ce0-a830-1d39983f159e.jpg)

## S2.3.4. E ect of flow-inlet misalignment GLYPH&lt;11&gt;

To test the e GLYPH&lt;11&gt; ect of inlet conditions on the experimental results, the inlet ports are manually punched onto the PDMS, therefore the tubing connection can be either aligned or misaligned

relative to the centreline of the microdevice. Here, Fig. S6 shows the evolution of PDD profiles under the misaligned configuration, corresponding to the channel in Fig. 1c of the main text.

Figure S6: PDDs monitored at x = 65 GLYPH&lt;6&gt; 65, 565 GLYPH&lt;6&gt; 65, 1220 GLYPH&lt;6&gt; 65, 2065 GLYPH&lt;6&gt; 65 µ maway from the entrance in the misaligned configuration, under a volume flow rate of Q = 0 2 : µ L min. /

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00250v2-figure-812182e3-e2f0-4fa1-a61d-ff6dcc866c47.jpg)

## S2.3.5. E ect of volume flow rate GLYPH&lt;11&gt;

The volume flow rate was found to have an impact on the initial PDD profile at x = 65 µ m. In general, a lower Q (e.g., Q = 0 2 : µ L min) yields a virtually flat density profile as Fig. 2a / of the main text shows, whereas a higher Q (e.g., Q = 4 0 : µ L min) gives rise to appreciable / density peaks near the walls as in Fig. S7b here. This dependence of initial RBC distribution on Q presumably arises because the cell deformation changes significantly within the investigated range of flow rates, thus non-trivially altering the lift forces imposed on the cell.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00250v2-figure-f02772d5-3a8a-4870-a0f3-14f4a26b93bb.jpg)

Figure S7: Experimental PDD profiles monitored at x = 65, 2065 µ mfor a H F = 1% RBC suspension in the aligned configuration (see Fig. 1b of the main text). The volume flow rates are (a) Q = 0 8 : µ L min; (b) / Q = 4 0 : µ L min. /

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00250v2-figure-cb02c91e-c15a-4113-818d-7464d3253de5.jpg)

## S2.4. Cell-free layer (CFL)

To measure the CFL at a ROI within the channel, a combined image of 300 frames via minimum intensity Z-projection is generated (Fig. S8a). Using this method, if a cell is passing through the ROI with lower intensity than the background, it will be registered; with a su GLYPH&lt;14&gt; cient number of image frames, the layer depleted of cells (CFL) close to the walls can be distinguished. The binary image of such a projection (Fig. S8b) presents a sharp boundary separating the RBC

core and the CFL. The average CFL thickness is then defined as half of the averaged di GLYPH&lt;11&gt; erence between the channel width and the RBC core along the ROI.

Figure S8: Determination of the CFL thickness. (a) Minimum intensity Z-projection image using a stack of 300 frames. (b) Binary image of the Z-projection for calculation of the CFL thickness.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00250v2-figure-c52fed74-612a-485b-b962-7df599048ef5.jpg)

Figure S9: (a) Thickness of the RBC core measured at selected windows of the channel. (b) Calculated CFL values (accounting for both walls) based on RBC core thickness and channel width. (a-b) are for experiment Q = 0 2 : µ L min. /

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00250v2-figure-ae2342da-504a-4cef-9f94-a332c72a2b4c.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00250v2-figure-348bd816-25f9-42ca-a6f4-a422bf4a8024.jpg)

## S3. Simulations

## S3.1. Model configuration

Fig. S10 shows the 3D computational model for simulating the RBC flow through a rectangular microchannel. Two geometries have been designed in the simulation to account for potential entrance e GLYPH&lt;11&gt; ects that exist in experiments due to inlet configuration, which primarily arise from the inflow contraction between the entry region and the main channel. The main channel part in both geometries is the same, which is 96 µ m wide, 30 µ m deep and 1.32 mm long. The only di GLYPH&lt;11&gt; erence is the front-end diameter of the tapered entry region (100 µ m long). The entry region for the first geometry has a front-end diameter of 96 µ m (hereafter referred to as CTRAC-S, meaning 'smaller (inflow) contraction'), and the one for the second geometry has a front-end diameter of 180 µ m (hereafter referred to as CTRAC-L, meaning 'larger (inflow) contraction'. All simulation data and plots presented in this paper are from CTRAC-L, unless otherwise specified). An identical entry region is also added to the channel outlet for removing cells from the simulation. The whole domain measures 1.52 mm in length, with a total of 26,201,764 lattice voxels (each GLYPH&lt;1&gt; x = 2 3 / µ m. All dimensions of the domain are labelled in Fig. S10, shown from three camera views.

Figure S10: Geometry of the computational model (CTRAC-L, larger contraction). (a) is the top view of the geometry, showing the widthwise span of the channel. In (a), the main channel (rectangular) is marked with two blue dotted lines. (b) shows cross-sectional views of the front-end and rear-end of the tapered entry region, the latter of which coincides with the channel entrance ( x = 0). (c) is the side view of the geometry, showing the depthwise span of the channel.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00250v2-figure-450f8a64-eedd-4c32-9421-dc7b93959d41.jpg)

Figure S11: A snapshot of the simulation at t = 0.37s, showing distribution of the 1% RBC suspension in the width and depth directions of the microchannel. Cells arriving at the end of the channel are removed from the system.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00250v2-figure-72bdbd00-5818-4d04-90a5-6071247386ee.jpg)

## S3.2. Simulation setup

Initially, a plasma flow (no cells) is driven from left to right along the x -axis direction with a fixed volume flow rate (Fig. S11). The flow is controlled by imposing a parabolic velocity profile (assuming Poiseuille flow) at the inlet and a reference pressure at the outlet. The no-slip condition on walls is implemented with the Bouzidi-Firdaouss-Lallemand method (2). Because it is computationally prohibitive to simulate the realistic experimental condition where the particle Reynolds number approaches zero, the fluid flow in our simulation is numerically accelerated to a state where the particle Reynolds number equals 0.03 (still within the low-inertia regime). After the plasma flow is converged, RBCs are randomly inserted from the cylindrical entry region in a continuous manner with fixed haematocrit ( H F = 1%). To make the experiments and simulations comparable, we establish the similitude law based on the capillary number Ca . For the target Ca = 0 6, we have fairly deformable RBCs in the channel assuming tank-treading motion, where : the cell membrane is fluidized and rolls around the cell interior.

Figure S12: Haematocrit ( Ht ) histograms (bars) at x = 3 3 : D h, 10 D h, 16 5 : D h, 23 D h downstream of the channel, showing (a,b,c) the widthwise and (b,c,d) the depthwise distribution of RBCs at H F = 1%. The inset in each panel shows the approximate Ht profiles enveloping the calculated histograms for visual clarity. (a-b) are for a simulation in a channel with AR = 3 2 and : Ca = 0 1. (c-d) are for a simulation in a channel with : AR = 1 7 and : Ca = 0 6. (e-f) are for a : simulation in a channel with AR = 1 0 and and : Ca = 0 6. :

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00250v2-figure-f78b7900-3ce1-4993-a626-d141ca2d565a.jpg)

## S3.3. Numerical analysis

## S3.3.1. Calculation of RBC velocities

The lateral or axial velocity for an RBC is calculated as the unit displacement of its center of mass along the specific direction over a unit time step in the simulation:

$$V _ { \mathrm l, y } = \Delta L _ { y } / \Delta t, \quad V _ { \mathrm l, z } = \Delta L _ { z } / \Delta t, \quad V _ { \mathrm x } = \Delta L _ { x } / \Delta t$$

Because the trajectories for all RBCs are recorded throughout the simulation, it is straightforward to group the cells into bins according to the location of their center of mass and conduct statistical analysis to estimate the time-average lift velocity (Fig. 4) or streamwise velocity (Fig. S15) of cells at a given position in the microchannel.

## S3.3.2. Haematocrit calculation

Once the RBCs are fully perfused inside the channel and the number of cells within becomes stable, snapshots of the simulation are extracted at given time intervals, which are then assembled to count the cells for the calculation of time-averaged haematocrit profiles. For this dilute suspension, the number of time steps is set as n t = 100 (with 2.22 ms interval between each) to ensure su GLYPH&lt;14&gt; cient cells existing in the sampling pool. Next, various ROIs alongside the channel length are selected to count cells locally, starting at a position 150 µ m away from the entrance of the main channel. In total, 10 target positions between 150 µ m and 1050 µ m are chosen, with a 100 µ mincrement. At each target position, we consider cells contained in a W GLYPH&lt;2&gt; H GLYPH&lt;2&gt; GLYPH&lt;1&gt; L sampling box. All cells within this box are counted based on their centre of mass and allocated into a number of subdivisions either across the widthwise ( W ) direction or the depthwise ( H ) direction. The local haematocrits Ht i ( ) in each subdivision are then calculated as

$$H t _ { ( i ) } = \frac { n _ { t b c ( i ) } V o l _ { t b c } } { n _ { t } V o l _ { b i n } }, \ \ V o l _ { b i n } = \frac { V o l _ { b o x } } { N _ { b i n } }, \ \ ( i = 1, 2, \dots, N _ { b i n } )$$

where n rbc( ) i is the number of cells detected in a given subdivision, Vol rbc = 100 fL is the standard volume of an RBC, Vol bin is the volume of each subdivision, Vol box = W GLYPH&lt;2&gt; H GLYPH&lt;2&gt; GLYPH&lt;1&gt; L ( GLYPH&lt;1&gt; L = 50 µ m)

is the volume of the sampling box, and N bin = 10 is the number of subdivisions in Vol box . Fig. S12 shows the obtained haematocrit histograms in both the width and depth directions of the channel from two simulations.

## S3.3.3. Calculation of RBC fluxes

Based on the tube haematocrits Ht i ( ) (subdivided into Nbin bins, Eq. S2) at a distance x downstream of the flow inlet, the cross-sectional RBC fluxes Q rbc( ) i allocated into the same number of bins can be calculated via the discharge haematocrit HD i ( ) and blood volume flow rate Q B( ) i :

$$Q _ { r b c ( i ) } = H _ { \text{D} ( i ) } Q _ { \text{B} ( i ) }, \quad ( i = 1, 2, \dots N _ { \text{bin} } )$$

$$H _ { \text{D} ( i ) } = H t _ { ( i ) } \left ( \frac { \bar { V } _ { \text{x} ( i ) } } { \bar { U } _ { \text{B} ( i ) } } \right ), \quad Q _ { \text{B} ( i ) } = \bar { U } _ { \text{B} ( i ) } A _ { \text{bin} }$$

where A bin = W GLYPH&lt;2&gt; H N = bin is the cross-sectional area of each bin, and ¯ V x( ) i and ¯ U B( ) i are the spatial-average RBC velocity and blood flow velocity within each bin, respectively. Combining Eqs. (S3) and (S4), we have the following equation to be used for RBC flux calculation:

$$Q _ { \text{rbc} ( i ) } & = H t _ { ( i ) } \bar { V } _ { \text{xi} ( i ) } A _ { \text{bin} }$$

## S3.3.4. Calculation of cell-free layer

To characterise the CFL in a rectangular channel with width W and depth H , we consider a cross-section at a distance x downstream of the flow inlet to be divided into four triangles by the diagonals and the wall edges. The heights of the four triangles, each perpendicular to a channel wall, are located at GLYPH&lt;18&gt; = 0 ; GLYPH&lt;25&gt;= 2 ; GLYPH&lt;25&gt;; 3 GLYPH&lt;25&gt;= 2, respectively. We define the RBC density GLYPH&lt;30&gt; ( h ; GLYPH&lt;18&gt; 0 ; x t ; ) at a point p x h ( ; ; GLYPH&lt;18&gt; 0 ) of longitudinal coordinate x , height coordinate h and angular coordinate GLYPH&lt;18&gt; 0 within a given triangle for a certain time instant t as 1 if that point is enclosed by at least one RBC membrane; otherwise, GLYPH&lt;30&gt; ( h ; GLYPH&lt;18&gt; 0 ; x t ; ) = 0. The time-mean RBC density within any of the fours triangles can then be obtained via integration over ( h ; GLYPH&lt;18&gt; 0 ) and average over t :

$$\bar { \phi } ( x, \theta ) = \frac { 1 } { N } \sum _ { i = 1 } ^ { N } \int _ { 0 } ^ { L _ { h } } \int _ { \theta - \Theta / 2 } ^ { \theta + \Theta / 2 } \phi ( h, \theta ^ { \prime }, x, t _ { i } ) h \tan | \theta ^ { \prime } - \theta | d \theta ^ { \prime } d h$$

wherein N is the number of time steps extracted from the simulation, L h is the height of the triangle and GLYPH&lt;2&gt; is its tip angle centred in the cross-section. For the left and right triangles:

$$L _ { \mathfrak { h } } = W / 2, \quad \Theta = 2 \arctan ( H / W )$$

and for the top and bottom triangles:

$$L _ { \mathfrak { h } } = H / 2, \quad \Theta = 2 \arctan ( W / H )$$

If we define the CFL as a thin layer in proximity to a given wall (base of the triangle) which contains only a negligible fraction of RBC density " , e.g., " &lt; 1 e GLYPH&lt;0&gt; 6, then the thickness of this layer GLYPH&lt;14&gt; ( x ; GLYPH&lt;18&gt; ) can be determined using numerical iterations:

$$\frac { 1 } { \bar { \phi } ( x, \theta ) N } \sum _ { i = 1 } ^ { N } \int _ { L _ { 0 } - \delta ( x, \theta ) } ^ { L _ { 0 } } \int _ { \theta - \Theta / 2 } ^ { \theta + \Theta / 2 } \phi ( h, \theta ^ { \prime }, x, t _ { i } ) h \tan | \theta ^ { \prime } - \theta | d \theta ^ { \prime } d h \geq \varepsilon \quad \text{(S9)}$$

By calculating the threshold thickness GLYPH&lt;14&gt; ( x ; GLYPH&lt;18&gt; ) satisfying the above criteria for each GLYPH&lt;18&gt; at consecutive x locations, we can obtain the CFL values for each wall along the channel-axis direction as in Fig. S13. The development of these CFLs is further plotted in log-log scale against a 1 3 slope / indicating power-law growth (Fig. S14).

Figure S13: Numerical CFLs calculated along each wall of the channel (Right Bottom Left Top-walls) from di / / / GLYPH&lt;11&gt; erent simulations. (a) H F = 1%, Ca = 0 6, AR : = 3.2; (b) H F = 1%, Ca = 0 6, AR : = 1.7; (c) H F = 1%, Ca = 0 6, AR : = 1.0; (d) H F = 1%, Ca = 0 1, AR : = 3.2; (e) H F = 2%, Ca = 0 6, AR : = 3.2; (f) H F = 5%, Ca = 0 6, AR : = 3.2. The inset in each panel shows a cross-sectional view of the RBC pattern at L = 23 D h away from the entrance, combining snapshots from 50 consecutive time steps from the corresponding simulation.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00250v2-figure-c4252039-8209-4ee9-b896-4c314171615c.jpg)

Figure S14: Log-log plots of the CFLs as in Fig. S13 against a 1 3 power-law trendline. The CFL values / GLYPH&lt;14&gt; cfl are normalised by the maximum CFL thickness GLYPH&lt;14&gt; max detected within the whole investigated range, and the development lengths L are normalised by the channel hydraulic diameter D h. The inset in each panel show a zoom-in view of the range L = 10 GLYPH&lt;24&gt; 28 Dh .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00250v2-figure-25b25b8b-5775-4138-974a-66f9197d62b3.jpg)

Figure S15: 3D flow in the rectangular channel. (a-b) (Simulation) Transient mid-plane velocity profiles near the entrance at positions labelled by white solid lines in the channel sketch. (c-d) (Simulation) Steady velocity profiles at x = 50 µ m. z and y values in the legends represent positions relative to the mid-plane in each direction, with z = 0 and y = 0 denoting the widthwise and depthwise mid-planes, respectively. (e-f) (Analytical) Shear rates (main frame) and shear-rate gradients (inset) calculated from asymptotic solutions in Ref. (3). The colors represent z or y values as in (c-d). The black vertical lines in (e-f) indicate the center of the outmost layer of cells at x = 50 µ mdetected in the RBC simulation along the (e) y -direction and (f) z -direction, respectively (see Fig. 4).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00250v2-figure-1d105621-02de-44a3-8784-35a2071ed1a6.jpg)

## S3.4. Geometric e GLYPH&lt;11&gt; ect and shear anisotropy

In our problem, where a laminar flow is pressure-driven in a rectangular microchannel of high aspect ratio ( W H = = 3 2), the flow field is complex because distinct velocity profiles co-: exist in the channel cross-section: a typical parabola characteristic of Poiseuille flow in the depth ( H ) direction and a blunted profile featuring plug-like flow in the width ( W ) direction (3). Additionally, the existence of entrance e GLYPH&lt;11&gt; ects further increases the complexity of the problem. These two aspects contribute to a plethora of scenarios for the velocity profile curvature (reflecting the shear-rate gradient), which on the one hand constitutes a mechanism for cell migration itself (by introducing symmetry-breaking deformation, see Ref. (4)) and on the other hand impacts the wall-induced lift (being a function of the particle shape, see Ref. (5)). We show the evolution of unperturbed mid-plane velocity profiles (in the absence of cells) near the entrance across both directions, where rich behaviour of the velocity profile can be found as the fluid accelerates through the contracted transition between the entry region and the channel itself (Fig. S15a-b). This geometry-induced velocity alteration potentially underpins the distinct initial distribution of RBCs between the W and H directions of the channel observed in Fig. 2a,d.

Fig. S15a-b reveal that the flow becomes fully-developed between x = 30 GLYPH&lt;24&gt; 50 µ m in our simulations. Fig. S15c-d show the converged velocity profiles velocity profiles for di GLYPH&lt;11&gt; erent zand y-planes at x = 50 µ m, respectively. Using the asymptotic solution of this microfluidic flow (3), we further calculate and plot the shear-rate profiles and shear-rate-gradient profiles. The cross-sectional pattern of shear rates (normalised by the time scale GLYPH&lt;28&gt; as in Eq. 2) turns out more intricate than often assumed (Fig. S15e-f), even in the narrow dimension ( H ) where linearity

Figure S16: Velocity profiles of the RBC-laden flow at x = 50 µ m. z and y values in the legends represent positions relative to the mid-plane in each direction, with z = 0 and y = 0 denoting the widthwise and depthwise mid-planes, respectively. The blood flow velocity U B is normalised by mean velocity of the unperturbed plasma flow ¯. u

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00250v2-figure-778dafa6-6663-489b-90fd-0deee1d0e6d0.jpg)

is commonly assumed for such a high aspect ratio ( W H = &gt; 3) (6). At close distances to the wall (e.g., y = 45 µ m in Fig. S15d), the H -direction velocity profile deviates significantly from a parabolic one (Fig. S15f), resulting in non-linear decrease of shear rate on approaching the channel centreline and varying shear gradients (inset of Fig. S15f). In the W direction where blunted velocity profiles prevail (Fig. S15c), the shear gradients vary both at close and large distances to the wall (inset of Fig. S15e), leaving the hydrodynamic lift even more unpredictable by the classic theory derived on the basis of ideal shear flows (4).

## S3.5. Disturbance of RBCs to the plasma flow

The velocity profiles from the RBC-laden flow at Ht 1% (Fig. S16) apparently do not deviate much from the unperturbed flow in the absence of RBCs (Fig. S15c,d). This implies that the sparse distribution of RBCs in the dilute limit only a GLYPH&lt;11&gt; ects the flow velocity to a limited amount and causes local disturbances without modifying the overall velocity profile significantly. Therefore, it is reasonable to analyse the migration of RBCs in the 1% suspension using the unperturbed plasma flow.

## S3.6. RBC lateral migration and axial acceleration

## S3.6.1. Cross-sectional RBC flux distribution

The Ht profile is co-determined by the lateral migration (toward channel centreline) and the streamwise acceleration (along channel axis) of RBCs, the former of which arises from the hydrodynamic lift e GLYPH&lt;11&gt; ect and the latter is associated with the Fåhræus e GLYPH&lt;11&gt; ect. Both e GLYPH&lt;11&gt; ects have an impact the distribution of local haematocrits and can a GLYPH&lt;11&gt; ect the dynamical evolution of the Ht profile. In this sense, the peaks emerging in the Ht profile (Fig. 2) may not be ideal markers of the 'pseudo-equilibrium' position for the inward migration of cells. Therefore, to ascertain the 'pseudo-equilibrium' positions suggested by Ht peaks, these two e GLYPH&lt;11&gt; ects need to be decoupled. To this end, we calculate the RBC fluxes, the sum of which strictly conserves along the channel axis and examine the cross-sectional distribution of the RBC fluxes. This allows us to exclude the e GLYPH&lt;11&gt; ect of cell acceleration and distinguish the amount of lateral migration of the cells (Fig. S17).

The evolution of RBC-flux distribution turns out qualitatively similar to that of the Ht profiles. It reveals a clear trend of inward RBC migration over time, with near-wall RBC fluxes moving away and subsequently accumulating at a 'pseudo-equilibrium' position in both directions of the channel. For the widthwise migration, this equilibrium position is found to be 0.3 times of half channel width away from the wall; for the depthwise migration, it is 0.7 times of half channel away. The inward shift of RBC fluxes is accompanied by the continuous expansion of a zero-flux zone in the vicinity of the wall, which corresponds to the gradually growing CFL.

Figure S17: (Simulation) Evolution of RBC fluxes Q rbc across the channel in (a) the width direction and (b) the depth direction. The four colors red, green, blue, black indicate a distance of x = 150 450 750 1050 ; ; ; µ maway from the entrance, respectively. All RBC fluxes are normalised by the volume flow rate of the unperturbed flow Q 0 and expressed in percentages.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00250v2-figure-d9fb577a-6414-4c01-a107-c9ebe4d370bc.jpg)

## S3.6.2. Axial acceleration of RBCs

Figure S18: (Simulation) Axial velocity Vx (streamwise, normalised by mean velocity of the unperturbed flow ¯) of u RBCs in the channel at x = 150 450 750 1050 ; ; ; µ maway from the entrance. (a) and (b) show the width and depth directions of the channel, respectively. For statistical analysis in both (a) and (b), the channel is fictitiously folded along the channel centreline and then divided into 10 bins, followed by the assignment of cells according to the position of their center of mass. Each error bar here shows the mean and standard deviation of the streamwise velocity.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00250v2-figure-4bc0b062-1d86-4d69-ab45-72fc8b2610b3.jpg)

## Supporting references

- 1. Losserand, S., G. Coupier, and T. Podgorski, 2019. Migration velocity of red blood cells in microchannels. Microvascular Research 124:30-36.
- 2. Bouzidi, M., M. Firdaouss, and P. Lallemand, 2001. Momentum transfer of a Boltzmannlattice fluid with boundaries. Physics of Fluids 13:3452-3459.
- 3. Bruus, H., 2008. Theoretical Microfluidics. Oxford Master Series in Physics. OUP Oxford.
- 4. Kaoui, B., G. H. Ristow, I. Cantat, C. Misbah, and W. Zimmermann, 2008. Lateral migration of a two-dimensional vesicle in unbounded Poiseuille flow. Physical Review E 77:021903.
- 5. Olla, P., 1997. The lift on a tank-treading ellipsoidal cell in a bounded shear flow. Journal de Physique II 7:1533-1540.
- 6. Coupier, G., B. Kaoui, T. Podgorski, and C. Misbah, 2008. Noninertial lateral migration of vesicles in bounded Poiseuille flow. Physics of Fluids 20:111702.
Qi Zhou Joana Fidalgo Lavinia Calvi Miguel O. Bernabeu Peter R. Hoskins Mónica S. N. Oliveira Timm Krüger
Spatio-temporal dynamics of dilute red blood cell suspensions in a microchannel flow at low Reynolds number
Microfluidic technologies are commonly used for the manipulation of red blood
cell (RBC) suspensions and analyses of flow-mediated biomechanics. To enhance
the performance of microfluidic devices, understanding the dynamics of the
suspensions processed within is crucial. We report novel aspects of the
spatio-temporal dynamics of RBC suspensions flowing through a typical
microchannel at low Reynolds number. Through experiments with dilute RBC
suspensions, we find an off-centre two-peak (OCTP) profile of cells contrary to
the centralised distribution commonly reported for low-inertia flows. This is
reminiscent of the well-known "tubular pinch effect" which arises from inertial
effects. However, given the conditions of negligible inertia in our
experiments, an alternative explanation is needed for this OCTP profile. Our
massively-parallel simulations of RBC flow in real-size microfluidic dimensions
using the immersed-boundary-lattice-Boltzmann method (IB-LBM) confirm the
experimental findings and elucidate the underlying mechanism for the
counterintuitive RBC pattern. By analysing the RBC migration and cell-free
layer (CFL) development within a high-aspect-ratio channel, we show that such a
distribution is co-determined by the spatial decay of hydrodynamic lift and the
global deficiency of cell dispersion in dilute suspensions. We find a CFL
development length greater than 46 and 28 hydraulic diameters in the experiment
and simulation, respectively, exceeding typical lengths of microfluidic
designs. Our work highlights the key role of transient cell distribution in
dilute suspensions, which may negatively affect the reliability of experimental
results if not taken into account.