Cell-free layer development and spatial organization of healthy and rigid red blood cells in a microfluidic bifurcation

Abstract

Bifurcations and branches in the microcirculation dramatically affect blood flow as they determine the spatiotemporal organization of red blood cells (RBCs). Such changes in vessel geometries can further influence the formation of a cell-free layer (CFL) close to the vessel walls. Biophysical cell properties, such as their deformability, which is impaired in various diseases, are often thought to impact blood flow and affect the distribution of flowing RBCs. This study investigates the flow behavior of healthy and artificially hardened RBCs in a bifurcating microfluidic T-junction. We determine the RBC distribution across the channel width at multiple positions before and after the bifurcation. Thus, we reveal distinct focusing profiles in the feeding mother channel for rigid and healthy RBCs that dramatically impact the cell organization in the successive daughter channels. Moreover, we experimentally show how the characteristic asymmetric CFLs in the daughter vessels develop along their flow direction. Complimentary numerical simulations indicate that the buildup of the CFL is faster for healthy than for rigid RBCs. Our results provide fundamental knowledge to understand the partitioning of rigid RBC as a model of cells with pathologically impaired deformability in complex in vitro networks.

---

1. Introduction

Blood circulation in the human body is a crucial process ensuring the delivery of oxygen and various nutrients to our organs through the circulatory system. The circulatory system is a complex network of bifurcating and branching vessels. Such bifurcating vessels can strongly affect the distribution and hematocrit partition of passing red blood cells (RBCs).^1,2 The rigidity of RBCs is often thought to impair blood flow and affect the spatiotemporal organization of flowing RBCs. Rigid RBCs with impaired deformability are found in multiple diseases, such as malaria, diabetes, sickle cell disease, or acanthocytosis.^3–8 However, the fundamental mechanisms through which rigid RBCs modify blood flow have not been characterized extensively.

One of the non-trivial features of RBC flow in the circulation is their partitioning through vessel bifurcations.⁹ Recently, many studies have focused on parameters affecting the deviation from the well-established empirical model from Pries et al.^10,11 Among them, several studies have pointed out the role of cell focusing and the development of the cell-free layer (CFL) along the daughter vessels.¹² The formation of a CFL is crucial for blood flow in vivo, as it reduces its hydrodynamic resistance.^13,14 Furthermore, it is important for various biomedical applications such as plasma separation, especially in complex geometries,^15–18 and highly relevant for disease diagnostics of cells with impaired deformability.^19–21 Such an impaired RBC deformability severely impacts the viscosity and the shear-thinning behavior of blood and RBC suspensions, thus affecting hemorheology, flow resistance, and microvascular perfusion.^22–25 Moreover, RBC deformability can influence the emergence of a CFL or cell-depleted zones in complex microfluidic geometries that are often used in lab-on-a-chip devices.^11,26,27

In general, the CFL in microfluidic devices depends on various factors, including hematocrit, channel dimensions, and flow rate.^16,18 Moreover, changes in the CFL in dilute suspensions can also arise from geometric features of the channel, such as confinement,^28–30 constrictions,^15,27 and bifurcations.^11,31 Therefore, the formation of a CFL of healthy RBC suspensions in complex geometries, such as vessel networks has received increasing attention in recent experimental and numerical investigations.^11,32–36 Although the biophysical RBC properties, such as their deformability, were found to influence the CFL as well,^27,37 detailed knowledge about the effects of RBCs with impaired deformability on the spatiotemporal RBC organization and the partitioning in complex geometries remains scarce.

In this study, we explore how the rigidity of RBCs modifies their focusing and subsequent formation of a CFL in a bifurcating microfluidic channel. Therefore, we artificially rigidify RBCs using glutaraldehyde and examine their spatial distribution at various positions in the microfluidic T-junction. We investigate the evolution of the RBC distribution across the channel width, as well as the formation of a CFL along the mother and the subsequent daughter channels. Our results on rigid RBC are compared with investigations of healthy cells at the same RBC concentrations (0.1–5% Ht). We observe that the RBC rigidity drastically alters the RBC organization at the bifurcation after the mother channel, resulting in distinct lateral distributions of RBCs in the daughter channels. The magnitude of this effect increases with increasing inertia. Furthermore, we show how the CFL develops in the daughter vessels of the bifurcation. Our numerical simulations demonstrate that the distance that is required to reach a steady CFL in the straight daughter vessel is shorter for healthy than for rigid RBCs. Our work offers further insights into the flow behavior of RBCs with impaired deformability and how rigid cells can impair blood circulation.

2. Materials and methods

2.1 Experimental

2.1.1 Sample preparation. Blood is taken with informed consent from healthy voluntary donors. It is suspended in phosphate-buffered saline (PBS) solution (Gibco PBS, Fisher Scientific, Schwerte, Germany) and centrifuged at 1500 g for five minutes to separate RBCs and plasma. Sedimented RBCs are subsequently resuspended in PBS and the centrifugation and washing steps are repeated three times. For the final suspensions, hematocrit concentrations of 0.1% Ht, 1% Ht, and 5% Ht are adjusted in a PBS solution that contains 1 g L⁻¹ bovine serum albumin (BSA, Sigma-Aldrich, Taufkirchen, Germany). Furthermore, samples with artificially rigidified RBCs are prepared. For this, washed RBCs are incubated in a 0.1% glutaraldehyde (GA, grade I solution, Merck KGaA, Darmstadt, Germany) solution for one hour according to Abay et al.³⁸ Subsequently, RBCs are washed with PBS to remove excess GA and are suspended in a PBS/BSA solution at the same concentrations as for healthy RBCs. The RBC treatment using GA fixes the RBC in stasis. In this study, we are using the term “rigid” to highlight the lack of deformability in GA-treated cells, as opposed to the deformable nature of healthy RBCs.

Blood withdrawal, sample preparation, and experiments were performed according to the guidelines of the Declaration of Helsinki and approved by the ethics committee of the ‘Ärztekammer des Saarlandes’ (approval number 51/18).

2.1.2 Microfluidic setup. The RBC suspensions are pumped through a microfluidic chip that contains a T-junction geometry. The microfluidic device is fabricated using polydimethylsiloxane (PDMS, RTV 615A/B, Momentive Performance Materials, Waterford, NY) through standard soft lithography.³⁹ The T-channel has a height of H = 54 ± 1 μm in the z-direction and consists of an inlet channel, referred to as the mother channel, with a width of W_M = 104 ± 2 μm, and a length L_M = 3 cm in the flow direction. At the T-bifurcation, the mother channel splits into two daughter outlet channels, each with a width of W_D = 55 ± 1 μm and a length of L_D = 1.75 cm. Fig. 1a shows a schematic representation of the microfluidic T-junction and the used coordinate systems. The mother and daughter channels are connected with rigid medical-grade polyethylene tubing (0.86 mm inner diameter, Scientific Commodities, Lake Havasu City, AZ) to the sample and waste containers, respectively. The microfluidic device is mounted on an inverted microscope (Eclipse TE2000-S, Nikon, Melville, New York), equipped with a LED illumination and a 40× air objective (Plan Fluor, Nikon, Melville, NY) with a numerical aperture NA = 0.6. A high-precision pressure device (OB1-MK3, Elveflow, Paris, France) is used to pump the suspensions through the channel at constant pressure drops of p = 200 mbar, 400 mbar, 600 mbar, and 800 mbar.

