In vitro investigations of red blood cell phase separation in a complex microchannel network

Living tissues rely on the supply of oxygen and nutrients as well as on the removal of carbon dioxide and other metabolites²¹ to maintain a physiological state. The microcirculation is the system where these mass exchange phenomena are happening. Microvascular networks feature a complex topology with highly interconnected capillaries bifurcating in an irregular fashion.^10,22 Blood flowing in the microcirculation cannot be treated as a homogeneous Newtonian fluid because it contains deformable red blood cells (RBCs), whose dimensions are comparable to the vessel diameter.⁵⁰ The RBC concentration is heterogeneous in space and time, and there are even RBC-free vessels.⁴⁷ Therefore, the particulate nature of blood and phase separation effects must be considered for a correct interpretation of fluid dynamics at the microscale. In such microvascular networks, the dynamics of RBCs flowing through sequentially diverging and converging bifurcations plays an important role in the local perfusion.

Nonuniform RBC partitioning at diverging bifurcations leads to the heterogeneous RBC distribution in the daughter branches and ultimately impacts the oxygen delivery to living tissues.⁴¹ The daughter branch with a higher flow rate tends to receive a disproportionally higher RBC fraction (the Zweifach-Fung effect¹⁷) at the cost of a RBC reduction in the branch with a lower flow rate. The extreme case is a daughter vessel with a very high flow rate drawing all RBCs and leaving the low-flow branch with only plasma. Fung¹⁷ reported that the RBC flow at microvascular bifurcations is self-regulated: RBCs entering the high-flow branch cause an increased RBC density in this branch, which leads to a higher local resistance. Eventually, this will reduce the flow in this branch to the point when RBCs start entering the other branch. This self-regulating mechanism leads to a reduction of flow disparities between daughter vessels. Schmid et al.⁴⁹ showed that such well-balanced bifurcations dominate cerebral microvascular networks in mice and support robust perfusion characteristics. The Zweifach-Fung effect is strongest in capillaries, where the vessel diameter is comparable to the RBC diameter (strong confinement).⁵⁰ In larger vessels, RBCs tend to follow the streamlines of the plasma flow.⁶ 

Given the high importance of the microcirculation in human physiology, blood flow distribution at microvascular bifurcations has been thoroughly studied over the past decades. Phase separation has been investigated with two-dimensional^6,7,12,55 and three-dimensional²⁵ numerical simulations in symmetric and asymmetric bifurcations. In particular, Barber et al.^6,7 and Doyeux et al.¹² described mechanisms related to the cell partitioning causing the individual cell trajectories to divert from fluid streamlines. In vitro experiments have also been reported where the flow partitioning of rigid particles,⁴⁶ lipid vesicles,¹² or red cells^{9,10,14,47,53} was studied in models of symmetric and asymmetric bifurcations (T-shaped or Y-shaped). Fenton et al.¹⁴ and Roman et al.⁴⁷ studied the effects of the feeding hematocrit and the local bifurcation geometry. A reduction of the phase separation for increasing hematocrit was observed in symmetric bifurcations.⁴⁷ Carr and Wickam,⁹ Roberts and Olbricht,⁴⁶ and Sherwood et al.⁵³ showed how the proximity of two diverging bifurcations increases the asymmetry of the phase separation. Sherwood et al.⁵³ demonstrated that in the case of weak confinement (particle to channel diameter ratio equal to 0.16), the phase separation asymmetry was due to a skewed hematocrit distribution in the parent vessel of the bifurcations. In an experimental study of RBC partitioning with an asymmetric bifurcation, Clavica et al.¹⁰ observed a reduction and even an inversion of the Zweifach-Fung effect with increasing flow rate.

A number of computational studies have investigated large-scale networks, where the geometry of the networks was either idealized^1,39,48 or obtained from the reconstruction of in vivo images.^{5,18,20,31,42,44} Numerical results of idealized³⁹ or realistic⁵ networks showed that the RBC distribution was heterogeneous, and this was even the case if the networks were symmetric and the boundary conditions at the inlet were kept constant.³⁹ Balogh and Bagchi⁵ developed a direct numerical simulation technique for blood flow modeling in microvascular networks with in vivo-like features. They showed how RBC partitioning could oscillate between classical partitioning (the Zweifach-Fung effect) and reverse partitioning (inversion of the Zweifach-Fung effect^10,52) over time. Different cellular-scale mechanisms were identified that could affect the local RBC distribution: (I) skewness of the hematocrit profile in the parent branch; (II) RBC lingering at branching points;⁴ (III) local self-regulation due to cell-cell and cell-vasculature interactions; and (IV) flow balancing due to transient capillary obstructions.⁵ 

In contrast to computational models, in vitro studies to investigate the local RBC partitioning for complex and anatomically plausible microvascular networks^15,28,34,54 are scarce. In particular, there is a lack of experimental research with channel sizes comparable to the RBC diameter (⁠$≤10μm$⁠), which can probably be attributed to technical challenges of microfabrication.

In the present study, a complex yet idealized capillary network was designed and fabricated. The network topology was inspired by studies previously reported in the literature.^34,35,44,48 All microchannels had a width of $10μm$ (particle to channel diameter ratio $∼0.6$⁠) to reproduce the typical confinement that RBCs experience in microvascular capillary networks. This allowed us to provide for the first time quantitative in vitro data on phase separation in a complex network with many sequential diverging and converging bifurcations for different inflow conditions. This also includes data for bifurcations with RBC-free daughter vessels and information on the flow threshold for RBC-free vessels. It will be shown that the observed phase separation is rather well described by the model of Pries et al.,⁴¹ although some bifurcations exhibited an inversion of the Zweifach-Fung effect. Further, we investigated the relation between the hematocrit distribution in the parent vessels of a series of diverging bifurcations and the RBC partitioning, with a particular interest on the consequences for the Zweifach-Fung effect.

