Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC3529161

Formats

Article sections

Authors

Related links

Microcirculation. Author manuscript; available in PMC 2012 December 23.

Published in final edited form as:

PMCID: PMC3529161

NIHMSID: NIHMS419430

DMITRY A. FEDOSOV,^{*}^{†} BRUCE CASWELL,^{‡} ALEKSANDER S. POPEL,^{§} and GEORGE EM KARNIADAKIS^{*}

Address for correspondence: George Em Karniadakis, Division of Applied Mathematics, Brown University, Providence, RI 02912, USA. Email: ude.nworb.mad@kg

The publisher's final edited version of this article is available at Microcirculation

See other articles in PMC that cite the published article.

Blood is modeled as a suspension of red blood cells using the dissipative particle dynamics method. The red blood cell membrane is coarse-grained for efficient simulations of multiple cells, yet accurately describes its viscoelastic properties. Blood flow in microtubes ranging from 10 to 40 *μ*m in diameter is simulated in three dimensions for values of hematocrit in the range of 0.15–0.45 and carefully compared with available experimental data. Velocity profiles for different hematocrit values show an increase in bluntness with an increase in hematocrit. Red blood cell center-of-mass distributions demonstrate cell migration away from the wall to the tube center. This results in the formation of a cell-free layer next to the tube wall corresponding to the experimentally observed Fahraeus and Fahraeus–Lindqvist effects. The predicted cell-free layer widths are in agreement with those found in *in vitro* experiments; the results are also in qualitative agreement with *in vivo* experiments. However, additional features have to be taken into account for simulating microvascular flow, e.g., the endothelial glycocalyx. The developed model is able to capture blood flow properties and provides a computational framework at the mesoscopic level for obtaining realistic predictions of blood flow in microcirculation under normal and pathological conditions.

Blood is a physiological fluid that consists of erythrocytes or RBCs, leukocytes or WBCs, platelets, and plasma containing various molecules and ions. Its main functions are the transport of oxygen and nutrients to cells of the body, removal of waste products such as carbon dioxide and urea, and circulation of molecules and cells that mediate the organism’s defense and immune response and play a fundamental role in the tissue repair process. RBCs constitute approximately 45% of the total blood volume, WBCs around 0.7%, and the rest is taken up by blood plasma and its substances. Under healthy conditions human RBCs have a biconcave shape of approximately 8 *μ*m in diameter and are highly deformable, which allows them to pass through narrow capillaries with a diameter several times smaller than the RBC size. Due to a high volume fraction of RBCs, the rheological properties of blood are mainly determined by the RBC properties.

*In vitro* experiments [11,28,30] of blood flow in glass tubes with diameters ranging from 3–4 to 1,000 *μ*m have shown a dependence of the apparent blood viscosity on the tube diameter, RBC volume fraction or hematocrit, cell aggregability, and flow rate. In tubes with diameters larger than 500–600 *μ*m blood can be assumed to be a nearly Newtonian fluid with a constant effective viscosity, whereas in smaller tubes it shows complex rheological behavior. The apparent blood viscosity reaches its minimum at the tube diameter of about 7 *μ*m and shows a significant increase with respect to this minimum as the tube diameter decreases or increases.

An increase of the apparent viscosity with increasing tube diameter is known as the Fahraeus–Lindqvist effect [11]. This effect is attributed to a cross-stream migration of RBCs in tube flow leading to the formation of two phases [3,16]: a flow core consisting mainly of RBCs and a CFL or a cell-depleted layer next to the tube wall devoid or depleted of cells. The CFL has a lower viscosity in comparison with the RBC core and serves as a lubrication layer reducing the effective blood viscosity [3,26]. The cross-stream migration of RBCs in tube flow is governed by cell-wall hydrodynamic interactions, which drive the cells away from the wall, and by cell–cell hydrodynamic interactions, which tend to disperse RBCs [15]. The CFL formed by cell migration depends on the tube diameter, the hematocrit value, and the flow rate, and is correlated with the apparent blood viscosity such that increasing the CFL thickness reduces the apparent viscosity of blood.

