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

**|**HHS Author Manuscripts**|**PMC3085988

Formats

Article sections

- Abstract
- 1. Introduction
- 2. RBC aggregation
- 3. Bifurcation geometry
- 4. Boundary conditions and Chimera grids
- 5. Results
- 6. Discussion
- References

Authors

Related links

Comput Methods Biomech Biomed Engin. Author manuscript; available in PMC 2011 May 3.

Published in final edited form as:

PMCID: PMC3085988

NIHMSID: NIHMS283084

Bong Chung,^{a,}^{*} Sangho Kim,^{b,}^{1} Paul C. Johnson,^{b,}^{2} and Aleksander S. Popel^{a,}^{3}

The publisher's final edited version of this article is available at Comput Methods Biomech Biomed Engin

See other articles in PMC that cite the published article.

Aggregate formation of red blood cells (RBCs) in a postcapillary venular bifurcation is investigated with three-dimensional computer simulations using the Chimera grid method. Interaction energy between the RBCs is modelled by a depletion interaction theory; RBCs are modelled as rigid oblate ellipsoids. The cell–cell interactions of RBCs are strongly dependent on vessel geometry and shear rates. The experimental data on vessel geometry, pseudoshear rates, and Dextran concentration obtained in our previous *in vivo* RBC aggregation study in postcapillary venules of the rat spinotrapezius muscle were used to simulate RBC aggregation. The computational results were compared to the experimental results from the *in vivo* study. The results show that cells have a larger tendency to form an aggregate under reduced flows. Aggregate formation also depends on the angle and location of the cells before they enter the bifurcation region. Comparisons with experimental data are discussed.

Red blood cell (RBC) aggregation is known to affect blood flow in microcirculation (Popel and Johnson 2005). Numerous *in vivo* and *in vitro* experiments showed an inverse relationship between the relative viscosity and shear rate at low shear rates (Bishop et al. 2001a, 2004). Earlier *in vivo* studies using rat spinotrapezius muscle revealed that blood velocity profiles become blunt at a low shear rate due to RBC aggregation (Bishop et al. 2001b). At a low pseudoshear rate (defined as the ratio of mean velocity to vessel diameter), aggregation forces between cells exceed hydrodynamic forces that tend to separate them, resulting in formation of RBC aggregates in microvessels, especially in venules, thus increasing effective viscosity. The level of the bluntness of velocity profiles due to RBC aggregates determines non-Newtonian characteristics of blood in microvessels. There have been many studies of constitutive relationships between viscosity and shear rate that describe non-Newtonian behaviour of blood such as power law and Casson models (Das et al. 2007). However, detailed mechanisms of RBC aggregation at the cellular level have not been clearly understood.

In a recent study (Kim et al. 2005), we examined the aggregation process in postcapillary venules of the rat spinotrapezius muscle, which is the site in the venous network where the aggregates first appear as the RBCs emerge from converging capillaries. Dextran (DEX) 500 (0.6% in plasma) was used to induce RBC aggregation *in vivo*, since rat RBCs do not have the innate ability to aggregate or aggregate only weakly. In that study, we showed that most aggregates (>90%) formed in a localised region of the postcapillary venule, approximately 15–30 μm from the vessel entrance, i.e., from the apex of the bifurcation. No aggregates formed in the region less than 15 μm from the entrance and few new aggregates formed in the region 30–50 μm from the entrance. It has been proposed that aggregates did not form in the initial region 0–15 μm from the entrance because the velocity profile was not fully developed (the entrance effect) and that in the distal region more than 30 μm from the entrance all the red cells whose properties were favourable to aggregate formation, e.g. the older cells (Rampling et al. 2004) had already formed aggregates and the remainder were younger cells that had a much lower aggregability. These experiments demonstrated that the level of RBC aggregation depends on flow rate, vascular geometry and molecular weight and concentration of DEX.

Different numerical methods have been used in computational analyses of blood flow with discrete RBCs. Such methods include a Boundary Integral Method (Pozrikidis 2001), an Immersed Boundary Method (IB; Eggleton and Popel 1998; Bagchi et al. 2005; Bagchi 2007), an Immersed Finite Element Method (Liu et al. 2004) and a Lattice Boltzmann Method (LBM; Sun et al. 2003; Sun and Munn 2005; Zhang et al. 2007). The IB method in combination with a finite element method for the RBC membrane has been used to simulate a single three-dimensional RBC in a simple shear flow (Eggleton and Popel 1998). The IB method was also used to simulate aggregation of two two-dimensional RBCs in a shear flow (Bagchi et al. 2005). The LBM method was used to simulate flow of multiple two-dimensional RBCs, modelled as rods with rounded ends, in a straight channel to reproduce Fahraeus and Fahraeus-Lindqvist effects (Sun and Munn 2005).