We used a microfluidic device comprising a complex honeycomb network of microchannels whose dimensions are inspired by in vivo length scales.³³ A similar design was used by Schmid et al.⁴⁸ to perform in silico studies and by Merlo^34,35 for in vitro experiments. Microchannels had rectangular cross sections (⁠$width=10μm$⁠, height$=8μ$m), and a length of $85μm$⁠. Wider drainage microchannels (⁠$width=100μm$⁠, $height=8μm$⁠) were placed upstream and downstream of the capillary network at the inlet and outlet. The typical range of RBC velocities in the microcirculation (⁠$0−2mm/s$⁠)^23,36 was reproduced by a constant pressure head at the inlet of the device.

The geometry of the microdevice was designed using DraftSight (Dassault Systèmes, Vélizy-Villacoublay, France) and transferred onto a chrome photomask (JD Photodata, Hitchin, UK). The master of the microdevice was made from a silicon wafer (Prolog Semicor Ltd, Kiev, Ukraine) and fabricated using conventional soft-lithography. The process involved the following steps: coating of the wafer with a SU-8 negative photoresist (⁠$thickness=8μm$⁠, GM1060 from Gersteltec Sarl., Pully, Switzerland) and patterning by exposure to UV light (Suss MicroTec Lithography GmbH, Garching, Germany) through the photomask. The microfluidic device was then made of polydimethylsiloxane (PDMS) via replica molding by pouring liquid PDMS (Sylgard 184, Dow Corning, Midland, MI, USA) with a monomer/curing agent weight ratio of 10:1 onto the master device. After degassing for approximately 25 min, the microfluidic chip was cured overnight at $60°C$⁠. In the final fabrication step, the casted PDMS was bonded to a flat PDMS layer by oxygen-activated plasma treatment (Harrick Plasma, Ithaca, NY, USA). The resulting device had a honeycomblike geometry with 176 microchannels and the actual dimensions (Fig. 1) were measured a posteriori using a profilometer and an Axioplan microscope (Carl Zeiss AG, Jena, Germany).

Fresh venous blood (10 ml) from specific pathogen-free Large White pigs was provided by the Institute of Virology and Immunology (IVI) in Mittelhäusern, Switzerland. The 10 ml blood samples were obtained exclusively in the context of 300–400 ml blood drawings from adult pigs performed for the purpose of in vitro studies at the IVI under the animal licenses BE 88/14 and BE131/17 delivered by the veterinary authorities of the Canton of Bern. The samples were centrifuged at $1800×g$ for $10min$ to separate cellular components and plasma. The supernatant was discarded, and 0.5 ml of RBCs was collected. A solution of Bovine Serum Albumin (BSA) (Sigma-Aldrich, St Louis, MO, USA) dissolved at $1%$ in Phosphate-Buffered Saline (PBS) (Sigma-Aldrich, St Louis, MO, USA) was used to wash RBCs to reduce echinocytosis.⁴⁵ After washing, the RBC samples were centrifuged for 5 min at $2500×g$⁠.

The suspending medium (hereafter referred to as plasma) was prepared to match the density of the RBCs to avoid sedimentation. As suggested by Roman et al.,⁴⁷ it contained 65% Glucose-Albumin-Sodium-Phosphate (GASP) buffer (PBS with 5.5 mM glucose and $4%$ BSA) with $35%$ of stock solution (⁠$90%$ Optiprep (Sigma-Aldrich, St Louis, MO, USA) + 10% GASP 10 times concentrated), resulting in a density of $1090kg/m3$ and a viscosity of $1.96×10−3Pas$ at $20°C$⁠, which is in the physiological range of plasma.¹⁷ RBCs were diluted in the medium such that the hematocrit was $10%$⁠. We chose a reservoir hematocrit $Hf=10%$⁠, which is within the physiological range for microcirculation.^24,40 Due to the centrifugation stages, no other biological components but RBCs were expected to be in the final suspension. Experiments were carried out at ambient temperature (⁠$T=21°C$⁠) within 12 h after the blood collection to ensure a healthy state of the RBCs for the whole duration of the experiment.

The microchannel network was degassed and prefilled as described by Clavica et al.¹⁰ The blood sample was placed in a 10 ml centrifuge tube, was mounted on a vertical linear-motion stage, and connected to the inlet of the microdevice. We performed experiments imposing a hydrostatic pressure difference between the fluid level in the reservoir and the device outlet (Fig. 2). Two different driving pressures were used: $ΔP1=35.3mbar$ and $ΔP2=47.1mbar$ (hereafter referred to as low velocity and high velocity experiments, respectively). Pressure-driven flows have short response times when the pressure is changed and they are stable over long periods, making the steady-state study of phase separation possible.⁴⁷ Attention was paid to maintain the pressure difference constant throughout the experiments such that pressure variations due to changing fluid levels in the reservoir were negligible. The microdevice was placed on an inverted microscope (Nikon, Japan) with a $10×$ air objective (⁠$NumericalAperture=0.45$⁠) and illuminated by a LED lamp. A $XYZ$ micrometer stage was used to focus the light on a region of interest (ROI) of $512×512pixels$ (Fig. 1). For each experiment, the hydrostatic pressure difference was set first. Once a steady-state condition was reached, videos of 25 s were recorded at 120 frames/s using a high-speed camera (ORCA-flash 4.0, Hamamatsu, Japan).