Fig. 1 Schematic representation of the microfluidic setup and the experimental data analysis routine. (a) T-junction geometry, consisting of an inlet mother channel with width W_M, height H, and length L_M and two outlet daughter channels each with width W_D, height H, and length L_D. Dashed magenta boxes correspond to different regions of interest (ROI) for data acquisition in the mother and daughter vessels. (b) Image processing routine to determine the CFL thickness. Scale bars represent 20 μm. 250 single images are stacked and the standard deviation and median of the intensity are derived from the inverted image. Along each vertical pixel line, the intensity's standard deviation and the intensity's median are calculated. Two representative plots of the two parameters are shown below the images, corresponding to a horizontal position marked by the dashed white lines. Red dots correspond to a threshold value of 0.5% of the measured maximum inverted intensity to detect the RBC flow and blue dots highlight the peak position of the median intensity, which corresponds to the position of the walls. The red and blue lines in the final image correspond to the determined channel borders and borders of the RBC flow, respectively.

Microfluidic experiments are performed both with symmetric as well as asymmetric flow rates in both daughter vessels. For the symmetric case, both daughters have the same length L_D. To introduce an asymmetric flow rate, the length of one daughter vessel is shortened to L_D* = 3L_D/4, L_D/2, or L_D/4, while the length of the other daughter channels is kept constant at L_D. Shortening the length of one daughter channel decreases the hydraulic resistance in the channel, hence, leading to an increase in the flow rate in the short daughter vessel. We define the asymmetry ratio as AR = L_D/L_D*, which is one in the case of a symmetric bifurcation.

2.1.3 Data acquisition and analysis. To determine the evolution of the RBC distribution across the channel width and to detect the CFL, the RBC flow is recorded at different regions of interest (ROIs) along the flow direction. Magenta boxes in Fig. 1a indicate ROIs corresponding to two positions at the beginning and end of the mother channel before the bifurcation, respectively, one at the T-bifurcation, as well as three positions in the daughter channels after the bifurcation.

At low hematocrit (≤ 1% Ht), a high-speed camera (MEMRECAM GX1, NAC Image Technology, Salem, MA) with a frame rate of 20 000 fps (frames per second) is used, which allows us to resolve the trajectories and determine the velocities of individual RBCs by means of particle-tracking velocimetry (PTV). Therefore, the positions of individual RBCs are determined in each frame and linked in consecutive frames to calculate the RBC trajectories using a self-written Matlab (the MathWorks, Natick, MA) PTV algorithm. Based on those trajectories, the velocity profile and the RBC distribution across the channel width in the y-direction are calculated.

To assess the relative influence of inertial forces, we calculate the Reynolds number Re, which relates the inertial to viscous forces in the system. Here, we define it as Re = u_maxD_hρ/η, with the fluid density ρ = 1 g cm⁻³, the maximum cell velocity u_max in the channel center, the fluid's dynamic viscosity η = 1.2 mPa s, and the hydraulic diameter of the rectangular microfluidic channel D_h = 2W_MH/(W_M + H). Representative plots of the velocity profiles in the mother and daughter channels are shown in Fig. S1 in the ESI.† Based on the applied pressure drop and the resulting maximum cell velocity u_max (see Fig. S2 in the ESI†), the Reynolds number is always greater than one (Re = 6–24), hence inertial forces cannot be neglected. Moreover, we use the PTV analysis to assess the flow in both daughter channels in the case of an asymmetric bifurcation. The fraction of the volumetric flow rate as a function of the asymmetric ratio is in good agreement with the analytical solution based on the length ratio L_D*/L_D between both channels (see Fig. S3 in the ESI†).

At higher RBC concentrations, individual cells cannot be discriminated and the cell distribution across the channel width and the CFL cannot be determined based on the RBC trajectories. Therefore, RBC flow is recorded with a high-speed camera (Fastec HiSpec 2G, FASTEC Imaging, San Diego, CA) at a lower frame rate of 1500 fps. A self-written Matlab program is used to determine the channel borders, the core RBC flow, and thus, the thickness of the CFL, as exemplified in Fig. 1b. First, 250 images of the recorded sequence are inverted and then stacked. Second, the pixel-wise standard deviation and median of the intensity in the stack are calculated. Regions with a high standard deviation correspond to locations in the channel where the image intensity fluctuates the most, hence the core RBC flow. The median of the intensity is most prominent at the channel borders, which do not change their position over the image sequence. Based on the peak positions of the median inverted intensity and a threshold of 0.5% of the measured maximum inverted intensity, the channel borders and the region of the RBC flow are determined, respectively. Subsequently, the thickness of the CFL is calculated for each vertical pixel line. A representative depiction of the resulting final image of the border and core RBC flow detection is shown in the bottom left of Fig. 1b. Representative single images and image stacks for all used RBC concentrations are shown in Fig. S4 in the ESI.† Experimental CFL results are shown as mean values with corresponding standard deviations between different measurements as error bars.

2.2 Simulations

We simulate the fluid motion in three dimensions using the standard D3Q19 lattice Boltzmann method (LBM)^40,41 with the Bhatnagar–Gross–Krook (BGK) collision operator.⁴² The RBCs are modeled as hyperelastic biconcave capsules endowed with shear and bending resistance. We use the Skalak strain model to account for the RBC membrane's shear elasticity and local area dilatation.^43,44 The resulting elastic force is computed using a linear finite element method (FEM).⁴⁵ The Skalak model has been initially tailored for RBCs and used successfully to recover their dynamics and morphologies under different flow conditions.^46–49 The bending is accounted for using the Helfrich free energy developed for bilayer lipid membranes.⁵⁰ The bending force is obtained through the functional derivative of the Helfrich free energy, and its terms are numerically evaluated using discrete differential geometry operators for triangulated meshes.^48,51 The impermeability and inextensibility of the RBC membrane are fulfilled by imposing penalty functions on both the global volume and area.⁴⁷ The coupling between the fluid and the particles is based on the immersed boundary method (IBM).⁵² Additional details on the numerical method are provided in Appendix A at the end of the manuscript.

Simulations are performed in a straight channel with the same dimensions as the daughter vessels in the microfluidic experiments to determine the steady-state behavior of the CFL. Since the magnitude of the CFL decreases with increasing RBC concentration, we limit our numerical simulations to the lowest hematocrit of 0.1% Ht that was tested experimentally. In the numerical simulations, the CFL is calculated based on the trajectories of the RBCs’ center of mass. Since the RBCs have a spatial extension between 2 μm and 4 μm around their center of mass, the CFL obtained from the numerical simulation should be reduced by a value between 2 μm and 4 μm to be compared with experimental measurements. Hence, numerical data of the CFL is shown as the cell's center of mass +3 μm and with an error of ±1 μm.