Our earlier computational study on RBC aggregation showed that the Chimera grid scheme, developed by Dougherty (Dougherty 1985) is applicable to simulations of three-dimensional motion and aggregation of RBCs, which are modelled as rigid ellipsoidal particles in narrow tubes; the methods can be expanded to describe the flow of multiple deformable particles. We employed the depletion interaction theory (Neu and Meiselman 2002) to model the interaction energy between the cells for aggregation. In the present study, we apply the Chimera grid scheme to investigate the aggregation of RBCs in the region of a postcapillary bifurcation, based on the depletion interaction theory. The RBC is deformable; its membrane that envelopes the hemoglobin solution with five times larger viscosity than plasma’s is composed of lipid bilayer, its associated proteins and underlying cytoskeleton (spectrin, actin, other proteins), which exhibit flexibility. Nevertheless, our interest in the present study is not in investigating detailed mechanisms of aggregation process between cells but to examine shear-dependent aggregation based on RBC’s initial position and angle as will be described below. While our computational method can handle deformable cells through the arbitrary Lagrangian–Eurerian technique, the RBCs are modelled as rigid oblate ellipsoidal particles with eccentricity *e* = 0.3 to reduce computational costs. The value of eccentricity mimics the dimensions of the rat RBC. Based on the Chimera grid scheme, we aim at (i) simulating the motion of two RBCs with and without aggregation in a converging bifurcation and (ii) to qualitatively and quantitatively comparing the computational results with our earlier *in vivo* experimental data.

The simulation was performed for RBC motion with and without aggregation under normal (*Re* = 0.1) and reduced (*Re* = 0.01) flows. The various cases by initial locations and orientations of two cells were considered. Initial orientations and locations of the cells determine the velocities of the cells and as a result, a distance between the cells at the region of bifurcation which plays a key role in aggregation. Because aggregation between the cells depends on the balance between hydrodynamic forces and attractive forces between the cells, we assume that initial cell orientation and location could play important roles on aggregation of two cells at the bifurcation region. Therefore in the present study, we show the effects of initial location and orientations of two cells on aggregation at normal and reduced flows with and without DEX.

We model RBCs as rigid oblate ellipsoids with eccentricity ratio of 0.3, mimicing the thickness to width ratio of rat RBC, as illustrated in Figure 1. The density of RBC is close to the density of plasma. Accordingly, RBCs are considered to be purely buoyant in plasma flow. The depletion interaction theory applied to RBC aggregation (Neu and Meiselman 2002) proposes that a depletion layer between two adjacent cells restricts the space available to polymer molecules (e.g. DEX) between cells by displacing solvent into the bulk phase, and this generates osmotic pressure between the depletion zone and the bulk phase. According to the theory, total interaction energy per unit surface *W*_{T} between two infinite surfaces due to polymers such as DEX and Poly Ethylene Glycol is the sum of depletion energy *W*_{D} and electrostatic repulsive energy *W*_{E} (with negligible van der Waals interaction)

$${W}_{\mathrm{T}}={W}_{\mathrm{D}}+{W}_{\mathrm{E}},$$

(1)

where

$$\begin{array}{c}\hfill {W}_{\mathrm{D}}=-2\Pi \left(\Delta -\frac{d}{2}+\delta -P\right)\hfill \\ \hfill \text{when}\left(\frac{d}{2}-\delta +P\right)<\Delta \hfill \\ \hfill {W}_{\mathrm{D}}=0\phantom{\rule{1em}{0ex}}\text{when}\left(\frac{d}{2}-\delta +P\right)>\Delta \hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill {W}_{\mathrm{E}}=& \frac{{\sigma}^{2}}{{\delta}^{2}\epsilon {\epsilon}_{\mathrm{o}}{k}^{3}}\mathrm{sinh}\left(k\delta \right)({\mathrm{e}}^{k\delta -kd}-{\mathrm{e}}^{-kd})\phantom{\rule{1em}{0ex}}\text{when}\phantom{\rule{thinmathspace}{0ex}}d\ge 2\delta \hfill \\ \hfill {W}_{\mathrm{E}}=& \frac{{\sigma}^{2}}{{\delta}^{2}\epsilon {\epsilon}_{\mathrm{o}}{k}^{3}}(2k\delta -kd)-({\mathrm{e}}^{-k\delta}+1)\mathrm{sinh}(k\delta -kd)\hfill \\ \hfill & -\mathrm{sinh}\left(k\delta \right){\mathrm{e}}^{-kd}\phantom{\rule{1em}{0ex}}\text{when}\phantom{\rule{thinmathspace}{0ex}}d<2\delta ,\hfill \end{array}$$

(3)

where π, Δ, *d, δ* and *P* are the osmotic pressure term, depletion thickness, intercellular distance, RBC glycocalyx thickness (5 nm) and penetration depth; *σ, ε, ε*_{o} and *k* are the surface charge density of RBC, the relative permittivity of the solvent, the permittivity of the vacuum and the Debye–Huckel length, respectively. The above formulation yields an optimum polymer concentration for interaction energy because the penetration depth, *P*, is an inverse exponential function of the polymer concentration. The interaction energy increases, reaches the maximum and decreases as the concentration increases. Furthermore, due to the electrostatic energy, interaction energy has an optimum with respect to the intercellular distance. The interaction energy gradually increases to a maximum and decreases to zero as two surfaces approach each other. A strong repulsive force results when the distance between the surfaces is close to the sum of the thicknesses of RBC glycocalyx.

Bhattacharjee and collaborators approximated Derjaguin’s formula to find intercellular energy between two curved surfaces using the surface element integration (Bhattacharjee et al. 1998). Based on the surface element integration, total interaction energy can be expressed as

$$U={\int}_{{S}_{1}}({\mathbf{n}}_{1}\cdot {\mathbf{k}}_{1})({\mathbf{n}}_{2}\cdot {\mathbf{k}}_{2}){W}_{\mathrm{T}}\left(h\right)\mathrm{d}{S}_{1},$$