The video sequences were imported and preprocessed in Matlab (Mathworks, Natick, MA, USA) to correct for background illumination differences, remove the background, and convert the color from gray scale to black and white. The processed frames were imported in an open-source program for Particle-Tracking Velocimetry (PTV)⁸ (PTVlab, Mathworks, Natick, MA, USA). The output data from PTVlab were postprocessed with Matlab to automatically derive the following quantities for each frame: average RBC velocity [mm/s], RBC line density $[RBCs/m]$ and tube hematocrit (⁠$Ht$⁠, i.e., the ratio of the total RBC volume to channel volume). The average RBC velocity is the mean of the velocities of all RBCs identified in the region of interest (ROI, Fig. 3) for a given frame. The RBC line density is the number of red cells present in the ROI for a given frame divided by the length of the channel segment included in the ROI. The raw data for the number of RBCs within the ROI were first smoothed with a second-order Savitzky-Golay filter in Matlab (Mathworks, Natick, MA, USA), using a window of 199 points corresponding to intervals of 1.66 s. These data were then used to compute RBC flux and blood flow rate (see Sec. II D 2). The PTV algorithm was already compared in a previous study¹⁰ with the line scan method (LSM)²⁷ to ensure accurate results. In the present study, we performed a more detailed validation (see Sec. III A). To this end, line scans were sampled from the microchannel and stacked in a temporal sequence. The RBCs were identified from line traces in these stacks and their velocity was computed from the slope of these line traces. Additionally, the data from the particle-tracking algorithm were checked by verifying the mass balance for the RBC flux at bifurcations.

The technique for the measurement of the tube hematocrit $Ht$ relied on the identification of the RBCs by PTVlab. The discharged hematocrit in the reservoir was $Hd=10%$ and, therefore, the tube hematocrit $Ht$ was $10%$ or less because of the Fåhraeus effect.¹³ This low hematocrit allowed counting of individual RBCs.⁴⁷ The tube hematocrit for each frame could then be calculated as

where $Nrbc$ is the number of RBCs in the ROI, $Vchannel$ is the volume of the microchannel segment in the ROI (⁠$Vchannel=LROI×WROI×H$⁠), and $MCVrbc$ is the RBC mean corpuscular volume. The average volume of a single porcine RBC^2,56 is $MCVrbc=56μm3$⁠. The mean hematocrit was computed as the temporal average of the tube hematocrit.

Data on the number and position of RBCs were used to compute line density profiles in the parent vessels of four diverging bifurcations (see Sec. III C 2) in the function of the lateral coordinate (y-coordinate). To this end, the ROI was divided parallel to the axial direction into 18 subregions (each 1 pixel wide). We assigned each RBC to one of those regions depending on its centroid position determined by particle tracking. This allowed us to count the number of RBCs at a given lateral position in the microchannel. Line density profiles were computed by dividing the mean RBC number at each lateral position by the length of the ROI [Fig. 9(b)], where the lateral coordinate was normalized by the width of the ROI (⁠$y∗=y/WROI$⁠). To quantify the skewness of the line density profiles, the hematocrit skewness index⁵³ $Sh$ for all divergent bifurcations considered in the study was computed as

where $Sh=0$ for a perfectly symmetric hematocrit profile along the width of the microchannel, while $Sh=0.5$ for a completely skewed hematocrit profile.

For further validation, the tube hematocrit measurements obtained from particle tracking were compared to measurements obtained from image sequences, which were postprocessed according to the photometric method by Roman et al.⁴⁷ More information on this validation can be found in  Appendix B.

Images of RBCs flowing in the microfluidic network were used to calculate the velocity in the parent and in the daughter branches of the diverging bifurcations. The PTV algorithm allowed to compute the velocity $Urbc(x,y)$ in the plane $XY$ for all the RBCs identified in the region of interest. The mean RBC velocity (⁠$Urbc$⁠) is the mean of the velocities of all RBCs identified in the ROI for a given frame averaged over all frames included in the image sequence. Assuming a mean velocity $Urbc$ in the ROI (⁠$width=WROI$⁠, $length=LROI$⁠, $height=H$⁠), we compute the mean RBC flux as

Starting from the mean RBC velocity $Urbc$⁠, the total blood flow rate could be calculated according to

where $χ$ is a coefficient taking into account the velocity difference between the plasma and RBCs. The plasma velocity is typically smaller than the RBC velocity due to the Fåhraeus effect.¹³ However, we assumed $χ=1$ (similar to Sherwood et al.⁵³) implying that the mean velocity of the total blood is equal to the mean RBC velocity. We will see below that one of the main results of this study (phase separation diagram, Fig. 8) is independent of the choice of $χ$⁠.

The RBC flux and blood flow rate data were corrected according to the procedure introduced by Pries et al.⁴¹ and used in other experimental studies.^47,53 This approach allowed to ensure mass conservation at each bifurcation, to reduce scatter in the experimental data and the uncertainty in the fitting.^41,53 The corrected flow values (marked by hats) for the generic uncorrected flow rates $QP$⁠, $Q1$⁠, and $Q2$ were calculated as

where the weighting factors were given by

The fractional RBC flux (⁠$Ψi$⁠) and the fractional blood flow (⁠$Φi$⁠) in the daughter branch $i$ are given by

where $i=1,2$ denotes the daughter branch 1 and 2 of a generic diverging bifurcation (Fig. 3). Note that the coefficient $χ$ defined in Eq. (4) cancels in the definition of $Φi$⁠, such that the phase separation results are independent of the choice of $χ$⁠.

The measured values for $Ψi$ and $Φi$ were compared with the phase separation law proposed by Pries et al.,⁴² which is defined as

where $A$⁠, $B$⁠, and $X0$ are dimensionless parameters. $A$ describes the asymmetry between the daughter branches and is given by $A=−15.47×[((D1h)2−(D2h)2)/((D1h)2+(D2h)2)]×[(1−Hd)/DPh]$⁠. In this study, $A=0$ because the bifurcations are symmetric (⁠$D1h=D2h$⁠). $B$ denotes the sigmoidal shape of the fitting function and is expressed by $B=1+8.13×(1−Hd)/DPh$⁠. $X0$ defines the minimal blood flow fraction $Φi$⁠, which is required to draw RBC into the daughter branch. It is computed as $X0=1.12×(1−Hd)/DPh$⁠. The numerical values for the parameters $A$⁠, $B$⁠, and $X0$ were finally multiplied by a correction factor⁴⁷ $MCVratio=56/913=0.85$ to account for the volume difference between porcine and human RBCs.