3. Results

3.1 RBC focusing and distribution at low concentrations in the microfluidic T-junction

At finite inertia (Re > 1), particles and cells in dilute suspensions can show ordering and flow focusing in straight vessels.⁵³ In the mother channel, we observe a focusing of RBCs along the flow direction, as shown in Fig. 2a for healthy and rigid RBC at 0.1% Ht and Re = 24. The top and bottom panels in Fig. 2a correspond to an x-position close to the beginning (x_M ≈ 0) and end (x_M = L_M) of the mother channel, respectively. At the beginning of the mother vessel (x_M ≈ 0, Fig. 2a top), cells are homogeneously distributed across the channel width for both healthy and rigid RBCs. Furthermore, RBCs flow in the close vicinity of the channel walls. Hence, we do not experimentally observe a CFL at x_M ≈ 0, indicated by the histogram bars at y_M/W_M = ±0.5. However, at the end of the mother vessel (x_M = L_M), we find pronounced differences in the distributions for healthy and rigid RBCs. Healthy cells are focused on three equilibrium positions, namely one in the channel center and two close to the borders. Rigid RBCs preferentially flow in the center of the channel. Furthermore, both healthy and rigid RBCs exhibit a pronounced CFL close to the channel walls at x_M = L_M. The observed focusing effect towards distinct positions for healthy and rigid cells emerges for Re ≥ 12 and its magnitude increases with increasing Re. Complementary results for other investigated Re are shown in Fig. S5 in the ESI.†

Fig. 2 Flow of a 0.1% Ht RBC suspension through the microfluidic T-junction. (a) Probability density distributions (pdf) for healthy (left) and rigid (right) RBCs at the beginning (top) and end (bottom) of the mother channel at Re = 24. Dashed vertical lines indicate the position of the channel borders. (b) CFL at the end of the mother channel as a function of Re at 0.1% Ht (top) and the RBC concentration at Re = 18 (bottom). Horizontal lines correspond to the optical resolution limit of two pixels. (c) Distribution of healthy (top) and rigid (bottom) RBCs at three positions along the flow direction in the daughter vessels for the symmetric partitioning at Re = 24. (d) Stacks of 500 images showing RBC partitioning at the bifurcation for different asymmetry ratios AR. The red and blue lines in the image stack correspond to the channel borders and the RBC flow, respectively. For AR = 4, representative trajectories are shown. Scale bars represent 50 μm. (e) Relative RBC flux in the shorter daughter vessel n*/n₀ (top) and the ratio of the upper CFL between the short and the long daughter CFL*/CFL as a function of AR.

Based on the RBC distributions across the channel width and the positions of the vessel borders, the CFL thickness is determined at the end of the mother vessel x_M = L_M. Fig. 2b shows the dependency of the CFL for both rigid and healthy cells as a function of Re at 0.1% Ht (top) and as a function of the RBC concentration at a constant Re = 18 (bottom). Within the investigated pressure drop range, the CFL only increases slightly for both rigid and healthy cells with increasing Re. However, we observe a strong decrease in the CFL thickness with increasing RBC concentration, as shown in the bottom panel of Fig. 2b. In general, rigid RBCs form a smaller CFL compared to healthy RBCs, in agreement with previous studies.¹⁶ This difference is most pronounced at low RBC concentrations.

At the end of the mother channel, the RBCs reach the T-bifurcation, where the flow splits into the two daughter branches. Similar to the mother channel, we perform RBC tracking at the beginning (D1: x_D ≈ 0), the middle (D2: x_D ≈ L_D/2), and the end (D3: x_D ≈ L_D) of the daughter branches. For the sake of simplicity and based on the top view projection of the T-junction, we refer to positions at −0.5 ≤ y_D/W_D < 0 as up and to positions at 0 ≤ y_D/W_D ≤ 0.5 as down. Hence, the upper CFL in the daughter branches corresponds to the continuation of the CFL in the mother channel. Fig. 2c shows the RBC distribution across the daughter channel width for healthy (top) and rigid (bottom) RBCs after the symmetric bifurcation. Right after the T-bifurcation in D1, we observe a strong accumulation of healthy RBCs at an off-centered position at y_D/W_D ≈ −0.28 close to the upper channel wall of the daughter vessel. With increasing y_D/W_D, the local cell concentration at this position continuously decreases. Note that in D1, RBC flow in close proximity to the down wall at y_D/W_D = 0.5, while a large upper CFL (−0.5 ≤ y_D/W_D ≤ −0.3) is formed. As the RBCs flow through the daughter vessels, the peak in the distribution diminished and a more uniform distribution is found in D2 and D3 for the healthy RBCs. Furthermore, a CFL appears at 0.3 ≤ y_D/W_D ≤ 0.5 as the RBCs flow in the daughter branches.

In contrast to the healthy cells, the distribution for rigid RBCs exhibits a broad off-centered peak around y_D/W_D ≈ 0.2 right after the bifurcation in D1, as shown in the bottom panel of Fig. 2c. While flowing along the daughter channel, more rigid RBCs migrate towards the upper half of the channel (−0.5 ≤ y_D/W_D ≤ 0). However, the peak in the distribution at y_D/W_D ≈ 0.2 persists for rigid RBCs until the exit of the daughter vessel at D3. With increasing Re, it seems that these distinct focusing effects for healthy and rigid RBCs increase, as shown for other investigated Re in Fig. S6 in the ESI.†

We further examine how an asymmetric bifurcating flow affects the RBC partitioning and the CFLs at the bifurcation. Fig. 2d shows representative image stacks for the bifurcating flow of a 0.1% Ht healthy RBC suspension at different asymmetry ratios AR. This asymmetry leads to different flow rates as well as to different CFLs in both daughter vessels. To assess the relative RBC partitioning and flux, individual RBCs are tracked at the T-junction, representatively shown for 90 trajectories at AR = 4 in Fig. 2d. The two colors represent trajectories that end in either the left or the right daughter channel. The fraction of RBCs that end in the shorter daughter vessel n* with respect to the total number of RBCs that enter the bifurcation n₀ is plotted as a function of the asymmetry ratio for healthy and rigid cells in the top panel of Fig. 2e. In the case of a symmetric bifurcation (AR = 1), we find n*/n₀ = 0.5. However, with increasing AR > 1, the number of RBCs that flow into the short daughter channel increases. This increase seems to be more pronounced for rigid than for healthy RBCs. Moreover, as visible from the snapshots in Fig. 2d, two distinct CFLs emerge close to the two channel borders in each daughter vessel. In Fig. 2c, we observe that the down CFL at the impacting wall is essentially zero, while the upper CFL closer to the mother channel is more pronounced. The ratio between the two upper CFLs of the short and long daughter CFL*/CFL is shown in the bottom panel of Fig. 2d. Here, we find that with increasing AR the CFL*/CFL ratio decreases.

3.2 CFL development in symmetric daughter vessels