(4)

where **n** _{1} and **n** _{2} are the outward unit normal vectors of the surface elements of particles 1 and 2, respectively, and *S*_{1} is the surface of particle 1. Here, *h* is the local distance between the two surface elements, **k** _{1} and **k** _{2} are the unit normal vectors of two parallel planes with the direction parallel to a line *L*, which connects the centres of two particles. Therefore, **k** _{1} = −**k** _{2}. The detailed description is given in our earlier study (Chung et al. 2007).

Based on the depletion interaction energy, from Equations (1) to (4), an interaction force and torque between two RBCs were derived in Appendix A (see Equations (A6) and (A9)). The interaction force based on (Neu and Meiselman 2002) was modified using the Morse potential adopted by Liu et al. (2004). A force per unit surface by Liu et al. (2004) is represented by

$${\mathbf{F}}^{A}=A\left({\mathrm{exp}}^{2\beta ({h}_{\mathrm{o}}-h)}-{\mathrm{exp}}^{\beta ({h}_{\mathrm{o}}-h)}\right),$$

(5)

where *h*_{o} is the intercellular distance between the two surface elements when the force becomes zero and *A, β* are the parameters explicitly obtained from Equation (A6) with Equation (5) by comparing the optimum force and zero force occurred at two points at a distance. Figure 2 shows the force vs. intercellular distance. We adopted the force given by Equation (5) because the constant attractive force over a range of intercellular distances shown in Figure 2 was regarded as non-physical; the repulsive component of the force was modified computationally to prevent RBCs overlap as explained below. Torque in Equation (A9) is also modified using the modified force in Equation (5).

Because of the attractive force described in Equation (5) (see also Figure 2), two spheroids modelled as RBCs begin to approach each other as they are close enough to be in the region where the interaction energy acts. Total interaction force over the entire surface of an RBC may increase until the closest surface of the cell overlaps the other RBC. As described earlier, the minimum cellular distance (10 nm) is modelled to be the sum of the thicknesses of glycocalyx covering the entire surfaces of two cells. Accordingly, the cells cannot approach within 10 nm in the model. Energy distribution near the minimum distance shows that the electrostatic energy increases infinitely as the minimum distance is approached. Physically, when two cells approach each other, the interaction forces (attractive and repulsive) between the entire surfaces of cells reach equilibrium and the cells form an aggregate. Due to the computational limitation of the grid size, the interaction energy was scaled up to have the interaction force between cells within the computational interaction range equal to the force generated in the physiological range using an intercellular distance scale factor 50, as in our previous study (Chung et al. 2007). After the minimum distance between the cell surfaces reaches within our minimum computational length (500 nm), if the cells overlap each other, we imposed a repulsive force which has the same magnitude with opposite direction as that of the interaction force between the cells at the time, so that the total force became zero to prevent the overlap. At this point, we assume that the two RBCs form an aggregate. Based on this aggregation criterion, we obtained distances from the apex of bifurcation and time required for two cells to form an aggregate after leaving the daughter branches.

A bifurcation, which represents the entrance to a postcapillary venule, consists of two converging daughter branches representing capillaries and a straight segment. Figure 3(a) shows a representative bifurcation with two RBCs approaching the region of converging flows. The construction of a bifurcation is made by overlapping two individual diverging curved tubes. The combination of two types of finite volume cells, 5- and 6-sided cells, generates a tube. Each cell has generalised coordinates, which are used to discretise Navier–Stokes equations (Chung et al. 2007). Figure 3(b) shows the angles, *α*_{L} and *α*_{R}, between the centre axis of the parent vessel and of the two daughter branches. *D, D*_{L} and *D*_{R} in the figure denote the diameter of the parent vessel, the diameter of the left daughter branch and the diameter of the right daughter branch, respectively. We use

$${\alpha}_{\mathrm{L}}={\alpha}_{\mathrm{R}}=26.5\xb0,\phantom{\rule{1em}{0ex}}\frac{{D}_{\mathrm{L}}}{D}=\frac{{D}_{\mathrm{R}}}{D}=0.5,$$

(6)

which approximately correspond to the geometry of the postcapillary venules investigated in *in vivo* study (Kim et al. 2005). We assume that the daughter capillaries lie in the same plane with the venule. Accordingly, the longitudinal cut of the bifurcation matches the side view of the bifurcation.

The governing equations and method of numerical solution are presented in Appendix A. For the simulation of motion of two RBCs with and without aggregation in the bifurcation, no-slip conditions at the tube wall, Γ_{w}, and a fully developed flow condition at the outlet, Γ_{o} (parent branch), are specified

$$\mathbf{v}\mid {\Gamma}_{\mathrm{w}}=0,$$

(7)

$$\mathbf{v}\mid {\Gamma}_{\mathrm{o}}=-{V}_{\mathrm{max}}\left(1-\frac{{r}^{2}}{{R}^{2}}\right),$$

(8)

where V_{max}, *r* and *R* represent the maximum fluid velocity at the outlet, radial position and radius of the parent branch, respectively. At the inlets of two daughter branches, Γ_{i}, we specify zero axial velocity gradient

$$\frac{\partial \mathbf{v}}{\partial \mathbf{z}}\mid {\Gamma}_{\mathrm{i}}=0,$$

(9)

where *z* is the axial position in the daughter branch (see Figure 3(a)).