Next to the definitions for $B$ and $X0$ given by Pries et al.,⁴² we tried to find the best set of parameters (⁠$B^$ and $X^0$⁠) by minimizing the Root Mean Square Error (RMSE) between the present in vitro results and the model of Pries et al.⁴² [Eq. (13)]. To this end, two parameters $α$ and $β$ were introduced as

To analyze the blood flow heterogeneity at individual bifurcations, the relative velocity difference between the daughter branches was computed according to Schmid et al.⁴⁹ as

where $Urbc1$ and $Urbc2$ are the velocities in the daughter vessels of a diverging bifurcation.

Likewise, we computed the average balance factor $Bdb¯$ for diverging bifurcations introduced by Schmid et al.⁴⁸ according to

where $Φk$ is the blood flow rate fraction for one of the daughter branches of the $k$th bifurcation and $n$ is the total number of diverging bifurcations. The average balance factor $Bdb¯$ ranges from 0 for a completely unbalanced network to 1 for a perfectly balanced network. A network with $Bdb¯>0.8$ was considered well balanced.⁴⁸ 

There were four channels (Fig. 3, blue arrows) that carried no RBCs at all. Therefore, they required a different method for the computation of the flow rate. To this end, we estimated the flow velocity $X$ in the RBC-free branch (Fig. 4) from mass balance equations at the adjacent bifurcations. To reduce the effect of measurement errors, the estimates from two neighboring bifurcations were averaged. This resulted in

where $U¯A,U¯B,U¯C,$ and $U¯D$ represent the mean RBC velocities according to Fig. 4 (data are reported in  Appendix A). Note that the mass balance is formulated here with velocities because the cross sections are equal for all branches of the network.

First, velocity measurements are validated against the Line Scan Method²⁷ (LSM) and the mass balance is verified at the bifurcations (Sec. III A). Second, results on the hematocrit distribution in the microfluidic network (Sec. III B) are shown. Finally, phase separation data for 21 diverging bifurcations and an analysis for the minimization of the root mean square error between present results and literature models are presented (Sec. III C). Because the microchannel network is symmetric with respect to the x axis, the results are only reported and discussed for to the upper half of the microfluidic device.

The measured RBC velocities varied between 0.027 mm/s and 0.400 mm/s for the low velocity experiment. For the high velocity experiment they ranged between 0.051 mm/s and 1.726 mm/s. We validated these measurements against the Line Scan Method (LSM), which can be considered a reliable method in the limit of low hematocrit,⁴⁷ because the line trace associated to each RBC is clearly discernible and the computation of the RBC velocity as the slope of this line is a straightforward process.²⁷ Our PTVlab results were validated by considering a sample of 11 parent branches from the high velocity experiment (⁠$ΔP2=47.1mbar$⁠) because it represented the worst case scenario in which the resolution of the PTVlab data was the potential limiting factor for the correct quantification of RBC velocity. The results presented in Fig. 5(a) show good agreement between the methods. The average of the relative difference between PTVlab and LSM is $7.4%±4.4$ (⁠$mean±SD$⁠). This justified the choice of the particle-tracking method for the RBC velocity for both experiments (⁠$ΔP1=35.3mbar$ and $ΔP2=47.1mbar$⁠).

The method for the velocity computation was additionally tested by verifying that the presented data were consistent with the RBC flux mass balance at diverging bifurcations. For all diverging bifurcations with RBCs flowing in both daughter branches, the RBC flux in the parent branch was compared with the sum of the flux measured in the daughter branches. Results in Fig. 5(b) show good agreement between inflow and outflow at diverging bifurcations. The presented results indicate that the velocity measurements were accurate and respected mass conservation.

Following Roman et al.,⁴⁷ an apparent deviation from mass conservation was computed for each diverging bifurcation as

where $QrbcP$ is the RBC flux in the parent branch and $Qrbc1$ and $Qrbc2$ are the RBC fluxes in the daughter branches. The mean deviation for the low velocity experiment is $E¯=7.18%±3.52$ (⁠$mean±SD$⁠, range 1.65%–12.90%). The mean deviation for the high velocity experiment is $E¯=7.51%±4.28$ (⁠$mean±SD$⁠, range 1.61%–14.64%). These results are in the range of other experimental studies^47,53 with a similar approach for the RBC flux computation. The deviation from mass conservation inherent to the experimental nature of the study explained a posteriori the need of correction factors (⁠$WP$⁠, $W1$⁠, and $W2$⁠, Sec. II D 2) introduced by Pries et al.⁴¹ for the correct quantification of phase separation at diverging bifurcations. The mean deviations computed in this study are in the same range as the results reported by Pries et al.⁴¹ 

The distribution of RBCs and the flow field were analyzed for both experiments, which provided the basis for the investigation of the phase separation. The first step was the RBC tracking to identify preferential pathways within the network. Even if the number of RBCs in microchannels was time-dependent, the direction of flow was always constant and no inversion of flow direction was observed during the experiments. Figure 3 shows the flow directions for each branch of the network.

The second step was the conversion of the time-averaged number of RBCs to the tube hematocrit (Sec. II D 1). The hematocrit distribution is heterogeneous in both experiments: there is an accumulation of RBCs in the center of the network, while there is a RBC reduction in the corners. Figure 6 shows the time-averaged normalized pixel intensity of the video frames for the low velocity experiment (a) and for the high velocity experiment (b), which provides a direct visual impression of the heterogeneous hematocrit distribution in the network with preferred RBC pathways and RBC-free branches.

The hematocrit distributions expressed in function of the y-coordinate [Fig. 7(a)] are similar for both experiments. The distributions show that the main RBC accumulation did not occur close to the axis of symmetry, but in the central region (4 ≤ y ≤ 5) of the half-network.

