Effects of erythrocyte deformability and aggregation on the cell free layer and apparent viscosity of microscopic blood flows

Abstract

Concentrated erythrocyte (i.e., red blood cell) suspensions flowing in microchannels have been simulated with an immersed-boundary lattice Boltzmann algorithm, to examine the cell layer development process and the effects of cell deformability and aggregation on hemodynamic and hemorheological behaviors. The cells are modeled as two-dimensional deformable biconcave capsules and experimentally measured cell properties have been utilized. The aggregation among cells is modeled by a Morse potential. The flow development process demonstrates how red blood cells migrate away from the boundary toward the channel center, while the suspending plasma fluid is displaced to the cell free layer regions left by the migrating cells. Several important characteristics of microscopic blood flows observed experimentally have been well reproduced in our model, including the cell free layer, blunt velocity profile, changes in apparent viscosity, and the Fahraeus effect. We found that the cell free layer thickness increases with both cell deformability and aggregation strength. Due to the opposing effects of the cell free layer lubrication and the high viscosity of cell-concentrated core, the influence of aggregation is complex but the lubrication effect appears to dominate, causing the relative apparent viscosity to decrease with aggregation. It appears therefore that the immersed-boundary lattice Boltzmann numerical model may be useful in providing valuable information on microscopic blood flows in various microcirculation situations.

FIG. 1

Representative snapshots of the flow developing process of normal RBCs with moderate aggregation. The red and blue circles are fluid tracers inside RBCs and in plasma, respectively; and the black circles are a membrane marker on one RBC. The four fluid tracers in plasma (blue circles) have also been labeled for the convenience of discussion. The simulation domain is 19.5×58.5 μm and the time step is 1.25×10⁻⁸ s.

FIG. 2

Trajectories of the four plasma tracers in Figure 1. The nine circles on each trajectory curve indicate the simulation time, which are from t =0 to t = 8 ×10⁶ with an interval of 10⁶ time steps (Δt = 1.25 × 10⁻⁸ s). The tracer positions are presented in lattice units (h = 0.195 μm).

FIG. 3

(a) The total and RBC flow rates, (b) the averaged θ distributions, and (c) the averaged velocity profiles during the developing process. In (b) and (c), the profiles are taken at t = 0.4, 1.0, 3.0, and 7.0 ×10⁶ in the arrow directions. The dashed lines are the initial θ distribution in (b) and the parabolic velocity profile of a pure plasma flow under the same pressure gradient in (c). Q_p and U_p, respectively, are the flow rate and maximum velocity of such a plasma flow. The time t and position y are presented in time steps (Δt = 1.25 ×10⁻⁸ s) and lattice units (h = 0.195 μm), respectively. Other quantities are dimensionless.

FIG. 4

The relative stable states of the blood flows with different deformabilities (left: normal; right: less deformable) and aggregation strengths (from top to bottom: none, moderate, and strong). The flow fields are also displayed by arrows. The simulation domain is 19.5 × 58.5 μm.

FIG. 5

(a) θ distributions and (b) velocity profiles at the relative stable states with different deformabilities and aggregation strengths. The dashed lines are the initial θ distribution in (a) and the parabolic velocity profile of a pure plasma flow under the same pressure gradient in (b). U_p is the maximum velocity of such a plasma flow. The line thickness indicates the cell deformability (thick: less deformable, and thin: normal), and the color is for the aggregation strength (black: none, blue: moderate, and red: strong). The position y is presented in lattice units (h = 0.195 μm), and other quantities are dimensionless.

FIG. A.1

A schematic of the immersed boundary method. The open circles are fluid nodes and the filled circles represent the membrane nodes. The membrane force calculated at node x_m is distributed to the fluid nodes x_f in the 2h × 2h square (dashed lines) through Eq.(A.8); and the position x_m is updated according to x_f velocities through Eq.(A.11).