The boundary condition at the surface of the RBC, Γ_{c}, is specified using the cell’s translational and rotational velocities (see also Figure 3(a))

$$\mathbf{v}\mid {\Gamma}_{\mathrm{c}}=({\mathbf{v}}_{\mathrm{c}}+\omega \times {\mathbf{r}}_{\mathrm{c}})\mid {\Gamma}_{\mathrm{c}}.$$

(10)

Chimera grid method can employ a number of individual overset grids to construct a complex geometry such as a bifurcation in which multiple particles are placed. The communication between the grids can be made by identifying the interpolation points (fringe points). We consider two computational domains to solve the flow past a spherical or ellipsoidal particle in a tube. The computational domains represent the flow in the tube and the flow near the particle (Figure 4(a)). We define the tube domain to be the ‘major grid’ representing the fluid in the tube and spherical or ellipsoidal domain as the ‘minor grid’, placed over the major grid, representing the fluid near the particle. The minor grid includes the sphere or ellipsoid. These two independent grids can communicate through nonconservative trilinear interpolation at the boundary of each grid. The figure also shows ‘fringe points’ (solid circles) that surround the particle in the major grid. Fringe points of the minor grid are marked as open circles. Pressure and velocity are interpolated at these fringe points. The grid points on the surfaces overlapped by the two diverging curved tubes, producing a bifurcation are used as fringe points for trilinear interpolation between the two grids, representing the two tubes as shown in Figure 4(b). The figure shows the fringe points of the left or the right tubes. The fringe points are represented by the solid circles in the red coloured regions. The locations of fringe points of the left and right tubes are identical due to geometrical symmetry. The way to identify the fringe points between the grids representing RBCs when they are close each other is the same as that given in our previous study (Chung et al. 2007); the interpolation method is also explained in detail in that study.

Based on the geometric parameters in Equation (6), we simulated four cases under normal (control) and reduced flows with and without DEX 500, respectively. For the normal and the reduced flows, the pseudoshear rates were 250 and 25 s^{−1}, respectively; the corresponding pseudoshear rates in the experimental study (Kim et al. 2005) ranged between 240 and 27 s^{−1}. The corresponding Reynolds number

$$\mathit{Re}=\frac{\rho VD}{\mu}$$

(11)

for the normal and reduced flows was 0.1 and 0.01, respectively; here *V* is the average velocity at the outlet, Γ_{o}, and *D* is the diameter of the parent branch. Polymer concentration *C*_{p} was fixed at 15 kg m^{−3} for aggregation.

For a microvessel (venule) with 20 μm diameter, the average plasma velocities in venule are 0.5 cm s^{−1} and 0.05 cm s^{−1} at *Re* = 0.1 and 0.01, respectively. The dimensions for the model are similar to those found in postcapillary venules where aggregates first form *in vivo* (Kim et al. 2005). Flow fields in a bifurcation without cells at *Re* = 0.1 and 0.01 are presented in Appendix B to show where the flow is nearly fully developed in the parent vessel.

The initial positions and angles of cells were chosen based on a strategy that accounts for distinct events that might occur along with flows in the capillary. For example, we account for three possible cases for the cell’s angle: one for the longest axis of cell parallel to the flow direction, one for the longest axis of the cell perpendicular to the flow direction, and one for the longest axis of the cell tilted with respect to the flow direction. Furthermore, there are three possible choices for the cell’s radial position: one for the centre of RBC on the centre axis of daughter branch, the others for the centre of RBC off the centre axis of daughter branch (left or right); and finally two possible choices for axial distance between two cells: one for RBCs positioned at the same axial distance, the other for RBCs positioned at different axial distances. Based on the strategy in conjunction with the consideration that RBCs are likely to form an aggregate only when they are close to each other in the bifurcation region, we chose the initial positions and angles of cells by the combination of the above possible events to consider two cases for RBC motion without aggregation and four cases for RBC motion with aggregation.

In this case, *Re* was fixed at 0.1. Each cell was initially positioned near the inner surface of the daughter branch as shown in Figure 5(a). The initial positions and tilted angles of two cells were symmetrical. The initial velocity profile is also shown. Due to the difference between the fluid velocities in the parent and daughter branches, the cell velocity decreases after the cell enters the bifurcation region. Figure 5(b) shows cell trajectories with time. The two cells start to approach each other after the cells leave the daughter branches because the fluid streams from the daughter branches converge at the centreline of parent branch near the apex of bifurcation. The distance between the cells then remains nearly constant when the cells reach the far downstream region where the flow is fully developed.

Cells were initially located asymmetrically in this case. Cell A was placed more upstream in the right daughter branch compared to the location of Cell B. In addition, Cell B was closer to the inner surface of the left daughter branch. The simulation was performed at *Re* = 0.1. Initial velocity profile is shown in Figure 6(a). The velocity of Cell A is larger than that of Cell B due to the parabolic velocity profiles in daughter branches. Thus, Cell A catches up with Cell B after entering the region of bifurcation. The trajectories of cells are given in Figure 6(b); the velocities of Cells A and B decrease after leaving the daughter branches and then, the velocity of Cell B becomes larger than that of Cell A where the flow is fully developed. This is because Cell B moves closer to the centreline of the parent branch due to its initial position.