Similar to the hematocrit distribution, the RBC flux [Fig. 7(b)] is also heterogeneous. Apart from the expected result that the mean RBC flux increases with the perfusion pressure and the flow velocities, the RBC flux is higher close to the axis of symmetry and decreases toward the periphery of the network. A RBC flux reduction of $80.2%$ and $87.1%$ was computed from the center (⁠$y=1$⁠) to the periphery of the network (⁠$y=8$⁠) for the low velocity and high velocity experiment, respectively.

In order to study the phase separation, the RBC flux fraction $Ψi$ in each daughter branch of the diverging bifurcations was plotted in the function of the respective blood flow fraction $Φi$ (Fig. 8). The results were compared with Pries’ separation law⁴¹ (empirical fitting of in vivo data with a sigmoid function).

Pries’ law appears to agree well for the experimental results obtained for both perfusion pressures. The majority of the data points follows the Zweifach-Fung effect and indicates that branches with lower blood flow fraction (⁠$Φi<0.5$⁠) receive even less RBC flux (⁠$Ψi<Φi$⁠) and vice versa. Exceptions to this general rule are represented by data in the blue shaded regions in Fig. 8. In this case, the bifurcation behavior is reversed: the branch with a lower blood flow fraction (⁠$Φi<0.5$⁠) is receiving more RBCs (⁠$Ψi>Φi$⁠) and vice versa. This phenomenon, known as inversion of the Zweifach-Fung effect^10,52 or reverse partitioning,⁵ has been previously reported in experimental studies for single asymmetric bifurcations¹⁰ and simplified rectangular networks.⁵² To the best of our knowledge, it is reported here for the first time for a larger and more realistic in vitro network. The experimental data shown here are in good agreement with the numerical results reported by Balogh and Bagchi⁵ on RBC partitioning in simulated microvascular networks with in vivo like features.

Data points in the gray regions in Fig. 8 represent diverging bifurcations with one RBC-free daughter branch. Following the procedure explained in Sec. II D 3, the plasma flow rate was computed in the RBC-free branch (data for all variables of interest are reported in  Appendix A). This result indicates that there exists a flow threshold below which no RBCs are entering a branch. This threshold corresponds to the parameter $X0$ in Pries’ law. For both experiments, a flow threshold $X0=0.19$ was observed, which means that for a blood flow fraction $Φi<0.19$⁠, we found a RBC flux fraction $Ψi=0$ and, correspondingly, $Ψi=1$ for $Φi>0.81$⁠.

The local phase separation at diverging bifurcations has consequences on the flow velocities at the global scale of the whole network. Table I shows that for the low velocity experiment, $41.2%$ of the bifurcations had a velocity difference $|ΔrUdb|<40%$ while only $11.8%$ of the bifurcations had $|ΔrUdb|>80%$ (Sec. II D 2). For the high velocity experiment, the blood flow distribution was more heterogeneous: $35.3%$ of bifurcations had a $|ΔrUdb|>80%$⁠. These results were supported by the averaged balance factor $Bdb¯$ [Eq. (16)] for the whole network. $Bdb¯=0.79$ and $Bdb¯=0.70$ were found for the low velocity experiment and for the high velocity experiment, respectively. These results indicate that the network is less balanced at higher velocities.

The path along the bifurcations I, II, III, and IV indicated in Fig. 3 and in Fig. 6 features a series of diverging and converging bifurcations. Figure 9(a) presents data extracted from Fig. 8 corresponding to the diverging bifurcations I, II, III, and IV. For the shown results, the data were sampled from the daughter branch, which is closer to the upper half of the parent vessel [$0<y∗<0.5$⁠, gray region in the inset of Fig. 9(a)].

In general, the phase separation behavior for bifurcations closer to the inlet of the microdevice (Bif. I and II) was reasonably well predicted by Pries’ separation law for both experiments. For the more distal bifurcations III and IV the deviation from Pries’ law increased and we observed an inversion of the Zweifach-Fung effect. This effect was stronger for the high velocity experiment.

A correlation was found between the phase separation results for the bifurcations I, II, III, and IV and the local hematocrit distribution in their respective parent vessels. Figure 9(b) shows that the normalized line density profiles underwent significant shape changes from the first to the last bifurcation (see also  Appendix B). For the low velocity experiment, the profile in the parent vessel of bifurcations I was symmetric and almost flat while it was skewed with asymmetric peaks for bifurcations III and IV. Results from the high velocity experiment showed a similar shape transition with an even more skewed profile at the last bifurcation.

To quantitatively demonstrate the relation between the skewness and RBC partitioning, all divergent bifurcations were sorted based on the type of partitioning they presented (i.e., classical or reversed) and the mean hematocrit skewness index (⁠$S¯h$⁠) was computed for each of the two groups. We found a $S¯h,classical=0.11±0.1$ and $S¯h,reverse=0.21±0.1$ for the bifurcations with classical and reverse partitioning, respectively. The Wilcoxon rank sum test was performed to compare the hematocrit skewness index for the two groups of bifurcations (significance for p $<0.05$⁠). $p=0.014$ was obtained, meaning that the bifurcations with reverse partitioning have a statistically higher hematocrit skewness index.

In conclusion, for flat and symmetric line density profiles the RBC partitioning is more likely to follow the Zweifach-Fung effect and Pries’ separation law. For skewed line density profiles, the RBCs rather adhere to the streamlines and the phase separation is less sensitive to blood flow fraction differences.^6,51 As a result, reverse partitioning is more likely.

The phase separation model proposed by Pries et al.⁴¹ was based on the empirical fitting of in vivo data from arterial bifurcations in which the flow rate of the daughter branches was manually imposed by occluding a downstream vessel. In this study, the flow rates of daughter branches at diverging bifurcations was a natural consequence of the downstream microchannel resistances and the discrete contributions to resistance due to RBCs.

