A two-phase core-plasma model for microvascular blood flow: Comparative analysis of hemodynamic models
Maya Salame, Marianne Fenech

TL;DR
This paper introduces a new two-phase model to better understand blood flow in small vessels by capturing the behavior of red blood cells and plasma.
Contribution
The Core-Plasma Model integrates Newtonian and non-Newtonian elements to accurately represent microvascular blood flow dynamics.
Findings
The model outperforms traditional approaches in capturing velocity and shear rate profiles in vitro.
It effectively represents the dynamic interplay between red blood cells and the cell-free layer near vessel walls.
Results are consistent across varying flow rates, hematocrit levels, and suspending media.
Abstract
Microcirculatory blood flow exhibits complex non-Newtonian behavior, including shear-thinning properties and the formation of a cell-free layer (CFL)—a plasma-rich region near vessel walls. While traditional rheological models such as Newtonian, Power Law, and Carreau describe certain flow characteristics, and empirical models like the double-parameter power fit have been used to capture velocity profiles, these approaches fall short in fully characterizing the dynamic interplay between red blood cells (RBCs) and plasma. This study introduces the Core-Plasma Model, a two-phase framework that integrates Newtonian and non-Newtonian elements to represent the RBC-rich core and surrounding CFL. In vitro experiments in 25 μm and 50 μm round channels across varying flow rates, hematocrit levels (5–20%), and suspending media (PBS and native plasma) demonstrate the model’s superior ability to…
Genes, proteins, chemicals, diseases, species, mutations and cell lines named across the full text — each resolved to its canonical identifier and authoritative record.
Click any figure to enlarge with its caption.
Fig 1
Fig 2
Fig 3
Fig 4
Fig 5
Fig 6
Fig 7
Fig 8
Fig 9
Fig 10
Fig 11- —http://dx.doi.org/10.13039/501100000038Natural Sciences and Engineering Research Council of Canada
- —http://dx.doi.org/10.13039/501100000038Natural Sciences and Engineering Research Council of Canada
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
TopicsBlood properties and coagulation · Cardiovascular Health and Disease Prevention · Trauma, Hemostasis, Coagulopathy, Resuscitation
Introduction
Blood is a complex suspension of cells in plasma, exhibiting rheological behavior that varies significantly with vessel size, flow conditions, and hematocrit (Ht). In larger vessels, such as arteries and veins, blood behaves approximately as a Newtonian fluid, since the influence of individual red blood cells (RBCs) is averaged out by the dominant plasma flow [1,2]. However, in the microcirculation—where vessel diameters drop below 100 m—the discrete nature of RBCs, along with their interactions with plasma and vessel walls, becomes much more pronounced. These microscale effects give rise to non-Newtonian behavior that cannot be ignored in accurate models of microvascular blood flow [3].
One of the most physiologically important phenomena in microvascular blood flow is the formation of the cell-free layer (CFL)—a plasma-rich region near the vessel walls depleted of RBCs. The CFL plays a critical role in reducing flow resistance, improving oxygen delivery, and influencing wall shear stress [2,4]. Its formation is governed by the interplay of hydrodynamic forces, RBC aggregation, and shear-induced migration [5–7]. Recent microfluidic investigations have revealed how RBC deformability, hematocrit, and channel geometry contribute to CFL dynamics under physiologically relevant shear conditions, particularly within lab-on-a-chip models that reproduce in vivo capillary environments [8,9]. As blood flows through narrow vessels, RBCs tend to migrate toward the center, creating a distinct core-plasma separation.
Blood viscosity, defined as the resistance to deformation under shear, is central to describing flow dynamics. Importantly, blood is shear-thinning—its viscosity decreases with increasing shear rate due to the breakup of RBC aggregates and alignment of cells with the flow direction [10,11]. This non-Newtonian behavior is further modulated by the Fåhraeus-Lindqvist effect, where apparent viscosity decreases with vessel diameter as the CFL becomes more pronounced [4,12].
Accurate modeling of blood rheology in the microcirculation is essential for predicting flow behavior and understanding microvascular physiology and pathology. Reports of altered microcirculatory perfusion in conditions such as sepsis underscore the clinical relevance of microscale hemorheology [13]. Classical models often fail to fully capture the complexity of blood flow at this scale, particularly under varying hematocrit levels and shear conditions. Contemporary studies emphasize the need for hybrid approaches combining experimental micro-PIV and computational fluid dynamics to quantify velocity fields and local viscosity distributions in confined geometries [8,9]. This study aims to address this gap by systematically evaluating commonly used blood rheology models and proposing an improved framework that explicitly incorporates the impact of CFL formation.
The following section introduces the key rheological and fitting models used to represent blood flow in microvessels. These models range from simple Newtonian approximations to more advanced non-Newtonian formulations, each with varying levels of complexity and physiological relevance.
Blood flow modeling
Several models have been developed to characterize blood rheology in the microcirculation. Each model offers distinct advantages and limitations based on its ability to capture shear-thinning behavior and microstructural effects such as RBC aggregation and CFL formation.
Newtonian model.
The Newtonian model assumes constant viscosity and a linear relationship between shear stress and shear rate. The velocity profile in a cylindrical tube is given by:
where is the pressure drop, is fluid viscosity, L is the vessel length, R is the vessel radius, and r is the radial position from the center of the vessel. Although simple, this model fails to capture the shear-thinning and non-Newtonian effects observed in the microcirculation [11].
Power law model.
The Power Law model introduces non-Newtonian behavior by relating viscosity to shear rate:
where Kp is the consistency index and np is the flow index. This model effectively captures shear-thinning at moderate shear rates, where blood viscosity decreases as shear rate increases. However, it does not account for the viscosity plateau observed at low shear rates, where RBC aggregation increases viscosity, or at high shear rates, where the minimum viscosity is reached. Instead, the viscosity tends to infinity as shear rate approaches zero and to zero at high shear rates—an unrealistic behavior in physiological conditions [11].
Carreau model.
The Carreau model addresses the limitations of the Power Law by modeling viscosity over a wide range of shear rates:
where and are viscosities at zero and infinite shear rates, is the time constant, and nc is the power index. This model is widely used in hemorheology for its ability to represent the entire shear-thinning spectrum of blood [11]. However, it does not explicitly account for microscale phenomena such as core–plasma phase separation or the formation of a CFL near vessel walls.
Double-parameter power fit.
To improve velocity profiling in microcirculatory flow, a Double-Parameter Power (DPP) Fit was introduced by Koutsiaris et al. (2009) [14]. Unlike traditional rheological models that derive velocity from fluid properties such as viscosity, this approach is a purely empirical fitting model designed to match experimentally observed velocity distributions.
where k1 controls the velocity gradient near the wall and k2 controls the bluntness of the core. Although it cannot predict viscosity, it effectively characterizes flow shapes, particularly in the presence of a CFL [15].
Core-plasma model.
Building on these foundations, this study introduces a Core-Plasma Model that represents blood flow as a two-phase system: a non-Newtonian, RBC-rich core surrounded by a Newtonian plasma layer near the vessel walls. By explicitly capturing phase separation and presence of a CFL, the model aims to improve predictions of shear-dependent viscosity changes and flow behavior in microvessels, offering greater physiological relevance than traditional single-phase models. The model is validated against experimental microfluidic data collected across varying channel sizes, hematocrits, and suspension mediums. Ultimately, this approach seeks to refine our understanding of microcirculatory blood flow and enhance predictive models for biomedical applications.
Materials and methods
Fluid preparation
Fresh human blood was collected from a single adult donor under ethics approval from the University of Ottawa (H-03-19-3441). Written informed consent was obtained from the donor prior to participation, as documented using a university-approved consent form. The sample used in this study was collected in March 2024. RBCs were suspended in either phosphate-buffered saline (PBS) or native plasma. The same donor was used, and samples were collected separately for testing the 25 m and 50 m channels. For each test, PBS suspensions were tested on day 1 and plasma suspensions were tested on day 2.
For PBS suspensions, whole blood was centrifuged to remove plasma and buffy coat. The isolated RBCs underwent three PBS washing cycles to eliminate residual plasma proteins, platelets, and white blood cells. RBCs were then re-suspended in PBS supplemented with 0.9 mg/mL glucose to support cell viability and 31.5% OptiPrep® (Sigma Aldrich, Ref. D1556) to stabilize suspension density and minimize sedimentation. Hematocrit (Ht) levels were adjusted to 5%, 10%, 15%, and 20%.
For plasma suspensions, RBCs were similarly washed and re-suspended in density-adjusted native plasma, modified by mixing with dried OptiPrep®. This preparation preserved plasma proteins to promote RBC aggregation. Fluorescent tracer particles (0.87 m, FluoroMax, Thermo Fisher) were added to each sample to enable PIV measurements.
Experimental setup
The experimental setup, shown schematically in Fig 1, was designed to simultaneously measure velocity profiles and CFL thickness in microchannels. Microchannel Fabrication: Microchannels were fabricated using borosilicate glass capillaries (25–50 m diameter) to model microvessels, following the method developed by Chartrand et al. (2023) [16]. The channels were coated with Poly(L-Lysine)-PEG methyl ether (PLL-PEG) to minimize red blood cell (RBC) adhesion to the glass surface. Measurement Acquisition: Pressure-driven flow was precisely controlled and monitored using a flow controller and sensor system. Velocity profiles were measured using micro-particle image velocimetry ( PIV) with fluorescent tracers, while CFL thickness was captured in parallel using high-speed imaging. Flow rates were adjusted to span a physiological yet broad shear-rate range (118–4200 s^−1^), using channel dimensions and a nominal reference viscosity of 3.5 cP only for preliminary setpoints [11]. For each hematocrit (5%, 10%, 15%, 20%) and suspending medium (PBS or plasma), data were acquired at ten pressure setpoints (20–200 mbar), yielding ten paired measurements (one PIV sequence and one CFL recording) per condition. The apparent viscosity was then computed for each condition from measured pressure–flow data via the laminar Poiseuille relation and used for all analyses (detailed in S5 File). Imaging and PIV were performed in the central section of the capillary after the pressure and volumetric flow rate reached a stable plateau. The CFL was recorded only under these steady conditions; no progressive drift in CFL thickness was observed during acquisition.
Experimental setup for simultaneous velocity and CFL measurements.Pressure-driven flow was regulated using a Fluigent Flow EZ system, with real-time flow monitoring via an in-line flow sensor. The microchannel, a 25–50 μm diameter borosilicate glass capillary, was mounted on an inverted microscope equipped with a dual-camera port. A 532 nm laser illuminated fluorescent tracer particles for velocity measurements using a LaVision FlowMaster micro-particle image velocimetry (μPIV) system and a CCD camera. Simultaneously, a high-speed camera (HSC) captured red blood cell (RBC) distributions at 100 fps with a 20× objective for cell-free layer (CFL) analysis.
Image post-processing: At each flow rate, high-speed images were analyzed with a MATLAB application implementing the gradient-based spatiotemporal method of Fenech et al. (2023) to detect the RBC core and compute the optical CFL thickness ( ) [17]. For velocity post-processing, image pairs were captured at each flow rate and analyzed using Nguyen’s cross-correlation algorithm with adaptive interrogation windows [18]. Velocity profiles were averaged, and noisy data were filtered based on standard deviation to generate reliable two-dimensional flow fields. The velocity acquisition process is illustrated in Fig 2.
μPIV velocity profile acquisition workflow.Step-by-step process used to extract velocity profiles from microfluidic blood flow experiments using μPIV. (1) Double-pulse images: Sequential images captured using a high-speed camera with fluorescent tracer particles suspended in the flow. (2) Cross-correlation: The Nguyen correlation algorithm [18] is applied to determine the displacement of tracer particles between image pairs. (3) Vector plot: Velocity vectors are generated based on particle displacement, visualizing flow direction and magnitude. (4) Post-processed image: The velocity field is refined by averaging multiple image pairs, with velocity data (red) and standard deviation (green) plotted. Regions where the standard deviation is zero indicate reliable velocity measurements. (5) Velocity profile: Truncated and extracted velocity profile is displayed, showing the parabolic distribution typical of microchannel flow, with a peak at the centerline (red line).
Mathematical analysis of velocity profiles
The developed core–plasma model.
The Core–Plasma Model captures the biphasic nature of blood by dividing the flow domain into two distinct regions: a Newtonian CFL and a non-Newtonian RBC-rich core. This approach enables the integration of both Newtonian and shear-thinning behavior in a single analytical framework.
The governing assumptions include steady, incompressible, laminar, axisymmetric, and fully developed flow in a cylindrical tube, with body forces neglected. The domain is partitioned as follows:
Newtonian CFL layer : This region represents the plasma-rich CFL and is modeled as a Newtonian fluid with constant viscosity . The velocity profile is derived by solving the simplified Navier–Stokes equation, applying the no-slip boundary condition at the vessel wall (r = R) and continuity of shear stress at the core–CFL interface ( ).Non-Newtonian core : The RBC-rich core is modeled as a power-law fluid, where viscosity depends on local shear rate. The solution enforces the symmetry condition at the centerline (r = 0) and ensures continuity of velocity at the interface ( ).
This two-region analytical model provides a simplified yet representative approximation of microvascular blood flow, enabling parameter extraction (e.g., flow indices) under physiologically relevant conditions:
where represents the CFL thickness.
The shear rate for the Core–Plasma model is obtained by differentiating the velocity profile u(r) with respect to the radial position r, and is given by:
where:
is the pressure drop, is the Newtonian viscosity of plasma,L is the vessel length,r is the radial distance from the channel center,R is the total vessel radius, is the CFL thickness,Kcp and ncp are the consistency and flow indices governing the non-Newtonian behavior in the Core–Plasma (cp) model.
A velocity distribution, shown in Fig 3, illustrates an example of the Core-Plasma Model overlaid on experimental data, highlighting distinct regions in the RBC-rich core and the surrounding CFL.
Velocity profile illustrating the Core-Plasma Model in a 50 μm channel.Visualization of the Core-Plasma Model for a 15% hematocrit suspension in plasma flowing through a 50 μm microchannel under 180 mbar pressure. The left panel shows a high-speed image of RBC distribution at high contrast with overlaid core and CFL boundaries. The right panel plots the velocity profile extracted from micro-particle image velocimetry (μPIV) data, with segmented fits for the RBC core (red dashed) and CFL region (yellow dashed). The vertical dashed line indicates the optically measured cell-free layer thickness (δo).
Direct measurement of plasma viscosity ( ) is not feasible in our experimental system. Instead, we estimate it using the apparent viscosity ( ) calculated from Poiseuille’s Law. This estimated value is used as an initial parameter for optimization in the Core-Plasma Model. The approach aligns with the Newtonian assumption applied to the CFL, where plasma is the dominant fluid component.
Poiseuille’s Law assumes laminar flow of a Newtonian fluid with constant viscosity and provides a practical method for estimating in microcirculatory settings due to its analytical simplicity [19]. In this study, is computed from experimentally measured flow rate and pressure drop using the expression:
As detailed in S1 Fig, the pressure drop across the upstream tubing is negligible. The resulting serves as an effective estimate for and provides a baseline for comparing Newtonian and non-Newtonian flow behavior in confined geometries.
Shear rate and hydrodynamic CFL determination.
The shear rate is identified by locating the boundary between the CFL and the RBC core, where an intersection occurs. This intersection is marked by an abrupt change in shear rate ( ), reflecting the distinct rheological properties of the plasma-rich CFL, which behaves more like a Newtonian fluid, and the RBC-dense core, which exhibits non-Newtonian characteristics. The shear rate profile is approximated by two slopes: one originating at the vessel center (r = 0) and extending towards an unknown intersection point (ip), and another slope up to the vessel wall (r = R). The goal is to determine the intersection point, so we can identify the hydrodynamic CFL ( ). Mathematically, the shear rate ( ) is estimated as:
where:
m1 is the slope of the first linear region (RBC core),m2 is the slope of the second region (CFL), is the shear intersection point, representing the hydrodynamic RBC core boundary.
To determine , we minimize the squared error between the measured shear rate data and the piecewise linear model:
where represents the fitted shear rate values and the additional term ensures that the transition point remains within the physical boundary . Once is optimized, the hydrodynamic CFL thickness is computed as:
Shear rate at the intersection: The shear rate at the intersection ( ) quantifies the rate of deformation at the boundary between the RBC core and the CFL. It is defined as:
This value represents the final shear rate within the RBC core before transitioning into the CFL, marking the point where the flow dynamics shift from a non-Newtonian to a Newtonian regime. Fig 4 presents a representative shear rate profile, illustrating how linear fits differentiate the RBC core and CFL regions. The plots correspond to half of the microchannel width.
Shear rate at the intersection of a 15% Ht suspension in plasma within 50 μm and 25 μm microchannels.Shear rate profiles for a 15% hematocrit suspension in plasma under an applied pressure of 180 mbar, shown for a 50 μm (left) and 25 μm (right) microchannel. Linear fits are applied to the core and cell-free layer (CFL) regions, with dashed vertical lines indicating the optically measured (δo) and hydrodynamically estimated (δh) CFL boundaries. Black triangles represent the numerically derived shear rate obtained from velocity profile data collected via μPIV.
Results
Velocity profiles and model performance
Figs 5 and 6 illustrate representative velocity and shear rate distributions for a 15% hematocrit (Ht) suspension in plasma, flowing through 50 m and 25 m microchannels under an applied pressure of 180 mbar. The velocity profiles (Fig 5A and Fig 6A) exhibit slightly blunted shapes, with experimental PIV measurements showing strong agreement with theoretical model fits. In both cases, the boundaries of the CFL, denoted as and , are well-defined and closely aligned, with a notably thinner CFL observed in the 25 m channel.
Velocity and shear rate characteristics of a 15% hematocrit suspension in plasma within a 50 μm microchannel.Flow behavior under an applied pressure of 180 mbar is presented: (A) Velocity profile across the channel width, comparing experimental μPIV data (black triangles) with model fits. (B) Shear rate profile derived from the velocity data, with model fits applied to the core and CFL regions. (C) Root mean square (RMS) error of the velocity model fits across the radial position. (D) RMS error of the shear rate model fits across the radial position. Dashed vertical lines or stars represent optically measured (δo) and hydrodynamically estimated (δh) CFL boundaries.
Velocity and shear rate characteristics of a 15% hematocrit suspension in plasma within a 25 μm microchannel.Flow behavior under an applied pressure of 180 mbar is presented: (A) Velocity profile across the channel width, comparing experimental μPIV data (black triangles) with model fits. (B) Shear rate profile derived from the velocity data, with model fits applied to the core and CFL regions. (C) Root mean square (RMS) error of the velocity model fits across the radial position. (D) RMS error of the shear rate model fits across the radial position. Dashed vertical lines or stars represent optically measured (δo) and hydrodynamically estimated (δh) CFL boundaries.
Model performance is evaluated using the root mean square (RMS) error between experimental data and fitted profiles. Velocity RMS errors (Figs 5C and 6C) remain minimal in the core region but increase near the CFL boundary, reflecting the greater sensitivity of the velocity gradient in these regions. Shear rate distributions (Figs 5B and 6B) reveal a distinct intersection between the RBC-rich core and the surrounding CFL, particularly pronounced in the 50 m channel. The RMS errors in shear rate (Figs 5D and 6D) are elevated near the CFL boundary in both channel sizes, with more pronounced discrepancies observed in the 25 m channel.
A summary of velocity and shear RMS errors is provided in Fig 7. The Core-Plasma model shows the most agreement with experimental data among the rheological models, especially in the 50 m channel, where phase separation effects such as the formation of a cell-free layer are more pronounced. Here, the Core-Plasma model demonstrates statistically significant reductions in both velocity and shear rate RMS errors compared to all other rheological models (p < 0.0001 in most comparisons).
Model comparison of RMS velocity and shear rate errors across microchannel sizes.Root mean square (RMS) values for velocity (left panels) and shear rate (right panels) are compared across five models—Carreau, Double-Parameter Power (DPP), Newtonian, Power Law, and Core-Plasma—using blood suspensions at hematocrit levels of 5%, 10%, 15%, and 20% in both PBS and plasma. Data are grouped by microchannel size: 50 μm (top row) and 25 μm (bottom row). Box plots show the distribution of RMS values across all tested conditions. Statistical comparisons between models were performed using one-way ANOVA with Sidak’s multiple comparisons test. Significance levels are indicated as follows: ns = not significant, * p<0.05, ** p<0.01, *** p<0.001, **** p<0.0001.
The Carreau model consistently exhibits the highest RMS errors, with significant differences observed against all other models in the 50 m channel for both velocity and shear rate data (p < 0.0001), confirming its limitations in accurately describing non-Newtonian blood behavior under confined geometries. The Newtonian and Power Law models yield intermediate performance, showing moderate error levels without significant differences between them in most comparisons, especially in the 25 m channel where performance across models tends to converge.
Interestingly, model performance diverges more sharply in the 50 m channel, suggesting that larger channel diameters amplify differences in each model’s ability to capture flow profiles, likely due to enhanced shear gradients and cell-free layer formation. In the 25 m channel, however, no significant differences were observed between models in shear rate RMS errors, indicating that at small scales, reduced separation between RBC-rich and plasma-rich regions may diminish the influence of complex rheological modeling.
As for the purely empirical fitting, the Double-Parameter Power (DPP) fit shows the lowest RMS velocity errors across most conditions. Statistically, the DPP model performs significantly better than Carreau and Newtonian models in the 25 m channel (with p < 0.05 and p < 0.01, respectively), and is not significantly different from Core-Plasma Model.
To further investigate the performance of the Core-Plasma model, we segmented the fitted domain into two distinct regions: the RBC-rich core and the surrounding CFL, and independently evaluated the model’s error in each. As shown in Fig 8, both velocity RMS and shear rate RMS errors were significantly lower in the core region compared to the CFL across both microchannel sizes (50 m and 25 m). In the 50 m channel, the velocity RMS in the core was markedly lower than in the CFL (p < 0.0001), and a similar trend was observed for the shear rate RMS (p < 0.001). This pattern held true in the more confined 25 m channel as well, where the core exhibited significantly lower velocity RMS (p < 0.001) and shear rate RMS (p < 0.01) values than the CFL.
*Regional error analysis of the Core-Plasma model in 50 μm and 25 μm microchannels.Velocity (left) and shear rate (right) root mean square (RMS) errors are compared between the core and cell-free layer (CFL) regions of the Core-Plasma model for both 50 μm (top row) and 25 μm (bottom row) channels. Errors were significantly lower in the core region across all cases, indicating better model agreement in the RBC-rich central zone compared to the CFL. Statistical analysis was performed using paired t-tests; significance levels are indicated as follows: **p<0.01, ***p<0.001, ***p<0.0001.
Characterization of non-Newtonian parameters
Non-Newtonian parameters were extracted by fitting experimental velocity profiles to several rheological models, as summarized in Table 1. Results are reported as mean (standard deviation), with standard deviations shown in parentheses, separately for PBS and plasma suspensions in both 25 m and 50 m microchannels, allowing direct comparison of channel size and suspending medium effects. Literature values are included to benchmark physiological relevance and model fidelity.
Table 1: Averaged Non-Newtonian parameters by channel size and suspension medium, with literature comparison.Values represent the mean (standard deviation) across all tested flow rates (20–200 mbar) and hematocrit levels (5–20%) for both plasma and PBS suspensions, separated by microchannel size. Individual data points for each condition are provided in the Supplementary Information.
Apparent viscosity, estimated from Poiseuille flow assumptions, was 1.98 (0.21) mPa s in PBS and 1.67 (0.38) mPa s in plasma suspensions (25 m channel), with similar trends observed in 50 m channels. Plasma-suspended fluid yielded lower viscosity values. All values fall within the physiological range of 1.50–2.00 mPa s reported by Wells et al. (1962) [20]. The pressure gradient remained consistent across conditions at approximately Pa/m. Fig 9 shows the dependence of apparent viscosity on pressure gradient fitted with a power-law model on log–log axes.
Apparent viscosity as a function of pressure gradient in microchannels.Apparent viscosity, derived from Poiseuille flow assumptions, is plotted as a function of pressure gradient for 50 μm (left) and 25 μm (right) microchannels. Data are shown across different hematocrit levels and suspension types (PBS and plasma). Fitted with a power-law model on log–log axes.
The infinite-shear viscosities ( ) range from 1.55 (0.42) to 1.86 (0.22) mPa s and closely aligned with previously reported ranges (1.6–2.3 mPa s) [10]. The variation of with experimental pressure gradient across different hematocrits and suspensions is provided in S4 Figs, where plasma suspensions have a lower , particularly in the 25 m channel. The time constant ( ) remained almost identical across all cases (3.31 s with negligible variation). In contrast, zero-shear viscosity ( ) showed considerable variability—particularly in PBS suspensions and the narrower 25 m channel—ranging from 7.27 (5.39) to 13.2 (18.4) mPa s. These lower values compared to Mehri et al. (26.9–118.6 mPa s) may reflect differences in hematocrit, channel confinement, and estimation method. This variability of with experimental pressure gradient across different hematocrits and suspensions is provided in S4 Figs. Fig 10 displays the Carreau-modeled local viscosity profiles as a function of shear rate. The 50 m channel shows a steeper initial drop in viscosity, while the 25 m channel displays an earlier plateau, indicating reduced shear variation due to geometric confinement.
Local viscosity from the Carreau Model as a function of shear rate.Local viscosity is plotted against local shear rate for a 15% hematocrit suspension in plasma flowing through 50 μm (left) and 25 μm (right) microchannels under an applied pressure of 180 mbar. Curves are fitted using the Carreau model with condition-specific rheological parameters.
The Carreau flow behavior index (nc) ranged from 0.354 (0.009) to 0.413 (0.177), with higher values observed in PBS and larger channels, indicating reduced shear-thinning under these conditions. These values are consistent with the reported physiological range (0.353–0.369) [10].
The Power Law model consistently underestimated non-Newtonian effects, producing flow indices (np) close to 1 across all conditions (0.931 (0.016) to 0.950 (0.014)), suggesting nearly Newtonian behavior. Its consistency index (Kp) was considerably lower than the literature range (9.9–62.2 mPa·s*^n^), with the lowest values observed in plasma suspensions. The behaviour of n_p_ and K_p_ with pressure gradient across different hematocrits and suspensions is provided in S4 Figs, where n_p_ remains near unity and K_p*_ shows no clear trend.
The Double-Parameter Power (DPP) model effectively captured the shape of velocity profiles, particularly the bluntness at the core and near the walls. Wall shape factors (k1) ranged from 0.374 (0.238) to 0.867 (0.494), while core shape factors (k2) reached 3.96 (2.56) in plasma and wide channels. These values fall within the range reported by Pitts et al. (1.74–3.66) [15], confirming the model’s sensitivity to both geometry and fluid phase. The behaviour of k1 and k2 with pressure gradient across different hematocrits and suspensions is provided in S4 Figs, where k1 generally increases with in narrow channels, while k2 shows a decreasing trend across most conditions.
Lastly, the Core-Plasma model produced flow indices (ncp) ranging from 0.976 (0.278) in plasma to 1.17 (0.362) in PBS (25 m), indicating nearly Newtonian behavior in the RBC-rich core. The consistency index (Kcp) was slightly elevated in PBS (0.00310 (0.00519)) compared to plasma (0.00263 (0.00388)), supporting the trend of higher effective viscosity in PBS suspensions. While this model aligns well with observed velocity profiles, the relatively large standard deviations in Kcp values—particularly for plasma—highlights the impact of the CFL that is often neglected in blood modeling. This behaviour of ncp and Kcp with pressure gradient across different hematocrits and suspensions is provided in S4 Figs, where variation is prominent.
The flow behavior index (ncp) varied between 0.588 and 2.58 across all conditions. For PBS suspensions, ncp ranged from 0.674 (5% Ht, 160 mbar) to 2.58 (10% Ht, 160 mbar). For plasma suspensions, values ranged from 0.447 (10% Ht, 120 mbar) to 2.09 (10% Ht, 60 mbar). The consistency index (Kcp) exhibited a wide distribution. In PBS, values ranged from (5% Ht, 100 mbar) to (5% Ht, 160 mbar). In plasma, Kcp ranged from (10% Ht, 60 mbar) to (10% Ht, 120 mbar). Across both suspensions, variation in ncp and Kcp was observed across pressure and hematocrit levels, with no consistent directional trend.
Cell-free layer
Fig 11 presents a statistical comparison of CFL thickness measurements, confirming that is significantly greater than under all experimental conditions. A one-way ANOVA followed by Sidak’s multiple comparison test is performed (as detailed in S2 Text). In the 50 m channel, the median values of are 7.89 m in PBS and 6.91 m in plasma, while corresponding values are 4.76 m and 5.09 m, respectively. In the 25 m channel, CFL thickness is markedly reduced, with median values of 1.12 m (PBS) and 1.15 m (plasma), and values of 0.524 m and 0.825 m, respectively.
*Comparison of cell-free layer thickness across suspensions, channel sizes, and measurement methods.Cell-free layer (CFL) thickness measurements—optical (δo) and hydrodynamic (δh)—are shown for both PBS (blue) and plasma (yellow) suspensions at hematocrit levels of 5%, 10%, 15%, and 20% in 50 μm and 25 μm microchannels. Box plots display the distribution of CFL values grouped by measurement method and channel size. Statistical analysis was performed using one-way ANOVA with Sidak’s multiple comparisons test. Significance levels: **p<0.01, ***p<0.0001.
All comparisons show strong statistical significance (p < 0.0001), with the exception of the 25 m plasma condition, which remains significant but with a lower confidence level (p = 0.0065). These findings underscore the influence of both channel size and suspending medium on CFL measurements, with the 50 m channel exhibiting more pronounced differences between and than the 25 m channel.
Discussion
Strengths and limitations of rheological models and fits
Among the physically grounded models, the Core-Plasma model offers the best overall agreement with experimental data. As shown in Fig 7, the Core-Plasma model consistently yields lower RMS errors than the Newtonian, Power Law, and Carreau models, particularly in the 50 m channel where phase separation is more pronounced. This strong performance underscores the model’s strength in capturing two-phase flow behavior in microvascular environments. Notably, its ability to account for shear rate intersection between the core and the CFL contributes to its accuracy in reproducing experimental flow fields.
To better understand the model’s internal behavior, we separately analyzed RMS error in the core and CFL regions of the Core-Plasma model (Fig 8). Across both channel sizes, RMS errors were significantly lower in the core region than in the CFL. This pattern reflects the model’s more consistent performance in the RBC-rich zone, where the rheological behavior is relatively stable and dominated by shear-dependent viscosity. Conversely, the higher error observed in the CFL likely arises from increased local variability near the wall, including shear gradients, RBC migration, and potential artifacts in boundary detection. These findings suggest that refining the model’s treatment of the CFL, more accurate plasma viscosity estimation, or improved intersection modeling—could further enhance its predictive power.
Interestingly, the core region error in the 25 m channel was even lower than that observed in the 50 m channel, suggesting more RBC packing in highly confined geometries. This increased confinement may lessen lateral migration and flow fluctuations, making the velocity profile within the core more uniform and thus easier to model. However, the overall performance of the Core-Plasma model in the 25 m channel is slightly diminished due to the reduced CFL thickness and more ambiguous separation between phases, which challenges the two-phase assumptions on which the model is built.
The Power Law and Newtonian models provide intermediate fits. The Power Law model, although better suited than Newtonian assumptions, oversimplifies the system by treating blood as a homogeneous fluid with continuously varying viscosity, thereby ignoring discrete core–CFL dynamics. The Newtonian approximation, assuming constant viscosity, predictably underperforms in shear-thinning-dominated conditions but remains useful as a computationally efficient baseline for benchmarking.
The Carreau model consistently exhibits the poorest fit across both channel sizes. Its limited performance likely results from a mismatch between the model’s target shear range and the experimentally observed shear rates in this study. While the Carreau model is designed to span a wide spectrum of shear rates—from near-zero to extremely high values—the flow regimes probed here are predominantly in the mid-to-high shear range. As a result, the Carreau model fails to appropriately capture the shear-thinning behavior that defines blood flow in confined microchannels under physiological conditions.
The Double-Parameter Power (DPP) fit demonstrates the highest accuracy in capturing velocity profiles across both 25 m and 50 m microchannels, as illustrated in Fig 7. It effectively models the characteristic profile bluntness near the centerline, reflecting experimental observations in both RBC-rich and plasma-dominant regions. However, while the DPP model exhibits superior empirical performance, it remains a purely curve-fitting approach lacking a rheological foundation. This absence of physical grounding limits its predictive value under altered physiological conditions or nonstandard geometries, where extrapolation beyond the fitted regime becomes unreliable.
Interpretation of non-Newtonian parameters
The non-Newtonian parameters extracted from experimental velocity profiles reveal consistent and physiologically relevant trends across flow conditions, hematocrit levels, and suspending media. Apparent viscosity, estimated from Poiseuille flow assumptions, remained within the expected physiological range of 1.50–2.00 mPa s [20]. Plasma suspensions consistently yielded lower values than PBS across both channel sizes (Fig 9), which is likely due to the enhanced development of the CFL in plasma, reducing near-wall viscosity and overall flow resistance. This effect was especially prominent in the 50 m channel, where phase separation is more fully developed and pressure gradients are higher.
Although shear-thinning behavior is expected to be weak at hematocrits below 20% and at the high shear rates examined here, the Carreau model was applied to parameterize the viscosity–shear-rate relationship for comparison with literature data. With the Carreau model, the infinite-shear viscosity ( ) remained consistently between 1.55–1.86 mPa s, in agreement with values reported by Mehri et al. [10]. The time constant also showed near-identical values (3.31 s) across all conditions, indicating robust model fitting and consistency. However, the zero-shear viscosity ( ) exhibited substantial variability: higher values were consistently observed in PBS (up to 13.2 mPa s), particularly in the narrower 25 m channel, while lower values (as low as 7.27 mPa s) were associated with plasma suspensions. These results suggest that confinement and suspending medium jointly influence RBC aggregation. In plasma-RBC suspensions, inducing higher shear rate can enhance RBC deformability and hydrodynamic lift, supporting near-wall dispersion and a lower effective aggregation locally—producing thicker CFLs, as described in Salame et al. (2025); conversely, at low shear rate, aggregated RBCs are less prone to disaggregation and tend to maintain their structural integrity, thus reducing CFL thickness [22]. Compared to Mehri et al. (who reported values of 26.9–118.6 mPa s using wider rectangular channels and lower hematocrits), our lower values reflect both our use of cylindrical geometry and our focus on high-shear, physiologically relevant conditions (118–4200 s^−1^ vs. 5–35 s^−1^). The Carreau flow behavior index (nc) increased in PBS and in the 50 m channel, indicating a transition toward more Newtonian-like flow, consistent with reduced aggregation and more uniform RBC distribution at high shear [22]. A detailed analysis of the Carreau parameters is provided in S4 Figs.
In contrast, the Power Law model showed limitations in confined microvascular conditions. The flow behavior index (np) remained between 0.931–0.950 for all conditions, indicating an overestimation of Newtonian behavior and insufficient representation of the non-linear viscosity gradients. Consistency indices (Kp) were also notably low, especially in plasma suspensions, and fell well below the range reported by Mehri et al. (9.9–62.2 mPa·s*^n^*). These outcomes highlight the model’s inability to accommodate distinct phase separation and shear-dependent structuring observed experimentally in microchannels. A detailed analysis of the Power Law parameters are outlined in S4 Figs.
The Double-Parameter Power (DPP) model revealed flow structure behaviour by separating bluntness contributions from the wall and core regions. The wall bluntness parameter (k1) decreased with increasing pressure in the 50 m channel, indicating steeper velocity gradients near the wall as RBC aggregation was reduced under higher shear. In contrast, in the 25 m channel, k1 increased slightly with pressure, suggesting reduced wall shear gradients potentially due to geometric confinement and limited CFL development. The core bluntness parameter (k2) exhibited clear pressure-dependent behavior: in the 50 m channel, k2 decreased at higher pressure, consistent with more parabolic velocity profiles resulting from RBC disaggregation. However, in the 25 m channel, k2 increased under elevated pressure, reflecting a flatter velocity core, likely due to confinement-driven clustering of RBCs. This effect was most pronounced in plasma suspensions, where k2 reached values as high as 3.96 at low pressure before decreasing toward 1.65 with increasing pressure, indicating a transition from blunt to more parabolic flow. These patterns are consistent with prior findings by Pitts et al. [15] and highlight the DPP model’s sensitivity to confinement, shear rate, and suspension composition. A detailed analysis of the Double-Parameter Power Fit parameters is provided in S4 Figs.
The Core-Plasma model demonstrated the highest robustness and physical relevance across conditions. Flow indices (ncp) remained consistent in PBS (1.17 0.36) and plasma (0.98 0.28) in the 25 m channel, showing with near-Newtonian behavior in the RBC-rich core under high shear rate. Similar trends were observed in the 50 m channel. Values of ncp were often greater than 1 in PBS, especially at intermediate hematocrits (10%–15%) and pressures between 100–160 mbar. In plasma, ncp was generally lower, and frequently dropped below 1, indicating stronger shear-thinning and alignment with native plasma behavior. The consistency index (Kcp) varied significantly, especially in plasma suspensions. Some plasma conditions yielded extremely low values (e.g., at 60 mbar, 10% Ht), while others reached peaks above 10^−2^ or even 10^−1^ (e.g., at 120 mbar, 10% Ht). This large spread reflects spatial heterogeneity within the CFL, likely due to wall effects, RBC margination, or local flow fluctuations not fully captured by the two-phase approximation. By contrast, PBS suspensions exhibited more stable Kcp values across pressures and hematocrit levels, consistent with the absence of protein-mediated aggregation and reduced variability at the wall. No monotonic trends in either ncp or Kcp were observed with respect to pressure or hematocrit alone. Instead, the results point to an interdependent influence of geometric confinement, suspension composition, and flow rate in shaping the apparent rheology. The increased variability in Kcp, particularly in plasma and at intermediate pressures, suggests that improved modeling of the CFL—potentially through spatially resolved viscosity profiles—could further enhance the Core-Plasma framework. A detailed analysis of the Core-Plasma parameters are outlined in S3 Table. Taken together, these findings highlight the strength of the Core-Plasma model in capturing core behavior and phase separation in microvascular systems, while also pointing to limitations in its treatment of wall-adjacent regions. In future work, region-specific parameterization and direct measurement of local viscosity—particularly within the CFL—may help refine model accuracy for predictive applications in microcirculatory modeling and in vitro diagnostics.
Reliability of CFL measurements
The comparison between optical ( ) and hydrodynamic ( ) measurements of CFL thickness, presented in Fig 11, underscores both the strengths and limitations inherent to each approach. Optical measurements, obtained via high-speed imaging coupled with image processing, provide a more consistent and reliable representation of CFL thickness across experimental conditions. A key advantage of this method is the ability to visually inspect and verify the CFL boundary, offering an additional layer of validation that enhances confidence in the results.
In contrast, the hydrodynamic approach, which derives CFL thickness from PIV velocity profiles, systematically underestimates the CFL. This discrepancy likely stems from methodological limitations associated with PIV, including depth-of-correlation effects and out-of-plane motion, both of which compromise the precision of velocity gradient detection near the CFL–core interface [23,24]. Furthermore, the interface between the plasma-rich layer and the RBC-rich core is not sharply delineated; rather, it constitutes a transitional zone where RBCs intermittently migrate in and out. This dynamic behavior introduces uncertainty into velocity-based measurements, further reducing the reliability of hydrodynamic estimates.
While PIV remains an invaluable technique for quantifying velocity fields and shear rate distributions, these findings highlight its limitations in accurately resolving sharp spatial transitions such as the CFL boundary. For applications requiring precise measurement of CFL thickness in microvascular flow studies, optical methods should be prioritized due to their higher spatial resolution and reduced susceptibility to artifacts.
Conclusion
Among the rheological models evaluated, the Core-Plasma model exhibited the strongest predictive capability, effectively capturing the two-phase characteristics of microvascular blood flow, including the separation between the red blood cell (RBC)-rich core and the plasma-rich cell-free layer (CFL). Its ability to reproduce shear rate discontinuities near the channel wall, particularly in the 50 m geometry where phase separation is more pronounced, reinforces the relevance of two-phase modeling in confined microcirculatory environments.
Despite its strengths, the Core-Plasma model is limited by simplifying assumptions—most notably, the approximation of plasma viscosity using apparent viscosity values and the use of linear fitting approaches that may overlook localized velocity gradient complexities. Additionally, discrepancies between hydrodynamic CFL thickness ( ), derived from PIV velocity profiles, and optical CFL thickness ( ), derived from high-speed imaging, emphasize the need for caution when interpreting near-wall flow behavior based on velocity data alone.
The Double-Parameter Power (DPP) fit, while purely empirical, consistently provided the most accurate fits across all conditions, effectively capturing velocity profile bluntness and wall gradients. However, its lack of physical interpretability limits its utility for extrapolative or mechanistic modeling.
Future work should prioritize direct measurement of plasma viscosity, improved spatial resolution near CFL boundaries, and integration of more physiologically grounded parameters into model formulations. These enhancements would improve model accuracy, strengthen physiological relevance, and support applications in microfluidic design and investigation of pathological microvascular disorders.
This work contributes to hemorheology by establishing a mechanistic two-phase Core–Plasma model that bridges experimental microfluidic data and physical interpretation of blood flow at the microscale. While single-phase or purely empirical models often fail to represent the spatial heterogeneity of confined blood flow, the present approach explicitly accounts for the RBC-rich core and the plasma-dominant CFL, providing parameters that can be experimentally validated through PIV and high-speed imaging. Such a framework aligns with recent efforts in microcirculatory research to couple image-based velocity field measurements with predictive modeling tools that can describe flow partitioning and transport phenomena in lab-on-a-chip systems and small vessels. By quantifying how confinement, hematocrit, and plasma composition influence microvascular resistance, the Core–Plasma model advances current methods for evaluating microcirculatory function and offers a platform adaptable to diverse applications—from biomimetic microfluidic diagnostics to the mechanistic study of pathological flow alterations observed in inflammation and critical illness.
Supporting information
S1 FigPressure drop analysis.Pressure drop ( ) is plotted against flow rate (Q) for 50 m (top) and 25 m (bottom) microchannels using PBS and plasma suspensions. Imposed and corrected values are compared to validate pressure bounds and confirm that upstream tubing losses are negligible.(PDF)
S2 TextStatistical analysis methods.Description of Root Mean Square (RMS) error, one-way ANOVA, and Sidak multiple-comparison tests used to evaluate model fits and statistical significance.(PDF)
S3 TableCore–Plasma model parameters (ncp and Kcp) across experimental conditions.Fitted parameters for the Core–Plasma model showing flow behavior and consistency indices across applied pressures, hematocrit levels, and suspending media (PBS or plasma).(PDF)
S4 FigsCharacterization of non-Newtonian and fitting parameters.Supplementary figures summarizing model-derived viscosity and rheological parameters from the Carreau, Power Law, Core–Plasma, and Double-Parameter Power (DPP) models across hematocrit levels, suspensions, and channel diameters.(PDF)
S5 FileViscosity and flow analysis.Measured volumetric flow rates and corresponding Reynolds numbers for 25 m and 50 m microchannels. Estimated wall shear rates demonstrate that experimental conditions reproduce physiologically realistic microvascular flow regimes. Estimated viscosity of base fluids and apparent viscosity are provided.(PDF)
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Mittal R, Iaccarino G. Immersed boundary methods. Annu Rev Fluid Mech. 2005;37(1):239–61. doi: 10.1146/annurev.fluid.37.061903.175743 · doi ↗
- 2Fedosov DA, Caswell B, Popel AS, Karniadakis GE. Blood flow and cell-free layer in microvessels. Microcirculation. 2010;17(8):615–28. doi: 10.1111/j.1549-8719.2010.00056.x 21044216 PMC 3529161 · doi ↗ · pubmed ↗
- 3Hendricks N. The microcirculation. South Afr J Anaesth Analg. 2020:S 62–5. doi: 10.36303/sajaa.2020.26.6.s 3.2540 · doi ↗
- 4Fåhræ us R, Lindqvist T. The viscosity of the blood in narrow capillary tubes. American Journal of Physiology-Legacy Content. 1931;96(3):562–8. doi: 10.1152/ajplegacy.1931.96.3.562 · doi ↗
- 5Ji B, Yang Z, Feng J. Oil-coated bubble formation from submerged coaxial orifices. Phys Rev Fluids. 2021;6(3):033602.doi: 10.1103/physrevfluids.6.033602 · doi ↗
- 6Leighton D, Acrivos A. The shear-induced migration of particles in concentrated suspensions. J Fluid Mech. 1987;181:415–39. doi: 10.1017/s 0022112087002155 · doi ↗
- 7Liu Z-C, Adrian RJ, Hanratty TJ. Reynolds number similarity of orthogonal decomposition of the outer layer of turbulent wall flow. Physics of Fluids. 1994;6(8):2815–9. doi: 10.1063/1.868169 · doi ↗
- 8Lai B-J, Zhu L-T, Chen Z, Ouyang B, Luo Z-H. Review on blood flow dynamics in Lab-on-a-Chip systems: an engineering perspective. Chem Bio Eng. 2024;1(1):26–43. doi: 10.1021/cbe.3c 00014 39973974 PMC 11835182 · doi ↗ · pubmed ↗