The pseudoshear rate was reduced from 250 to 25 s^{−1}. The corresponding *Re* was 0.01. Interaction energy between two cells was not considered. The initial positions of the cells were identical to those of Case 1 for normal flow without DEX. The behaviour of cells is very similar to that of the cells in the previous case (Figure 7(a)). The velocities of the cells are identical due to the symmetry. The velocities decrease in the bifurcation region and become constant.

Using the pseudoshear rate 25 s^{−1}, asymmetrical initial positions of two cells were considered. The positions of the cells were same as those of Case 2 in normal flow without DEX (Figure 7(b)). Like the previous case, Cell A catches up with Cell B in the region of the bifurcation, and Cell B moves downstream faster than Cell A since Cell B moves close to the centreline of the parent branch in the fully developed region as shown in the figure.

Considering an interaction energy between the cells, in the presence of DEX 500, we investigated aggregation of two cells initially positioned close to the inner surfaces of the left and right daughter branches at *Re* = 0.1 for Case 1. The cells enter the region of bifurcation symmetrically as shown in Figure 8(a). When the cells reach the range where the interaction forces are effective, the aggregation starts. Figure 9 demonstrates how aggregation forces act between the two cells to form an RBC aggregate over time for the four cases that we considered for the normal flow. Based on the aggregate formation criterion introduced above, the cells converge into the bifurcation region due to the flow and start to form an equilibrium configuration after the initiation of aggregation process. The process ends at a distance about 8 μm from the apex at *t* = 8 × 10^{−4} s (Case 1 in Figure 9). Due to the continuity of the fluid flow, the velocities of the cells decrease when they reach the bifurcation and then become constant in the fully developed region. Note that there is no difference between the velocities of the cells at all times because of the symmetry of the initial position of cells.

In Case 2, two cells were initially located asymmetrically as shown in Figure 8(b), similar to Case 2 of normal flow without DEX. Like the behaviour of the cells in the previous case, Cell A catches up with Cell B in the region of the bifurcation. The velocity of Cell A is higher than that of Cell B until they approach within the range of interaction energy as in the previous case. However, once the cells form an aggregate, the velocities of the cells become almost equal as in Figure 8(b). The aggregate was formed off the centreline of the parent branch because of the asymmetrical initial positions of the cells in this case. The force curve as a function of time for Case 2 in Figure 9 illustrates that the cells begin to aggregate and form an aggregate at a distance about 8 μm from the apex at *t* = 9.6 × 10^{−4} s in this case.

In Case 3, Cell A was initially located at the centreline of the right branch, and Cell B was positioned close to the inner surface of the left branch so that the velocity of Cell B can be initially lower than that of Cell A as shown in Figure 8(c). The cells move into the bifurcation region with different velocities, and their velocity becomes nearly constant after the initiation of the aggregation process. The force curve representing Case 3 in Figure 9 shows that the aggregation process initiates at *t* = 7.2 × 10^{−4} s and ends at the distance about 15 μm from the apex at *t* = 16.8 × 10^{−4} s in this case.

In Case 4 for the normal flow, the cells were positioned asymmetrically as shown in Figure 10(b), which illustrates cell trajectories. Cell A was located further downstream than Cell B, and both cells were initially close to the centrelines of the daughter branches. Unlike the previous cases, the angle of the position of Cell A was tilted so that the larger surface of the cell faces the flow. The cells move into the bifurcation, and converge through the centreline of the parent bifurcation as shown in Figure 10(b). Figure 10(a) shows a velocity profile at the intermediate time. As shown in the force curve of Case 4 in Figure 9, the cells initiate an aggregation process when they are sufficiently close to each other at *t* = 7.2 × 10^{−4} s. However, unlike the previous cases, the cells do not aggregate but are separated. As in Figure 9, unlike the other curves, the force curve representing Case 4 shows that the force increases and then, slowly decreases to zero, which demonstrates that a mechanical force exceeds an attractive force between the cells; thus, the cells do not aggregate.

Based on the interaction energy corresponding to the presence of DEX 500, the aggregation behaviour of two cells was investigated at *Re* = 0.01. The initial locations of two cells for four cases (Cases 1–4) were identical to those of the cells in the cases of normal flow with DEX. For Case 1, the cells approach each other in the bifurcation region, and then an aggregation force between the cells begins to act when the surfaces of two cells are within the interaction range. Figure 11(a) shows the trajectories of the cells. Like Case 1 of normal flow with DEX, the cells form an aggregate to leave the bifurcation region. The force curve for Case 1 is shown in Figure 12. The cells form an aggregate at a distance 8 μm from the apex of the bifurcation at *t* = 7.2 × 10^{−3} s.

For Case 2, the cells were initially located at the same initial positions as Case 2 in the normal flow with DEX 500. As shown in Figure 11(b), Cell A in the right branch moves into the bifurcation faster than Cell B, and they come closer to each other in the region where an attractive force acts. They aggregate at a distance 8 μm from the apex at *t* = 8.8 × 10^{−3} s as shown in the force curve for Case 2 in Figure 12 and leave the bifurcation region with the same velocity.

In Case 3, Cell A was initially located close to the centre of the right branch, and Cell B was positioned closer to the inner surface of the left branch. As in Case 3 in the normal flow with DEX 500, due to the flow patterns, Cell B moves slower than Cell A. The cells approach each other due to the difference between cell velocities and the flow pattern in the bifurcation so that they begin to aggregate in the interaction range. The force curve of Case 3 in Figure 12 shows that the two cells aggregate at a distance about 14 μm from the apex at *t* = 16.8 × 10^{−3} s.