A quantification of the Root Mean Square Error (RMSE) was performed in order to understand to which extent the phase separation law of Pries⁴¹ is able to fit the present experimental data and to find the best set of parameters (⁠$X^0$ and $B^$⁠) to describe the in vitro data. According to the methodology presented in Sec. II D 2, colormaps were produced to quantify the RMSE in the function of two parameters $α$ and $β$ (Fig. 10). The RMSE for the original Pries’ curve (⁠$α=β=1$⁠) is $RMSElow=0.101$ and $RMSEhigh=0.166$ for the low and high velocity experiment, respectively.

The RMSE was also quantified for an updated version of the Pries’ separation law proposed by Rasmussen et al.⁴³ By exploiting a Bayesian data analysis framework (posthoc analysis), a new set of parameters (⁠$X0$ and $B$⁠) was proposed to better fit in vivo data with the well-known sigmoid function of Pries. Using the conventions of this study, the set of parameters introduced by Rasmussen et al.⁴³ corresponds to $X^0=α×X0$ and $B^=β×B$ with $α=0.08$ and $β=0.9$ resulting in $RMSElow=0.110$ and $RMSEhigh=0.171$⁠.

In a recent review, the need for more experimental data to validate computational models of RBC flow in capillary networks¹⁶ was highlighted. To help fill this gap, the present study featured a microvascular network model with complex yet idealized topology with microchannels $10μm$ wide or less. The impact of RBCs on the global dynamics of the network as well as on the local hemodynamics at diverging bifurcations was investigated. To the best of our knowledge, in vitro quantitative phase separation data from such networks have not been reported yet.

On the scale of the whole network, the hematocrit (Fig. 6) and the RBC flux distribution [Fig. 7(b)] showed that the RBC flow was heterogeneously distributed throughout the network and preferential pathways were present. The basic mechanisms leading to preferential pathways are not specific to this network and they are a consequence of the effective hydraulic resistance of the microchannels and the phase separation law.³⁹ It should be highlighted that these pathways developed freely in an in vitro microvascular network model and they were not a consequence of locally imposed boundary conditions on single bifurcations. The results on the relative velocity difference (Table I) and the average balance factor $Bdb¯$ of 0.79 (low velocity) and 0.70 (high velocity) (Sec. III C 1) suggest that the self-regulation process described in Sec. I was less effective for higher velocities. In their numerical study with a hexagonal network similar to ours, Schmid et al.⁴⁸ reported that the average balance factor ranged between $0.7<Bdb¯<0.8$ for an inlet tube hematocrit $Ht,inlet<0.10$⁠, which is comparable to our experimental data. The maximum in the balance factor (⁠$Bdb¯≈0.84$⁠) was obtained for $Ht,inlet=0.25$⁠. As suggested by Schmid et al.,⁴⁸ a sufficient number of RBCs in the network is needed to balance the downstream bifurcations. Considering that the maximum inlet tube hematocrit calculated in the present experimental study was $Ht,inlet<0.10$⁠, it is reasonable to believe that the efficiency of self-regulation mechanisms was limited.

We found that the RBC splitting at diverging bifurcations is not proportional to the blood flow fraction in the daughter branches, which agrees with the consensus in the literature that phase separation at diverging bifurcations is dominated by the Zweifach-Fung effect.^{5,41,47,52,53} The data from the low velocity experiment can be well described by the empirical fitting of in vivo measurements proposed by Pries et al.,⁴¹ the RMSE is consistent with similar studies found in the literature.⁴⁷ The data from the high velocity experiment are in reasonable agreement with Pries’ separation law for $Φ<0.4$ or $Φ>0.6$ (see Fig. 8). A deviation is present for $0.4<Φ<0.6$ that might originate from the inversion of Zweifach-Fung effect^10,52 and the history effect³⁴ as explained further below. Results suggest that at most branching points the majority of RBCs is more likely to enter the high-flow branch (classical partitioning). In some bifurcations, RBCs preferred to enter the low-flow branch (blue shaded regions of Fig. 8). This phenomenon is a representation of the reverse partitioning (i.e., inversion of the Zweifach-Fung effect) and we found it to be more likely for higher flow velocities. This dependence on the flow velocity is consistent with the findings of an in vitro study with asymmetric single bifurcations.¹⁰ 

A complete description of phase separation is not trivial and in the limit of low hematocrit, two competing phenomena affect the RBC partitioning: either RBCs enter the branch with a higher flow rate or they follow the streamlines.⁵⁰ It was demonstrated in a computational model⁶ that RBCs flowing along the centerline of the channel adhere to the Zweifach-Fung effect, while RBCs located closer to the walls more likely follow the streamlines of the plasma flow. In our experiments, both phenomena were likely to be present because the particle to channel diameter ratio was equal to 0.6, meaning cell margination or cell lateral focusing were possible.

A path featuring multiple diverging bifurcations was considered in Fig. 9 and the relation between phase separation and RBC spatial distribution in the parent vessel was investigated. We found that the RBCs are more likely to split according to the Zweifach-Fung effect if their distribution in the parent vessel is symmetric and homogeneous [Fig. 9(b), Bif. I]. By contrast, if the RBC distribution is skewed in the parent vessel [Fig. 9(b), Bif. III] or if the RBCs are marginated [Fig. 9(b), Bif. IV], it is more likely that they follow the streamlines and enter the closer branch. The latter situation supports the hypothesis¹⁰ that a depletion in the central region of the parent vessel may affect the phase separation by reducing the number of RBCs entering the high-flow branch (i.e., reducing the Zweifach-Fung effect). As an example, the upper daughter branch of bifurcation IV [blue rectangle in Fig. 9(a)] received a higher RBC flux fraction than the blood flow fraction because in the parent vessel, the RBC line density was much higher on the upper side of the parent vessel [$0<y∗<0.5$⁠, Fig. 9(b)]. The peculiar shape of the line density profiles for bifurcations III and IV has to be attributed to the fact that the RBC distributions in these parent vessels are in turn the result of two separate RBC distributions from converging branches. The skewed profile in the parent of bifurcation III is the result of the confluence of a channel with very low hematocrit and a channel with very high hematocrit [Fig. 6(c)]. The bimodal distribution in the parent of bifurcation IV is probably due to a lack of mixing of the converging RBC fluxes at the inflow to this vessel [Fig. 6(d)]. This interpretation is consistent with numerical and experimental results reported for a Y-shaped converging bifurcation, where a cell-depleted layer of RBC was found in the center of a microchannel after a confluence.^3,29,34 Moreover, we observed a statistically significant higher hematocrit skewness index for the reverse partitioning and a greater deviation from Pries’ law (i.e., reverse partitioning) for the high velocity. Thus, the hematocrit skewness is higher in the high velocity experiment. The local mechanisms that may explain the origin of the skewness and its velocity dependence are probably related to the dynamics of RBCs, which is governed by the complex interaction between lift force, shear rate gradient experienced by the cell and, in the case of concentrated RBC suspension, nonlinear anisotropic shear-induced diffusion.^19,32,37 It is the resultant of these mechanisms at RBC scale that may lead to skewed hematocrit profiles.

