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 (
![[x with macron]](/corehtml/pmc/pmcents/x0078x0304.gif)
,

) 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 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/
μm
2) 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 ). 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. 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 () and convex shapes at high interaction strengths (). At very high interaction strengths (), 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 - 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/
μm
2. For the case with no intercellular interaction in , 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 (). With the application of the same shear flow with γ = 20s−1, the aggregates rotate about 270° and deform largely in the flow (). However, it appears that this interaction is not large enough to hold the cells together, and at γt = 43.75 (), the two end cells are separated, leaving a two-cell aggregate rotating in the domain center for four rounds (). At last, the two center cells are also separated by the viscous forces and the four cells become completely detached (). 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 . 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 () and 1.3 × 10−6 μJ/μm2 (). At this higher shear rate, the rouleau with a weak interaction De = 5.2 × 10−8 rotates only about 90° () and is quickly broken into two two-cell aggregates (), which are not stable either and are finally split into four individual cells (). Note that compared to the case with same weak interaction at γ = 20 s−1 in , 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 (). The compact rouleau is separated into two two-cell aggregates after rotating about 270° at γt = 40.625 (). 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.