In Case 4, two cells were initially at the same positions as those of the cells in Case 4 in the normal flow with DEX 500. As shown in Figure 13(b), the cells move into the bifurcation region, and then come closer each other to be in the interaction range. Figure 12 shows that they start to aggregate at *t* = 7.2 × 10^{−3} s and form an aggregate at a distance about 10 μm from the apex at *t* = 12 × 10^{−3} s. Unlike the behaviours of two cells in the normal flow with DEX, the cells aggregate, because the attractive force is larger than the mechanical force. The aggregability of the cells in the most cases (Cases 1–3) were nearly similar but, the aggregability was dependent on the initial locations and angles of the cells as shown in the results of Case 4 in both the normal and reduced flows with DEX. Therefore, under reduced flows, cells likely have a higher possibility to form an aggregate in the bifurcation region, depending on both cells locations and orientations before cells enter the parent branch.

Given the venule geometry modelled from the previous experimental data, the RBC flow from the capillaries to the venule was simulated. The RBC was modelled as an ellipsoidal oblate spheroid with an eccentricity 0.3 that mimics the shape of a rat RBC. Two spheroids representing RBCs were initially located in the daughter branches (capillaries) and suspended in plasma with and without DEX 500.

The computational results show that RBCs have higher tendencies to form an aggregate under the reduced flow. This phenomenon occurs because the ratio of the mechanical forces that tends to separate the cells, to aggregation force is approximately ten times smaller than that under the normal flow. However, our results show that aggregation also depends on the angle and location of the cells before the cells enter the bifurcation region, and in the most cases under the normal and reduced flows, the cells had a tendency to form an aggregate. Since cells experience higher shear forces near the vessel wall due to velocity profile, the aggregation of cells in the region near the centreline may not be highly affected by Reynolds number. For this reason, the effect of initial angles of the cells must be one of critical factors that affect aggregation of the cells because of the difference in mechanical interaction forces between the cells. The differences in the forces due to the initial angles of the cells may result in non-aggregation between cells under high *Re* flows even if cells are suspended near the centreline of the tube. Therefore, cells suspended under normal flows have a smaller tendency to aggregate. This is a qualitative agreement with the experimental data. Our earlier *in vivo* study showed that cells aggregated approximately 15–30 μm from the apex of the bifurcation, which has been proposed to be the distances where the flow becomes fully developed in the bifurcation. However, the computational results showed that in most cases for aggregation, the cells formed an aggregate at distances approximately 8–15 μm from the apex, where the flow is still developing. The difference in the results may occur because of the bifurcation geometry. The distance where the diameter of the venule *in vivo* becomes nearly constant from the apex of the bifurcation was longer than that of the computational geometry whereas the diameter of the computational parent branch was nearly constant from the apex. Hence, cells in *in vivo* geometries in the experiment had likely longer distances to come closer to each other to be within an interaction range between cells.

Because we considered three-dimensional geometry with associated computational expense and a limit on the grid size, scaled interaction energy was used. This scaling factor may overestimate *in vivo* intercellular energy between RBCs because in reality attractive forces between cells act only in a local region within about 50 nm intercellular distance between the surfaces of cells. The scaling energy or force then covers larger surfaces of an RBC so that it causes overestimation in total force acting on the cell. Nevertheless, the present simulation agrees qualitatively with our earlier *in vivo* experimental data. The role of mechanical forces which may separate the cells is found to depend upon the location and the angles of the cells from this study. The additional consideration of effects of interactions between more than two cells in a bifurcation and the deformability of RBCs should provide further information on aggregation mechanisms *in vivo*. However, within the limitations of our computational study, the three-dimensional simulation of RBC aggregation in a postcapillary venular bifurcation using an idealised computational geometry based on the Chimera grid method is in qualitative agreement with our *in vivo* experimental study. The Chimera grid method could be extended and modified to take into account RBC deformability; thus, the present study provides a methodological framework for further developments.

This work was supported by NIH grants NHLBI RO1 HL52684 and HL18292. P.C. Johnson is also a Senior Scientist at La Jolla Bioengineering Institute. We thank Dr Junfeng Zhang for reading the manuscripts and critical comments.

The details of the formulation and numerical analysis are presented in Chung et al. (2007); here we briefly outline the methodology. The fluid representing the blood plasma can be considered to be incompressible and Newtonian, so that the equations governing the flow can be written in the integral forms as

$$\int \mathbf{v}\cdot \mathbf{n}\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}A=0,$$

(A1)

$$\int \rho \frac{\partial \mathbf{v}}{\partial t}\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}V+\int \rho (\mathbf{v}\cdot \nabla )\mathbf{v}\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}V=-\int p\mathbf{n}\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}A+\int \tau \cdot \mathbf{n}\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}A,$$

(A2)

where **v**, *p*, *ρ* and **τ** are the velocity, pressure, density of the fluid and stress tensor, respectively. The stress tensor is given by

$$\tau =\mu (\nabla \mathbf{v}+\nabla {\mathbf{v}}^{T}),$$

(A3)

where *μ* is the viscosity of the plasma.

RBCs suspended in blood flow experience a drag force due to the flow and an interaction force between the two cells when they are close to each other. By ignoring body force, the equation of motion for the particle is