In vitro studies of phase separation are mostly performed in single bifurcations in which the velocity and hematocrit profiles in the parent vessel are fully developed and symmetric.^{12,14,26,47,53} Unfortunately, this configuration may not necessarily be a good representation of the hemodynamics in vivo. In the case of a complex network with diverging/converging bifurcations, the hematocrit profile disturbed by the flow field at a first bifurcation may not be able to regain symmetry before the next bifurcation.^9,32,53 The length required to regain complete symmetry was estimated to be even greater than 11 diameters for disturbed flow conditions.³⁸ The low radial RBC dispersivity⁵³ and the short time a RBC spends within the channel provide evidence that the flow history can be important if the bifurcations are close to each other or if the flow rates are sufficiently high.⁹ The possibility to compute RBC line density profiles along a path with multiple diverging bifurcations allowed us to study this phenomenon and it indicated that a history effect, as described by Merlo,³⁴ is indeed present in our model of a microvascular network. A skewed RBC distribution across the parent vessel does not have enough time to regain symmetry and leads to reversed cell partitioning, which is unlikely to be seen in a model of single bifurcation with homogeneous RBC distribution in the parent vessel. A detailed analysis can be found in  Appendix C to demonstrate that an equilibration of the hematocrit profile between two successive bifurcations was unlikely to be achieved in the present microchannel network.

The agreement of the present in vitro results on the phase separation with the computational results of Balogh and Bagchi^4,5 is striking. Both studies highlighted the simultaneous occurrence of classical and reverse partitioning in a microvascular network and related these phenomena to the upstream hematocrit distribution and to the influence of sequential bifurcations. Our experimental approach allowed to overcome a limitation of their numerical study, namely, the modeling of the RBC membrane. In their model, Balogh and Bagchi did not consider the membrane viscosity, which has an impact on the time scale of shape recovery when the RBCs are deformed. This may have a negative consequence on the study of cellular mechanisms at the bifurcations. However, the consistency between experiments and simulations resolved this potential issue and validates the numerical model of Balogh and Bagchi.^4,5

In vitro experiments involve simplifications with respect to in vivo investigations. Capillaries are cylindrical and nonuniform over their length while in our study we had straight microfluidic channels with rectangular cross sections. Our aim was to produce microchannels with a characteristic width $≤10μm$ to mimic the anatomical dimensions of the capillaries. Due to technological limitations, fabricating rectangular channels was the only available option (networks of microchannels with circular cross section are limited today to channel size $30−40μm$⁠).^10,30 We are aware that rectangular cross sections may alter the local fluid dynamics: RBCs do not occupy the whole lumen of the channels leaving open regions at the corners for the plasma to flow. We tried to limit these effects by reproducing the confinement that RBCs experiences in the capillaries using channel sizes slightly bigger than the RBCs. Further, the surface properties of the internal walls of capillaries were not reproduced. In vivo, there is an endothelial surface layer with glycocalyx while in vitro, there is only a nonspecific adsorption of bovine serum albumin.⁴⁷ The good agreement between our in vitro results and the empirical fitting of in vivo data indicated that these simplifications are nondominant determinants of phase separation as long as averaged variables are considered and cell-cell interactions are weak.⁴⁷ 

Another limitation is related to the blood mimicking fluid used in the experiments. In the RBC solution, there were no other biological components than RBCs. This prevented the study of the interactions between blood components that may affect the cell screening and phase separation at bifurcations in vivo.

As reported by Roman et al.,⁴⁷ the quite large experimental variability has to be accounted as another limitation. The experimental variability might be attributed to technical difficulties experienced with such small microchannels. For the same reason, we chose to report experiments with only a reservoir hematocrit $Hd=10%$ because it was challenging (if not impossible) to perfuse the network in a controlled manner with more concentrated RBC suspensions. In fact, microchannel obstructions were more likely to happen at high hematocrit determining a perturbation in the network and biasing the study of the phase separation.

In the present study, we developed an in vitro microvascular network model to study the RBC phase separation at diverging bifurcations. For the first time, we provided quantitative data on the RBC distribution and separation in a complex in vitro network model. These data confirmed computational results on RBC partitioning in realistic microvascular networks.⁵ 

The data reproduced the classical result that branches receiving the higher blood flow fractions typically received an overproportional RBC flux fraction. An inversion of this classical phase separation behavior was observed in the case of a skewed hematocrit profile in the parent vessel of the bifurcations, especially for high-flow rates. This provided evidence that the flow history has to be considered to properly describe the RBC phase separation in complex networks with multiple diverging bifurcations. Finally, we were able to identify a critical blood flow fraction of $Φ=0.19$ below which only plasma enters a downstream branch.

In conclusion, the nonuniform RBC distribution and partitioning underline the importance of treating blood as a biphasic fluid and the need of considering the discrete and particulate nature of blood to properly describe fluid dynamics at the microscale, especially in microvascular networks where the RBC dynamics becomes complex and often counterintuitive.³⁹ 

