Abstract
Erythrocyte aggregation and dissociation play an important role in the determination of hemodynamical properties of blood flow in microcirculation. This paper intends to investigate the adhesion and dissociation kinetics of erythrocytes through computational modeling. The technique of immersed boundary-fictitious domain method has been applied to the study of erythrocyte aggregates traversing modeled stenotic microchannels. The effects of stenosis geometry, cell membrane stiffness, and intercellular interaction strength on aggregate hemodynamics including transit velocity are studied. It is found that the width of the stenosis throat and shape of stenosis have a significant influence on the dissociation of the aggregates. Moreover, horizontally orientated erythrocyte aggregates are observed to dissociate much easier than their vertical counterparts under the same simulation conditions. Results from this study contribute to the fundamental understanding and knowledge on the biophysical characteristics of erythrocyte aggregates in microscopic blood flow, which will provide pathological insights into some human diseases, such as malaria.
1. Introduction
Erythrocyte (red blood cell or RBC) aggregation and dissociation are a common and complex biophysical process in vivo [1, 2]. The tendency of erythrocytes to aggregate increases significantly in a number of human diseases, such as malaria, cancer, and some hemorheological disorders [3–7] due to the alternations in cell membrane mechanical and adhesive properties. For example, when people are infected by malaria in which erythrocytes are parasitized by Plasmodium falciparum, erythrocytes tend to form aggregates in the shape of stack of coin because the surface of the cell becomes stiff and adhesive after infection. In fact, The malaria-infected erythrocytes could be more than ten times rigid [8–10] and adhesive [11, 12] in comparison with healthy ones. These rouleaux structures break into smaller aggregates or even become individually dispersed as the shear stress in blood increases. This aggregation phenomenon of erythrocytes can lead to drastically altered blood flow dynamics and is also responsible for the non-Newtonian behavior of blood, which is of great scientifical and clinical interest [13–15].
Arterioles are the primary location of blood flow in microcirculations where erythrocytes are the essential constituents. Due to their relatively small size in diameter (20–50 μm), the effect of vessel wall thickening and luminal narrowing on the blood flow in arterioles may be significant. Evidence shows that retinal arteriolar stenosis is closely related to coronary heart disease in woman [16]. Recently, significant advances have been made on experimental design for microfluidic systems. For example, micropipette aspiration, optical tweezers, and microchannel are recently developed techniques that have been applied to the study of the behavior of vesicles and erythrocytes. Due to the development in fabrication of novel microfluidic device for separating and manipulating erythrocytes, there has been a growing interest in investigation of blood flow dynamics at microscale [17–20]. Experimental in vitro studies of individual cancer cell [21] and erythrocyte [19, 22] have been done in microfluidic devices, to name a few. In particular, single erythrocyte dynamics has been studied in flow along sinusoidal patterned microchannels [22], and both orientational angle and shape oscillations have been observed.
However, blood flow exhibits a complex and rich behavior due to the cell-cell, cell-structure, and cell-fluid interactions in microcirculation, and the current experimental techniques have several limitations for this topic, especially when the aggregation of erythrocytes is involved. For example, it is difficult to obtain three-dimensional flow information using these methods. Moreover, microvasculature in human body consists of very complex network of circular channels which is impossible to be fabricated by today’s technology. On the other hand, it is critical to understand the underlying biophysical mechanisms of physiological and pathological phenomena for the design of advanced microfluidic devices. Numerical simulation is an alternative to overcome these problems and has made remarkable success. Fenech et al. [23] managed to simulate erythrocyte aggregation phenomena with a large number of cells by a particle dynamic model. In their model, the erythrocytes were designed as spheres instead of biconcave, deformable objects. Recent numerical studies on the erythrocyte aggregation have taken into account the rheological aspect and dynamic motion of the cells in blood flow [24–27]. It has been shown numerically that aggregability is an important determinant factor of the hemodynamics and rheological behavior of blood in microcirculation [28, 29]. Studies also showed that erythrocyte aggregability was closely related to the cell deformability and cell-cell interactions [30]. However, these numerical simulations have been conducted in regular straight channels or tubes. There have not been many microscopic simulations performed with geometrically irregular boundaries, which is a common situation in human cardiovascular system.
The present study aimed to investigate the dynamics of erythrocytes aggregates in microchannel with stenoses by an immersed boundary-fictitious domain scheme. The goal of this study focused on the structure-induced aggregate dissociation taking place in the irregular-shaped vessel. To achieve this, a previously developed membrane model [30, 31] was used to simulate the dynamics of erythrocytes; that is, the erythrocytes were modeled as membrane particles connected by springs with stretch/compression resistance and bending rigidity. By varying the stretch/compression and bending constants, healthy and malaria-infected erythrocytes were modeled. The immersed boundary method has been coupled with the fictitious domain method to deal with the complex flow behavior in this irregular domain geometry. The effect of geometry of the channel, mechanical property of the cell membrane, intercellular strength, and initial aggregate orientation were also analyzed.
2. Mathematical Formulation and Numerical Method
The blood flow region is a two-dimensional microchannel with two symmetrically formed stenoses, and erythrocyte aggregates move in the stenotic channel under the influence of hydrodynamical force. The fluid (blood plasma) is assumed to be incompressible, Newtonian with constant density and constant viscosity so that the Navier-Stokes equations can be applied. The flow region is expressed as . In order to solve the fluid flow and the fluid-cell interactions in this irregular-shaped domain, the fictitious domain method (FDM) was combined with the immersed boundary method (IBM).
2.1. Immersed Boundary-Fictitious Domain Scheme
The fictitious domain method (FDM) and its applications to fluid flow problems have been discussed elsewhere [32, 33]. To employ the FDM, the flow region is embedded in the smallest rectangular domain, which is denoted by . Then the fluid flow containing erythrocytes is solved in the bigger domain , and the no-flow condition in the solid region is treated as constrains. Through this method, the irregular-shaped domain is extended to regular shaped so that a simple structured computation grid instead of unstructured mesh can be used, which substantially reduces computational complexity of the algorithm. The method can be described by the following extended Navier-Stokes equations: where and are the fluid velocity and pressure anywhere in the flow; is the fluid density; and is the fluid viscosity. The body force term is introduced to account for the force acting on the fluid/cell interface. The boundary conditions are such that on no-slip condition is applied and at the inlet and outlet of the channel, periodic flow condition is enforced. A detailed description of the solution method of (1) and (2) can be found in [32, 33].
In this study, the fluid-cell interaction dealt with the immersed boundary method (IBM) developed by Peskin [34]. Based on this method, the boundary of the deformable object is easily calculated by the following scheme: first, the force located at the immersed boundary node affects the nearby fluid mesh nodes through a discrete function: where is the uniform finite element mesh size and with the 1D discrete -functions being The force generated by the deformation of the membrane is then substituted into the external force term of (1) only for membrane particles; next, the movement of the immersed boundary node is affected by all the nearby fluid mesh nodes through the same discrete -function: Finally, after each time step , the position of the immersed boundary node is updated by
2.2. Erythrocyte Model
In this paper, we adopted the spring model introduced in [31] and modeled individual erythrocyte as cytoplasm enclosed by a membrane, which is represented by the two-dimensional network consisting of finite number of membrane particles connected by springs. The springs can change length with respect to its reference length and rotate with respect to neighboring springs. By doing so, elastic energy is stored in the spring, and the shape of the enclosed area changes. The total energy of the erythrocyte membrane stored, , is the sum of the total elastic energy for stretch/compression , the total elastic energy for bending , and the penalty function for area : In the above equations, is the length of the spring; is the angle between two neighboring springs; is the total number of the spring elements; and are spring constants for changes in length and bending angle, respectively; and is the constant for area reduction. The equilibrium enclosed area is denoted by , while is the time-dependent area of the erythrocyte.
Based on the principle of virtual work, the elastic spring force acting on the membrane particle is then with being the position of the th membrane particle. The force generated by the deformation of the membrane is substituted into the external force term of (1).
Then the shape of the erythrocyte is obtained using the elastic spring model based on minimum energy principle. Initially, the cell is assumed to be a circle with radius formed by the springs. When the area of the circle is reduced, each membrane particle moves according to the following equation of motion: Here, denotes the time derivative and and represent the mass and the viscosity of the membrane. The position of the th membrane particle is solved by a discrete analogue of (10) via a second-order finite difference method. The total elastic energy stored in the membrane decreases as the time elapses. The final shape of the erythrocyte is obtained when the total elastic energy is minimized.
It is important to note that the bending constant is closely related to the rigidity of the membrane. The higher the value of the bending constants, the more rigid the cell membrane. Based on this property, the malaria-infected cell and the normal healthy cell can be modeled by changing the values of the bending constants in the spring model.
2.3. Cell-Cell Interactions
The intercellular energy is modeled by the Morse potential function [26] as follows: where is the distance between the cell surface; is the equilibrium distance at zero force, which is set to be 0.049 m; is the intercellular interaction strength; and is the scaling factor. Thus, the intercellular force has the form
It is important to note that this simple model does not illustrate the underlying intercellular interaction mechanism which still remains unclear nowadays.
3. Numerical Results and Discussion
The simulation parameters are listed in Table 1. A constant pressure gradient was imposed at the inlet and the outlet of the channel so that a fluid flow is established from left to right. We have chosen the pressure gradient based on the velocity profile for tube flow , although our simulations have been done in two-dimensional channels. The maximum flow velocity for the tube flow was 12.5 cm/s unless otherwise stated. In addition, periodic boundary conditions were assumed at the left and right boundaries of the domain.
The shape of erythrocytes has been chosen as in the simulations, which represents the typical biconcave shape of the erythrocyte cell observed in human blood. In particular, the length of the cell with is about 7.6 μm and the thickness is about 2.1 μm.
3.1. Formation of Aggregates
The erythrocyte model was verified by forming four-cell aggregates under no-flow condition prior to introducing the cells in the vessel. The simulations have been done in a 15 μm × 15 μm rectangular domain. The distance between the centers of the two neighboring cells is 3 μm so that the intercellular force is only attractive initially. As time elapses, erythrocytes attract and move toward each other until the balance of the attractive force and the repulsive force is achieved and the equilibrium configuration is obtained. The results are demonstrated in Figure 1 for the equilibrium configuration. It can be seen that, at equilibrium position, the cells deform themselves to form a loose or compact aggregate that depends on the strength of the intercellular force and the deformability of the cell membrane. The rouleaux formed under strong intercellular interaction for the more deformable cells (Figure 1(b)) are more compact with the biconcave shape of the two middle cells being totally lost. The obtained results have been compared with previous simulation results [27, 30] and good agreement has been found.
(a)
(b)
(c)
(d)
(e)
(f)
3.2. Dissociation of Aggregates
Erythrocyte motions and dissociation in stenotic microvessels were analyzed in this section. Two types of channels were used in the simulations, namely, the nonstenotic channels and the stenotic channels. For both types, the total length of the channel is fixed at 100 μm to ensure the fully development of the flow. For the microchannels with two trapezoidal-shaped stenoses, the stenoses formed at the location so that the length of the stenosis throat is 12 μm for all cases. The aggregates formed in Figure 1 are placed near the left outlet initially and then move with the fluid when the flow starts. The effect of several important factors, such as severity of the stenosis, membrane stiffness, intercellular strength, initial aggregate orientation, and shape of stenosis, is studied in this section. The aggregate orientation is defined based on the orientation of the erythrocyte cells assembly in the aggregate. The aggregate is defined as horizontally orientated if the cells are horizontal, and the aggregate is called vertical if the cells are vertically located.
3.2.1. Nonstenotic and Stenotic Channels
The first case considered was that of a four-cell aggregate horizontally located in the channel with the cell center coaxial with the channel. Snapshots of the numerical simulations for the aggregate traversing a microchannel without or with stenosis are presented in Figures 2, 3, and 4. In these figures, four time instances are presented with both velocity vector field and axial velocity magnitude shown in pairs.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
The results in Figure 2 show that, in a nonstenotic channel, the aggregate flows with the blood plasma. The cells deform under the hydrodynamic force with the velocity of the cells being slower than the velocity of the plasma. The two outside cells in the aggregate are being peeled off along the flow process while the two middle ones still form an aggregate. At some stage, the two outside cells are dissociated almost completely. The flow field around the cells is not parabolic any more but of some distortion. The distortion becomes larger when the aggregate becomes more dispersed.
To investigate the effect of stenosis severity to the dissociation of the aggregates in the simulation channel, we have also performed the same simulations with width of stenosis throat and 8 μm. It can been seen from Figure 3 that when the channel is constricted the flow field altered from parabolic profile. At the throat of the channel, the velocity of the blood plasma is much higher. When the aggregate enters the constriction, the flow field is disturbed. Again the velocity of the cells is slower than that of pure plasma. It can also be observed that the two outside cells have been peeled off earlier than in the straight channel and the two middle cells are more deformed upon exiting the first constriction. After the second constriction, the two peeled off cells are lagging further behind. At the same time, the two middle cells are driven further apart by the flow. The results for the channel with 8 μm stenosis are shown in Figure 4. The flow is more blocked and the maximum velocity at the stenosis throat decreases sharply. Moreover, the aggregate is totally dissociated to individual cells after traversing the two constrictions. Similarly, disturbance of velocity magnitude contour has been observed around erythrocyte cells.
3.2.2. Membrane Stiffness
The effect of the stiffness of the cell membrane was investigated by changing the membrane constants. The snapshots of the aggregate formed in Figure 1(e) traversing through a 12 μm stenosis are illustrated in Figure 5. This figure indicates that when the cells are more rigid the cells experience less deformation. The increasing membrane stiffness prevents the cells from large deformation and formation of sharp edges. Therefore, normal biconcave shape of erythrocytes is more conserved as the cells become less deformable. Similar to Figure 3, the two outside cells are peeled off by the viscous force and the two middle cells remain a smaller aggregate by the intercellular force.
(a)
(b)
(c)
(d)
3.2.3. Intercellular Strength
In this simulation, the effect of intercellular strength was studied. The simulation pressure gradient was decreased by a factor of five compared to the previous simulations. The intercellular strength in Figure 7 was hundred times greater than that in Figure 6. It can be noticed from the slight discrepancies between these two figures that the stronger the intercellular strength, the more difficult to dissociate the aggregates.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
Another noticeable point is that for all the previous simulations the two individual cells move away from the aggregate after the first stenosis because of the divergent characteristic of the streamline. It is also observed in this section that when the cells approach and traverse the stenosis, the velocity magnitude contours are completely disturbed. The standard symmetric velocity profile is recovered some time after the cells passing through the stenosis.
3.2.4. Initial Aggregate Orientation
In this section, we conducted the simulations of a four-cell aggregate traversing a stenotic channel while the cells were placed vertically with the centers of the erythrocytes at the axis of the channel. The results are presented in Figures 8, 9, 10, and 11 which correspond to the configurations in Figures 1(c)–1(f), respectively, for different membrane stiffness and intercellular strength. Comparing these figure with the previous ones, different behaviors have been observed for the entire period of simulations. In summary, when the cells are vertically located initially, the aggregates are not completely dissociated by the fluid viscous force after passing through the two constrictions in a 8 μm stenotic channel. The cells deform themselves into nonsymmetric slipper-like shapes (Figure 9) or symmetric parachute-like shapes (Figures 10 and 11) at the throat of the channel. Due to the higher fluid velocity and deformation, the aggregates in fact become more disperse in the stenotic part while forming a compact one upon exiting the stenosis. The cells exhibit a periodic shape transition from stretched to deformed shapes as they flow with the blood plasma. Another noticeable point in these simulations is that although complete breakage did not happen for the aggregates, they become less compact compared to the initial configurations, especially for the case of strong intercellular strength and less rigid cells.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
In addition, velocity plots show that the velocity magnitude of the fluid flow is smaller at the erythrocyte vicinity. It can be seen from the velocity vector plots that the velocity profile changes from parabolic to a blunt form. This result is also observed in previous numerical studies conducted in straight channels [35].
3.2.5. Shape of Stenosis
The effect of the shape of stenosis on the dissociation of the erythrocyte aggregates has been studied. We considered stenosis with radial symmetrically sinusoidal shape instead of trapezoidal shape in this section. The length of the base and the maximum of the height of the sinusoidal-shaped stenosis have been chosen the same as in Figure 4 which are 24 m and 6 m, respectively. Compared with Figures 4 and 8, a similar dissociation process is observed in sinusoidal-shaped stenosis (Figures 12 and 13). However, the profile of velocity is less disturbed and the velocity of the flow in sinusoidal shaped stenosis is much higher than in trapezoidal-shaped one.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
3.3. Transit Velocity
It is important to note that the transit time of erythrocytes in the microvessel depends on both the cell deformability and the friction encountered by the cell during its entry into the microchannel. Hence, we define the transit velocity of a single cell being the average velocity of the cell as it passes through the microchannel and travels the same distance. Similarly, the transit velocity of an aggregate is defined as the averaged transit velocity of all cells in the aggregate. Transit velocities of erythrocyte aggregates in the modeled stenotic channel were calculated and analyzed in this section for various cases, and the results are presented in Figure 14.
First of all, the effect of severity of the stenosis on the transit velocity was considered by varying the width of the throat of the vessel for the erythrocyte aggregates. Figure 14 shows a scatter plot of transit velocity against stenosis severity, which is defined as the percentage of the channel that is blocked by the stenosis; for example, the stenosis severity of a 12 μm stenosis is 40% by the definition.
From the graph, it can be seen that the transit velocity decreases relatively slowly when the severity increases from 0% to 40% and sharply when the severity increases from 40% to 60% for the same simulation parameters chosen. The reason for this is that since the pressure difference at the inlet and the outlet remains the same, when the vessel is severely stenotic, the flow is more blocked, which slows down the transit velocity of the aggregates.
The effect of intercellular strength is demonstrated by comparing Figure 14(a) with Figure 14(b) for horizontal aggregates. The small discrepancy between the two plots indicates that the effect of intercellular strength on the transit velocity is less significant. The plots in Figures 14(c) and 14(d) show that when the aggregate is vertically orientated the transit velocity again follows the same trend as the horizontal case; however the transit velocities have a significant increase for all the simulations. The aggregates consist of more rigid cells, and the transit velocity decreases slightly compared to the more deformable one for all the stenosis severity. It has been observed experimentally [9] that for individual erythrocyte the transit velocity decreases with the increase of the membrane stiffness. Our results show that erythrocyte aggregates have the similar behavior as nonaggregated cells.
Figures 14(e) and 14(f) reveal the transit velocity of aggregates in sinusoidal stenotic channels. Because the length of the stenosis throat is significantly shorter, the flow is less blocked and the transit velocities are much higher for the sinusoidal-shaped stenotic channel than for the trapezoidal-shaped one.
4. Conclusions
We have applied the immersed boundary-fictitious domain method to the erythrocyte aggregates in microscale blood flow through modeled stenotic arterioles. A spring model was adopted to describe the deformability of the erythrocyte membrane. The intercellular interaction between the cells was characterized by the Morse potential function. We focused on the effect of the stenotic structure on the dissociation of the aggregates. We also investigated the dependence of the dissociation on the cell membrane stiffness, intercellular strength, and initial orientation of the cell aggregates. The transit velocity of the erythrocyte aggregates traversing the stenotic channel was also studied in this paper.
The results demonstrate that the velocity of the blood flow in stenotic channel is decreased by the stenosis. The velocity at the stenosis throat is significantly higher than that at the nonstenotic part of the channel. At vicinity of the aggregates, the parabolic profile of the flow is disturbed and flattened by the existence of the cells. At stenosis throat, the cells undergo larger deformation than at the nonstenotic part. The aggregates may be dissociated completely by increasing the severity of the stenosis; that is, increasing shear rate at the narrow part of the vessel facilitates the dissociation of aggregates. Moreover, the initial orientation prior to entering the stenosis has great influence on the dissociation of the erythrocyte aggregates. Horizontally orientated aggregates have a tendency to be dissociated easily by the fluid viscous force and the cell-structure interactive force than vertically orientated aggregates. We also simulated with erythrocytes of various stiffness and intercellular interactive strength to mimic the healthy and the abnormal cells, and the obtained results agree well with previous simulation findings. In addition, we have explored the dependence of the transit velocity of the aggregates in the stenotic channel on simulation parameters, and the results are presented in this paper.
The study shows that the structure of the vessel has a significant effect on the rheological behavior of erythrocyte aggregates. It is of medical interest since the erythrocyte aggregation critically affects blood hemodynamics in microvessels and is therefore a nonnegligible factor for capillary blockage and other serious malfunctions of microcirculation. Although the simulations are performed in two dimensions, some similar behavior is expected for the more realistic three-dimensional situations. However, it should be noted that there is crucial difference between two-dimensional and three-dimensional cases. The simulations on the three-dimensional microscale blood flow are undertaken and will be presented in the future studies. Furthermore, the simulation technique will be applied to address many important questions, such as erythrocyte aggregates in a bifurcated channel or in compliant blood vessels.
Acknowledgment
The authors acknowledge the support of the National Natural Science Foundation of China (11074109).