*In vivo* experiments [20,22,27] of blood perfused through microvessels have shown that the Fahraeus–Lindqvist effect is mediated by a CFL formed next to the vessel walls. However, blood flow resistance in microvessels *in vivo* is markedly higher than that in microtubes [29,36]. Microvessels in comparison with glass tubes are elastic, lined with endothelium, relatively short, and may be irregular in shape. These differences are likely to attenuate RBC migration resulting in a decreased CFL and an increased effective viscosity. In addition, special features of the vessel walls (e.g., geometry, glycocalyx layer) may contribute to spatial variations of the CFL that tend to diminish the effect of the CFL on the apparent blood viscosity [33].

Computational modeling of blood flow in microtubes or microvessels is challenging as a homogeneous continuum description is not adequate due to the motion of individual deformable cells. However, there exist several inhomogeneous continuum approaches to model blood flow in microvessels including the inhomogeneous continuum model of [24] and the two-phase fluid model [33], where two fluids with different viscosities represent the RBC flow core and the CFL, respectively. Several computational approaches to model blood flow as a suspension of cells have recently been developed. These include two-dimensional models [1,34] and three-dimensional simulations [7,21,23]. Even though two-dimensional models may qualitatively capture the behavior of an RBC suspension, their quantitative predictions have to be considered cautiously, as RBC motion and deformation in a flow are inherently three-dimensional. Three-dimensional modeling of an RBC suspension is limited due to computational expense such that the number of cells modeled in Refs. [21,23] was on the order of *O*(10). Dupin *et al.* [7] were able to simulate blood flow in microchannels with a characteristic size up to 30 *μ*m employing *O*(10^{2}) RBCs. We have developed an efficient parallel code that employs the multiscale RBC model described in Refs. [13,14]. This code allows us to simulate the total number of RBCs on the order of *O*(10^{4}) at reasonable computational cost.

The study is organized as follows. The next section presents details of the multiscale RBC model. The following section contains results of simulations of blood flow in tubes with diameters ranging from 10 to 40 *μ*m for different hematocrit values. Flow velocity profiles, the Fahraeus– Lindqvist effect, flow resistance, and the corresponding CFLs will be examined and compared with available experiments. Finally, we conclude with a brief discussion.

RBCs and their enclosing and surrounding fluids are modeled within the DPD method framework [19]. DPD is a mesoscopic particle-based simulation technique, where each particle represents a *cluster* of atoms or molecules rather than an individual atom. A detailed description of the DPD method can be found in Appendix S1.