This study has been financially supported by the Swiss National Science Foundation under Grant No. CR2312_166707. We thank the work of Alexandra Homsy and Laure Jeandupeux (Haute Ecole Arc Ingénierie, University of Applied Sciences Western Switzerland, La Chaux-de-Fonds, Switzerland) for the fabrication of the master for the microfluidic device. We gratefully acknowledge the Institute of Virology and Immunology (IVI, University of Bern) for providing us the blood samples. We acknowledge the Organs-on-chip Technologies Group (ARTORG Center, University of Bern) for sharing their laboratory equipment.

Considering a pair of bifurcations sharing a RBC-free channel, the plasma flow was computed by mass conservation as

where $U¯A,U¯B,U¯C,$ and $U¯D$ represent the mean RBC velocities according to Fig. 4. The raw data from the particle-tracking velocimetry that have been utilized to compute the mass conservation are reported in Tables II and III for the low and high velocity experiments, respectively. These data have been ultimately used to compute the blood flow fraction and the RBC flux fraction in the resulting diverging bifurcations with a RBC-free daughter branch. During the postprocessing of the image sequences, we noted that a spurious white blood cell flowed through the bifurcation named “Q4, Bif. III” during the high velocity experiment. For this reason, the velocity data were biased and this bifurcation was excluded from further analysis.

A photometric method has been reported by Roman et al.⁴⁷ to compute the tube hematocrit for high concentrations of RBCs. The same method has been applied here to qualitatively validate the line density profiles presented in Sec. III C 2. For each frame of the acquisition, the pixel intensity along the lateral direction $Y$ at a fixed longitudinal position $xref$ was computed [$I(xref,y)$]. The mean pixel intensity was obtained as the time-averaged intensity over all the frames of the acquisition. As a reference, the pixel intensity at the same longitudinal position was computed with only plasma in the channel [$I0(xref,y)$]. The light attenuation at any point (⁠$xref,y$⁠) was defined as the ratio of $I0(xref,y)/I(xref,y)$⁠. The lateral optical density profile (⁠$OD(xref)$⁠) was obtained by averaging along the lateral direction $Y$ as follows:

where $⟨⟩y$ denotes the averaging along the transversal direction. Figure 11 shows the optical density profiles obtained for the bifurcations I, II, III, and IV for the low velocity experiment.

The hematocrit is linearly proportional to the optical density; therefore, the higher the OD, the higher the tube hematocrit at a specific lateral position. The optical density profiles are in reasonable qualitative agreement with the line density profiles presented in Fig. 9(b). These results demonstrate the validity of the method developed to obtain the line density profiles and the equivalence to another method found in the literature for the hematocrit computation.^34,47

A detailed analysis is performed to demonstrate that the “history effect” was present in the capillary network and the relaxation of the hematocrit profile between two successive bifurcations was not likely to be achieved due to the short length of the microchannels.

First, velocimetry data obtained with the particle tracking can be useful to quantitatively assess the contribution of the lateral RBC velocity (i.e., the velocity component in the $y$-direction) to the magnitude of the RBC velocity. Figure 12 shows the ratio between the lateral RBC velocity and the magnitude of the RBC velocity $Uy,rbc/Urbc$ for all diverging bifurcations considered in our study as a function of the magnitude of the RBC velocity. The median ratio is $5.4×10−3$ and $5.6×10−3$ for the low and high velocity experiments, respectively. Given that the RBC velocity in the $Uy,rbc$ is on average three orders of magnitude smaller than $Urbc$⁠, it is reasonable affirm that the displacement in the lateral direction is very limited.

Second, a dimensional analysis can provide insight on the contribution of the RBC advection in the main flow direction with respect to the RBC dispersion in the lateral direction. Considering the length $L$ of a microchannel and the mean RBC velocity $Urbc,mean$ determined by particle-tracking, the advection time scale (i.e., the time for a RBC to pass through a channel of length $L$⁠) is given by

Using the value of radial dispersivity reported by Cokelet¹¹ (⁠$D=10−7mm2/s$⁠), it is possible to compute the dispersion time scale (i.e., the time required by RBCs to disperse over the width $W$ of the microchannel) as:

The fact that the related Peclet number $Pe=Tadv/Tdisp≪1$ (which holds even if we had considered the minimum RBC velocity) indicates that the RBC flow is dominated by advection and the radial dispersion is of minor relevance in this configuration.

Third, by looking at single-cell migration models,³² we estimated the typical time a RBC requires to migrate from the wall (⁠$y∗=W/2$⁠) through the Cell-Free Layer (CFL) (⁠$y∗=(W/2)−d$⁠) as

Considering the optimized parameters provided by Losserand et al.³² and the geometrical dimensions of our microchannels (⁠$δS=1.3$⁠, $ξS=1.1×10−2$⁠, $R0=2.37μm$ for porcine RBCs, $d=2μm$⁠, $Urbc,max=1.726mm/s$⁠, $W/2=5μm$⁠), we obtained $τL=0.178s$⁠.

Similarly, we estimated the equilibration time due to the shear-induced diffusion as

Considering $f=0.06$ from Grandchamp et al.¹⁹ and $Hd=10%$⁠, we obtained $τD=2.149s$⁠.

Following Losserand et al.,³² the length a RBC requires to pass through the CFL is $τL×Urbc,max≈307.2μm$ and the one to reach equilibration is $τD×Urbc,max≈3.72mm$⁠, which is longer than the length of our microchannels (⁠$L=85μm$⁠). This is consistent with the results provided by Losserand et al.³² where lengths of $292.5μm$ and $930μm$ were needed to study the RBC transversal migration in microchannels $3.9μm$ and $6.2μm$ wide, respectively. Therefore, the hematocrit distribution relaxation is not likely to happen over the short distance between two successive bifurcations and an history effect is present in the microchannel network.