To quantify the development of the CFL in the daughter vessels after the bifurcation in-depth, we employ the image processing routine shown in Fig. 1b. Therefore, 50 images are stacked at each ROI and the position of the RBC core flow and the channel borders are determined. Representative images for the three positions along the daughter vessels are shown in Fig. 3a. Similar to the histogram-based analysis in Fig. 2, the notations up and down correspond to the CFL in negative and positive y_D -direction, respectively.

Fig. 3 CFL development along the flow direction in the daughter vessels. (a) Representative superimposed images for a 0.1% Ht RBC suspension at three positions in the daughter channel at Re = 18 for (top) rigid and (bottom) healthy RBCs. The red and blue lines in the images correspond to the channel borders and the RBC flow, respectively. Scale bars represent 20 μm. Up and down CFL correspond to the CFL in negative and positive y_D -direction, respectively. (b) Up and down CFL for a 0.1% Ht RBC suspension as a function of x_D for different Re.

Fig. 3b shows the development of the asymmetric up and down CFLs along the flow direction in the daughter for different Re. At Re = 6, the upper CFL is constant in the first two regions and slightly decreases at the end of the daughter vessel. Rigid RBCs show similar behavior as healthy RBCs at Re = 6. However, they generate a smaller CFL compared to healthy cells. At x_D = 0, the down CFL is essentially zero, since cells that flow in the center of the mother channel are pushed against the wall in the stagnation area of the bifurcation. As the cells travel along the daughter, the down CFL consistently builds up for both rigid and healthy RBCs. Similar to the upper CFL, the down CFL for healthy cells is larger than for rigid cells at the end of the daughter channel. However, for Re = 6, the thickness of the down CFL does not seem to have reached the thickness of its upper counterpart, thus forming an asymmetric CFL at the end of the daughter branch at x_D = L_D.

Increasing Re affects the CFL development in the daughter vessel as follows:

(i) At the beginning of the daughter vessel at x_D = 0, increasing Re leads to an increase in the upper CFL for both healthy and rigid RBCs. Due to the increased inertia, RBCs that initially flow in a central position of the mother channel are pushed into the T-junction wall opposite the mother channel. Therefore, the down CFL is still zero. Moreover, RBCs that flow at streamlines close to the channel walls in the mother channel are pushed away from the wall to positions closer to the center of the daughter channel, thus forming a larger CFL directly after entering the daughter branch at x_D = 0, reminiscent of cell-depleted zones generated by lip vortices in other complex flow fields of RBC suspensions.⁵⁴

(ii) In the second region of interest at x_D = L_D/2, the thickness of the upper CFL decreases for Re > 12 compared to the upper CFL at the entrance of the daughter vessel. Opposite to this trend, the down CFL increases.

(iii) At the end of the daughter vessel at x_D = L_D, experiments at Re ≥ 12 result in the formation of an equal thickness in the up and down CFL. This symmetric CFL at x_D = L_D is independent of Re within the investigated regime.

Similar behavior regarding the CFL development in the daughter branches is also observed for higher RBC concentration, as shown in Fig. S7 in the ESI.† However, increasing the RBC concentration overall reduces the thickness of the CFL and results in a symmetric CFL already at Re = 6 for RBC suspension with ≥1% Ht.

3.3 Assessment of the steady-state CFL through numerical simulations

We further employ numerical simulations to determine the steady-state CFL of dilute healthy and rigid RBC suspensions at Re > 1. Fig. 4 highlights the results of the numerical simulations regarding the development of the CFL at Re > 1. Representative developments of the CFL along the flow direction for healthy and rigid RBCs at Re = 12 are shown in Fig. 4a. Here, the x-axes are normalized by the length of the daughter branches L_D in the microfluidic experiments. At a given Re, healthy RBCs form a larger steady-state CFL than rigid cells, as indicated by the black dashed horizontal lines in Fig. 4a. The same behavior is observed for the other investigated Re, as shown in Fig. S8 in the ESI.† The distance x_s and the time t_s = x_s/u_max that are required to reach the steady-state behavior in the numerical simulations are shown in Fig. 4b as a function of Re. Here, we define x_s when 95% of the value for the mean steady-state CFL is reached. We find that at Re = 12, healthy RBCs reach the plateau CFL already around 0.2 × L_D, while rigid cells require a longer length of approximately 0.25 × L_D to achieve their steady-state CFL. Moreover, x_s does not systematically depend on Re for both healthy and rigid RBCs. The time t_s to reach the CFL plateau decreases with increasing Re, as shown in the lower panel of Fig. 4b.

Fig. 4 Numerical simulations of healthy and rigid RBC suspension at 0.1% Ht in a straight channel. (a) CFL development along the flow direction at Re = 12. The x-axes in (a) are normalized by the length of the daughter branches L_D in the microfluidic experiments. The enveloping shaded areas correspond to ±1 μm, as discussed in the Method section. Black, dashed horizontal lines indicate the mean steady-state CFL and the black circles represent the distance x_s that the RBCs travel until the steady-state CFL is reached. (b) Length x_s (top) and corresponding time t_s (bottom) that are required to achieve a steady-state CFL in the simulations as a function of Re. (c) Comparison between the experimental and numerical (hatched bars) results for different Re.

The CFL for the experiment and the steady-state CFL in the simulations are shown in Fig. 4c for different Re. Here, we plot the upper CFL for the microfluidic experiments at the end of the daughter vessel (see Fig. 3), which corresponds to the formation of a symmetric CFL. In Fig. 4c, we find that the CFL in the simulations does not significantly depend on Re within the investigated range for both healthy and rigid RBCs. In the daughter vessel, we observe smaller CFLs for rigid than for healthy cells in the simulations, in agreement with our microfluidic experiments. Further, the behavior of CFL with Re is similar to the experimental observations in the mother channel (see Fig. 2b). Overall, the CFL determined in the numerical simulations is slightly higher by roughly 1–2 μm compared to the experiments, similar to other numerical and microfluidic studies.⁵⁵ These differences could be due to the correction that is applied to determine the RBC surface based on the center of mass from the numerics. Another factor that could contribute to these differences could arise from a slight misalignment of the x–y-plane of the microfluidic chip with the optical axes of the microscope. Due to the relatively deep channel (H = 54 μm), having the chip slightly tilted can lead to an underestimation of the CFL.

4. Discussion

Flowing through the inlet mother channel, RBCs are focused to specific positions along the channel width for healthy and rigid RBCs, as shown in Fig. 2a. For healthy RBCs, the emergence of two pronounced peaks close to the walls of the channel at the end of the mother channels is in contrast to the idea of an RBC core flow near the channel center, accompanied by a CFL near the channel walls.^2,56 Non-uniform RBC distributions in straight channels have previously been reported.^55,57 Zhou et al.⁵⁵ observed so-called off-center two-peak profiles for healthy RBCs in straight channels with dimensions on the same order of magnitude as the mother channel in this study. The authors argue that the emergence of these off-center two-peak profiles is co-determined by the spatial decay of hydrodynamic lift and the global deficiency of RBC dispersion in dilute suspensions.⁵⁵ The positions of the peaks move towards the channel centerline along the flow direction until they stabilized around y/W ≈ ±0.4, similar to the peak positions observed in Fig. 2a at the end of the mother channel. However, in their study, these peculiar RBC distributions emerged at negligible inertia (Re < 2 × 10⁻⁴).⁵⁵ For dilute suspensions with hematocrits <1% shown in Fig. 2 at Re > 1, focusing effects can severely impact the spatiotemporal distribution of particles and cells in tubes and microchannels.

In general, lateral migration and focusing of particles in circular pipes due to inertial lift forces were first reported by Segré and Silberberg.⁵⁸ Inertial focusing is mainly driven by two dominating opposing forces. On one hand, the wall-induced lift force, due to the interaction between the particle and the channel wall, pushes the particle away from the wall. On the other hand, the shear gradient-induced lift force, due to the curvature of the velocity profile in the channel, directs the particle away from the channel center.^37,53 These forces drive particles and cells towards specific equilibrium positions in the channel cross-section, which depend on the channel geometry, particle size and rigidity, and flow velocity. In square or rectangular microchannels, focusing of particles or cells is achieved for particle Reynolds number Re_p of the order of equal or larger than one, where Re_p = Re(a/D_h)² with the particle or cell size a. Further, a certain length L_f is required for particles to reach stable equilibrium positions along the flow direction.³⁷ For rectangular channels, different stable focusing positions located near the faces of the channel exist. Depending on the aspect ratio of the channel cross-section, the distance from the inlet, and Re, particles and cells can be focused towards the wide channel faces or both the short and wide faces.^59–65 For the experiments shown in Fig. 2, and with the RBC diameter a = 8 μm, we find Re_p ≈ 0.24, L_f ≈ 17 mm, and Re = 24. Although the Re_p is smaller than one, we expect inertial lift forces to partially affect RBC ordering and influence the observed distributions of healthy RBCs at Re > 1. The lift force acting on the RBCs and hence, the focusing and separation of cells, further depend on their shape.^24,66–68 In our study, the RBCs are rigidified in stasis where the cells exhibit a biconcave disk-like shape. While the healthy RBCs are able to adapt their shape according to the flow conditions in the channel, rigid RBCs retain their shape fixed in stasis. Hence, the difference in the distribution of rigid RBCs shown in Fig. 2a emerges as a consequence of their impaired deformability, which results in different shapes in flow, in agreement with studies on particles and ellipsoids⁶⁶ and rigid RBCs.⁶⁷ Furthermore, Shen et al.¹¹ investigated the effect of RBC deformability on the inversion of hematocrit partition at microfluidic bifurcations. At low RBC concentrations, they observed an inversion of the hematocrit profiles that we observe for healthy and rigid cells. However, instead of using GA to irreversibly fix the RBC shape, the authors used different dextran concentrations to tune the viscosity contrast between the inner cytosol and surrounding fluid. This changes the dynamics of lift and hydrodynamic interactions of cells, as well as their transient shape dynamics,⁶⁹ and might therefore lead to contradictory behavior with respect to the results reported in Fig. 2a.

At the bifurcation at the end of the mother channel, we observe RBC partitioning, which is affected by the asymmetry ratio of the T-channel, as shown in Fig. 2c and d, as well as distinct distributions for healthy and rigid RBC, shown in Fig. 2e. Early theoretical model predictions⁷⁰ and experiments⁷¹ on microvascular bifurcations revealed a nonlinear relation between the fractional RBC flux and the fractional bulk flow rate at Re ≪ 1. The classical Zweifach-Fung effect states that the daughter vessel with the highest flow rate will collect a higher hematocrit.⁹ One simple rationalization is that the RBCs focus on the streamlines with a higher flow rate, while the CFL is formed and contains the streamlines with lower flow rates. The daughter vessel with the highest flow rate then collects blood enriched in RBCs, while the other collects a bigger proportion of blood depleted from RBCs. The focusing and characteristics of the CFL are important parameters influencing the partitioning of RBCs at bifurcations. Recent studies highlighted that the asymmetry of cell focusing observed after a bifurcation can significantly modify the partitioning of RBCs in subsequent bifurcations.^11,12,36,72 Additionally, it has been shown that at low hematocrit, RBC flow tends to deviate from the Zweifach-Fung effect and that deformability plays an important role in this reverse behavior of partitioning.¹¹

Additionally to our results on the RBC focusing phenomenon, we provide further insight into the CFL formation along the mother and daughter channels and their dependency on the applied pressure drop, hematocrit, and cell rigidity, as highlighted in Fig. 2 and 3. The formation of a CFL is attributed to the tendency of RBCs to laterally migrate away from the vessel wall. The main driving forces for this migration include size exclusion effects due to the finite size of the RBC, boundary interactions, and macromolecular layer exclusion effects due to the inner lining of the vessel, and the curvature of the velocity profile in the vessel that produces a tendency for migration toward the vessel centerline.^1,2,16,73 As shown in Fig. 2a and 3, we do not observe a CFL at the beginning of the mother channel and similarly, at the beginning of the daughter vessel for the down CFL.

For straight microfluidic channels, the CFL growth along the flow direction was investigated before. Zhou et al.⁵⁵ recently investigated the CFL growth in a straight rectangular channel comparable to the dimensions of the mother channel used in this study, however, at negligible inertia Re ≪ 1 using numerical simulations and microfluidic experiments. At low RBC concentrations, they observed a build-up of the CFL that followed a power-law behavior with exponents between 0.26 and 0.4. The CFL in their simulations increased over a length of 28 × L/D_h without saturation. In their microfluidic experiments, the CFL growth occurred over more than 46D_h without reaching an equilibrium for RBC concentrations of 1%. In our study at Re > 1, the RBC travel over roughly 350 × L/D_h at the end of the daughter vessel showing both saturation and recovery of symmetry. Based on the evolution of the lower CFL in the daughter vessel in the experiments, shown in Fig. 3, we consider the CFL to have reached a steady state at x_D = L_D/2 (≈ 175D_h) for 0.1% Ht. Our numerical simulations further show that the length to reach the steady-state CFL is not significantly influenced by Re. However, the simulations also demonstrate that healthy RBCs achieve their equilibrium CFL under flow faster than rigid cells. For higher RBC concentrations in the experiments, the CFL saturates even faster, as shown in Fig. S7 in the ESI.† This is due to enhanced RBC interactions and collisions in the RBC core flow that increase with the hematocrit and that push the cells toward the channel walls.

