Blood is simulated with a number of RBCs suspended in a solvent. A single fluid for both external and internal solvents was used, so that particle reflections at the membrane surface were omitted. Turning off reflections of solvent particles at the RBC membranes results in a reduction of the computational cost by approximately two to three times. To approximate membrane impenetrability by solvent particles, the DPD repulsive force between them and the membrane vertices of RBCs was employed. This does not fully guarantee separation of the internal and external fluids by the RBC membrane, but it nearly eliminates membrane crossing by solvent particles. To identify the effect of such a simplification, we have run several simulations of blood flow in small tubes (D = 10 μm) with the separation of the internal and the external fluids (through reflections), and no significant differences in the simulated blood properties with those obtained using the simplified model were found.
Suspended RBCs assumed the diameter D0
= 7.82 μ
m, the Young’s modulus Y0
= 18.9 μ
N/m, the bending rigidity 2.4 × 10−19
J, and the membrane viscosity ηm
= 0.022 Pa·s, which were verified through a number of mechanical, rheological, and dynamics tests on a single RBC in References [13
]. We employed the RBC model with the following parameters:
= 2.2, n
= 2, ka
= 4,900, kd
= 100, and kv
= 5,000. The DPD interactions among different particle types (solvent [S] and wall [W] particles, RBC vertices [V]) are presented in .
Table 1 DPD simulation parameters in blood flow simulations. ‘S’ and ‘W’ denote solvent and wall particles, respectively, while ‘V’ corresponds to RBC vertices. The parameters a and γ are the DPD conservative (more ...)
The random force coefficients for different interactions were obtained using the energy unit kBT
= 0.0945 calculated according to the energy scale in equation (8)
. The RBC volume was assumed to be slightly larger than that of the triangulated network (V0
= 92.45), as there exist repulsive interactions between solvent particles and RBC vertices resulting in a thin solvent-free layer next to the RBC membrane. The effective RBC volume was obtained from an analysis of the three-dimensional solvent distribution next to the membrane in an equilibrium simulation, and was set to V′
= 105. Thus, the cell volume fraction or tube hematocrit was calculated as follows
is the number of RBCs in the volume Vt
is the tube radius, and L
is the tube length.
RBCs are represented by 500 DPD particles forming a triangulated network with the characteristic length of each edge (spring) to be approximately 0.54 μ
m. The solvent is simulated by a collection of free DPD particles that correspond to small fluid volumes of blood plasma. The fluid particles are “soft” volumes due to soft interparticle interactions, and their effective radius can be estimated to be approximately 100–150 nm [25
]. The simulated solvent provides hydrodynamic interactions among RBCs and between RBCs and the vessel walls. However, it does not apply any restriction on the simulated spatial resolution, as two RBCs are able to come closer than the effective diameter of the solvent particles. Therefore, this method can be successfully applied for blood flow simulations with relatively high hematocrit values such as Ht
Two types of EV interactions among cells were considered. The first case is shown in , where the repulsive force coefficient between membrane vertices aV–V is set to 100. This method introduces a non-zero screening length between two membrane surfaces governed by the cutoff radius of the repulsive interactions rc = 0.5. Hence, the choice of a smaller cutoff radius may result in overlapping of cells, whereas a larger one would increase the screening distance between cells, which may be unphysical and may strongly affect the simulation results at high volume fractions of RBCs. The second method to enforce EV interactions among cells employed reflections of RBC vertices on the membrane surfaces of other cells. The repulsive force coefficient in this case was set to aV–V = 2 yielding the screening length between two RBC surfaces to be virtually zero. These two methods of EV interactions will be called “repulsion” and “reflection,” respectively. In addition, several simulations employed a net repulsion of RBCs from the wall by setting the repulsive force coefficient between the wall particles and the cell vertices to aW–V = 30. They employ “reflection” EV interactions and will be denoted as “wall force” further in text. This was performed to minimize the effect of near wall density fluctuations of the suspending fluid, which influences simulation results of blood flow in small tubes (e.g., D = 10 μm). The relative performance of the three methods is discussed below.
The blood flow is driven by a uniform body force per unit mass (f) applied to both solvent and RBC particles in the flow direction. Such a force application is equivalent to the pressure gradient ΔP/L = ρf, where ΔP is the pressure drop over the tube length L and ρ is the suspension’s mass density. The wall is modeled by frozen DPD particles with the same radial distribution function as the fluid particles in combination with a specular reflection of the solvent and RBC particles at the fluid–wall interface. The number densities of both solvent and wall were set to nS = nW = 3. In all simulations, we assume periodic boundary conditions along the flow.
Blood Velocity Profiles
As an initial condition, RBCs were placed randomly within the tube. After the flow was turned on, the simulations were run for a time that was long enough to achieve steady state, which was characterized by a time independent fully developed velocity profile. The time to reach steady state is equal to approximately one second in physical units. After this time, the RBC migration is considered to be fully achieved and blood flow measurements are performed.
RBCs in Poiseuille flow migrate toward the tube center forming a core in the flow. shows several sample trajectories of individual RBCs (left) as flow develops and a snapshot of RBCs (right) flowing in a tube of diameter D = 20 μm after steady state is achieved. An RBC core formation is clearly observed with a thin plasma layer or CFL (not shown) next to the tube walls. The cell core of the flow results in plug-like velocity profiles shown in for tubes of diameters D = 10 μm (left) and D = 40 μm (right) for different Ht values.
Sample trajectories of individual RBCs (left) as flow develops and a snapshot of RBCs (right) in Poiseuille flow in a tube of a diameter D = 20 μm after steady state is achieved. Ht = 0.45. x is the coordinate along the cylinder axis.
Figure 2 Typical velocity profiles of blood flow in tubes of diameters D = 10 μm (left) and D = 40 μm (right) for different Ht values employing “repulsion” EV interactions. Dashed lines show the corresponding parabolic profiles (more ...)
Velocity profiles and subsequent number density distributions are shown over the half tube because of flow axisymmetry. Dashed lines in correspond to the parabolic profiles of the Newtonian plasma in absence of the cells for the same pressure gradients. Velocity profiles are plotted in physical units, where the time scale was defined in equation (7)
= 0.0012 Pa·s being the plasma viscosity. The velocity curves are averaged over the tube cross-section and over 2 × 105
time steps, which corresponds to the total time of 0.1 seconds. The pressure gradients employed were 2.633 × 105
and 6.582 × 104
Pa/m for tubes of diameters 10 and 40 μ
m, respectively. In the case of low Ht
(e.g., 0.15), the velocity profiles closely follow parabolic curves in the near-wall region. In the central region of the tube, a substantial reduction in velocity is found for all volume fractions in comparison with the parabolic profiles indicating a decrease in the flow rate given by
is the cross-sectional area. As the volume fraction of RBCs is increased, the flow rate decreases for the same pressure gradient. For small tube diameters (e.g., D
= 10 μ
m), we observe plug-flow profiles with a nearly flat velocity in the center, whereas for larger tubes the velocity profiles only slightly resemble plug-flow. Moreover, flat velocity profiles are extended to the tube walls for larger Ht
values indicating a wider RBC core and a smaller CFL. presents blood flow characteristics for different tube diametes and Ht
values. The table includes mean flow velocities defined as
, mean shear rates
, pressure gradients ΔP
, centerline flow velocities vc
, and the corresponding centerline velocities vp
of the Newtonian plasma in the absence of the cells. The mean flow velocities increase for tubes having larger diameters with comparable shear rates.
The typical shear rates employed in in vitro
and in vivo
experiments are above c
= 30 per second. The influence of shear rate on the apparent blood viscosity appears to be not significant for the average shear rates above c
. At shear rates considerably lower than c
, strong effects of shear rate on viscosity are expected [31
], which is attributed to RBC sedimentation and aggregation properties. The shear rates employed in simulations are high enough to allow quantitative comparison of the blood flow properties with the available experimental data.
The Fahraeus effect, discovered in in vitro
experiments of blood flow in glass tubes [10
], shows an increased value of discharge hematocrit (Hd
) measured at the tube exit in comparison with that of the cell suspension before the tube entrance. Hd
is defined as the volume fraction of RBCs exiting a tube per unit time. This effect is directly correlated with cell migration to the tube centerline. Hence, the RBC core is moving faster than the average blood velocity resulting in an increased Hd
value measured at the tube discharge. We define Hd
in simulations as follows
is the average cell velocity calculated by averaging over velocities of all cell vertices and also averaged in time after the stationary state is reached. Pries et al.
] have compiled a number of experiments on blood flow through glass tubes of different diameters at various Hd
values to obtain an empirical relation between Ht
where the tube diameter D
is in micormeters. compares Hd
values in simulations with those from experiments [28
] for different tube diameters and cell volume fractions. The left figure corresponds to the “repulsion” EV interactions among the cells with the parameters in that result in a small screening distance between two RBC membrane surfaces discussed in the Modeling parameters section. The right plot contains simulation results obtained with “reflection” EV interactions yielding essentially no screening length among RBCs. Agreement between simulations and empirical values is excellent for low Ht
values, whereas discrepancies are found for large Ht
and small tube diameters (D
= 10 μ
m). EV interactions appear to be of importance at Ht
= 0.45 as the cells are closely packed. The introduced screening length among the cells with “repulsion” EV interactions results in a wider RBC core compared with that using “reflection” EV interactions. This yields a smaller CFL and flow rate (shown in the next section) and consequently a smaller Hd
value as seen in .
Figure 5 Discharge hematocrits for various RBC volume fractions and tube diameters in comparison with the approximation in equation (12). “Repulsion” (left) and “reflection” (right) EV interactions are employed. Ht is the tube hematocrit. (more ...)
The remaining discrepancies for small tube diameters are due to solvent density fluctuations next to the wall as shown in . presents Hd values from simulations using the “wall force” method (see Modeling parameters section). Additional wall repulsion of RBCs reduces the effect of near-wall density fluctuations of blood plasma in tubes of small diameters yielding good agreement of the simulated Hd values with those in in vitro experiments. Note that the repulsive wall does not affect the simulation results for low Ht values and large tube diameters as the CFL is wider than the distance of the effective cell– wall interactions.
Hd for different RBC volume fractions and tube diameters obtained by the “wall force” method that utilize a net repulsion of cells from the wall. Ht is the tube hematocrit.
In conclusion, the “wall force” method yields superior results for Hd in comparison with the other two considered approaches. The “repulsion” method for EV interactions introduces a nonzero screening distance among RBCs, which may be unphysical. Cell–cell hydrodynamic interactions mediated by the simulated blood plasma are likely to be affected by this screening layer among the cells. The “reflection” approach eliminates the effect of the artificial screening layer to yield better results compared with the “repulsion” method. However, neither “repulsion” nor “reflection” method apply any treatments of the existing solvent density fluctuations near the walls. The “wall force” method extends the “reflection” approach to minimize the effect of density fluctuations and proves to be the best approach.
The CFL is a near-wall layer of blood plasma absent of RBCs (see ) since they are subject to migration to the tube center in Poiseuille flow. The fluid viscosity of the CFL region is much smaller than that of the tube core populated with RBCs providing an effective lubrication for the core to flow. The thickness of the CFL is directly related to the Fahraeus and the Fahraeus–Lindqvist effects. Thus, in small tubes, the CFL thickness is significant with respect to the tube diameter resulting in a smaller relative apparent viscosity and larger Hd in comparison with those in larger tubes, where the CFL thickness becomes negligible with respect to the tube diameter.
To determine the CFL thickness δ
, we measured the outer edge of the RBC core (see ), which is similar to CFL measurements in experiments [20
]. shows a sample CFL edge (left) from simulations and CFL thickness distribution (right) for Ht
= 0.45 and D
= 20 μ
m. The cell edge was measured by projecting the cell vertices of a snapshot of RBCs on the x
plane, where curves of the RBC core minimum and maximum were fitted. Discrete samples of δ
from the obtained curves were taken every 0.5 μ
m along the flow direction x
. In addition, several cell edge curves were extracted from a single snapshot by rotating the RBC core before projection on the x
plane, and a number of core snapshots corresponding to different times was averaged. A steep beginning of the CFL thickness distribution in (right) characterizes the global minimum of the RBC core, while the slowly decaying tail corresponds to a number of cavities in the RBC core formed by two or more adjacent RBCs as shown in (left). To characterize variations in CFL thickness the SD value σδ
was calculated. Also, the persistence of spatial variation was obtained by calculation of the correlation length of the cell edge pattern given by
An example of a CFL edge (left) and CFL thickness distribution (right) for Ht = 0.45 and D = 20 μm. x is the coordinate along the cylinder axis and δ is the CFL thickness.
) is the covariance of (δ
) and (δ
is the average CFL, and d
is the separation length. The covariance decreases to zero with increasing d
, and if C
) ≈ 0, then d′
is the correlation length of the cell edge pattern. presents the simulated correlation lengths for different Ht
values and tube diameters. The correlation length of the cell edge pattern increases for larger Ht
as the cell core is denser. Note that in case of Ht
= 0.15, the correlation length is comparable with the cell diameter indicating small correlations among separate cells. The length d′
appears to be independent of the tube diameter, in agreement with in vivo
] for the vessel diameter range of 20–50 μ
m. Furthermore, d′
values found in [20
] for Ht
= 0.42 are in the range of 25–35 μ
m, which is consistent with the simulated values of 27–28 μ
m for Ht
= 0.45 in .
Correlation lengths of the cell edge pattern for different tube diameters and Ht values. Ht is the tube hematocrit and D is the tube diameter.
shows CFLs for different tube diameters and cell volume fractions obtained using the “repulsion” and “reflection” EV interactions (left) and the “wall force” method in comparison with in vitro
] and in vivo
] (right). The simulated CFLs are consistent with the described Fahraeus and Fahraeus– Lindvuist effects. CFLs are wider for lower Ht
values and larger tube diameters indicating migration of RBCs to the tube centerline. CFL thicknesses for different EV interactions in (left) seem to have insignificant discrepancies. This appears to be true for larger tube diameters and lower cell volume fractions in agreement with the analogous independence of the relative apparent viscosity of an RBC suspension on EV interactions found in the Fahraeus–Lindqvist effect section. However, in case of small tube diameters and high Ht
values, even a small change in the CFL thickness may result in a considerable increase of the relative apparent viscosity. As an example, for the case of D
= 10 μ
m and Ht
= 0.45, the CFLs are 1.19, 1.3, and 1.45 μ
m for the “repulsion,” “reflection” EV interactions, and the “wall force” method, respectively, and the corresponding ηrel
values are 2.4, 2.002, and 1.6 ( and ) showing a significant decrease.
Figure 10 CFLs obtained in blood flow simulations employing the “repulsion” and “reflection” EV interactions (left) and the “wall force” setup in comparison with experimental data [20,22,30] (right) for various H (more ...)
A comparison of the simulated CFLs and those obtained in experiments in (right) shows partial agreement. Note that our simulations mimic blood flow in rigid and long tubes. Thus, the CFL measurements are carried out after the stationary state is reached and cell migration can be neglected. Bugliarello and Sevilla [2
] performed experiments of blood flow in glass tubes at various Hd
values shown by the “x” symbols in (right). The CFLs obtained in these experiments appear to be close to the simulated values if we calculate the corresponding Ht
= 0.14 and Ht
= 0.31) from equation (12)
. Reinke et al.
] conducted in vitro
experiments of blood flow at Hd
= 0.45 in glass tubes. The obtained CFL value for D
= 31 μ
= 0.45 suitable for comparison (diamond symbol in (right)) shows good agreement with simulations. Note that the simulated CFLs lie within the range of a scatter in the available in vitro
experimental data on CFLs. We are not aware of any other in vitro
experimental measurements of CFLs suitable for comparison.
CFLs from in vivo
] plotted in (right) show satisfactory agreement with simulations for high Ht
= 0.42–0.45), whereas for low Ht
, the correspondence is rather poor. In fact, available in vivo
measurements of CFLs show a scatter in the results. Yamaguchi et al.
] found the CFL in cat cerebral microvessels to be approximately 4 μ
m (not shown in [right]) and independent of vessel diameter. Maeda et al.
] reported CFLs in the range of 1–1.8 μ
m for vessel diameters 10–40 μ
m at Ht
= 0.45, where a rabbit mesentery was perfused. Kim et al.
] obtained CFLs in the range of 0.5– 3 μ
m increasing from the diameter 10–50 μ
m at Ht
= 0.42 where a rat cremaster muscle was used. The variability of in vivo
measurements of the CFL and their discrepancies with simulations shown in (right) may be due to several reasons such as the existence of the glycocalyx layer, variations in vessel width, use of a short vessel, close proximity of the site of CFL measurements to vessel bifurcations, vessel elasticity, and spatial resolution of the measurements.
The thickness of the glycocalyx layer is estimated to be approximately of 0.4–0.5 μ
]. The experimental data of Kim et al.
] shown in (right) include the glycocalyx thickness and yield satisfactory agreement with simulations. Subtracting its average thickness from the experimental CFL values would result in very thin CFLs especially for small vessel diameters (D
= 10 μ
m) and yield further discrepancies in comparison with simulations. The other referenced experimental sources do not explicitly state whether the thickness of the glycocalyx layer was included into CFL measurements. The largest discrepancies between simulated CFLs and those measured in in vivo
experiments are comparable with the glycocalyx thickness for Ht
= 0.42–0.45, while for lower Ht
values discrepancies appear to be more significant.
Variations in vessel diameter may locally alter blood flow and impact CFL measurements. Maeda et al.
] reported spatial variations in the vessel diameter on the order of several percent, while Kim et al.
] estimated a potential error in CFL measurements to be approximately 10% due to the change in the vessel diameter. Even though variations in the vessel diameter seem to have a rather insignificant effect on CFL measurements, it may influence migration of RBCs toward the center of the vessel, which is the main mechanism for CFL development. Moreover, RBC migration strongly depends on the vessel length and flow nonuniformity due to nearby vessel bifurcations. The perfused vessel length was not specified in experiments and it may not be long enough to allow RBCs to fully migrate to the vessel center and to form a steady core of the flow. In simulations, periodic boundary conditions in the direction of the tube axes were employed, which is essentially equivalent to the infinite length of the tube, and complete migration of RBCs can be considered to be fully achieved. Incomplete migration of RBCs in short vessels would result in the smaller CFLs found in experiments [22
] in comparison with those in simulations ( [right]).
Upstream flow conditions such as vessel bifurcations affect the width of the CFL [8
]. Kim et al.
] reported a nonaxisymmetric nature of CFL variations on opposite sides of the arteriole with a difference of 0.5 μ
m in the mean width of CFLs on the two sides. This difference in CFLs was attributed to a possible effect of upstream bifurcations. Note that the simulated CFLs here are found to be axisymmetric. Pries et al.
] reported that an axisymmetric cell profile is recovered after at least ten vessel diameters downstream of a bifurcation. However, this distance may not be long enough to account for full RBC migration. At high Ht
values, CFLs may not be very sensitive to the discussed geometrical uncertainties present in the in vivo
experiments, whereas at low Ht
they may significantly affect measured CFLs.
The magnitude of the CFLs reported by Maeda et al.
] plotted in (right) are smaller than those obtained in the simulations. In accordance, the reported values of the relative apparent viscosity for the hardened vessels of diameters in the range of 25–35 μ
m were 1.63, 2.0, and 2.83 for Ht
= 0.16, 0.3, and 0.45, respectively. These viscosities are considerably higher than those found in in vitro
] (1.29, 1.65, and 2.19 for D
= 30 μ
m and Ht
= 0.15, 0.3, and 0.45), which provides an additional explanation for the discrepancies found. Finally, spatial resolution of the experimental measurements may contribute to the discrepancies between the simulated CFLs and those obtained in in vivo
experiments. Kim et al.
] reported that their spatial resolution was 0.4 μ
m, which is on the order of the thickness of the glycocalyx layer, whereas the spatial resolution was not specified in Ref. [22
(left) presents SD values σδ
of the CFL thickness. Spatial variations of the CFL thickness characterized by σδ
appear to increase with tube diameter in agreement with experiments [20
]. The observed discrepancies between simulated and experimental SD values may be due to the differences discussed in the previous paragraph. The obtained SD values of the CFL thickness appear to be independent of Ht
and a model selected for EV interactions among cells (not shown here). (right) shows CFLs at Ht
= 0.45 for different tube diameters and shear rates
. The CFL thickness is slightly decreasing with increasing shear rate or equivalently flow rate. This is consistent with experiments [20
], where the same trend was found for the range of shear rates 100–500 per second. A decrease in the CFL thickness for higher flow rates is attributed to intensified hydrodynamic interactions among cells at higher shear rates that tend to expand the RBC core of the flow. Note that a decrease in the CFL thickness in this case may not necessarily result in an increase of the relative apparent viscosity because it corresponds to a characteristic viscosity of an RBC suspension averaged over the tube cross-section. A suspension of RBCs shows a shearthinning behavior with shear-dependent viscosity.
Figure 11 Spatial variations of the CFL thickness (SD) (left) and CFLs for different shear rates at Ht = 0.45 (right) for various Ht values and tube diameters. In vivo experimental data  for Ht = 0.42 are included in the left plot for comparison. Ht is the (more ...)