Modeling Parameters

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

*D*_{0} = 7.82

*μ*m, the Young’s modulus

*Y*_{0} = 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,

14]. We employed the RBC model with the following parameters:

*, x*_{0} = 2.2,

*n* = 2,

*k*_{a} = 4,900,

*k*_{d} = 100, and

*k*_{v} = 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

*k*_{B}*T* = 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 (

*V*_{0} = 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

where

*N*_{c} is the number of RBCs in the volume

*V*_{t} =

*πR*^{2}*L*,

*R* 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

*H*_{t} = 0.45.

Two types of EV interactions among cells were considered. The first case is shown in , where the repulsive force coefficient between membrane vertices *a*_{V–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 *r*_{c} = 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 *a*_{V–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 *a*_{W–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 *n*_{S} = *n*_{W} = 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 *H*_{t} values.

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) with

*η* = 0.0012 Pa·s being the plasma viscosity. The velocity curves are averaged over the tube cross-section and over 2 × 10

^{5} time steps, which corresponds to the total time of 0.1 seconds. The pressure gradients employed were 2.633 × 10

^{5} and 6.582 × 10

^{4} Pa/m for tubes of diameters 10 and 40

*μ*m, respectively. In the case of low

*H*_{t} (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

where

*A* 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

*H*_{t} values indicating a wider RBC core and a smaller CFL. presents blood flow characteristics for different tube diametes and

*H*_{t} values. The table includes mean flow velocities defined as

=

*Q/A*, mean shear rates

, pressure gradients Δ

*P*/

*L*, centerline flow velocities

*v*_{c}, and the corresponding centerline velocities

*v*_{p} of the Newtonian plasma in the absence of the cells. The mean flow velocities increase for tubes having larger diameters with comparable shear rates.

| **Table 2**Blood flow characteristics for different tube diameters and *H*_{t} values employing “repulsion” EV interactions. H_{t} is the tube hematocrit, D is the tube diameter, is the mean flow velocity, is the mean (more ...) |

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.

Fahraeus Effect

The Fahraeus effect, discovered in

*in vitro* experiments of blood flow in glass tubes [

10], shows an increased value of discharge hematocrit (

*H*_{d}) measured at the tube exit in comparison with that of the cell suspension before the tube entrance.

*H*_{d} 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

*H*_{d} value measured at the tube discharge. We define

*H*_{d} in simulations as follows

where

=

*Q/A* and

_{c} 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.* [

28] have compiled a number of experiments on blood flow through glass tubes of different diameters at various

*H*_{d} values to obtain an empirical relation between

*H*_{t} and

*H*_{d} given by

where the tube diameter

*D* is in micormeters. compares

*H*_{d} 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

*H*_{t} values, whereas discrepancies are found for large

*H*_{t} and small tube diameters (

*D* = 10

*μ*m). EV interactions appear to be of importance at

*H*_{t} = 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

*H*_{d} value as seen in .

The remaining discrepancies for small tube diameters are due to solvent density fluctuations next to the wall as shown in . presents *H*_{d} 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 *H*_{d} values with those in *in vitro* experiments. Note that the repulsive wall does not affect the simulation results for low *H*_{t} values and large tube diameters as the CFL is wider than the distance of the effective cell– wall interactions.

In conclusion, the “wall force” method yields superior results for *H*_{d} 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.

Cell-Free Layer

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 *H*_{d} 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,

22]. shows a sample CFL edge (left) from simulations and CFL thickness distribution (right) for

*H*_{t} = 0.45 and

*D* = 20

*μ*m. The cell edge was measured by projecting the cell vertices of a snapshot of RBCs on the

*x*–

*y* 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*–

*y* 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

where

*C*(

*d*) is the covariance of (

*δ*(

*x*) −

) and (

*δ*(

*x*+

*d*) −

),

is the average CFL, and

*d* is the separation length. The covariance decreases to zero with increasing

*d*, and if

*C*(

*d′*) ≈ 0, then

*d′* is the correlation length of the cell edge pattern. presents the simulated correlation lengths for different

*H*_{t} values and tube diameters. The correlation length of the cell edge pattern increases for larger

*H*_{t} as the cell core is denser. Note that in case of

*H*_{t} = 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* experiments [

20] for the vessel diameter range of 20–50

*μ*m. Furthermore,

*d′* values found in [

20] for

*H*_{t} = 0.42 are in the range of 25–35

*μ*m, which is consistent with the simulated values of 27–28

*μ*m for

*H*_{t} = 0.45 in .

| **Table 3**Correlation lengths of the cell edge pattern for different tube diameters and *H*_{t} values. *H*_{t} 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* experiments [

30] and

*in vivo* experiments [

20,

22] (right). The simulated CFLs are consistent with the described Fahraeus and Fahraeus– Lindvuist effects. CFLs are wider for lower

*H*_{t} 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

*H*_{t} 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

*H*_{t} = 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.

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

*H*_{d} 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

*H*_{t} values (

*H*_{t} = 0.14 and

*H*_{t} = 0.31) from

equation (12). Reinke

*et al.* [

30] conducted

*in vitro* experiments of blood flow at

*H*_{d} = 0.45 in glass tubes. The obtained CFL value for

*D* = 31

*μ*m,

*H*_{d} = 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* experiments [

20,

22] plotted in (right) show satisfactory agreement with simulations for high

*H*_{t} values (

*H*_{t} = 0.42–0.45), whereas for low

*H*_{t}, the correspondence is rather poor. In fact, available

*in vivo* measurements of CFLs show a scatter in the results. Yamaguchi

*et al.* [

37] found the CFL in cat cerebral microvessels to be approximately 4

*μ*m (not shown in [right]) and independent of vessel diameter. Maeda

*et al.* [

22] reported CFLs in the range of 1–1.8

*μ*m for vessel diameters 10–40

*μ*m at

*H*_{t} = 0.45, where a rabbit mesentery was perfused. Kim

*et al.* [

20] obtained CFLs in the range of 0.5– 3

*μ*m increasing from the diameter 10–50

*μ*m at

*H*_{t} = 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

*μ*m [

4,

35]. The experimental data of Kim

*et al.* [

20] 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

*H*_{t} = 0.42–0.45, while for lower

*H*_{t} values discrepancies appear to be more significant.

Variations in vessel diameter may locally alter blood flow and impact CFL measurements. Maeda

*et al.* [

22] reported spatial variations in the vessel diameter on the order of several percent, while Kim

*et al.* [

20] 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,

27]. Kim

*et al.* [

20] 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.* [

27] 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

*H*_{t} values, CFLs may not be very sensitive to the discussed geometrical uncertainties present in the

*in vivo* experiments, whereas at low

*H*_{t} they may significantly affect measured CFLs.

The magnitude of the CFLs reported by Maeda

*et al.* [

22] 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

*H*_{t} = 0.16, 0.3, and 0.45, respectively. These viscosities are considerably higher than those found in

*in vitro* experiments [

28] (1.29, 1.65, and 2.19 for

*D* = 30

*μ*m and

*H*_{t} = 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.* [

20] 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

*H*_{t} and a model selected for EV interactions among cells (not shown here). (right) shows CFLs at

*H*_{t} = 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.