Besides this dependency of the CFL growth on the RBC concentration, the magnitude of the CFL decreases with increasing hematocrit, similar to previous works that demonstrated the strong role of hematocrit in determining the RBC organization and the CFL thickness of healthy RBCs in straight channels.^57,74 Additionally, we find that the CFL thickness is always smaller for rigid RBCs compared to healthy RBCs under the same experimental conditions (Re and hematocrit) and at the same position (x_D or x_M), in agreement with previous studies.^26,27,33,34 In bifurcating microfluidic channels, the partitioning and flow of RBC suspensions and the emergent CFL revealed heterogeneous RBC distributions, skewed and blunt velocity profiles, and an enhancement of the thickness of the CFL at higher hematocrit.^{31,32,75–78} Yamamoto et al.³⁵ studied the partitioning of RBCs through asymmetric bifurcating microchannels at 0.5% Ht. They observed that with an increasing fractional flow rate in the daughter vessels, the CFL ratio between the mother and the daughter vessel decreases, similar to our results shown in Fig. 2d. While other numerical studies already showed that rigidification of RBCs tends to decrease the CFL width,^79,80 those previous studies do not provide a comparison with experimental data as our study does. This is particularly relevant since most of the combined studies showed some discrepancies between simulations and experiments.^20,55

Zhou et al.³⁶ found that changing the fractional flow rate between the daughter vessels leads to a complete depletion of RBCs in one branch of the successive bifurcation. Similarly, the distinct RBC distributions of healthy and rigid RBCs in the different ROIs in the daughter vessels (see Fig. 2e) would lead to dramatically different RBC organizations in subsequent branching channels. Namely, the formation of a pronounced peak of rigid cell concentration close to the lower daughter walls is a potential mechanism through which rigid cells can perturb the healthy flow of RBCs in the circulatory network.

Understanding the distribution profile and CFL development in bifurcation channels is paramount to understanding the flow and RBC in vascular networks. In such in vitro network structures, which often consist of multiple series-connected bifurcations, the CFL in successive bifurcations is crucially impacted by the RBC distribution in the previous bifurcation.^{33,34,72,81,82} Merlo et al.⁷² recently experimentally demonstrated that upstream bifurcations drive the spatial distribution of healthy RBCs in model microfluidic networks. The authors showed that persistent perturbations in the hematocrit profile after diverging bifurcations affect RBC partitioning at the next bifurcations. These perturbations were shown to propagate from bifurcation to bifurcation.⁷² Based on our results regarding the hematocrit profiles of healthy and rigid RBCs in the daughter branches (see Fig. 2c), we would expect distinct RBC transport through such model in vitro networks. Besides these microfluidic investigations, the three-dimensionality of the CFL in vascular networks, and the underlying hydrodynamic mechanisms, have been highlighted in recent numerical investigations.^82–84 Here, our results provide experimental inside into the effect of RBC deformability on hemodynamic alterations such as hematocrit skewness and disproportionate partitioning along sequential bifurcations within a network. Moreover, the reported differences in the CFL development and spatial organization can have implications in the context of microfluidic designs. Particularly, it has been predicted in numerical works how the shape and properties of RBCs, droplets, and capsules affect their flow paths and migration in confined microchannel flows.^85–87

5. Conclusions

Studying the effect of RBC rigidity in microfluidic flows advances our knowledge of blood flow in vivo and is crucial to understand the impact of pathological RBC changes on their flow properties, especially in complex geometries, such as bifurcating vessels and networks. In this study, we performed microfluidic measurements at various positions along the channel flow direction in a bifurcating microchannel to understand the RBC flow behavior, their lateral organization across the channel width, and the CFL phenomenon covering a broad hematocrit and Re range of 0.1–5% Ht and Re = 6–24, respectively.

In the mother channel, we observe different RBC focusing patterns that result in the emergence of two peaks close to the channel walls for healthy cells, while rigid RBCs predominantly flow along the channel centerline. Arriving at the bifurcation, these differences lead to distinct RBC distributions in the daughter vessels for rigid and healthy RBCs, which persist until the end of the channel. Our microfluidic in vitro experiments on artificially hardened RBC with GA show that the partition at the level of the bifurcation depends strongly on RBC deformability. Further, we reveal how the bifurcation affects the development of the different CFLs in the daughter branches. Our numerical simulations further demonstrate that the length to reach the steady-state CFL depends on cell rigidity. Since the distribution of RBCs after a bifurcation has been shown to be a determining factor for the distribution of RBCs in the microcirculatory network,³⁶ the higher concentration close to the lower daughter walls demonstrated by our study for rigid cells is a potential mechanism through which blood flow can be modified by pathologically stiffened cells. Our findings highlight the importance of understanding the influence of RBC rigidity on partitioning in complex flow fields, and thus oxygen delivery in the microcirculatory network.

Data availability statement

All study data are included in this article and/or the Supplementary material.

Author contributions

Conceptualization, Y. R. and S. M. R.; methodology, data acquisition and analysis, Y. R., O. A., A. D., T. J. and S. M. R.; writing—original draft preparation, Y. R., A. D., O. A. and S. M. R.; writing—review and editing, Y. R., O. A., A. D., T. J., J. H., C. W. and S. M. R.; supervision, J. H., C. W. and S. M. R.; funding acquisition, J. H., C. W., and S. M. R.; all authors have read and agreed to the published version of the manuscript.

Conflicts of interest

The authors declare no conflict of interest. The funders had no role in the design of the study, in the collection, analyses, or interpretation of data, in the writing of the manuscript, or in the decision to publish the results.

Appendix A

The lattice Boltzmann method (LBM) is based on the discretized Boltzmann transport equation in time and phase space. We consider in this work the standard D3Q19 LBM-BGK model^40–42 that reads as

f_i(x + e_iΔt, t + Δt) − f_i(x, t) = Ω_i(x, t) + S_i(x, t)Δt, (1)

where f_i(x, t) is the single-particle distribution function at a position x and time t with a microscopic velocity e_i, where i = 1…19. Ω_i is the Bhatnagar–Gross–Krook collision operator (BGK) and is defined asΩ_i describes the relaxation of f_i toward its equilibrium f^eq_i within the relaxation time τ. f^eq_i corresponds to a truncated expansion of the Maxwell–Boltzmann distribution and reads asω_i are the lattice weights, which for a D3Q19 lattice read as 1/3, 1/18, and 1/36 for i = 1, i = 2…7, and i = 8…19, respectively. The kinematic viscosity of the fluid is related to the relaxation time with

ν = c_s²[τ − Δt/2], (4)

where is the lattice speed of sound. The macroscopic fluid density and velocity are calculated from the moments of the distribution functions such thatExternal and internal forces are added to eqn (1)via the source term S_i through the Shan and Chen coupling scheme.⁸⁸

The red blood cells (RBC) are modeled as hyperelastic biconcave capsules endowed with shear and bending resistance. We use the Skalak strain model to account for the RBC membrane's shear elasticity and local area dilatation.^43,44 The strain energy over the surface of the RBC (Ω_s) reads as

with I₁ = λ₁² + λ₂² − 2 and I₂ = λ₁²λ₂² − 1 are the two deformation invariants, λ₁ and λ₂ are the principal stretching ratios, and G_s is the shear elastic modulus. The area dilatation modulus G_A is related to G_S according to

G_A = G_S(1 + 2C), (7)