$$m\frac{\mathrm{d}{\mathbf{v}}_{\mathrm{c}}}{\mathrm{d}t}={\mathbf{F}}^{D}+{\mathbf{F}}^{A}=\mathbf{F},$$

(A4)

where *m* and **v**_{c} are the mass and translational velocity of the RBC, respectively, the drag force **F*** ^{D}* exerted by the fluid on the particle is in the form

$${\mathbf{F}}^{D}=-\int p\mathbf{n}\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}A+\int \tau \cdot \mathbf{n}\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}A$$

(A5)

and **F*** ^{A}* is the interaction force based on the depletion interaction model (Neu and Meiselman 2002). The force can be obtained by differentiating the interaction energy (Equation (4)) with respect to

$${\mathbf{F}}_{1}^{A}=-\frac{\partial U}{\partial \mathbf{h}}=-{\int}_{{Q}_{1}}({\mathbf{n}}_{2}\cdot {\mathbf{k}}_{2})\frac{\partial {W}_{\mathrm{T}}\left(h\right)}{\partial \mathbf{h}}\mathrm{d}{Q}_{1},$$

(A6)

where d*Q*_{1} is the projected surface of the surface element d*S*_{1} in Equation (4).

The moment of momentum equation for the particle is written in the form

$$I\frac{\mathrm{d}\omega}{\mathrm{d}t}+\omega \times I\cdot \omega ={\mathbf{T}}^{D}+{\mathbf{T}}^{A}=\mathbf{T},$$

(A7)

where *I* and **ω** are the moment of inertia and angular velocity of the cell, respectively, and **T*** ^{A}* is the torque due to the interaction force, and the torque

$${\mathbf{T}}^{D}=-\int {\mathbf{r}}_{\mathrm{c}}\times p\mathbf{n}\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}A+\int {\mathbf{r}}_{\mathrm{c}}\times \tau \cdot \mathbf{n}\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}A,$$

(A8)

where **r**_{c} is the position vector from the centre of mass and the surface of the cell. Torque **T*** ^{A}* based on the depletion interaction model can be obtained using Equation (A6)

$${\mathbf{T}}_{1}^{A}=-{\int}_{{Q}_{1}}{\mathbf{r}}_{\mathrm{c}1}\times ({\mathbf{n}}_{2}\cdot {\mathbf{k}}_{2})\frac{\partial {W}_{\mathrm{T}}\left(h\right)}{\partial \mathbf{h}}\mathrm{d}{Q}_{1},$$

(A9)

where **r**_{c1} is the distance vector from the centre of cell 1 to d*Q*_{1}.

Governing equations for fluid (Equations (A1)–(A3)) are solved by Chorin’s pressure correction method after discretising the equations based on the generalised coordinates. The nonlinear discretised equations are then linearised by Newton-Raphson method. The detailed explanation is given in Chung et al. (2007).

The translational and angular velocities of each cell are determined by the resultant force and torque in Equations (A4) and (A7). Numerically, we discretise these two equations based on the algorithm to obtain the negligible force and torque on a particle at a given time using a larger time step so that the particle’s acceleration becomes zero during the time step, which was proposed in our previous numerical study (Chung et al. 2007). The current particle translational velocities (Equation (A4)) are obtained by

$${\mathbf{v}}_{\mathrm{c}}^{n+1,m+1}={\mathbf{v}}_{\mathrm{c}}^{n+1,m}+\frac{\Delta {t}^{\prime}}{m}\frac{\mathbf{F}({\mathbf{v}}^{n+1,m+1},{p}^{n+1,m+1})+\mathbf{F}({\mathbf{v}}^{n},{p}^{n})}{2},$$

(A10)

where the superscripts *n* and *n* + 1 stand for old and new time, *m* and *m* + 1 stand for the previous and current iterations, respectively, and Δ *t* ′ is a time step whose value can be arbitrarily chosen. Likewise, Equation (A7) is solved by a fourth-order Runge Kutta method to obtain the current **ω**^{m+1} using the previous **ω*** ^{m}*, Δ

Navier–Stokes equations and the equations of motion of cells with including the interaction energy between cells are solved based on the iterative steps given in Chung et al. (2007). Tolerance given to solve the overall equations in the present study was 10^{−3}.

Flow fields in the bifurcation at *Re* = 0.1 and 0.01 are shown in Figures FiguresB1B1 and andB2,B2, respectively. As expected, the flows are similar and the velocity magnitude differ by a factor of 10. The flow rate *Q* at the outlet Γ_{o} must be equal to the sum of *Q*_{L} and *Q*_{R} by the continuity so that it satisfies

$$VA={V}_{\mathrm{L}}{A}_{\mathrm{L}}+{V}_{\mathrm{R}}{A}_{\mathrm{R}},$$

(B1)

where

$${A}_{\mathrm{L}}={A}_{\mathrm{R}}=\frac{A}{4}$$

(B2)

and *V, V*_{L} and *V*_{R} are the average velocity in parent, left and right daughter branches, respectively, and *A, A*_{L}, *A*_{R} are the area of parent, left and right daughter branches, respectively. Thus

$$V=\frac{{V}_{\mathrm{L}}+{V}_{\mathrm{R}}}{4}$$

(B3)

Figures FiguresB1B1 and andB2B2 also show the magnitude of average velocity in the bifurcation at *Re* = 0.1 and 0.01. In the downstream region of the parent branch where the flow is fully developed,

