Red blood cell lingering modulates hematocrit distribution in the microcirculation

Article Red blood cell lingering modulates hematocrit distribution in the microcirculation Yazdan Rashidi, 1, *Greta Simionato, 1,2 Qi Zhou, 3 Thomas John, 1 Alexander Kihm, 1 Mohammed Bendaoud, 1,4,5 Timm Kr€ uger, 3 Miguel O. Bernabeu, 6,7 Lars Kaestner, 1,8 Matthias W. Laschke, 2 Michael D. Menger, 2 Christian Wagner, 1,9 and Alexis Darras 1, * 1 Experimental Physics, Saarland University, Saarbruecken, Germany; 2 Institute for Clinical and Experimental Surgery, Saarland University, Homburg, Germany; 3 School of Engineering, Institute for Multiscale Thermofluids, University of Edinburgh, Edinburgh, United Kingdom; 4 Universit e Grenoble Alpes, CNRS, LIPhy, Grenoble, France; 5 LaMCScI, Faculty of Sciences, Mohammed V University of Rabat, Rabat, Morocco; 6 Centre for Medical Informatics, Usher Institute, University of Edinburgh, Edinburgh, United Kingdom; 7 The Bayes Centre, University of Edinburgh, Edinburgh, United Kingdom; 8 Theoretical Medicine and Biosciences, Saarland University, Homburg, Germany; and 9 Physics and Materials Science Research Unit, University of Luxembourg, Luxembourg, Luxembourg ABSTRACT The distribution of red blood cells (RBCs) in the microcirculation determines the oxygen delivery and solute transport to tissues. This process relies on the partitioning of RBCs at successive bifurcations throughout the microvascular network, and it has been known since the last century that RBCs partition disproportionately to the fractional blood flow rate, therefore leading to heterogeneity of the hematocrit (i.e., volume fraction of RBCs in blood) in microvessels. Usually, downstream of a microvascular bifurcation, the vessel branch with a higher fraction of blood flow receives an even higher fraction of RBC flux. However, both temporal and time-average deviations from this phase-separation law have been observed in recent studies. Here, we quantify how the microscopic behavior of RBC lingering (i.e., RBCs temporarily residing near the bifurcation apex with diminished velocity) influences their partitioning, through combined in vivo experiments and in silico simulations. We developed an approach to quantify the cell lingering at highly confined capillary-level bifurcations and demonstrate that it correlates with deviations of the phase-separation process from established empirical predictions by Pries et al. Furthermore, we shed light on how the bifurcation geometry and cell membrane rigidity can affect the lingering behavior of RBCs; e.g., rigid cells tend to linger less than softer ones. Taken together, RBC lingering is an important mechanism that should be considered when studying how abnormal RBC rigidity in diseases such as malaria and sickle-cell disease could hinder the microcirculatory blood flow or how the vascular networks are altered under pathological conditions (e.g., thrombosis, tumors, aneurysm). INTRODUCTION Partitioning of red blood cells (RBCs) through bifurcations of blood vessels determines how the oxygen is delivered to tissues and organs. Early studies have shown that this partitioning isa complex phenomenon(1–4) and, since then, it has been an area of intensive research, attracting a wide range of experimental (5–13) and numerical (14–20)studies.Moreaccurately, vessels with a higher flow rate tend to collect an even higher proportion of RBCs. This fact is known as the Zweifach-Fung effect and leads to heterogeneity in the hematocrit between various vessels (2). Uncovering the mechanisms for such a partitioning behavior is crucial to understand not only the transport of oxygen and solutes across vascular networks but also the development or remodeling of the networks themselves, as RBC dynamics were recently discovered to Submitted August 31, 2022, and accepted for publication March 13, 2023. *Correspondence: yazdan.rashid[email protected] or alexis.darras@ uni-saarland.de Yazdan Rashidi, Greta Simionato, and Qi Zhou contributed equally to this work. Editor: Timo Betz. SIGNIFICANCE We demonstrate in vivo that the lingering behavior of red blood cells (RBCs) at the apex of bifurcations modulates their partitioning in the microcirculation. The influence of this mechanism opens new ways of understanding how altered RBC properties in pathologies can hinder the microvascular blood flow, or influence the development or remodeling of the vascular networks. For instance, we show that rigid RBCs linger less at typical Y-shaped bifurcations and manifest a distribution pattern closer to the well-known Zweifach-Fung effect. We also highlight certain properties of the bifurcation geometry that contribute to the lingering intensity of RBCs. 1526 Biophysical Journal 122, 1526–1537, April 18, 2023 https://doi.org/10.1016/j.bpj.2023.03.020 Ó2023 Biophysical Society. This is an open access article under the CC BY-NC-ND license (http:// creativecommons.org/licenses/by-nc-nd/4.0/). underpin vessel remodeling, presumably by mediating the wall shear stress difference between neighboring branches via the route of effective blood viscosity (21). An empirical model that effectively recapitulates the time-average behavior of the Zweifach-Fung effect was developed by Pries et al. in the 1990s (22–24) and has been widely employed by the research community for studying microvascular blood flow. However, recent in vitro and in silico studies demonstrated that notable deviations from this empirical model can arise for various reasons (5,7,8,10,25,26). In particular, it has been shown that RBC partitioning at bifurcations in the smallest vessels of the microvascular network tend to deviate from the empirical predictions more severely (11,12,21,27,28). On the other hand, first established through numerical simulations (25) and subsequently validated by in vivo experiments (29), an intriguing behavior of RBCs in the microvasculature was highlighted: they linger at the apex of bifurcations (i.e., the cells temporarily residing near the branching point of vessels with diminished velocity) for certain periods of time, consequently modifying the dynamics of cells entering downstream vessels and the characteristic inter-cell distances featuring intermittent voids. Previously, Bagchi et al. (25) extensively characterized the time-dependent correlation between single-cell lingering times and geometry of the network, but the causal relationship between the lingering pattern and the average behavior of RBC distribution through the bifurcation was not explicitly investigated. Their following work (28) further combined a continuum model with the cellular simulation to distinguish the effects of plasma skimming and cell screening. More recently, Bagchi et al. reported how changes in the fraction of lingering population between rigidified and healthy cells are related to deviations from identity of the RBC flow rate and overall blood flow rate (30). Notwithstanding the progress above, a mechanistic understanding of the deviation from identity due to the Zweifach-Fung effect as well as the baseline behavior of how such deviations are associated with the proportion or absolute residence time of lingering cells is yet to be obtained. In the present study, we systematically analyze how the lingering of RBCs influences their partitioning at highly confined capillary-level bifurcations in the microcirculation (where the cells literally squeeze themselves through and negligible cell-free layer exists contrary to prior studies; cf. (22–24)). Through carefully designed in vivo experiments corroborated by complementary in silico simulations on timescales longer than the lingering time and relevant for physiological conditions, we show that global deviations from the empirical model by Pries et al. (24) are strongly correlated with normalized average duration of cell lingering events. We further reveal how lingering modulates RBCs partitioning, with a primary focus on mechanistic understanding of the deviation from the baseline behavior of Zweifach-Fung effect as postulated by the empirical model (24). The asymmetric drive from RBC lingering is found to either reinforce or revert the Zweifach-Fung effect. In addition, given that RBCs usually exhibit large deformations when lingering, we also investigate numerically how cells with increased stiffness differ from healthy ones in their lingering and partitioning behavior. Our results demonstrate that RBC lingering is a cause of statistical deviation from the Zweifach-Fung partitioning and provide further evidence that pathologically rigidified cells can impair the distribution of RBCs across the microcirculation (31–38) due to altered lingering properties. MATERIALS AND METHODS In vivo experiments Permissions Experiments were performed in Syrian golden hamsters according to the German legislation on protection of animals and were approved by the local governmental animal protection committee (permission number 25/2018). Hamsters were maintained on a standard 12/12 h day/night cycle and water and food were provided ad libitum. Detailed methods are described below for each type of investigation. Animal preparation and microscopy Hamsters with an age of 5–7 weeks, weighing 55–70 g, were used for the implantation of a dorsal skinfold chamber (39). The surgery was performed under deep anesthesia (150 mg/kg ketamin (Serumwerke Bernburg AG)/ 0.25 mg/kg domitor (Orion Pharma) intraperitoneal), with intraoperative pain medication by carprofen (10 mg/kg, Zoetis, subcutaneous). A previously described protocol (39) was followed. Briefly, the back of the hamsters was shaved and a titanium chamber consisting of two frames was implanted on the lifted dorsal skinfold. On the frame with a circular observation window (diameter of 10 mm) the cutis, subcutis, and retractor muscle were removed to expose the striated skin muscle for later observation of the microcirculation. The window was closed with a cover glass that was fixed with a snap ring. An implanted dorsal skinfold chamber is depicted in Fig. 1 a. Animals were allowed to recover for 72 h after the procedure. Hamsters were anesthetized as described above before intravital microscopy. A volume of 100 mL of the fluorescent plasma marker fluorescein isothiocyanate (FITC)-labeled (5%, 150 kDa, Sigma-Aldrich) was injected retro-orbitally and the animals were fixed on a Plexiglas stage, as illustrated in Fig. 1 b. Several capillary bifurcations in different areas of the chamber window were chosen for epifluorescence microscopy (Axio Examiner A1, Zeiss). FITC-labeled dextran was excited with a 480-nm light-emitting diode (LED) (Colibri 7, Zeiss), and image contrast was additionally enhanced by simultaneously using transmitted blue light, which is absorbed by hemoglobin, making RBCs appear darker. Imaging was performed with 20(LD A-Plan, NA ¼0.35, Zeiss), 50(LD EC Epiplan-Neofluar, NA ¼ 0.55, Zeiss), or 100(LD C Epiplan-Neofluar 100,NA¼0.75, Zeiss) long-distance objectives. Video acquisition was carried out with a digital camera (Hamamatsu Orca Flash 4.0, C13440) using the software ZEN 3.1 Blue (Zeiss). For each video, 2000 frames were recorded at different rates according to flow velocity (100–170 Hz). A 2 2 binning was applied during acquisition. An example of an obtained image is given in Fig. 1 c. Image analysis and lingering quantification Despite the animal being fixed on the Plexiglas stage, its breathing and/or muscles movement caused slight translations of the microscopic field of view in the image sequence (Fig. S1 a). To determine the motion of the RBCs in the vessels, we first correct those translations by translating each image by the maximum of its 2D correlation with the first image in RBC lingering and hematocrit distribution Biophysical Journal 122, 1526–1537, April 18, 2023 1527 the sequence. This process creates a movie where the vessels and surrounding tissues are still fixed (see Fig. S1 b). A Gaussian filter (with standard deviation of two pixels) is then applied to despeckel the images, and a mask is drawn around the vessels. The creation of this mask is facilitated by also considering the image showing the standard deviation of each pixel intensity along time, which is higher in the vessels with transient (dark) RBCs and offers a guide to the eye for the mask drawing (see Fig. S1 c). The diameters Dof the various vessels are computed as twice the average distance from the skeleton pixels of the masks to the border of the vessels, using standard Matlab functions bwmorph and bwdist (see Fig. S1 c). The image series is then inverted and binarized inside the mask area with a threshold based on a user-defined fraction of the Otsu threshold (40). The Otsu thresholding method was used because the average intensity of images is usually slightly flickering. However, the automated threshold algorithm overestimated the threshold required to detect the whole RBCs; hence, we implemented a user-defined ratio to adjust it (see Fig. S1 d). The hematocrit Hwas calculated based on the binarized image: H¼ 4s VRBC ARBCpRv, where 4sis the measured average area fraction of RBCs in the vessel, VRBC ¼56:5mm3and ARBC are the cell volume and the surface area (measured empirically in each vessel for a selection of single cells) of the RBCs respectively, and Rvis the vessel radius. This correction, similar to the approach in previous work (7), takes into account the deformation of the cells in the plane while avoiding any counting mistake that close RBCs (or the segmentation algorithm) could induce (see Supporting material for more details). This correction is validated using 2D projection slices of the numerical simulation, where the hematocrit value is known as simulation input; the error for the validation is below 1%. For detected RBCs, the standard Matlab function WeightedCentroid was used to calculate their center of mass. In this function, the center of mass is calculated based on the location and intensity ~ R¼PN i¼1PN j¼1½xi;yjIðxi;yjÞ= PN i¼1PN j¼1Iðxi;yjÞ, where xand yare the coordinates of each pixel and Iðx;yÞrefers to the intensity of the pixel at position ½x;y. A tracking algorithm was then used to determine the velocity of the detected RBCs (41). The variations of light intensity along the vessels, in conjunction with the collisions and transient aggregations of RBCs, did not allow us to reliably follow all RBCs along their entire trajectories. However, as we focused on bifurcations where the mother vessel ðMÞand the two daughter vessels have a hematocrit low enough to distinct single RBCs, this tracking allowed us to extract the spatial distribution of the RBC velocity (i.e., we performed particle tracking velocimetry (PTV)), on which we performed our further analyses. In this article, we will refer to the two daughter vessels as main daughter ðMDÞand secondary daughter ðSDÞ. By definition, MD is the daughter vessel with higher fractional blood flow rate and SD is the one with lower value. To detect the lingering of RBCs, we applied a method as in our prior work (29) by defining a circle with a diameter of 6mm (i.e., the main diameter of a hamster RBC) surrounding the bifurcation apex (Fig. 3 a). The cells were considered to linger if their center of mass was located in this circle and at the same time had a velocity lower than the local minimum detected in the probability density function (PDF) of a collective dataset of RBC velocities (Fig. 3 b). The velocity PDF displayed in Fig. 3 bwas obtained by considering the M vessel and the bifurcation area. In this study, flow and cell statistics were characterized within specific time intervals (10–20 s), longer than the transit time of individual RBCs but still short enough so that the flow is quasi-steady (i.e., no systematic trend and/or significant change in the volume flow rates can be identified). Within such intervals, self-consistent observables could be well defined. Numerical simulations Complementary simulations for time-dependent RBC flow in the representative capillary bifurcations were run with HemeLB (https://github.com/ hemelb-codes/hemelb) using the immersed-boundary-lattice-Boltzmann method following our previous approach (21). First, three-dimensional (3D) luminal surface models were reconstructed from binary masks of the capillary bifurcations using open-source software PolNet (42), assuming circular cross sections of varying diameter along each vessel. Then the flow domain enclosed by the luminal surface was uniformly discretized into cubic lattices of fine voxel size Dx¼0:25 mm, aimed at resolving cell dynamics with high resolution at the bifurcation apex. To initialize the simulation, inflow/outflow boundary conditions based on experimental time-average volume flow rates were imposed at the end of the M vessel and the MD branch (similarly defined as above for the experiments) of each bifurcation, and a reference pressure for the SD branch. No-slip boundary conditions were imposed at the vessel wall. The physical length corresponding to each simulation time step was Dt¼1:04 106s. The simulation for each bifurcation was initialized with (RBC-free) plasma flow. Once the plasma flow became converged, RBCs were randomly inserted from the Mvessel through a cylindrical flow inlet with constant feeding hematocrit ðHFÞas measured from experiments (see Supporting material section ‘‘simulation hematocrit and RBC initialization’’). When RBCs reached the end of MD or SD, they were removed from the simulation domain. Each RBC was modeled as a cytosol-filled capsule with an isotropic and hyperelastic membrane consisting of 5120 triangular facets. The mechanical properties of the RBCs were governed by elastic moduli governing different energy contributions, such as shearing and bending of the membrane. The cytosol was treated as a Newtonian fluid with the same viscosity as the suspending medium (i.e., blood plasma). The viscosity of the RBC membrane itself was not considered. To probe the effect of confinement, two different RBC diameters were considered (Drbc ¼6mmcorresponding to normal hamster RBC size and Drbc ¼8mmcorresponding to larger cells of similar size to a human RBC). Two different membrane stiffness levels (stain modulus ks¼5106;5105N=m) were considered for representing a normal and a hardened cell, respectively. For cell-cell and cellwall interactions, repulsive potentials inversely decaying with the distance were implemented. High-enough resolution was used to ensure the cellcell/cell-wall distance did not fall below one grid point. RESULTS AND DISCUSSION Determination of flow rates and baseline Zweifach-Fung partitioning To determine how the lingering of RBCs influences their partitioning, we used the empirical model of Pries et al. abc FIGURE 1 Experimental methods. (a) Dorsal skinfold chamber implanted on the back of a hamster. (b) Anesthetized hamster placed underneath the objective of an epifluorescence microscope. (c) An exemplar microvascular network imaged by fluorescence microscopy. The dyed plasma appears bright; RBCs in capillaries (see arrows) and the surrounding tissues are dark. To see this figure in color, go online. Rashidi et al. 1528 Biophysical Journal 122, 1526–1537, April 18, 2023 for the Zweifach-Fung effect (24) as a baseline. For this purpose, we first determined the blood flow rate Qin each vessel. Since the vessels were smaller than the main radius of the RBCs, we assumed a plug flow. We also assumed that the difference between the tube and discharge hematocrits was negligible in our observations. The flow rates within all bifurcation segments were therefore initially assessed from the velocities (V) and diameter (D) measurements of the individual segments CQD¼CVDpD 22 ;(1) similarly to previous studies (22,29,43). Mass conservation stipulates that the total volume flow rate from the M vessel QMis divided into the daughter branches and satisfy QM¼ QMD þQSD. The subscripts MD and SD refer to the main daughter and secondary daughter, respectively. To obey this conservation law and minimize experimental errors, values obtained from Eq. 1 were corrected following the normalization process described by Pries et al. in previous similar measurements (22). Similarly, the measured erythrocyte flow rates Ei¼HiQi,Hibeing the local hematocrit, were also corrected to obey EM¼EMD þESD. The empirical law for the Zweifach-Fung effect developed by Pries et al. (22–24) describes the relation between fractional blood flow, FQb, and the fractional erythrocyte flux, FQe, in the microvascular bifurcation. The fractional blood flows FQbfor both daughter vessels were calculated by dividing the flow rate Qin the respective daughter by the flow rate in the M vessel. The fractional erythrocyte flow rates FQe were computed as the ratio of the corrected erythrocytes flow rates (e.g., EMD=EMfor the MD). With these definitions, the model from Pries et al. (22–24) states where logitðxÞ¼lnðx=ð1xÞÞ and A,B,X0are related to the bifurcation geometry and feeding hematocrit in the M: where Hrefers to the measured hematocrit, calculated from the average area fraction of RBCs in the mask during the experiment (23,24). The subscripts indicate the vessel for which the quantity is considered. Due to the heterogeneous blood flow and diverse geometries in the microvasculature, the FQeand FQbdata extracted in vivo are naturally scattered (22,23), which may also apply for regular in vitro vascular networks (10,11). However, Pries et al. derived their empirical model for robust single-bifurcation predictability (see Fig. 4 in Ref. (22) and Fig. 2 in Ref. (23)) with in vivo data obtained through a classic protocol that enabled continual variation of the hematocrits in individual bifurcations. In the present study, we set out to investigate why individual bifurcation cases (a widely adopted approach (5–8, 16–18) to study the RBC partitioning behavior) would nicely scatter around this empirical limit as defined by Eq. 2, bearing in mind that the effects of flow ratio, vessel geometry, and feeding hematocrit were already accounted for in the original model. To do this experimentally, we consider a set of bifurcations selected in vivo to examine how the RBC partitioning in them deviates from the baseline split (as postulated by Eq. 2) and further correlate the deviation with experimental factors not considered by the empirical model. In particular, we focus on the lingering behavior, which appears a potential mechanism for such deviations in highly confined capillaries. Indeed, the initial study of Pries et al. analyzed arteriolar bifurcations with characteristic diameters about 8mm or larger (22); i.e., larger than the main diameter of hamster RBCs. In their studies of bifurcations with lower confinement level, the lingering behavior (if any) was probably not playing a significant role in the partitioning of cells. In the case presented here, studied vessel diameters range between 3.5 and 5.6 mm; i.e., considerably smaller than the diameter of the hamster RBCs. Such high confinement would favor AMD ¼13:29 D2 MDD2 SD 1 D2 MDD2 SD þ1ð1HMÞDM;ASD ¼AMD;B¼1þ6:981HM DM;X0¼0:9641HM DM; (3) FQe¼ 8 > > > > < > > > > : 0;if FQb<X 0 1;if FQb>1X0 logit1AþBlogitFQbX0 12X0;otherwise (2) RBC lingering and hematocrit distribution Biophysical Journal 122, 1526–1537, April 18, 2023 1529 the lingering of RBCs as the cells cannot easily avoid interacting with the vessel wall at the bifurcation. Lingering and partitioning Global description Experiments. Fig. 2 a–cshows three representative bifurcations (two Y shaped and one T shaped) where all required parameters could be extracted from the image sequences. A characteristic movie is provided as Video S1. We calculated the empirical predictions by Eq. 2 for each bifurcation (nine in total) and defined the deviations from the model as dZF ¼DFQeEX DFQeZF. Alternative definitions of the deviation (e.g., in signed form DZF) are discussed in the Supporting material (Fig. S7). We noted DFQe¼ FQeðMDÞFQeðSDÞ, the difference of fractional erythrocyte flow rate, and the expressions EX and ZF refer to the experimental data and the empirical prediction, respectively. The fractional blood flow rate FQbis given by the experimental measurements, and the predicted fractional erythrocyte flow rate FQZF eis computed through Eq. 2 (see Fig. 4 afor a schematic). To quantify the lingering of RBCs at each bifurcation, we first determined the average duration of the lingering events CtD. The total lingering time was considered as the sum of the duration of all detected lingering events. We then determined the ratio of this total lingering time to the number of the RBCs as the characteristic average lingering time CtD. Since this time is highly correlated with the flow rate and the transit time of an RBC through the bifurcation (i.e., the red circle defined in Fig. 3 a), we normalized it by ts¼ R=CVMD. The radius R¼3mm is the radius of the red circle defined in Fig. 3 aand VMthe average velocity in the M vessel. abc fed FIGURE 2 Representative bifurcations from experiments that are simulated. (a–c) Three bifurcations selected from the experiments hereafter referred to as BIF-a, BIF-b, and BIF-c, respectively. The M vessel, the MD, and SD are annotated with arrows indicating the flow direction. The plasma was stained with fluorescent dye (bright in the images). Dark areas in the vessels indicate the RBCs. The border of the masks used to analyze the bifurcations are depictedin colored symbols. Consistent colors and symbols are used for these bifurcations throughout the article. (d–f) Simulated RBC flow in the reconstructed domain of cropped bifurcations, respectively (a–c) as characterized experimentally. To see this figure in color, go online. FIGURE 3 RBC velocity distributions in the experiment and simulation. (a) Spatial distribution of RBC velocities obtained through particle tracking in experiment for BIF-b. Lingering is considered to occur if a cell is inside the red circle (6mm, size of a typical hamster RBC) with a velocity lower than the threshold velocity, which is determined as a local minimum in the PDF of velocities obtained from the M vessel and the bifurcation area (M þbf), as shown in (b). For comparison, the velocity PDFs for the MD and SD branches and the M vessel are also shown. (c) Corresponding color map of RBC velocities measured along simulated cell trajectories in BIF-b. (d) Velocity PDFs measured in different vessel branches of the simulated BIF-b as for the experimental data in (b). To see this figure in color, go online. Rashidi et al. 1530 Biophysical Journal 122, 1526–1537, April 18, 2023 The time tsis then a characteristic advection time that a cell would need to go through the lingering detection area (where the cell velocities diminish; Fig. 3 b). We defined the ratio of these two characteristic times as Pel¼CtD=ts. The ratio Pel defines the equivalent of a P eclet number, since it compares two characteristic times of different transport modes at the bifurcation (lingering and advection). As shown in Fig. 4 c for nine different bifurcations (see experimental details in Table S1), the deviation dZF from the Zweifach-Fung prediction Eq. 2 strongly correlates with Pel(Pearson correlation coefficient r¼0:74,p<0:05). Interestingly, we observed in some case that the deviation from the Zweifach-Fung effect caused by the lingering can revert the partitioning. Indeed, in the case of the (red) characteristic bifurcation highlighted in Fig. 2 a, the MD (with higher FQb) is the branch with lower FQe.Thisis often referred to in the literature as reverse partitioning (7,10,11,26). In our data, we only observed such reverse partitioning for the highest Pelz1:75. Simulations. The three representative bifurcations from experiments in Fig. 2 a–cwere also simulated under equivalent flow conditions and RBC volume fractions (see details in Table S1), with the lingering phenomenon of RBCs successfully reproduced (Fig. 2 d–f). These bifurcations were selected for numerical scrutiny because their experimental counterparts roughly cover the entire range of primary observables (e.g., Pel,dZF,dpt). A compilation of animation videos generated from numerical simulations are provided as Video S2.To quantify the lingering behavior, an equivalent procedure to the experimental data analysis was applied to analyze the simulation data too. Exemplar trajectories of RBCs near the bifurcation apex from simulations of varying feeding hematocrits are shown in Fig. S3. In line with the experimental PTV measurements, the cell velocities diminished as the RBCs approached the bifurcation apex (Fig. 3 c). After entering one of the daughter branches (MD or SD), the instantaneous velocity of RBCs would increase again. These patterns were well captured by the velocity PDFs compiled for individual vessel branches in each bifurcation (Fig. 3 d). Quantification of lingering events based on the simulated RBC velocities revealed that the percentages of lingering cells for the three bifurcations were 63%, 85%, and 4%, respectively. For hardened RBCs with a membrane strain modulus 10 times as large, the ratios became 23%, 77%, and 42%. Note that the proportion of lingering cells does not unequivocally determine Pel; rather, the average lingering time dominates, considered as the sum of time duration (monitored from all detected lingering events) divided by the total number of cell transits. In other words, a higher proportion of lingering cells does not necessarily lead to a higher P eclet number. For instance, a lower portion of cells lingers in BIF-a than in BIF-b, but they linger for longer time compared with their advection time, therefore causing a larger Pel. The correlation between the lingering intensity of RBCs and their partitioning at the bifurcations, namely Peland dZF, was evaluated similarly to the experimental data (Fig. 4 d). A strong association was found, indicating a linear increase of dZF against Pel(Pearson-rcorrelation coefficient 0.71, p<0:05). Nevertheless, the absolute magnitude of dZF and Pelwere substantially smaller (although in a proportional manner) than their experimental counterpart in Fig. 4 c.The FIGURE 4 Experimental and simulation results of RBC partitioning versus lingering. (a) Comparison of experimental RBC partitioning against the empirical predictions of the Zweifach-Fung effect (for BIF-a as in Fig. 2 a). The axes FQeand FQb are for the fractional RBC flux and fractional blood flow, respectively. The deviation from the empirical prediction was calculated as dZF ¼DFQeEX  DFQeZF.(b) The bars show experimental results against the lines, representing empirical predictions by Eq. 2. The data here are for the three characteristic bifurcations in Fig. 2 a–c, with consistent colors and symbols. Dashed lines and hollow symbols here refer to the SD vessel with lower fractional flow rate FQb. Note the inversion of the Zweifach-Fung effect in BIF-a. (c) Experimental deviation dZF from Eq. 2 as a function of the lingering P eclet number Pel. The points with error bars are experimental data, and the solid line shows linear regression fitting of the data points, with the shaded area indicating 95% confidence interval prediction from the fit. All data gathered from nine bifurcations (including the three bifurcations BIF-a, BIF-b, and BIF-c in Fig. 2 a–c) are included in this graph. The Pearson-rcorrelation coefficient (n ¼9) is r¼0.74 (p¼0.02). Error bars in (a–c) are computed by error propagation from the underlying experimental measurements. (d) Correlation between dZF and Pelin simulations (n ¼12). The crosses are simulation data and the solid line shows linear regression fitting of the data points. The Pearson-rcorrelation coefficient is r¼0.71 (p¼0.014). To see this figure in color, go online. RBC lingering and hematocrit distribution Biophysical Journal 122, 1526–1537, April 18, 2023 1531 highest P eclet number was only Pel¼0:65 in the simulations compared with a maximum value of Pel¼1:75 in the experiments.This discrepancy may have arisen from simulation configurations that are different from in vivo experiments; e.g., absence of endothelial surface layer (ESL; roughly 0.4–0.5 mm in thickness) and RBC glycocalyx in the numerical model. Indeed, although a repulsive cell-wall potential was implemented, it is possible that the ESL coated by ciliated structures interacts with the RBCs in more complex ways (20,44,45). For instance, it may contribute to enhanced confinement or give rise to an attraction force. Because these advanced cell-wall interactions have not been definitively characterized in the literature, without a tangible physical model, we could not capture such effects in our simulations. In any case, the correlation between the lingering P eclet number Peland the deviation from the Pries model was qualitatively similar between the experiment and the simulation. We note that no reverse partitioning (or inversion of the Zweifach-Fung effect) wasobservedin our simulations, which is in agreement with the fact that reverse partitioning in experiments only occurred for bifurcations with Pel>1:75. Lingering asymmetry in experiments For insights into how the lingering behavior could modify RBC partitioning at the bifurcation, we further examined its symmetry between the two daughter branches. Indeed, as can be seen from Video S1, the RBCs tend to enter the lower daughter branch ðMDÞalmost only if there is another cell already lingering at the bifurcation apex. Due to cell collisions at the apex, it was not possible to precisely determine the proportion of lingering cells entering each daughter merely by analyzing their trajectories. Therefore, we calculated the average advection time taof the cells between the (blue) mask in bifurcation and the (magenta) mask in the daughter (Fig. 5 a). This advection time was determined by calculating the convolution of the temporal hematocrit in the daughter mask and that in the bifurcation mask (see Fig. 5 aand Fig. S2 b). We then checked if there was a lingering event at the bifurcation (also monitored in a thin mask) at a time taearlier for each cell detection. We realized that, for such sorting, some lingering events at the bifurcation could not be attributed to any advected cell into the daughters. Indeed, the advection time of each cell can vary slightly, depending on its interactions with other cells near the apex. Therefore, after being temporally translated by ta, the cell detection can be located close to the lingering event but not exactly within its duration. For more accurate characterization, we therefore measured the proximity time tibetween the closest lingering event and the transit time (detection time in the daughter minus ta) for each cell (Fig. 5 b). Cells with a proximity time ti¼0 can then be determined as lingering. We further defined dpt¼jpSD pMDjas the difference between the proportion of lingering cells in the two daughter branches (Fig. 5 f). FIGURE 5 Asymmetry in RBC lingering at the bifurcation. (a) Example of mask set used to determine the advection time taof RBCs from the bifurcation (blue mask) to the SD branch (magenta mask) via convolution of the hematocrit in both masks (see main text and Fig. S2). Equivalent analysis is performed independently for the MD branch. (b) Extraction of the proximity time tiof each cell. The hematocrit variation over time, extracted from the daughter’s mask, is first translated by the average advection time ta(see Fig. S2). Then, for each peak of this hematocrit along time, a proximity time tiis obtained as the temporal distance from the closest lingering event. This process is performed independently for the two daughters. (c–e) CDF of the proximity time to lingering events for both daughters of the three characteristic bifurcations. When a significant difference (considerable dpt) is observed in the experiments, higher proportion of lingering cells (and a closer proximity with lingering events) is usually obtained in the daughter branch with the lower fractional RBC flux FQe.(f) The correlation between experimental dpt(difference in the proportion of lingering cells in SD and MD calculated from CDFs at ti¼0) and Pel. The Pearson-rcorrelation coefficient is 0.59 ðp¼0:09Þ. To see this figure in color, go online. Error bars are computed by error propagation from the experimental uncertainties of pSD and pMD. Rashidi et al. 1532 Biophysical Journal 122, 1526–1537, April 18, 2023 Because cells with a small proximity time timight also have been lingering or might have interacted with another cell lingering at the apex, we used the difference between the first two points of the cumulative density function (CDF) distribution of ti(Fig. 5 c–e) as error bars for pSD and pMD. As can be seen in Fig. 5 c–e, the proximity time tidistribution and the proportion of cells identified as lingering are distinct for the two daughters under large Pel. For the experimentally observed cases as highlighted in Fig. 5 f, the asymmetry in the proportion of lingering cells is weakly correlated with the lingering intensity Pel. In simulation, both dptand Pel for the simulated bifurcations (three in total) are smaller than their experimental counterpart but in a proportional manner with qualitative agreement (see Fig. S2 c). Origin of lingering in experiments To understand what causes the cells to linger, we identified geometrical features of the bifurcations correlating with the lingering P eclet number Pel. Due to the definition of the lingering events as a significant reduction of the RBC velocity at the bifurcation, it is likely that cells linger at the stagnation point; i.e., the intersection of the bifurcation’s wall and the flow divider plane. Note that a stagnation point is considered because the 3D vessels are projected in the mid-vessel plane. At this position, the cells will also experience the smallest drag force, likely not to push the cell toward any daughter vessel. The position of the stagnation point depends on the fractional flow rate of the daughter branches (46). The distance of the stagnation point from the center line Ls, which is depicted in the Fig. 6 a, can be approximated through QMD QSD ¼ 1þ2 parcsinLs RMþ2Ls pR2 Mffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi R2 ML2 s q 12 parcsinLs RM2Ls pR2 Mffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi R2 ML2 s q;(4) where RMis the radius of the M branch (see Supporting material for further justification). For the situations where numerical simulations were compared with the experimental data, positions of the stagnation point determined this way were located within one pixel of distance from the position indicated by the streamlines (see Fig. S4). This approximation can then be considered to be consistent with the available numerical data. Numerical simulations also showed that the shift of the stagnation point due to studied RBC transits were negligible compared with other experimental uncertainties (see Fig. S5). Since the cells at the stagnation point experience the smallest drag force, one can assume that the dominating forces on the cells arise from their interaction with the endothelial layer. In turn, this interaction probably depends on the geometry of the wall at this location, whichcan mostly be described through its curvature. To determine the curvature of each bifurcation, a circle was fitted on the neighborhood of the stagnation point FIGURE 6 Potential origin of RBC lingering. (a) Schematic highlighting the flow split and the stagnation point at bifurcations. Flow in the mother branch splits into flows QMD and QSD in the MD and SD branches, respectively. The stagnation point is indicated by the green point, whose distance from the apex Ls is highlighted with the arrow and scale symbol. The filled circle shows an idealized cross section of the mother vessel with simplified and flat flow separatrix. (b) Correlation of the lingering P eclet number Pelwith the curvature zof the bifurcation at the stagnation point. (c) Correlation between Peland the flow rate ratio between the SD and MD QSD=QMD.(d) Correlation between Peland the distance between the stagnation point and the bifurcation apex LS.(e) Correlation between LSand the curvature at the stagnation point z.In(b)–(e), the points are experimental data, the solid lines show linear regression fitting of the data points, and the shaded areas indicate the 95% confidence interval prediction from each fit. To see this figure in color, go online. Error bars in (b–e) are computed by error propagation from the underlying experimental measurements. RBC lingering and hematocrit distribution Biophysical Journal 122, 1526–1537, April 18, 2023 1533 (see Fig. S4). The curvature is defined as the inverse of the radius of this circle. To estimate the measurement uncertainty, we used multiple fittedcircles, each calculated using a different length of the neighborhood (lengths between 4and 8mm, centered on the stagnation point). The meanvalue of the curvatures is taken as the best approximation, whereas the standard deviation gives the error bars. The results showed that the curvature of the bifurcation at the stagnation point correlates with Pel(r¼0:59 and p¼0:09) (see Fig. 6 b). Interestingly, when testing correlations with further possibly influencing parameters, we found that the lingering P eclet number Pelalso correlates with the ratio of the flow rates between the SD and MD QSD=QMD (Fig. 6 c). Since this parameter determines the distance LSbetween the apex of the bifurcation and the stagnation point through Eq. 4, this parameter LSalso correlates with Pel(Fig. 6 d). Furthermore, the bifurcations where the stagnation point is further from the apex generally have a lower curvature at the stagnation point, as shown by Fig. 6 e. We therefore hypothesize that the flow rate ratio determines where the cells can linger, whereas the curvature of the endothelial layer at this lingering position determines the intensity of the lingering. The lingering P eclet number would then be determined by an interplay between global parameters of the bifurcation ðQSD =QMDÞand local geometries of the endothelial layer (the curvature at the stagnation point, z). When it comes to the experimental correlations, the result that a lower coefficient is found between Peland zthan Peland QSD=QMD likely arises from the fact that the determination of zalso relies on the approximation for LSand then the measurement of QSD=QMD. The error propagation on these successive quantities then likely decreased artificially the correlation between Peland the successively determined parameters LSand z. More geometrical parameters, which failed to demonstrate a significant correlation with Pel, are included in the Supporting material (see Fig. S6). Effect of cell rigidity and hematocrit in simulations Besides the curvature of the stagnation point at the apex, other factors may also affect the interaction between the cells and the bifurcation, leading to either weakened or enhanced RBC lingering. Taking a typical Y-type bifurcation for example (Fig. 7 a), three physiologically relevant aspects were numerically explored: the cell rigidity ks, the feeding hematocrit HF, and the cell size DRBC (Drbc ¼6mm unless otherwise specified). Increased cell rigidity (imitating diseased RBCs with higher-level stiffness) at the Y-type bifurcation is found to introduce strong cell-wall and cell-cell steric repulsion at the apex (see enlarged gaps near the branching point in Fig. 7 a), thus reducing the overall lingering frequency and average duration of individual lingering events as indicated by a smaller lingering P eclet number Pel(BIF-a in Fig. 7 b). The same is found for BIF-b too (Fig. 7 b). Contrarily, for a T-type bifurcation (less common bifurcation type in microvascular networks), stiffer RBCs tend to get stuck at the apex due to cell collisions and can lead to more intense lingering instead (see BIF-c in Fig. 7 b). The distinct effect of RBC deformability on its lingering behavior at the Y-type and T-type bifurcations observed here is in line with a recent report on capillary vascular network (30). On the other hand, for identical flow conditions, an enriched hematocrit (HF¼33%, relative to the baseline HF¼12%basedonexperimentalmeasurement)in the M branch has a weakening effect on RBC lingering at the bifurcation, whereas a reduced hematocrit ðHF¼10%Þhas a strengthening effect compared with the experimentally measured value HF¼12%for the same bifurcation (Fig. 7 c). The probable reason behind the decreased Pelin this case is that crowded cell traffic reduces the possibility of individual RBCs residing on the apex for extended periods of time as they are more likely to be pushed forward by following cells (‘‘herding’’ effect as reported in (47)). We also considered a different RBC size given the variability of cell morphology in vivo, which is modeled through increasing the default cell diameter DRBC ¼6mm from DRBC ¼6mmto DRBC ¼8mm while maintaining the same HF(Fig. 7 c). The increased RBC size leads to higher level of cell confinement in the capillaries and reduces the fluctuation amplitude of RBC lingering caused by hematocrit alteration. FIGURE 7 Effect of cell rigidity and feeding hematocrit HFon the lingering of RBCs in simulations. (a) Steric repulsion for (top) normal (N) and (bottom) hardened (H) RBCs at the apex of a Y-type bifurcation corresponding to BIF-a in Fig. 2.(b)Pelversus cell rigidity for normal and hardened RBCs at BIF-a, BIF-b, and BIF-c’’ (see Fig. 2). Solid lines are for simulations with HF matching corresponding experimental condition, whereas dashed lines are simulations with elevated hematocrit levels (twice as high for BIF-a and BIF-b, and three times as high for BIF-c). (c)Pel versus feeding hematocrit HFfor BIF-a with normal RBCs (in deformability) of two different sizes Drbc ¼6;8mm. To see this figure in color, go online. Rashidi et al. 1534 Biophysical Journal 122, 1526–1537, April 18, 2023