where C is a dimensionless number quantifying the local inextensibility of the membrane. The resulting elastic force is computed using a linear finite element method.⁴⁵ The Skalak model has been initially tailored for RBCs and used successfully to recover their dynamics and morphologies under different flow conditions.^46–49 The bending is accounted for using the Helfrich free energy developed for bilayer lipid membranes⁵⁰where κ₀, , and κ_G = Π_i=1²κ_i are the spontaneous, mean and Gaussian curvatures respectively. κ₁ and κ₂ are the principal curvatures, and κ_H and κ_G are the bending and Gaussian curvatures moduli. The bending force is obtained through the functional derivative of Helfrich free energy. The numerical evaluation of the different terms in the bending force is based on discrete differential geometry operators for triangulated meshes.^48,51 The impermeability and inextensibility of the RBC membrane are fulfilled by imposing penalty functions on both the global volume (V) and the area (A)Here K_V and K_A are free parameters; V₀ and A₀ are the reference volume and area of the undeformed RBC. The simulation parameters are based on the experimental ones.

Acknowledgements

This research was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – project number 349558021 (WA 1336/13-1, HA 4382/8-1, and RE 5025/1-2). Y. R. acknowledges funding by the Marie Skłodowska-Curie grant agreement no. 860436—EVIDENCE. A. D. acknowledges funding from the Young Investigator Grant of Saarland University.

Notes and references

