|Home | About | Journals | Submit | Contact Us | Français|
In this paper we develop a lattice Boltzmann algorithm to simulate red blood cell (RBC) behavior in shear flows. The immersed boundary method is employed to incorporate the fluid-membrane interaction between the flow field and deformable cells. The cell membrane is treated as a neo-Hookean viscoelastic material and a Morse potential is adopted to model the intercellular interaction. Utilizing the available mechanical properties of RBCs, multiple cells have been studied in shear flows using a two-dimensional approximation. These cells aggregate and form a rouleau under the action of intercellular interaction. The equilibrium configuration is related to the interaction strength. The end cells exhibit concave shapes under weak interaction and convex shapes under strong interaction. In shear flows, such a rouleau-like aggregate will rotate or be separated, depending on the relative strengths of the intercellular interaction and hydrodynamic viscous forces. These behaviors are qualitatively similar to experimental observations and show the potential of this numerical scheme for future studies of blood flow in microvessels.
Red blood cells (RBCs) are an important component in blood because of their large number density (~ 5×106/mm3) and their crucial role in oxygen transport. Typically, a human RBC has a biconcave shape of ~ 8 μm in diameter and ~ 2 μm in thickness. The interior fluid (cytoplasm) has a viscosity of 6 cP, which is about 5 times of that of the suspending plasma (~ 1.2 cP). The cell membrane is highly deformable so that RBCs can pass through capillaries of as small as 4 μm inner diameter with large deformation (Popel and Johnson, 2005; Mchedlishvili and Maeda, 2001). RBCs can also aggregate and form one-dimensional stacks-of-coins-like rouleaux or three-dimensional aggregates (Popel and Johnson, 2005; Stoltz et al., 1999; Baumler et al., 1999). The process is reversible and the rouleaux and aggregates can be broken by, for example, increasing the flow shear rate. This phenomenon is particularly important in microcirculation, since such rouleaux or aggregates can dramatically influence blood flow in microvessels. However, the underlying mechanism of RBC aggregation is not yet clear. Currently, there exist two theoretical descriptions of this process: the bridging model and the depletion model (Popel and Johnson, 2005; Baumler et al., 1999). The former assumes that macromolecules, such as fibrinogen or dextran, can adhere onto the adjacent RBC surfaces and bridge them together (Merill et al., 1966; Brooks, 1973; Chien and Jan, 1973a). The depletion model attributes the RBC aggregation to a polymer depletion layer between RBC surfaces, which is accompanied by a decrease of the osmotic pressure (Baumler and Donath, 1987; Evans and Needham, 1988). Detailed discussions of these two models can be found elsewhere [see, for example, a review by Baumler et al. (1999)].
The nature of blood flow changes greatly with the vessel diameter. In vessels larger than 200 μm, the blood flow can be accurately modeled as a homogeneous fluid. However, in arterioles and venules smaller than 25 μm, and also in capillaries with diameter of 4–10 μm, the RBCs have to be treated as discrete fluid capsules suspended in the plasma. Significant efforts have been devoted to numerically study the RBC behaviors in various flow situations. For example, Pozrikidis and coworkers (1995; 2001; 2003) have employed the boundary integral method for Stokes flows to investigate RBC deformation and motion in shear and channel flows; Eggleton and Popel (1998) have combined the immersed boundary method (IBM) (Peskin, 1977) with a finite element treatment of the RBC membrane to simulate large three-dimensional RBC deformations in shear flow. Recently, a lattice Boltzmann approach has also been adopted for RBC flows in microvessels, where the RBCs were represented as rigid rods in two dimensions (Migliorini et al., 2002; Sun et al., 2003; Sun and Munn, 2005). Bagchi (2007) has simulated a large RBC population in vessels of size 20-300 μm. However, RBC aggregation was not considered in these studies. Liu et al. (2004) modeled the intercellular interaction through a Morse potential, thus accounting for RBC aggregation. The cell membrane was represented by elastic elements. Bagchi et al. (2005) have extended the IBM approach of Eggleton and Popel (1998) to two-cell systems and introduced the intercellular interaction according to a legand-receptor binding model. Chung et al. (2006) have utilized the theoretical formulation of depletion energy proposed by Neu and Meiselman (2002) to study two rigid elliptical particles in a channel flow. Sun and Munn (2006) have also improved their lattice Boltzmann model by including an interaction force between rigid RBCs.
In this paper, we develop a two-dimensional lattice Boltzmann scheme to simulate multiple deformable RBCs in shear flow. The lattice Boltzmann method (LBM) was chosen to solve the incompressible fluid field for its ability to deal with complex boundary conditions and its advantage for parallel computation (Succi, 2001), both of which could be valuable in our future studies. Different from previous LBM studies of RBC flows (Migliorini et al., 2002; Sun et al., 2003; Sun and Munn, 2005), here we modeled the cells as deformable fluid capsules and hence the membrane mechanics and cytoplasm viscosity could be considered. IBM was also utilized to incorporate the fluid-membrane interaction, and a Morse potential to describe the intercellular interaction. Detailed theory and formulations will be given in the next Section. The algorithm and program have been validated by comparing simulation results with theoretical predictions, and excellent agreement was found. Finally, simulations of multiple deformable RBCs were conducted and the results demonstrated the effects of intercellular interactions and shear rate on the RBC rheological behaviors, which are also in qualitative agreement with experimental observations.
Here we employ the lattice Boltzmann model (LBM) (Zhang and Kwok, 2005; Zhang et al., 2004) to solve the flow field and the immersed boundary method (Peskin, 1977) to incorporate the fluid-membrane interaction. The LBM approach is advantageous for parallel computations due to its local dynamics and can be relatively easily applied to systems with complex boundaries (Succi, 2001). Detailed descriptions of these methods are available elsewhere (also see the Supplemental Materials).
Experiments have shown that the RBC membrane is a neo-Hookean, highly deformable viscoelastic material (Hochmuth and Waugh, 1987). Also it exhibits finite bending resistance that becomes more profound in regions with large curvature (Evans, 1983). For the two-dimensional RBC model in this work, according to Bagchi et al. (2005), the neo-Hookean elastic component of membrane stress can be expressed as
where Es is the membrane shear modulus and ε is the stretch ratio. To incorporate the membrane viscous effect, an extra term should be added to the stress expression as (Evans and Hochmuth, 1976; Mills et al., 2004; Bagchi et al., 2005)
where μm is the membrane viscosity. In addition, the bending resistance can also be represented by relating the membrane curvature change to the membrane stress with (Pozrikidis, 2001; Bagchi et al., 2005)
where Eb is the bending modulus and k and k0 are, respectively, the instantaneous and initial stress-free membrane curvatures. l is a measure of the arc length along the membrane surface. Therefore, the total membrane stress T induced due to the cell deformation is a sum of the three terms discussed above:
Here t and n are the local tangential and normal directions on the membrane.
The physiological and pathological importance of RBC aggregation has been realized and extensive experimental investigations have been performed (Popel and Johnson, 2005; Stoltz et al., 1999; Baumler et al., 1999; Kim et al., 2006; Kounov and Petrov, 1999; Chien et al., 1977); however, the underlying mechanisms of the RBC aggregation are still subjects of investigation. Both the bridging and the depletion models can describe certain aggregation phenomena, however, they fail to explain some specific observations (Baumler et al., 1999; Armstrong et al., 1999). In general, it can be assumed that the attractive interaction between RBC surfaces would occur when they are close, and repulsive interaction would occur when the separation distance is sufficiently small. The repulsive interaction includes the steric forces due to the glycocalyx and electrostatic repulsion from the negative charges on RBC surfaces (Liu et al., 2004). In previous studies, to model such intercellular interactions, Bagchi et al. (2005) adopted the formalism of a ligand-receptor dynamics according to the bridging model. However, this description involves a number of parameters whose values are not available from experiments. On the other hand, Chung et al. (2006) employed a mathematical description of the depletion model proposed by Neu and Meiselman (2002). It was noticed that such a description results in a constant (instead of a decaying) attractive force at large separation, which is not physically realistic.
In this study, we follow the approach recently proposed by Liu et al. (2004) to model the intercellular interaction energy ϕ as a Morse potential:
where r is the surface separation, r0 and De are, respectively, the zero force separation and surface energy, and β is a scaling factor controlling the interaction decay behavior. The interaction force from such a potential is its negative derivative, i.e., f(r) = −ϕ/r. This Morse potential is selected for its simplicity and realistic physical description of the cell-cell interaction; however, it does not imply any underlying mechanism. Nevertheless, other interaction models can be readily incorporated into our algorithm.
Before applying the methods described above to RBC flows, several validation simuations have been carried out, and excellent agreement with known theoretical results was observed. All simulations in this paper are performed with a D2Q9 lattice structure (Zhang and Kwok, 2006). The biconcave shape of an RBC is described by the ollowing equation suggested by Evans and Fung (1972):
where c0=0.207, c1=2.002, and c2=1.122. The non-dimensional coordinates (, ) are normalized by the radius of a human RBC (3.91 μm). The cell membrane is represented by 100 segments and the simulation domain is a 256×128 rectangle in lattice units (50 μm×25 μm, 0.196 μm per lattice unit). For these simulations, we applied periodic boundary condition in the horizontal (x) direction and modified bounce-back boundary conditions (Nguyen and Ladd, 2005) on the top and bottom to generate a shear flow field. Other simulation parameters employed in the present work are also listed in Table 1 with references. The intercellular interactions are represented as a Morse potential and evaluated through the Derjaguin approximation (Israelachvili, 1991) for cells with arbitrary shapes. For the parameters in the Morse potential Eq.(5) for the intercellular interaction, experimental values are not available. According to Neu and Meiselman (2002), the zero force distance r0 is ~14 nm. In our simulations, r0 is chosen as 2.5 lattice units (0.49 μm) due to the limitation of the grid size and the immersed boundary treatment employed in the computational scheme. This distance is relatively large compared to experimental data (Chien and Jan, 1973b) and the scheme could be improved by adopting a smaller lattice grid and/or the adaptive mesh refinement technology (Toolke et al., 2006). To improve the computational efficiency, as is typical in molecular dynamics simulations, a cut-off distance is adopted. The scaling factor β is adjusted so that the attractive force beyond the cut-off distance rc=3.52 μm is neglected (β = 3.84 μm−1). Since no experimental value is available, a range of interaction strengths De has been examined; and results with two specific values (5.2×10−8 and 1.3×10−6 μJ/μm2) will be presented to demonstrate the effect of intercellular interaction on the cell rheological behaviors in shear flows. As can be seen below, these two specific values are selected for the representative cell behaviors at the two shear rates γ = 20 and 50 s−1, which are in the range of typical hemorheologic experiments (Stoltz et al., 1999). A typical simulation of four cells in a shear flow takes ~10 MB memory and runs 0.036 s per time step on a 3 GHz Pentium IV desktop computer.
It has long been recognized that RBCs can aggregate to form stacks-of-coins-like rouleaux in the presence of fibrinogen or other macromolecules (Chien and Jan, 1973b; Skalak and Chien, 1987). In this section, we investigate four RBCs with different interaction strengths De under a no-flow condition. Similar to those employed in previous work (Bagchi et al., 2005; Liu et al., 2004), the cells are initially placed in the domain with a center-to-center distance of 3.52 μm (see Figure 2a). The closest gap between the rims of two adjacent cells is 0.95 μm so that the interaction between them is attractive. Under this intercellular interaction, the cells aggregate slowly and an equilibrium configuration is reached when the interaction is balanced with the membrane forces after cell deformation. Figure 1 shows the rouleaux shapes as a function of the interaction strength De. The relatively uniform intercellular distances between adjacent cells are in accordance with the experimental observation by Chien and Jan (1973b). The end RBCs display concave shapes at low interaction strengths (Figure 1a-c) and convex shapes at high interaction strengths (Figure 1d-f). At very high interaction strengths (Figure 1f), the rouleau contracts to a nearly spherical clump to maximize the contact area. Also we observe flat contact surfaces at low interaction strengths and curved contact surfaces at high interaction strengths. These findings are in qualitative agreement with previous experimental studies (Skalak and Chien, 1987), and hence demonstrate that the Morse potential model is a reasonable choice for the intercellular interaction.
The complex RBC behaviors in flow fields result from the highly coupled interplay of the fluid viscous stresses due to the flow, the membrane forces due to cell deformation, and intercellular interactions. In this section, we will quantitatively investigate these relationships by simulating a system of four biconcave cells in shear flow, which is similar to that observed in hemorheological experiments (Stoltz et al., 1999). The simulation strategy is similar to those employed in previous work (Bagchi et al., 2005; Liu et al., 2004). Rouleaux, formed under intercellular interactions in no-flow condition, as described above, are placed at the center of the simulation domain. A shear flow is applied by specifying the top and bottom boundary velocities, which have the same value but in opposite horizontal directions, in the modified bounce-back boundary conditions (Nguyen and Ladd, 2005). Different shear rates can be obtained by adjusting the boundary velocity.
Experiments show that RBC aggregates can be broken into smaller parts or even individual cells at high shear rate (Popel and Johnson, 2005; Stoltz et al., 1999). Figures Figures22--44 display representative snapshots of the RBC deformation and movement in a shear flow with shear rate γ=20 s−1. The slight asymmetry in streamlines is due to the plotting software. The strength of intercellular interaction in the Morse potential Eq.(5) are, respectively, De = 0, 5.2×10−8, and 1.3×10−6 μJ/μm2. For the case with no intercellular interaction in Figure 2, the cells maintain their initial position and shape in the no-flow condition, until the shear flow is imposed. Under the influence of the flow field, these RBCs move with the local fluid with some shape deformation and rotation. Since there is no interaction between these cells, they disperse and no aggregates are formed. Except for their different translational velocities due to individual locations in the shear flow, the deformation and rotation behaviors of these fours cells are almost identical to each other. During a rotation period, the cell shape alternates between the biconcave and a shortened and folded (resembling a reversed S) shapes. Such behaviors are very similar to those observed in previous experimental and numerical studies (Fischer and Schmid–Schonbein, 1977; Pozrikidis, 2003; Eggleton and Popel, 1998).
When a weak intercellular interaction is introduced (De = 5.2×10−8 μJ/μm2), the cells adhere to each other and a rouleau-like aggregate with smaller gaps is formed under no-flow condition (Figure 3a). With the application of the same shear flow with γ = 20s−1, the aggregates rotate about 270° and deform largely in the flow (Figures 3a-c). However, it appears that this interaction is not large enough to hold the cells together, and at γt = 43.75 (Figures 3d), the two end cells are separated, leaving a two-cell aggregate rotating in the domain center for four rounds (Figures 3e). At last, the two center cells are also separated by the viscous forces and the four cells become completely detached (Figures 3f). To overcome such shearing viscous force, a stronger intercellular interaction (De = 1.3 × 10−6 μJ/μm2) was introduced, and the corresponding rouleau behaviors are shown in Figure 4. The rouleau simply rotates in the shear flow with some configuration change and cannot be broken into individual cells at γ = 20 s−1 due to this strong interaction.
To examine the hydrodynamic effect on cell behavior, we increase the shear rate to γ = 50 s−1 for the above cases with weak and strong interaction strengths De = 5.2 × 10−8 (Figure 5) and 1.3 × 10−6 μJ/μm2 (Figure 6). At this higher shear rate, the rouleau with a weak interaction De = 5.2 × 10−8 rotates only about 90° (Figures 5a-c) and is quickly broken into two two-cell aggregates (Figure 5d), which are not stable either and are finally split into four individual cells (Figure 5f). Note that compared to the case with same weak interaction at γ = 20 s−1 in Figure 3, here at γ = 50 s−1 the rouleau is more easily separated into individual cells in a shorter time. Similar pictures can be seen for the case with strong interaction (Figure 6). The compact rouleau is separated into two two-cell aggregates after rotating about 270° at γt = 40.625 (Figures 6a-d). After that, the two-cell aggregates rotate and cannot be separated further even at this high shear rate γ = 50 s−1 with the strong intercellular interaction between cells.
The above simulations show that the dynamic RBC behaviors in shear flows are closely related to the applied shear rate and the intercellular interaction. A larger shear rate produces larger hydrodynamic forces and the RBC rouleaux are less stable and easier split into smaller aggregates or individual cells. On the other hand, a stronger interaction can generate more compact rouleaux and tends to hold the cells together in shear flows. As a result, rouleaux with a stronger interaction can survive a larger shear rate without being broken. These findings are in qualitative agreement with existing experimental observations and therefore demonstrate that our model could be useful in further microscopic hemodynamic and hemorheologic studies. As found by Bagchi et al. (2005), the membrane viscosity appears have no major influences on the cell behaviors in the situations studied here. Effects of other RBC properties, such as deformability and suspending viscosity, will be examined in a future study. We have also observed various cell shapes in these simulations, some of which are similar to previous three-dimensional simulations (Eggleton and Popel, 1998; Pozrikidis, 2001). Goldsmith and Marlow (1979) had also presented RBC deformations in concentrated ghost cell suspensions; however, a direct comparison is difficult due to the different conditions. Liu et al. (2004) had studied a very similar system of four RBCs in shear flow. Our results, in both the flow field and cell movement, are similar to those of Liu et al. (2004) when there is no intercellular interaction involving aggregation forces. However, the aggregation process was not investigated and the cell deformation was less evident in Liu et al. (2004). The small cell deformation could be due to the finite element modeling of the RBC membrane as a thin shell, in which the experimentally measured RBC membrane properties cannot be utilized directly.
We have presented a numerical algorithm to investigate microscopic hemodynamic and hemorheological behaviors of discrete RBCs in shear flow. The immersed boundary treatment has been combined with the lattice Boltzmann method to study deformable fluid capsules in flows. This algorithm has been applied to study the red cell behaviors in shear flow. Crucial factors, including the RBC membrane mechanics, plasma/cytoplasm viscosity difference, intercellular interaction, and the hydrodynamic viscous forces, have been considered in this numerical model. Simulations of four RBCs in shear flows have been performed, and the effects of the intercellular interaction and hydrodynamic forces on the cell aggregating and disaggregating processes have been investigated. Without the intercellular interaction, the cells remain separated and no aggregates are formed. By introducing the interaction, cells are attracted to each other and rouleau-like aggregates can be generated. The aggregate configuration is closely related to the interaction strength. Such an aggregate could rotate as a unit with certain deformation, or be separated under a shear flow, depending on the intercellular interaction and shear strength. These findings are qualitatively consistent with experimental observations, and therefore demonstrate that the proposed algorithm could be useful in microscopic blood flow studies. Extending this algorithm to three-dimensional situations is straightforward, as well as parallelizing it for better computational efficiency. Different intercellular interaction models can also be readily incorporated in the numerical scheme.
It should be mentioned that the Morse potential adopted in this study is only a qualitative description of the still unclear microscopic intercellular interactions. In principle, parameters involved in the Morse potential Eq.(5) could be obtained through comparing simulation results to observable macroscopic behaviors, such as rouleaux configurations, aggregation and dissociation processes, and apparent viscosities at different shear rates. For the sake of simulation efficiency, the zero force distance r0 in this study was chosen to be relatively large compared to experimental values. Similar assumptions and simplifications were made in previous studies [for example, see Fig. 6 in Liu et al. (2004)]. In addition, since the fluid trapped in a gap between two close RBC surfaces is immobilized, from a viewpoint of hydrodynamics, a difference in the gap sizes might not have significant impact on the global behaviors. Nevertheless, with the present knowledge of the RBC intercellular interactions and with presently used computational technology, the Morse potential could be a reasonable choice to model the RBC aggregation.
This work was supported by NIH grant HL52684. J. Zhang acknowledges the financial support from the Natural Science and Engineering Research Council of Canada (NSERC) through a postdoctoral fellowship at the Johns Hopkins University. P. C. Johnson is also a Senior Scientist at La Jolla Bioengineering Institute.
Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Simulation Movies Available as Supplemental Materials.