Characterization of Blood Flow in Capillaries by Numerical Simulation

This paper presents a numerical investigation of the axisymmetric, pressure driven motion of single file erythrocyte (i.e., red blood cell) suspensions flowing in capillaries of diameter 8-11 μm. Our study successfully recreates several important in vivo hemodynamic and hemorheological properties of microscopic blood flow, such as parachute shape of the cells, blunt velocity profile, and the Fahraeus effect, and they have been shown to have strong dependence on cell deformability, hematocrit and vessel size.


Introduction
Investigation of blood flow dynamics and erythrocyte (i.e., red blood cell) rheology in capillaries is of great importance not only because they are the major site of oxygen and nutrient exchange, but also because the proper microcirculatory function is primarily determined by the rheological behavior of red blood cells (RBCs) in these vessels.In capillaries, blood can be viewed as a suspension of RBCs in plasma due to the high volume fraction (about 99%) of RBCs in the blood cells.It has been shown that capillary flow dynamics is significantly affected by arrangement, orientation and deformability of RBCs in plasma suspension [1].
Measuring 5-11 µm in diameter, capillaries are the smallest blood vessels in the human body.The average size of normal red cells is approximately 7.7 µm in diameter and varies in thickness from ~2.8 µm at the rim to ~1.4 µm at the center.They assume a biconcave-discshape in the absence of flow.At the level of capillaries, RBCs tend to travel in single file, separated by gaps of plasma [2].An important property of the RBC is that the cell membrane is highly deformable.In capillaries, they deform and assume an axisymmetric parachute-like shape (coaxial with the flow axis) in order to reduce the flow resistance, leaving a cell-free plasma layer between the RBC and the vessel wall [3,4].It is also observed that under pathological conditions, such as malaria, infected RBCs become more rigid than healthy ones.The deformability of RBC is a critical determinant of blood flow in capillaries, and is the combined result of several mechanical and physiological properties.
Another major determinant of blood property is the hematocrit (H ct ), which is defined as the volume fraction of the RBCs to the total blood volume and varies from 29-41% for a one-year-old child to 38-46% for adults.When the term hematocrit is applied to microvessels, it can include two different estimates of red cell distribution, the tube hematocrit (H T ), defined as the instantaneous volume fraction of RBCs in a specified volume of capillary and the discharge hematocrit (H D ), defined as the proportion of RBCs flowing out the capillary.The Fahraeus effect [5] is characterized as the reduction of H T below H D which occurs in microvessels.The effect is due to the fact that the mean velocity of the RBCs (U c ) is higher than the mean bulk flow velocity (U m ) [6].The ratio of these velocities decreases with increasing hematocrit.
Indeed, the blood flow dynamics in capillaries has not been extensively explored.Within this context, the application of mathematical and numerical models has been shown to provide useful information on dynamic characteristics of blood flow under complex flow conditions that is difficult to obtain through in vivo or in vitro experiments.The nature of blood flow in capillaries requires a different model, if not more complex, from those for blood flow in big and medium sized vessels in which the fluid is well approximated by a continuum medium.At the same time, in order to obtain a better understanding of blood flow in capillaries, two other aspects must be considered in the numerical simulation as well.One aspect is that multiple cells must be considered to account for cell-cell hydrodynamic interaction because under high H ct conditions, the interaction between RBCs becomes significant.The other is that deformation of the cell must also be considered in the model, as it is a major determinant of many physiologically significant phenomena, such as formation of the cell-free layer, and the Fahraeus effect.Tsubota et al. [7] carried out a numerical study on effects of H ct on blood flow properties in a two-dimensional channel using particle method and suggested that shape of the red cells and the flow in capillaries is significantly affected by the hematocrit level.However, the study did not include the effect of the cell deformability and vessel size.Pozrikidis [8] investigated the axisymmetric motion of RBCs in cylindrical capillaries using a boundary-integral method.The effect of capillary radius and cell spacing on the discharge hematocrit and apparent viscosity was studied.Zhang et al. [9] developed an immersed boundary lattice Boltzmann approach to investigate the blood flow in microvessels.However, the deformability of the red cell was only considered to a limited extent in [8] and [9].In addition, there was a lack of quantitative analysis of the RBC and the flow behaviors in the study of [9].
The objective of this study is to predict the blood flow properties in capillaries using numerical methods.We present two-dimensional computational simulations of blood flow in vessels of diameter 8-11 µm at H T of 10-41%, taking into consideration the particulate nature of blood and cell deformation.The simulation is based on the numerical solution of the Navier-Stokes equations, and the red blood cell membrane is modeled as membrane particles connecting by springs.Also presented in this paper are parametric simulation studies on the effect of hematocrit, deformability, and size of the vessels on the shape change of the cells, plug-flow velocity profile, cell-free layer thickness, and the Fahraeus effect.A qualitative/quantitative comparison between the simulation results and experimental data is also presented in this paper.