$${V}_{\mathrm{max}}^{\mathit{fd}}=2V$$

(B4)

so that

$$\begin{array}{cc}\hfill & {V}_{\mathrm{max}}^{\mathit{fd}}=1\phantom{\rule{thinmathspace}{0ex}}\mathrm{cm}\phantom{\rule{thinmathspace}{0ex}}{\mathrm{s}}^{-1},\phantom{\rule{1em}{0ex}}{V}_{\mathrm{max}}^{\mathit{fd}}=0.1\phantom{\rule{thinmathspace}{0ex}}\mathrm{cm}\phantom{\rule{thinmathspace}{0ex}}{\mathrm{s}}^{-1}\hfill \\ \hfill & \phantom{\rule{thinmathspace}{0ex}}\mathit{Re}=0.1,0.01,\hfill \end{array}$$

(B5)

where ${V}_{\mathrm{max}}^{\mathit{fd}}$ is the maximum velocity in the region where the flow is nearly fully developed. The fluid from the left and right daughter branches flows into the bifurcation region. The flow in the bifurcation region was considered to be nearly fully developed at the distance between 17.2 and 20.2 μm from the apex of the bifurcation based on the mesh density used in this study, at both *Re* = 0.1 and 0.01 using the condition

$$\mid \frac{{V}_{\mathrm{max}}^{\mathit{fd}}-{V}_{\mathrm{c}}}{{V}_{\mathrm{max}}^{\mathit{fd}}}\mid \le 0.01\phantom{\rule{1em}{0ex}}\text{for fully developed region},$$

(B6)

where *V*_{c} is the centreline velocity of the parent branch.

- Bagchi P. Mesoscale simulation of blood flow in small vessels. Biophys J. 2007;92(6):1858–1877. [PubMed]
- Bagchi P, Johnson PC, Popel AS. Computational fluid dynamic simulation of aggregation of deformable cells in a shear flow. J Biomech Eng. 2005;127(6):1070–1080. [PubMed]
- Bhattacharjee S, Elimelech M, Borkovecb M. DLVO Interaction between colloidal particles: beyond Derjaguin’s approximation. Croat Chem Acta CCACAA. 1998;71(4):883–903.
- Bishop JJ, Popel AS, Intaglietta M, Johnson PC. Rheological effects of red blood cell aggregation in the venous network: a review of recent studies. Biorheology. 2001;38(2–3):263–274. [PubMed]
- Bishop JJ, Popel AS, Intaglietta M, Johnson PC. Effects of erythrocyte aggregation and venous network geometry on red blood cell axial migration. Am J Physiol Heart Circ Physiol. 2001;281(2):H939–H950. [PubMed]
- Bishop JJ, Nance PR, Popel AS, Intaglietta M, Johnson PC. Relationship between erythrocyte aggregate size and flow rate in skeletal muscle venules. Am J Physiol Heart Circ Physiol. 2004;286(1):H113–H120. [PubMed]
- Chung B, Johnson PC, Popel AS. Application of Chimera grid to modeling cell motion and aggregation in a narrow tube. Int J Numer Methods Fluids. 2007;53:105–128.
- Das B, Bishop JJ, Kim S, Meiselman HJ, Johnson PC, Popel AS. Red blood cell velocity profiles in skeletal muscle venules at low flow rates are described by the Casson model. Clin Hemorheol Microcirc. 2007;36(3):217–233. [PubMed]
- Dougherty FC. Development of a Chimera grid scheme with applications to unsteady problems [PhD thesis. Stanford University; Stanford (CA): 1985.
- Eggleton CD, Popel AS. Large deformation of red blood cell ghosts in a simple shear flow. Phys Fluids. 1998;10(8):1834–1845.
- Kim S, Popel AS, Intaglietta M, Johnson PC. Aggregate formation of erythrocytes in postcapillary venules. Am J Physiol Heart Circ Physiol. 2005;288(2):H584–H590. [PubMed]
- Liu Y, Zhang L, Wang X, Liu WK. Coupling of Navier–Stokes equations with protein molecular dynamics and its application to hemodynamics. Int J Numer Methods Fluids. 2004;46:1237–1252.
- Neu B, Meiselman HJ. Depletion-mediated red blood cell aggregation in polymer solutions. Biophys J. 2002;83(5):2482–2490. [PubMed]
- Popel AS, Johnson PC. Microcirculation and hemorheology. Annu Rev Fluid Mech. 2005;37:43–69. [PMC free article] [PubMed]
- Pozrikidis C. Effect of membrane bending stiffness on the deformation of capsules in simple shear flow. J Fluid Mech. 2001;440:269–291.
- Rampling MW, Meiselman HJ, Neu B, Baskurt OK. Influence of cell-specific factors on red blood cell aggregation. Biorheology. 2004;41(2):91–112. [PubMed]
- Sun CH, Munn LL. Particulate nature of blood determines macroscopic rheology: a 2-D lattice Boltzmann analysis. Biophys J. 2005;88(3):1635–1645. [PubMed]
- Sun CH, Migliorini C, Munn LL. Red blood cells initiate leukocyte rolling in postcapillary expansions: a lattice Boltzmann analysis. Biophys J. 2003;85(1):208–222. [PubMed]
- Zhang J, Johnson PC, Popel AS. Lattice Boltzmann simulations of multiple deformable red blood cells in shear flows. J Biomech. 2007;41:47–55. [PMC free article] [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. |