1. A. R. Pries and T. W. Secomb, Microcirculation, Elsevier, 2008, pp. 3–36.
2. T. W. Secomb, Annu. Rev. Fluid Mech., 2017, 49, 443–461.
3. J. Stuart and G. Nash, Blood Rev., 1990, 4, 141–147.
4. A. Dondorp, M. Nyanoti, P. Kager, S. Mithwani, J. Vreeken and K. Marsh, Trans. R. Soc. Trop. Med. Hyg., 2002, 96, 282–286.
5. A. Symeonidis, G. Athanassiou, A. Psiroyannis, V. Kyriazopoulou, K. Kapatais-Zoumbos, Y. Missirlis and N. Zoumbos, Clin. Lab. Haematol., 2001, 23, 103–109.
6. R. Mannino, D. R. Myers, Y. Sakurai, R. E. Ware, G. Barabino and W. Lam, Blood, 2012, 120, 818.
7. F. Reichel, M. Kräter, K. Peikert, H. Glaß, P. Rosendahl, M. Herbig, A. Rivera Prieto, A. Kihm, G. Bosman, L. Kaestner, A. Hermann and J. Guck, Front. Physiol., 2022, 13, 1–11.
8. A. Rabe, A. Kihm, A. Darras, K. Peikert, G. Simionato, A. K. Dasanna, H. Glaß, J. Geisel, S. Quint, A. Danek, C. Wagner, D. A. Fedosov, A. Hermann and L. Kaestner, Biomolecules, 2021, 11, 727.
9. Y. C. Fung and B. W. Zweifach, Annu. Rev. Fluid Mech., 1971, 3, 189–210.
10. A. Pries, K. Ley, M. Claassen and P. Gaehtgens, Microvasc. Res., 1989, 38, 81–101.
11. Z. Shen, G. Coupier, B. Kaoui, B. Polack, J. Harting, C. Misbah and T. Podgorski, Microvasc. Res., 2016, 105, 40–46.
12. A. Mantegazza, F. Clavica and D. Obrist, Biomicrofluidics, 2020, 14, 014101.
13. R. Fåhraeus, Physiol. Rev., 1929, 9, 241–274.
14. R. Fåhræus and T. Lindqvist, Am. J. Physiol. Content, 1931, 96, 562–568.
15. M. Faivre, M. Abkarian, K. Bickraj and H. A. Stone, Biorheology, 2006, 43, 147–159.
16. S. Kim, P. K. Ong, O. Yalcin, M. Intaglietta and P. C. Johnson, Biorheology, 2009, 46, 181–189.
17. N. Xiang and Z. Ni, Biomed. Microdevices, 2015, 17, 1–11.
18. S. Tripathi, Y. V. B. Varun Kumar, A. Prabhakar, S. S. Joshi and A. Agrawal, J. Micromech. Microeng., 2015, 25, 083001.
19. V. Clavera, P. Connes, L. Lanotte, C. Renoux, P. Joly, R. Fort, A. Gauthier, C. Wagner and M. Abkarian, Front. Phys., 2021, 712.
20. D. A. Fedosov, B. Caswell, A. S. Popel and G. E. Karniadakis, Microcirculation, 2010, 17, 615–628.
21. C. Makena Hightower, B. Y. Salazar Vázquez, S. Woo Park, K. Sriram, J. Martini, O. Yalcin, A. G. Tsai, P. Cabrales, D. M. Tartakovsky and P. C. Johnson, et al. , Wiley Interdiscip. Rev.: Syst. Biol. Med., 2011, 3, 458–470.
22. S. Chien, S. Usami, R. J. Dellenback and M. I. Gregersen, Science, 1967, 157, 827–829.
23. S. Chien, Annu. Rev. Physiol., 1987, Vol. 49, 177–192.
24. L. Lanotte, J. Mauer, S. Mendez, D. A. Fedosov, J.-M. M. Fromental, V. Claveria, F. Nicoud, G. Gompper, M. Abkarian, J. Mauer, D. A. Fedosov, G. Gompper, S. Mendez, F. Nicoud and J.-M. M. Fromental, Proc. Natl. Acad. Sci. U. S. A., 2016, 113, E8207.
25. A. Passos, J. M. Sherwood, E. Kaliviotis, R. Agrawal, C. Pavesio and S. Balabani, Phys. Fluids, 2019, 31, 091903.
26. H. Fujiwara, T. Ishikawa, R. Lima, N. Matsuki, Y. Imai, H. Kaji, M. Nishizawa and T. Yamaguchi, J. Biomech., 2009, 42, 838–843.
27. A. Abay, S. M. Recktenwald, T. John, L. Kaestner and C. Wagner, Soft Matter, 2020, 16, 534–543.
28. B. Kaoui, T. Krüger and J. Harting, Soft Matter, 2012, 8, 9246.
29. G. Tomaiuolo, L. Lanotte, G. Ghigliotti, C. Misbah and S. Guido, Phys. Fluids, 2012, 24, 051903.
30. C. Iss, D. Midou, A. Moreau, D. Held, A. Charrier, S. Mendez, A. Viallat and E. Helfer, Soft Matter, 2019, 15, 2971–2980.
31. J. M. Sherwood, E. Kaliviotis, J. Dusting and S. Balabani, Biomech. Model. Mechanobiol., 2014, 13, 259–273.
32. E. Kaliviotis, J. M. Sherwood and S. Balabani, Sci. Rep., 2017, 7, 44563.
33. D. Bento, A. I. Pereira, J. Lima, J. M. Miranda and R. Lima, Comput. Methods Biomech. Biomed. Eng. Imaging Vis., 2017, 6, 1–9.
34. D. Bento, C. Fernandes, J. Miranda and R. Lima, Exp. Therm. Fluid Sci., 2019, 109, 109847.
35. K. Yamamoto, H. Abe, C. Miyoshi, H. Ogura and T. Hyakutake, J. Med. Biol. Eng., 2020, 40, 53–61.
36. Q. Zhou, J. Fidalgo, M. O. Bernabeu, M. S. Oliveira and T. Krüger, Soft Matter, 2021, 17, 3619–3633.
37. H. Amini, W. Lee and D. Di Carlo, Lab Chip, 2014, 14, 2739.
38. A. Abay, G. Simionato, R. Chachanidze, A. Bogdanova, L. Hertz, P. Bianchi, E. van den Akker, M. von Lindern, M. Leonetti, G. Minetti, C. Wagner and L. Kaestner, Front. Physiol., 2019, 10, 1–14.
39. J. Friend and L. Yeo, Biomicrofluidics, 2010, 4, 026502.
40. R. Benzi, S. Succi and M. Vergassola, Phys. Rep., 1992, 222, 145–197.
41. S. Succi, The lattice Boltzmann equation for fluid dynamics and beyond, Oxford University Press, 2001.
42. P. L. Bhatnagar, E. P. Gross and M. Krook, Phys. Rev., 1954, 94, 511.
43. R. Skalak, A. Tozeren, R. P. Zarda and S. Chien, Biophys. J., 1973, 13, 245–264.
44. R. Skalak, Biorheology, 1973, 10, 229–238.
45. T. Krüger, F. Varnik and D. Raabe, Comput. Math. Appl., 2011, 61, 3485–3505.
46. T. Omori, H. Hosaka, Y. Imai, T. Yamaguchi and T. Ishikawa, Phys. Rev. E, 2014, 89, 013008.
47. M. Gross, T. Krüger and F. Varnik, Soft Matter, 2014, 10, 4360–4372.
48. K. Sinha and M. D. Graham, Phys. Rev. E, 2015, 92, 042710.
49. O. Aouane, A. Scagliarini and J. Harting, J. Fluid Mech., 2021, 911, A11.
50. W. Helfrich, Z. Naturforsch., C: J. Biosci., 1973, 28, 693–703.
51. M. Meyer, M. Desbrun, P. Schröder and A. H. Barr, Visualization and mathematics III, Springer, 2003, pp. 35–57.
52. C. S. Peskin, Acta Numerica, 2002, 11, 479–517.
53. D. Di Carlo, Lab Chip, 2009, 9, 3038.
54. M. Brust, C. Schaefer, R. Doerr, L. Pan, M. Garcia, P. E. Arratia and C. Wagner, Phys. Rev. Lett., 2013, 110, 078305.
55. Q. Zhou, J. Fidalgo, L. Calvi, M. O. Bernabeu, P. R. Hoskins, M. S. Oliveira and T. Krüger, Biophys. J., 2020, 118, 2561–2573.
56. G. R. Cokelet and H. L. Goldsmith, Circ. Res., 1991, 68, 1–17.
57. S. M. Recktenwald, K. Graessel, Y. Rashidi, J. N. Steuer, T. John, S. Gekle and C. Wagner, Phys. Rev. Fluids, 2023, 8, 074202.
58. G. Segré and A. Silberberg, Nature, 1961, 189, 209–210.
59. S. C. Hur, H. T. K. Tse and D. Di Carlo, Lab Chip, 2010, 10, 274–280.
60. A. T. Ciftlik, M. Ettori and M. A. M. Gijs, Small, 2013, 9, 2764–2773.
61. J. Zhou and I. Papautsky, Lab Chip, 2013, 13, 1121.
62. T. M. Geislinger and T. Franke, Adv. Colloid Interface Sci., 2014, 208, 161–176.
63. S. Tanaka and M. Sugihara-Seki, J. Phys. Soc. Jpn., 2022, 91, 083401.
64. M. Sugihara-Seki and N. Takinouchi, Micromachines, 2021, 12, 1175.
65. D. Oh, S. Ii and S. Takagi, Fluids, 2022, 7, 96.
66. M. Masaeli, E. Sollier, H. Amini, W. Mao, K. Camacho, N. Doshi, S. Mitragotri, A. Alexeev and D. Di Carlo, Phys. Rev. X, 2012, 2, 031017.
67. Y. Chen, D. Li, Y. Li, J. Wan, J. Li and H. Chen, Sci. Rep., 2017, 7, 1–9.
68. S. R. Bazaz, A. Mashhadian, A. Ehsani, S. C. Saha, T. Krüger and M. E. Warkiani, Lab Chip, 2020, 20, 1023–1048.
69. S. M. Recktenwald, K. Graessel, F. M. Maurer, T. John, S. Gekle and C. Wagner, Biophys. J., 2022, 121, 23–36.
70. S. Chien, C. D. Tvetenstrand, M. A. Epstein and G. W. Schmid-Schonbein, Am. J. Physiol. Circ. Physiol., 1985, 248, H568–H576.
71. B. M. Fenton, R. T. Carr and G. R. Cokelet, Microvasc. Res., 1985, 29, 103–126.
72. A. Merlo, M. Berg, P. Duru, F. Risso, Y. Davit and S. Lorthois, Soft Matter, 2022, 18, 1463–1478.
73. D. Katanov, G. Gompper and D. A. Fedosov, Microvasc. Res., 2015, 99, 57–66.
74. Q. M. Qi and E. S. G. Shaqfeh, Phys. Rev. Fluids, 2017, 2, 093102.
75. J. M. Sherwood, J. Dusting, E. Kaliviotis and S. Balabani, Biomicrofluidics, 2012, 6, 024119.
76. J. M. Sherwood, D. Holmes, E. Kaliviotis and S. Balabani, PLoS One, 2014, 9, e100473.
77. E. Kaliviotis, D. Pasias, J. Sherwood and S. Balabani, Med. Eng. Phys., 2017, 48, 23–30.
78. E. Kaliviotis, J. M. Sherwood and S. Balabani, Phys. Fluids, 2018, 30, 030706.
79. X. Yin, T. Thomas and J. Zhang, Microvasc. Res., 2013, 89, 47–56.
80. J. Zhang, P. C. Johnson and A. S. Popel, Microvasc. Res., 2009, 77, 265–272.
81. G. Li, T. Ye, B. Yang, S. Wang and X. Li, Phys. Fluids, 2023, 35, 021906.
82. A. Pskowski, P. Bagchi and J. D. Zahn, Biomicrofluidics, 2022, 16, 064104.
83. P. Balogh and P. Bagchi, J. Fluid Mech., 2019, 864, 768–806.
84. S. Ebrahimi and P. Bagchi, Sci. Rep., 2022, 12, 4304.
85. S. H. Bryngelson and J. B. Freund, Rheol. Acta, 2016, 55, 451–464.
86. R. Lu, Z. Wang, A.-V. Salsac, D. Barthès-Biesel, W. Wang and Y. Sui, J. Fluid Mech., 2021, 923, A11.
87. Y. Pang, S. Jiao, Y. Lu, Q. Yang, X. Wang and Z. Liu, Phys. Fluids, 2023, 35, 072013.
88. X. Shan and H. Chen, Phys. Rev. E, 1993, 47, 1815.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d3sm00517h

---