Numerical Simulations and Discussions
In capillary vessels, blood is considered as a multiphase fluid because the size of the RBC is comparable to the size of the vessel.It is anticipated that capillary flow dynamics should be similar to that of a suspension of deformable particles.Based on the fact that blood plasma, the liquid component of blood in which the blood cells are suspended, is composed of mostly water (90% by volume), the plasma flow in capillaries is assumed to be governed by the Navier-Stokes equations for the incompressible, Newtonian fluid where u and p are the fluid velocity and pressure anywhere in the flow; ρ is the fluid density; µ is the fluid viscosity; f accounts for the external body force.Both ρ and µ are assumed to be constant for the entire fluid.
The deformable shape of the RBCs is modeled by the elastic spring model described in Appendix A. The model has been validated and applied to a number of blood flow studies [10,11].It has been shown to have the ability to capture the deformation of the RBCs under various flow conditions [10,11], which is fundamental to the successful simulation of the flow behavior.It takes into consideration the structure of the RBC membrane skeleton.The membrane of the cell has strong resistance to changes in area/volume and shear deformation which is consistent with other people's findings [12].It also shows the ability of recovering to the initial biconcave shape after the removal of the external flow field [11].
The motion of the deformable cells and fluid domain are coupled together by the immersed boundary method [13] described in Appendix B. This method is an innovative approach to deal with the problem of modeling fluid flows interacting with a flexible, elastic boundary.The applications of immersed boundary method to the simulation of deformable cells can be found in, e.g., [10,11,[14][15][16].
For the simulation, the RBCs are suspended in blood plasma which is assumed to be incompressible, Newtonian and has a density ρ = 1.00 g/cm 3 and a dynamic viscosity µ = 0.012 g/(cm•s).The viscosity ratio which describes the viscosity contrast of the fluid inside and outside the RBC membrane is fixed at 1.0.The fluid domain is a two dimensional channel of horizontal length 20-25 µm.The width of the channel is varied from 8 µm to 11 µm.A single file of RBCs coaxial with the flow axis is placed vertically with uniform center to center distance, which can be adjusted to correspond varies of hematocrit.
The parameters in the simulation of the shape change of the RBCs are set as follows: the membrane mass m = 2.0 × 10 -4 g and the membrane viscosity γ = 8.8 × 10 -7 Ns/m.The penalty coefficient k s = k b × 10 4 .The spring constants are set as k l = k b .The bending constant is closely related to the rigidity of the membrane.A higher k b value results a less deformable cell.The biconcave shape obtained for s* = 0.42 resembles the normal physiological shape of the RBC very well and is used for the simulations in the current study.
For all computations, the grid resolution for the computational domain is 80 grid points per unit length with the unit length equal to 10 µm.To obtain a Poiseuille flow, a constant pressure gradient is prescribed as a body force.Periodic boundary conditions are imposed at the left and right boundary of the domain.The evolution of the periodic file of cells from a specified initial configuration with uniform cell spacing is computed.Simulations from the unstressed biconcave shape (Figure 1(a)) are performed until steadily translating deformed shapes (e.g., Figure 1(b)) are obtained in the flow.

RBC Deformation
The deformation of RBCs in plasma flow with varying hematocrit, cell deformability and vessel diameter is studied.Simulations are conducted in the channel described above.As the flow starts, the RBCs change quickly from its static biconcave shape to parachute shape.At the mean time, RBCs move down stream with an increasing velocity.The change of shape stops when the hydrodynamic force is balanced by the elastic force of the cell membrane and the cell velocity become constant.The equilibrium shape strongly depends on the hematocrit (Figure 2(a) and 2(b)) and the deformability of the cell membrane (Figure 2(c)).It is also depends on the width of the flow channel (Figure 2(d)).These results are in qualitative agreement with experimental [17] and simulation results [7,8].The vector fields of the flow velocity are also shown in Figure 2.
We also wish to quantitatively study how the shape of the red cells is affected by the factors such as hematocrit, vessel diameter and the deformability of the cell.To this purpose, equilibrium length L of the cell and the deformation index DI, defined as DI = w/L, where L and w are shown in Figure 1(b), are investigated in this paper.Firstly, equilibrium length L is plotted versus hematocrit for the channel width D = 8, 9, 10, and 11 µm in Figure 3(a).The length L increases almost linearly with the increase of H T for all the vessels in the range of H T = 10-41%.Figure 3(b) shows deformation index DI decreases almost linearly with the increase of H T over the same range for these vessels.In Figure 4, the cell-free layer thickness H is shown to decrease with the increase of H T for the 11 µm-capillary and the results are in good agreement with the experimental results provided by Albrecht et al. [18].We also wish to investigate the effect of the deformability of the membrane on the size of

Blunt Velocity Profile
Comparing to the parabolic profile of the Poiseuille flow for the pure plasma, the velocity profile is flat topped in the center region for the blood flow in capillaries containing RBCs under the same pressure gradient.The effect of hematocrit, deformability of the cell, and vessel size on the velocity profile has been investigated in this paper.
The effect of increasing hematocrit on the velocity profile is show in Figure 6(a) for a 10 µm-capillary.The profile for H T = 10% appears to be parabolic, but with significantly reduced centerline velocity.The flow shows a more and more blunt profile with increasing tube hematocrit H T .The maximum velocity at the centerline decreases rapidly when more RBCs are present.When the tube hematocrit reaches 41%, the maximum velocity at the centerline reduces to about 45% of that of the pure plasma.
As shown in Figure 6(b), we adopted six different bending constant k b values and compare the mean velocity in blood flow of RBC suspension of same tube hematocrit (H T = 21%).It is observed that the flow profile is less blunt with the increase of the RBC deformability.However, the distribution of axial velocity in the capillary seems less sensitive to the change of deformabilty than the hematocrit of the blood.Since a decrease in blunt flow radius causes a reduction in resistance to flow, we can conclude that RBC flexibility plays a major role in reducing the flow resistance of blood in capillaries.
The velocity profile are also studied for an 8 µm-capillary for H T = 10, 21, 31, and 41% as well and the results are shown in Figure 6(c).Similar to the 10 µmcapillary, the flow becomes blunter for higher H T .It is interesting to note that, comparing to the flow in the 10 µm-capillary, the flow in the 8 µm-capillary is less deviated from pure plasma Poiseuille flow for Ht = 10% and more blunt for higher hematocrit values.The reason is that in a narrower capillary with the same hematocrit, the intercellular space between two neighboring cells is larger, which allows the flow to develop more fully than in  a wider capillary.However, this effect is only significant when the hematocrit is low.As the hematocrit becomes higher, the flow is severely blunted by the increasing number of cells in the vessel.

The Fahraeus Effect
In capillary vessels, the red cells speed up relative to the plasma as they squeeze through the capillary.Since they must travel faster than the plasma, there must be fewer of them present to maintain the same proportions of cells and plasma as blood exits the capillary.This is the so-called Fahraeus effect.The reduction of the tube hematocrit H T to the discharge hematocrit H D due to the Fahraeus effect is related to the mean flow velocity U m and the cell velocity U c [6] In the laboratory study [18], the feed hematocrit H F instead of tube hematocrit H T was used as a control parameter.It has been shown by Yen and Fung [19] that H T < H F for the capillaries used in the study and therefore, our numerical results shown in Figure 7 agree well with the experimental results in [18] and simulation results of [8] which illustrated that the hematocrit ratio increased as the tube hematocrit increased.
The hematocrit ratio is also strongly related to the deformability of the red cells.Figure 8(a) shows the ratio increases with the increase of k b rapidly when the cell is more deformable and the increase become slower as the cell become more rigid.On the other hand, the hematocrit ratio is found to decrease linearly with the increase of the deformation index (Figure 8(b)).

Conclusion
In summary, we have simulated the dynamics of the blood flow and RBCs in capillary using a numerical approach.The results show that RBCs in narrow capillaries change to parachute shape in the flow field.The profile of the capillary flow was markedly blunt in comparison to the parabolic one which characterizes the pure plasma flow.The hematocrit ratio reduces from the value of unity (the Fahraeus effect) in these capillaries.Our study reveals that the RBC shape, bluntness of the flow profile, and the reduction of the hematocrit ratio are strongly depend on the tube hematocrit, deformation of the cell, and the size of the vessel.These findings are consistent qualitatively or quantitatively with other people's experimental and numerical results.We also find that the distribution of axial velocity in the capillary is more sen-sitive to the change of hematocrit than the deformability of the cells.The results indicate that the pressure difference in the blood flow has to increase in the capillary vessels in order to sustain the same flux rate of the red blood cells when the hematocrit or the rigidity of the cell increases.The study yields useful insights into understanding the dynamic characteristics of blood flow in capillaries.The potential applications of this study include the analysis of some pathological conditions, such as sickle cell disease (SCD), in which the RBCs become hard, pointed and sticky and shaped like crescents or sickles.This model could also be used to predict microscopic hemodynamic and hemorheological behaviors in more complex microcirculation situations, such as curved and stenotic microvessels, branches and post-capillary expansions.

Appendices Appendix A. RBC Model
In two-dimensional simulations, the biconcave shape of the RBC is approximated by the characteristic cross section in the plane that is parallel to the flow direction if the cell were in shear flow.This paper adopts the elastic spring model that has been proposed by Wada and Kobayashi [20] and used by Tsubota et al. [7] and Wang et al. [10] to obtain the shape of RBC in the absence of external force.Based on this model, the RBC membrane is approximated by a group of membrane particles connecting with the neighboring membrane particles by springs (Figure A1).The total elastic energy of the RBC membrane is defined as E = E l + E b + Γ s .In particular, E l is the energy for stretch/compression which is induced by the change of the length l of the spring with respected to its reference length l 0 and E b is the energy for bending due to the change in angle θ between two neighboring springs The shape change is stimulated by reducing the total area of the circle s 0 through a penalty function In Equations (A1) and (A2), N is the total number of the spring elements; k l and k b are spring constants for changes in length and bending angle, respectively.In Equation (A3), s and s e are the time dependent area and the equilibrium area of the RBC, respectively.Based on the principle of virtual work, the elastic spring force acting on the membrane particle i is with r i being the position of the i-th membrane particle.Initially, the RBC is assumed to be a circle which is discretized into a group of membrane particles so that springs are formed by connecting the neighboring particles.When the area is reduced, each RBC membrane particle moves according to the following equation of motion Here, (˙) denotes the time derivative; m and γ represent the mass and the viscosity of the membrane particle.The position r i of the i-th membrane particle is solved by a discrete analogue of Equation (A5) via a second-order finite difference method.The total elastic energy stored in the membrane decreases as the time elapse.The final shape of the RBC for a given area ratio s* = s e /s 0 is obtained as the total elastic energy is minimized.

Appendix B. Immersed Boundary Method
The boundary of the deformable structure is discretized spatially into a set of boundary nodes.The force located at the immersed boundary node X affects the nearby fluid mesh nodes x through a 2-D discrete δ-function where h is the uniform finite element mesh size and with the 1-D discrete δ-functions being The movement of the immersed boundary node X is also affected by the surrounding fluid and therefore is enforced by summing the velocities at the nearby fluid mesh nodes x weighted by the same discrete δ-function After each time step Δt, the position of the immersed boundary node is updated by Equations ( 1) and ( 2) are numerically solved by a finite element technique combined with the immersed boundary method described in this section.

Figure 1 .
Figure 1.RBC shape obtained using the spring model: (a) Equilibrium shape of RBC under no-flow condition; (b) Parachute shape of RBC in Poiseuille flow.red cell in the blood flow.Figure 5(a) and Figure 5(b) show the dependence of the cell length L and deformation index DI on membrane bending constant k b , respectively.The length L (respectively deformation index DI) increases (respectively decreases) with the increase of k b and the rate of change is less severe for rigid cells.

Figure 5 (
a) and Figure 5(b) show the dependence of the cell length L and deformation index DI on membrane bending constant k b , respectively.The length L (respectively deformation index DI) increases (respectively decreases) with the increase of k b and the rate of change is less severe for rigid cells.

Figure 7 .
Figure 7. Dependence of hematocrit ratio on vessel size for various hematocrit levels.The spring bending constant K b = 5 × 10 -15 Nm.

Figure A1 .
Figure A1.The elastic spring model of the RBC membrane.