The RBC membrane is represented by *N _{v}* DPD particles with coordinates {

$${V}_{\text{s}}=\sum _{j\in 1,\dots ,{N}_{\text{s}}}\left[\frac{{k}_{\text{B}}{Tl}_{\text{m}}(3{x}_{j}^{2}-2{x}_{j}^{3})}{4p(1-{x}_{j})}+\frac{{k}_{p}}{(n-1){l}_{j}^{n-1}}\right],$$

(1)

where *l _{j}* is the length of the spring

To incorporate the membrane viscosity into the RBC model, a dissipative force is introduced for each spring. Following the general framework of the fluid particle model [9], we can define dissipative ${\mathbf{F}}_{ij}^{\text{D}}$ and random ${\mathbf{F}}_{ij}^{\text{R}}$ forces for each spring, which satisfy the fluctuation–dissipation balance providing consistent temperature of the RBC membrane in equilibrium. The forces are given by

$${\mathbf{F}}_{ij}^{\text{D}}=-{\gamma}^{\text{T}}{\mathbf{v}}_{ij}-{\gamma}^{\text{C}}({\mathbf{v}}_{ij}\xb7{\mathbf{e}}_{ij}){\mathbf{e}}_{ij},$$

(2)

$${\mathbf{F}}_{ij}^{\text{R}}\text{d}t=\sqrt{2{k}_{\text{B}}T}\left(\sqrt{2{\gamma}^{\text{T}}}\text{d}\overline{{\mathbf{W}}_{ij}^{\text{S}}}+\sqrt{3{\gamma}^{\text{C}}-{\gamma}^{\text{T}}}\frac{\text{tr}[\text{d}{\mathbf{W}}_{ij}]}{3}\mathbf{1}\right)\xb7{\mathbf{e}}_{ij},$$

(3)

where *γ*^{T} and *γ*^{C} are dissipative parameters, *v** _{ij}* is the relative velocity of spring ends,

The bending resistance of the RBC membrane is modeled by

$${V}_{\text{b}}=\sum _{j\in 1,\dots ,{N}_{s}}{k}_{\text{b}}[1-cos({\theta}_{j}-{\theta}_{0})],$$

(4)

where *k*_{b} is the bending constant, θ* _{j}* is the instantaneous angle between two adjacent triangles having the common edge

In addition, the RBC model includes the area and volume conservation constraints, which mimic area-incompressibility of the lipid bilayer and incompressibility of the interior fluid, respectively. The corresponding energy is given by

$${V}_{\text{a}+\text{v}}=\sum _{j\in 1,\dots ,{N}_{t}}\frac{{k}_{\text{d}}{({A}_{j}-{A}_{0})}^{2}}{2{A}_{0}}+\frac{{k}_{\text{a}}{(A-{A}_{0}^{\text{tot}})}^{2}}{2{A}_{0}^{\text{tot}}}+\frac{{k}_{\text{v}}{(V-{V}_{0}^{\text{tot}})}^{2}}{2{V}_{0}^{\mathit{tot}}},$$

(5)

where *N _{t}* is the number of triangles in the membrane network,

The relation of the parameters of triangulated surface models of membranes to the macroscopic parameters (shear, area-compression, and Young’s moduli) of the elasticity theory has been previously discussed in Refs. [5,17,32]. Extension of the linear analysis of [5] for a regular hexagonal network allows us to uniquely relate the model parameters and the network macroscopic elastic properties, see References [12,13] for details. The derived shear modulus of the membrane is given by

$${\mu}_{0}=\frac{\sqrt{3}{k}_{\text{B}}T}{4p{l}_{\text{m}}{x}_{0}}\left(\frac{{x}_{0}}{2{(1-{x}_{0})}^{3}}-\frac{1}{4{(1-{x}_{0})}^{2}}+\frac{1}{4}\right)+\frac{\sqrt{3}{k}_{p}(n+1)}{4{l}_{0}^{n+1}},$$

(6)

where *l*_{0} is the equilibrium spring length and *x*_{0} = *l*_{0}/*l*_{m}. The area-compression *K* and Young’s *Y* moduli are equal to 2*μ*_{0} + *k*_{a} + *k*_{d} and 4*Kμ*_{0}/(*K* + *μ*_{0}), respectively.

The relation between the model bending coefficient *k*_{b} and the macroscopic bending rigidity *k*_{c} of the Helfrich model [18] can be derived as
${k}_{\text{b}}=2{k}_{\text{c}}/\sqrt{3}$ for a spherical membrane [12,13]. This expression describes bending contribution of the energy in Equation (4), but may not fully represent actual bending resistance of the RBC membrane as membrane bending may also result in local in-plane deformations.

The membrane shear viscosity *η*_{m} is related to the dissipative parameters *γ*^{T}, *γ*^{C} as
${\eta}_{\text{m}}=\sqrt{3}({\gamma}^{\text{T}}+{\gamma}^{\text{C}}/4)$. Here *γ*^{T} accounts for a large portion of viscous contribution, and hence *γ*^{C} is set to *γ*^{T}/3 in all simulations.

In practice, the given macroscopic RBC properties serve as an input to be used to calculate the necessary mesoscopic model parameters from the equations above without any manual adjustment. A simulation of an RBC in equilibrium shows that the membrane may develop local bumps due to stress anomalies in a membrane triangulation as a network on a closed surface cannot consist of triangles whose edges have the same lengths. Such local stress artifacts depend on the network regularity and the ratio of the membrane elastic and bending contributions given by the Föppl–von Kármán number
$\kappa ={YR}_{0}^{2}/{k}_{\text{c}}$, where
${R}_{0}=\sqrt{{A}_{0}^{\text{tot}}/(4\pi )}$. To eliminate the stress artifacts, we employed a “stress-free” model obtained by computational annealing. Thus, the equilibrium length
${l}_{0}^{i}$ of each spring is set to the edge length after triangulation for *i* = 1, …, *N*_{s} and the individual maximum spring extension is set to
${l}_{\text{m}}^{i}={l}_{0}^{i}\times {x}_{0}$ (*x*_{0} is a constant). This modification provides a network free of local stress anomalies. The necessary spring parameters (*p* and *k _{p}*) in Equation (1) are calculated individually for each spring for given

Both internal and external fluids are simulated by a collection of free DPD particles and are separated by the RBC membrane through bounce-back reflections at the membrane surface. Moreover, a dissipative force between fluid particles and membrane vertices is set properly to account for the no-slip boundary conditions at the membrane surface. More details on boundary conditions can be found in Appendix S1.

Scaling between DPD model units (M) and physical units (P) adopts the following length and time scales

$${r}^{\text{M}}=\frac{{D}_{0}^{\text{P}}}{{D}_{0}^{\text{M}}}(\text{meters}),\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\tau =\frac{{D}_{0}^{\text{P}}}{{D}_{0}^{\text{M}}}\frac{{\eta}^{\text{P}}}{{\eta}^{\text{M}}}\frac{{Y}^{\text{M}}}{{Y}^{\text{P}}}(\text{seconds}),$$

(7)

where *r*^{M} is the model unit of length, *D*_{0} is the cell diameter, and *η* is the characteristic viscosity (e.g., internal/external fluids or membrane). Additionally, we can define scaling of the energy per unit mass (*k*_{B}*T*) as follows

$${({k}_{\text{B}}T)}^{\text{M}}=\frac{{Y}^{\text{P}}}{{Y}^{\text{M}}}{\left(\frac{{D}_{0}^{\text{P}}}{{D}_{0}^{\text{M}}}\right)}^{2}{({k}_{\text{B}}T)}^{\text{P}}.$$

(8)

This section presents results on the modeling of blood flow in microtubes. Blood was modeled as a suspension of RBCs using the model described above. Blood flow was simulated in tubes of diameters ranging from 10 to 40 *μ*m. The properties of the modeled cell suspensions flowing through tubes and corresponding CFLs were evaluated for different flow rates and tube hematocrit values.

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:
${D}_{0}^{\text{M}}=8.06,\phantom{\rule{0.16667em}{0ex}}{\mu}_{0}^{\text{M}}=100$*, 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 **...**

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

$${H}_{t}=\frac{{N}_{\text{c}}{V}^{\prime}}{{V}_{t}},$$

(9)

where *N*_{c} is the number of RBCs in the volume *V _{t}* =

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 Table 1, 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.

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. Figure 1 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 Figure 2 for tubes of diameters *D* = 10 *μ*m (left) and *D* = 40 *μ*m (right) for different *H _{t}* 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. *H*_{t} = 0.45. x is the coordinate along the cylinder axis.

Typical velocity profiles of blood flow in tubes of diameters *D* = 10 *μ*m (left) and *D* = 40 *μ*m (right) for different *H*_{t} values employing “repulsion” EV interactions. Dashed lines show the corresponding parabolic profiles **...**

Velocity profiles and subsequent number density distributions are shown over the half tube because of flow axisymmetry. Dashed lines in Figure 2 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

$$Q={\int}_{A}v(r)\text{d}A,$$

(10)

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. Table 2 presents blood flow characteristics for different tube diametes and

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

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.

Figure 3 shows center-of-mass distributions of RBCs in tubes of *D* = 10 *μ*m (left) and *D* = 40 *μ*m (right) for different cell volume fractions. The center-of-mass distributions for lower *H _{t}* values appear to be more confined around the tube centerline. This results in a larger CFL that can be roughly estimated as

Center-of-mass distributions of RBCs in tubes of *D* = 10 *μ*m (left) and *D* = 40 *μ*m (right) for different *H*_{t} values with “repulsion” EV interactions. Dotted lines denote the corresponding CFL thicknesses. *H*_{t} is the tube hematocrit **...**

Figure 4 presents the number densities of the solvent normalized by their average values in tubes of diameters 10 *μ*m (left) and 40 *μ*m (right) for different *H _{t}* values. The average number density of solvent is given by
${n}_{\text{s}}^{\text{avg}}={N}_{\text{s}}/{V}_{t}$, where

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

$${H}_{\text{d}}=\frac{{\overline{\nu}}_{\text{c}}}{\overline{\nu}}{H}_{t},$$

(11)

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

$$\frac{{H}_{t}}{{H}_{\text{d}}}={H}_{\text{d}}+(1-{H}_{\text{d}})(1+1.7{\text{e}}^{-0.35D}-0.6{\text{e}}^{-0.01D}),$$

(12)

where the tube diameter *D* is in micormeters. Figure 5 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 Table 1 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

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. *H*_{t} is the tube hematocrit. **...**

The remaining discrepancies for small tube diameters are due to solvent density fluctuations next to the wall as shown in Figure 4. Figure 6 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.

The Fahraeus–Lindqvist effect [11] corresponds to a decrease in the apparent blood viscosity with decreasing tube diameter found in experiments of blood flow in glass tubes [28]. The apparent viscosity is defined as follows

$${\eta}_{\text{app}}-\frac{\pi \mathrm{\Delta}{PD}^{4}}{128QL}=\frac{\mathrm{\Delta}{PD}^{2}}{32\overline{\nu}L}.$$

(13)

The apparent viscosity increases for higher *H _{t}* values since higher cell crowding yields larger flow resistance. It is more convenient to consider the relative apparent viscosity defined as

$${\eta}_{\text{rel}}=\frac{{\eta}_{\text{app}}}{{\eta}_{\text{s}}},$$

(14)

where *η _{s}* is the solvent viscosity. Figure 7 shows the simulated

The CFL is a near-wall layer of blood plasma absent of RBCs (see Figure 1) 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 Figure 1), which is similar to CFL measurements in experiments [20,22]. Figure 9 shows a sample CFL edge (left) from simulations and CFL thickness distribution (right) for *H _{t}* = 0.45 and

An example of a CFL edge (left) and CFL thickness distribution (right) for *H*_{t} = 0.45 and *D* = 20 *μ*m. x is the coordinate along the cylinder axis and *δ* is the CFL thickness.

$$C(d)=\langle (\delta (x)-\overline{\delta})(\delta (x+d)-\overline{\delta})\rangle ,$$

(15)

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. Table 3 presents the simulated correlation lengths for different *H _{t}* values and tube diameters. The correlation length of the cell edge pattern increases for larger

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.

Figure 10 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 Figure 10 (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

A comparison of the simulated CFLs and those obtained in experiments in Figure 10 (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 Figure 10 (right). The CFLs obtained in these experiments appear to be close to the simulated values if we calculate the corresponding *H _{t}* values (

CFLs from *in vivo* experiments [20,22] plotted in Figure 10 (right) show satisfactory agreement with simulations for high *H _{t}* values (

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 Figure 10 (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

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 (Figure 10 [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

The magnitude of the CFLs reported by Maeda *et al.* [22] plotted in Figure 10 (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

Figure 11 (left) presents SD values *σ _{δ}* of the CFL thickness. Spatial variations of the CFL thickness characterized by

Blood is modeled as a suspension of RBCs in a Newtonian fluid. Blood flow in microtubes under healthy conditions is simulated for different *H _{t}* values and tube diameters. Velocity profiles for different

The Fahraeus effect is characterized by a decrease in the discharge hematocrit when the tube diameter is increased, whereas the Fahraeus–Lindqvist effect describes a decrease in the blood apparent viscosity with decreasing tube diameter. The developed blood flow model is able to closely capture changes in the discharge hematocrit and the relative apparent viscosity for *H _{t}* values in the range 0.15–0.45 and tube diameters between 10 and 40

Fahraeus and Fahraeus–Lindqvist effects are directly related to the CFL next to the wall by providing an effective lubrication for the viscous RBC core to flow at much less friction. The simulated CFLs are in agreement with those found in *in vitro* and *in vivo* experiments; however, experimental CFL values show a considerable scatter that does not allow a very rigorous quantitative comparison. Discrepancies between the simulated CFLs and *in vitro* data are likely to be related to a potential RBC sedimentation, RBC aggregability, and differences in flow rates employed, whereas blood flow conditions *in vivo* are far more complex. In particular, CFLs *in vivo* are additionally affected by the presence of the glycocalyx layer, variations in vessel diameter, upstream flow conditions such as vessel bifurcations, and the length of a straight vessel segment where measurements are taken.

In conclusion, the developed model is able to capture blood flow properties and hence it may aid to obtain realistic predictions of blood flow in microcirculation.

This work was supported by NIH Grant Number R01Hl094270 and simulations were performed on the Cray XT5 at NSF/NICS.

- DPD
- Dissipative particle dynamics
- RBC
- Red blood cell
- WBC
- White blood cell
- CFL
- Cell-free layer
- EV
- Excluded volume
- SD
- Standard deviation

Additional Supporting Information may be found in the online version of this article (http://www.microcirculationjournal.com):

Appendix S1. Blood flow and cell-free layer in microvessels.

Please note: Wiley-Blackwell are not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

1. Bagchi P. Mesoscale simulation of blood flow in small vessels. Biophys J. 2007;92:1858–1877. [PubMed]

2. Bugliarello G, Sevilla J. Velocity distribution and other characteristics of steady and pulsatile blood flow in fine glass tubes. Biorheology. 1970;7:85–107. [PubMed]

3. Cokelet GR, Goldsmith HL. Decreased hydrodynamic resistance in the two-phase flow of blood through small vertical tubes at low flow rates. Circ Res. 1991;68:1–17. [PubMed]

4. Damiano ER. The effect of the endothelial-cell glycocalyx on the motion of red blood cells through capillaries. Microvasc Res. 1998;55:77–91. [PubMed]

5. Dao M, Li J, Suresh S. Molecularly based analysis of deformation of spectrin network and human erythrocyte. Mater Sci Eng C. 2006;26:1232–1244.

6. Discher DE, Boal DH, Boey SK. Simulations of the erythrocyte cyto-skeleton at large deformation. II. Micropipette aspiration. Biophys J. 1998;75(3):1584–1597. [PubMed]

7. Dupin MM, Halliday I, Care CM, Munn LL. Lattice Boltzmann modeling of blood cell dynamics. Int J Comput Fluid Dyn. 2008;22(7):481–492.

8. Enden G, Popel AS. A numerical study of plasma skimming in small vascular bifurcations. J Biomech Eng. 1994;116:79–88. [PubMed]

9. Espanol P. Fluid particle model. Phys Rev E. 1998;57(3):2930–2948.

10. Fahraeus R. The suspension stability of the blood. Physiol Rev. 1929;9:241–274.

11. Fahraeus R, Lindqvist T. Viscosity of blood in narrow capillary tubes. Am J Phys. 1931;96:562–568.

12. Fedosov DA. PhD thesis. Brown University; Providence, RI, USA: 2010. Multiscale modeling of blood flow and soft matter.

13. Fedosov DA, Caswell B, Karniadakis GE. A multiscale red blood cell model with accurate mechanics, rheology, and dynamics. Biophys J. 2010;98(10):2215–2225. [PubMed]

14. Fedosov DA, Caswell B, Karniadakis GE. Systematic coarse-graining of spectrin-level red blood cell models. Comp Meth Appl Mech Eng. 2010;199:1937–1948. [PMC free article] [PubMed]

15. Goldsmith HL. Red cell motions and wall interactions in tube flow. Fed Proc. 1971;30:1578–1590. [PubMed]

16. Goldsmith HL, Cokelet GR, Gaehtgens PR. Robin Fahraeus: evolution of his concepts in cardiovascular physiology. Am J Physiol. 1989;257:H1005–H1015. [PubMed]

17. Gompper G, Kroll DM. Random surface discretizations and the renormalization of the bending rigidity. J Phys I France. 1996;6(10):1305–1320.

18. Helfrich W. Elastic properties of lipid bilayers: theory and possible experiments. Z Naturforsch C. 1973;28:693–703. [PubMed]

19. Hoogerbrugge PJ, Koelman JMVA. Simulating microscopic hydrodynamic phenomena with dissipative particle dynamics. Europhys Lett. 1992;19(3):155–160.

20. Kim S, Long LR, Popel AS, Intaglietta M, Johnson PC. Temporal and spatial variations of cell-free layer width in arterioles. Am J Physiol. 2007;293:H1526–H1535. [PubMed]

21. Liu Y, Liu WK. Rheology of red blood cell aggregation by computer simulation. J Comput Phys. 2006;220:139–154.

22. Maeda N, Suzuki Y, Tanaka J, Tateishi N. Erythrocyte flow and elasticity of microvessels evaluated by marginal cell-free layer and flow resistance. Am J Physiol. 1996;271(6):H2454–H2461. [PubMed]

23. McWhirter JL, Noguchi H, Gompper G. Flow-induced clustering and alignment of vesicles and red blood cells in microcapillaries. Proc Natl Acad Sci USA. 2009;106(15):6039–6043. [PubMed]

24. Moyers-Gonzalez MA, Owens RG. Mathematical modelling of the cell-depleted peripheral layer in the steady flow of blood in a tube. Biorheology. 2010;47:39–71. [PubMed]

25. Pan W, Fedosov DA, Caswell B, Karniadakis GE. Hydrodynamic interactions for single dissipative-particle-dynamics particles and their clusters and filaments. Phys Rev E. 2008;78(4):046706. [PubMed]

26. Popel AS, Johnson PC. Microcirculation and hemorheology. Annu Rev Fluid Mech. 2005;37:43–69. [PMC free article] [PubMed]

27. Pries AR, Ley K, Claassen M, Gaehtgens P. Red cell distribution at microvascular bifurcations. Microvasc Res. 1989;38:81–101. [PubMed]

28. Pries AR, Neuhaus D, Gaehtgens P. Blood viscosity in tube flow: dependence on diameter and hematocrit. Am J Physiol. 1992;263(6):H1770–H1778. [PubMed]

29. Pries AR, Secomb TW, Gessner T, Sperandio MB, Gross JF, Gaehtgens P. Resistance to blood flow in microvessels *in vivo*. Circ Res. 1994;75:904–915. [PubMed]

30. Reinke W, Gaehtgens P, Johnson PC. Blood viscosity in small tubes: effect of shear rate, aggregation, and sedimentation. Am J Physiol. 1987;253:H540–H547. [PubMed]

31. Reinke W, Johnson PC, Gaehtgens P. Effect of shear rate variation on apparent viscosity of human blood in tubes of 29 to 90 microns diameter. Circ Res. 1986;59:124–132. [PubMed]

32. Seung HS, Nelson DR. Defects in flexible membranes with crystalline order. Phys Rev A. 1988;38(2):1005–1018. [PubMed]

33. Sharan M, Popel AS. A two-phase model for flow of blood in narrow tubes with increased effective viscosity near the wall. Biorheology. 2001;38:415–428. [PubMed]

34. Sun C, Munn LL. Particulate nature of blood determines macroscopic rheology: a 2D lattice-Boltzmann analysis. Biophys J. 2005;88:1635–1645. [PubMed]

35. Vink H, Duling BR. Identification of distinct luminal domains for macromolecules, erythrocytes, and leukocytes within mammalian capillaries. Circ Res. 1996;79:581–589. [PubMed]

36. Weinbaum S, Tarbell JM, Damiano ER. The structure and function of the endothelial glycocalyx layer. Annu Rev Biomed Eng. 2007;9:121–167. [PubMed]

37. Yamaguchi S, Yamakawa T, Niimi H. Cell-free plasma layer in cerebral microvessels. Biorheology. 1992;29:251–260. [PubMed]

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |