Depletion-model-based numerical simulation of the kinetics of red blood cell aggregation under sinusoidal pulsatile flow

Abstract

BACKGROUND:

Previous numerical modeling studies on red blood cell (RBC) aggregation have elucidated the inverse relationship between shear rate and RBC aggregation under steady flow. However, information on the cyclic variation in RBC aggregation under pulsatile flow remains lacking.

OBJECTIVE:

RBC aggregation was simulated to investigate the complex interrelationships among the parameters of RBC motion under pulsatile flow.

METHODS:

A two-dimensional particle model was used to simulate RBC motion driven by hydrodynamic, aggregation, and elastic forces in a sinusoidal pulsatile flow field. The kinetics of RBCs motion was simulated on the basis of the depletion model.

RESULTS:

The simulation results corresponded with previously obtained experimental results for the formation and destruction of RBC aggregates with a parabolic radial distribution during a pulsatile cycle. In addition, the results demonstrated that the cyclic variation in the mean aggregate size of RBCs increased as velocity amplitude increased from 1 cm/s to 3 cm/s under a mean steady flow of 2 cm/s, as mean steady flow velocity decreased from 6 cm/s to 2 cm/s under a velocity amplitude of 1.5 cm/s, and as stroke rate decreased from 180 beats per minute (bpm) to 60 bpm.

CONCLUSIONS:

The present simulation results verified previous experimental results and improved the current understanding of the complex spatiotemporal changes experienced by RBC aggregates during a sinusoidal pulsatile cycle.

Keywords

Red blood cell aggregation rouleaux mean aggregate size hemodynamics sinusoidal pulsatile cycle Poiseuille flow

1. Introduction

Blood consists mainly of erythrocytes, leukocytes, platelets, and plasma. Erythrocytes (red blood cells, RBCs) are the most prevalent type of blood cell and are the principal determinant of the effective viscosity of blood. In plasma, RBCs form an aggregate or network mediated by macromolecules, such as fibrinogen, in a process known as RBC aggregation. Excessive RBC aggregation plays an important role in the development of cardiovascular pathologies, such as microcirculatory flow disorders and vascular thrombosis [1]. RBC aggregation is an important rheological factor that contributes to elevated viscosity under low shear rates. Blood viscosity and shear rate are inversely related under steady Couette and Poiseuille flows [2–4]. However, studies on RBC aggregation under pulsatile flow remain limited, and the mechanisms underlying RBC aggregation under pulsatile flow have not been fully investigated.

In vitro and in vivo experiments have been conducted to study RBC aggregation in whole blood under pulsatile flow [5–9]. An experimental study on porcine whole blood under pulsatile flow in a Couette flow system showed that cyclic variation in RBC aggregation occurs at 40% hematocrit [10]. However, shear stress drives Couette flow, whereas pressure gradients drive Poiseuille flow and physiological blood flow as well. Pressure gradients may cause variations in the kinetics and aggregation of RBCs in blood flow. The shear rate of arterial blood flow inherently changes radially in the blood vessel and temporally during a cardiac cycle. Some in vitro and in vivo studies have demonstrated that cyclic and radial variations in blood echogenicity are mainly attributed to RBC aggregation [5,11–15]. An in vitro tube experiment with porcine whole blood showed that cyclic and radial variations in echogenicity are dependent on peak blood flow speed at systole, stroke rate, and mean steady flow, suggesting that under pulsatile flow, RBC aggregation is affected by the combined effects of shear rate and flow acceleration [16,17]. Despite all these data, to date, the mechanisms involved in RBC aggregation have not been fully elucidated. In particular, current knowledge on RBC behavior under pulsatile flow remains lacking.

To elucidate the relationship between microscopic RBC interactions and macroscopic rheological behavior, RBC aggregation in Couette flow systems under constant and pulsatile shear rates have been numerically simulated on the basis of the depletion model [18]. Simulation results showed that under constant shear rate, the mean aggregate size (MAS) and shear rate are modally related and that maximum aggregation occurs at approximately 0.1 s⁻¹. The cyclic variation in MAS under pulsatile shear rate is a function of shear rate, shear rate history, and initial aggregate size [18]. However, computer simulations can only account for the relationship between shear rate and RBC aggregation under Couette flow. The driving forces and flow dynamics in Couette flow may be different from those in physiological blood flow. For example, the Couette flow system provides a constant shear rate normal to the flow direction whereas in a tube or blood vessel a parabolic velocity profile ensues with shear rate that varies in the radial direction. Therefore, the present study focused on the depletion-model-based simulation of the microscopic interactions of RBCs during RBC aggregation in a two-dimensional tube under sinusoidal pulsatile flow.

In this mechanical model, the motion of RBCs in a two-dimensional tube under sinusoidal pulsatile flow is driven by interaction and hydrodynamic forces, and the spatial and temporal behavior and aggregation of RBCs were observed. Hemodynamic parameters, such as sinusoidal flow velocity amplitude, mean flow, and stroke rate, were controlled to investigate the effects of hemodynamic parameters on RBC aggregation in a tube. The calculation of MAS at different flow conditions enabled the investigation of the effects of each hemodynamic parameter on RBC aggregation under sinusoidal pulsatile flow.

2. Methods and theory

2.1. Simulation approach

In this simulation, each RBC particle was assumed to be a sphere of an equivalent radius of 4 μm and mass of 2.94 ×10⁻¹³ kg. For the initial condition, 596 RBC particles were randomly distributed without overlapping throughout a rigid two-dimensional tubular space with a diameter of 0.1 mm and a length of 1 mm and a corresponding hematocrit of 30%. Each particle traveled along the streamline. Particles leaving the end of the tube were positioned at the entrance of the tube at the same y-coordinate with the same acceleration, velocity, and force vectors. The time step utilized at each iteration was 0.1 ms, which was smaller than that set by the condition that the maximal displacements for all RBCs at the highest velocity should be smaller than half of the tube length. In our model, RBC interaction forces were adopted from the RBC particle aggregation model proposed by Fenech et al. on the basis of the depletion model [18]. In accordance with Newton’s second law, particle acceleration was derived from a net force composed of interactional and hydrodynamic forces. The position and velocity of a particle in the same streamline in the y-axis can be calculated from uniform motion acceleration. Therefore, the movement of the particles was determined by Newton’s second law with the interaction forces and hydrodynamic force. Single RBCs were marked with red circles. Two or more RBCs were considered aggregates if the distances between their centers were within 2% of the RBC diameter, as proposed by Fenech et al [18], and these RBCs were marked with black circles [18]. The code for obtaining RBC positions was implemented in MATLAB. The aggregation kinetics of RBC particles were represented by MAS as a function of time. In addition, changes in RBC hydrodynamic parameters (e.g., velocity amplitude, mean velocity, stroke rate) were analyzed through aggregation kinetics.

A flowchart describing the simulation is shown in Fig. 1. Initially, RBCs were randomly distributed without aggregation. The initial velocity field under sinusoidal pulsatile flow was given as a function of space and time. A new velocity field was calculated after the time elapsed by dt (0.1 ms) from the initial time. Then, the new position and distance of RBCs were calculated, and RBC aggregates were analyzed in accordance with distance. In case of aggregation, interaction and hydrodynamic forces were calculated. Otherwise, hydrodynamic force without interaction forces was applied. Then, the net force was calculated on the basis of Newton’s second law to obtain the acceleration, velocity, and distance of the particles at another elapsed time. The simulation process loop was repeated for 1 second.

2.2. Numerical modeling of blood and RBC particles

The simulation model considers blood as a collection of particles that interact with each other in plasma. In the depletion model used in this study, the particles are driven by hydrodynamic force and two interaction forces, as shown in Fig. 2. The acceleration of a particle i over time t is given in Eq. (1) in accordance with Newton’s second law. The two interaction forces (f _ij), an aggregation force ( $f_{ij}^{a}$ ), and an elastic force ( $f_{ij}^{e}$ ) are introduced in Eq. (2) to describe RBC aggregation [18]. $\begin{eqnarray}\displaystyle m_{i}\frac{dv_{i}}{dt} & = & \displaystyle \mathop{\sum }_{i\neq j}f_{ij}+f_{i}^{h}\end{eqnarray}$ 1 $\begin{eqnarray}\displaystyle f_{ij} & = & \displaystyle f_{ij}^{a}+f_{ij}^{e},\end{eqnarray}$ 2 where m _i is the mass of a red blood cell i, and v _i is its velocity vector. The vector f _ij represents the force of blood cell j that acts on blood cell i, and $f_{i}^{h}$ is the estimated hydrodynamic force. f _ij is the sum of aggregation force $f_{ij}^{a}$ and elastic force $f_{ij}^{e}$ . Through an iterative process, the velocities and positions of all particles at time t + dt are deduced from the acceleration calculated from total force at each instant t, where dt is 0.1 ms.

2.3. Elastic force

The calculation of elastic force was based on the granular interaction model and that of aggregation force was based on the depletion model as shown in Eq. (3) [18]. $\begin{eqnarray}\displaystyle f_{ij}^{e} & = & \displaystyle k(2R-d)^{\frac{3}{2}}\boldsymbol{n}\quad \text{if }d<2R\nonumber\\ \displaystyle & = & \displaystyle 0\quad \text{otherwise}\end{eqnarray}$ 3 where d represents the distance between the centers of mass of two RBCs, R is the radius of an RBC particle, and n is the normalized vector pointing from one RBC to the other. RBCs are not considered rigid and may be compressed by a factor modulated by the elastic modulus k, which is fixed at 3 ×10⁻⁶ N/m [18]. This force happens when the distance between the centers of two RBCs is smaller than the cell diameter 2R.

2.4. Aggregation force

The depletion model argues that a depletion layer will develop near the RBC surface if it is in contact with a solution containing macromolecules or proteins and the loss of the configuration entropy of surface molecules is not balanced by adsorption energy [19,20]. According to Neu and Meiselman’s theoretical model employed herein, only depletion and electrostatic interactions need to be considered, and steric interactions caused by the overlapping glycocalyx can be ignored [21]. Given the existence of high electrostatic repulsion, cell–cell distances at which the minimal interaction energy occurs are always greater than twice the thickness of the cell’s glycocalyx layer [21]. Therefore, the main behavior of the aggregation force Eq. (4) between two RBCs is simply illustrated as electrostatic repulsion force due to RBC surface charge and osmotic attractive force due to polymer depletion near the RBC surface. The interaction energy between RBCs in polymer solutions is calculated using Morse potential in Eq. (5) [18,21,22]. $\begin{eqnarray}\displaystyle f_{ij}^{a} & = & \displaystyle 2DAB(e^{\text{2}B({\delta}_{0}-{\delta})}-e^{B({\delta}_{0}-{\delta})})\boldsymbol{n}\end{eqnarray}$ 4 $\begin{eqnarray}\displaystyle {\Phi}_{ij} & = & \displaystyle D(e^{2B({\delta}_{0}-{\delta})}-2e^{B({\delta}_{0}-{\delta})}),\end{eqnarray}$ 5 where A is the area of RBC surface in face-to-face opposition with another RBC. D is the coefficient of surface energy. DA (=10⁻²⁵ J) represents energy. B (=10⁷ m⁻¹) is a scaling factor that is inversely related with depletion layer thickness. 𝛿 is the cell membrane distance between two RBCs (𝛿 = d −2 R). 𝛿₀ (=11 nm) is the 𝛿 for which the aggregation force is zero.

2.5. Hydrodynamic force

It is found that hydrodynamic force acts on a spherical particle undergoing arbitrary time-dependent motion in an otherwise quiescent fluid [23]. Equation (6) includes the pseudo-steady Strokes drag, which depends on the past history of particle motion; the Basset memory integral, which includes viscous force and inertial force; and a purely inertial contribution called added mass force. The hydrodynamic force was expressed as the summation of these three forces [23]: $\begin{eqnarray}\displaystyle f_{i}^{h}(t) & = & \displaystyle -6{\pi}{\mu}RU_{p}(t)-6{\pi}{\mu}R\left(\frac{R^{2}}{{\pi}v}\right)^{\frac{1}{2}}\int _{-\infty }^{t}\frac{\dot{U}_{p}(s)ds}{(t-s)^{\frac{1}{2}}}-\frac{2}{3}{\pi}{\rho}R^{3}\dot{U}_{p}(t)\end{eqnarray}$ 6 where U _p = v ₀ − v _i, v ₀ is the velocity of the fluid in the absence of particles, and v _i is particle velocity. 𝜇 and v are the dynamic and kinematic viscosities of blood, respectively. In this work, the dynamic viscosity is 2.7 ×10⁻³ Pa ⋅ s. 𝜌 is the density of blood cell particle (1096 kg/m³). For flow with a low Reynolds number, inertial forces in the Navier–Stokes equations may be neglected [24]. In our simulation, the maximum shear rate at the inner wall of a Newtonian fluid flowing within a tube was 1500 s⁻¹. Then, the Reynolds number Re = a ² S𝜌∕𝜇 is considerably less than unity, where S is the characteristic macroscopic shear rate; a is the characteristic linear dimension of the suspended particle; and 𝜌 and 𝜇 are the density and dynamic viscosity of the fluid, respectively. Thus, the Basset force and the added mass force can be neglected, and the fluid force (Eq. 6) is given only by the instantaneous Stokes drag. Equation (7) represents the velocity field of sinusoidal pulsatile flow. It combines with Poiseuille flow in the y-axis and the sinusoidal velocity profile to travel at a certain speed in the x-axis in addition to steady flow. $\begin{eqnarray}\displaystyle v(t,x,y) & = & \displaystyle \{A_{f}\sin ({\omega}t-k_{f}x+{\varphi})+B_{f}\}\left(1-\frac{4}{h^{2}}y^{2}\right)\end{eqnarray}$ 7 A _f and B _f are the amplitudes of velocity variation and mean velocity, respectively. k _f is the wave number (k _f = 𝜔/c [rad/mm]), and 𝜔 is the angular frequency [rad/s], and its traveling speed is c =1 mm/s. h is the tube diameter (h =0.1 mm), and 𝜑 is the phase (𝜑 =0). The sinusoidal pulsatile flow velocity profile with velocity amplitude (A _f =1.5 cm/s) summed with the mean speed of the blood (B _f =2 cm/s) is shown in Fig. 3(b). This flow is not volume-conserving and is incompatible with the incompressibility condition but is idealized to illustrate time-dependent variations in spatial plots to examine the interaction of three forces under sinusoidal flow during RBC aggregation.

3. Results

The motion of all RBC particles correspond to the flow velocity profile at a stroke rate of 60 bpm, and the snapshots of particle positions and velocity profile at a certain time are shown in Figs 3(a) and (b), respectively. The black particles in the region of interest (ROI) (0.1 mm × 0.1 mm) represent RBC aggregation. MAS for the ROI of Fig. 3(a) was calculated. Figure 4 shows the formation of RBC aggregates, which developed a parabolic shape. At initial time, the RBCs were randomly distributed without overlapping, as shown in Fig. 4(a). At t =0.01 s, the RBCs began to aggregate and form large rouleaux within a very short period, as shown in Fig. 4(b). Then, the rouleaux formed a parabolic shape at t =0.05 s, as shown in Fig. 4(c). The RBCs and rouleaux located near the center of the tube moved faster than those located in the boundaries. The central rouleaux enlarged and the parabolic shape was broken at t =0.1 s, as shown in Fig. 4(d). Three different velocity amplitudes of 0.5, 1, and 1.5 cm/s with a mean velocity of 2 cm/s and their corresponding MAS are plotted in Fig. 5. MAS increased with time and reached the maximal value under the lowest velocity range and decreased again. The maximal MAS appeared near the lowest velocity range between 0.5 and 0.7 s under three different amplitudes. As velocity amplitude increased, the peak MAS increased, and mean MAS over a cycle also slightly increased. These behaviors resulted in drastic cyclic variations. As shown in Fig. 6, three processes of broken parabolic shapes appeared under three different velocity amplitudes of 0.5, 1, and 1.5 cm/s. The cyclic and radial variations in RBC aggregation exhibited a broken parabolic shape at the velocity amplitude of 1.5 cm/s. Cyclic and radial variations weakened as velocity amplitude decreased. The mean flow velocities were changed to 2, 3, and 6 cm/s with the same velocity amplitude of 1.5 cm/s, as shown in Fig. 7(a), and (b), (c), and (d) show the MAS variations with the corresponding mean flow velocity profiles. As mean steady flow velocity increased from 2 cm/s to 6 cm/s, the cyclic variation in MAS weakened and almost disappeared at 6 cm/s. The peak MAS appeared near the lowest velocity range but shifted from 0.9 s to 0.5 s as mean velocity increased from 2 cm/s to 6 cm/s. MAS at three different stroke rates (60, 120, and 180 bpm) is shown as a function of time under the constant velocity amplitude of 1.5 cm/s and mean velocity of 2 cm/s in Fig. 8. Figure 8(a) represents the corresponding flow velocity profiles at the central line of the tube for three stroke rates. Figures 8(b), (c), and (d) show the variation in MAS with three corresponding velocity profiles for three stroke rates. As the stroke rate increased, the peak MAS and cyclic variation decreased.

4. Discussion

4.1. Principal findings

To our knowledge, this is the first study to use the depletion model for the numerical simulation of RBC aggregation under sinusoidal pulsatile flow. The simulation results illustrated cyclic and radial variation in RBC aggregation and the formation and breakage of a parabolic rouleaux. Similar cyclic and radial variation patterns of blood echogenicity were also observed from in vitro tube experiments using porcine whole blood [16,17] and harmonic images of the human carotid artery [12]. We observed that variations in MAS increased with the velocity amplitude but decreased with mean velocity and stroke rate. These trends corresponded with those revealed by experimental results for porcine blood echogenicity. However, the motion of each RBC motion at a certain radial position and time is difficult to investigate through both in vitro and in vivo experiments given the limitations of physical conditions and measurement systems. Simulation studies enable the analysis of the formation and spatiotemporal changes of RBC aggregates and the development of parabolic rouleaux. Therefore, this simulation study on the kinetics of RBC aggregation helps improve our understanding of the complicated mechanism of RBC aggregation and its cyclic and radial variation under sinusoidal pulsatile flow.

4.2. Cyclic and radial variation of rouleuax under sinusoidal pulsatile flow

As shown in Fig. 5, at the initial time, particles were located at random positions in the whole region that corresponded to a hematocrit of 30%. Under sinusoidal pulsatile flow, RBCs were more likely to form rouleaux at the lowest velocity region with time. Then, the accumulated rouleaux increased in size, and the central rouleaux moved faster to form a parabolic shape, which broke as a result of fast rouleaux movement at the central part of the tube. The simulation results were in agreement with the cyclic and radial variation pattern observed from in vitro experiments on porcine blood echogenicity [16,17]. The formation and breakdown of a parabolic shape are shown as the bright collapsing ring (BRCR) in cross-sectional B-mode images. Under certain hemodynamic conditions, the BRCR appears as a bright hyperechoic ring that converges from the periphery to the center of the tube and eventually collapses during a pulsatile cycle [16,17]. Meanwhile, longitudinal ultrasound images have been obtained to further describe the processes of RBC aggregation, and a parabolic shape was observed in B-mode longitudinal images [25,26].

Some papers reported the angular dependence of rouleaux under in vitro steady and oscillatory flow experiments [17,27,28]. For porcine and equine blood, an angular dependence was observed at certain specific shear rate conditions, suggesting the cone-shaped aggregated RBCs and orientation at an angle of about 25^∘ with respect to the tube axis [27,28]. Under oscillatory flow, the large rouleaux were formed and aligned at an angle of about 25^∘ near the tube center during flow acceleration from porcine whole blood [17]. This angular dependence may be from partial parabolic shape of rouleaux distribution near the tube center from the simulation results in Figs 4(c) and (d) and 6(c). However, ultrasonic investigations of RBC motion at a certain radial position and time and parabolic rouleaux formation are difficult given the limitations of the physical and measurement conditions of in vitro porcine blood experiments. By contrast, numerical simulation enabled the direct observation of the dynamic processes of RBC aggregation and the formation of a parabolic shape. The simulation facilitated the analytical observation of the formation and breakdown of a parabolic shape from the BRCR phenomenon. Previous experiments and the present numerical simulation provided similar RBC aggregation trends and parabolic profiles despite their different hemodynamic conditions.

4.3. Velocity amplitude dependence on RBC aggregation

We studied the effects of different hydrodynamic parameters, such as velocity amplitude, mean velocity, and stroke rate, on RBC aggregation under sinusoidal pulsatile flow. As shown in Fig. 5, peak MAS and its mean value increased as velocity amplitude increased. When the velocity amplitude icreased, additional rouleaux formed as a result of the high velocity difference. RBCs were more likely to interact to form large rouleaux in the low velocity region. In our previous experiments, as the peak speed increased from 10 cm/s to 25 cm/s, the BRCR phenomenon became apparent at the high speed levels. At a peak speed of 10 cm/s, cyclic and radial variations in echogenicity were minimal. The variation in echogenicity is observed as a BRCR of high contrast at a peak speed of 25 cm/s [16]. In contrast to the simulation results, however, high echogenicity is observed during the systolic phase of the flow. This difference may be attributed to low viscosity, which may cause the phase shift of the peak MAS. Moreover, hysteresis may account for this difference because MAS is a function of shear rate and shear rate history. In contrast to the simulation results, however, high echogenicity is observed during the systolic phase of the flow. Finally, the differences in velocity profiles between pulsatile flow in the previous experiment and sinusoidal pulsatile flow with a mean velocity in the simulation may contribute to this difference. Further systematic studies will be required to further understand this difference.

Under three different velocity amplitudes of 0.5, 1, and 1.5 cm/s, RBC aggregation intensified at 0.03, 0.05, and 0.1 s, respectively, as shown in Fig. 6. These results indicated that high velocity amplitude results in drastic cyclic and radial variations in RBC aggregation and the formation of a parabolic rouleaux. In porcine blood experiments, temporal variation became apparent as peak speed increases during systole; this behavior coincides with that demonstrated by the simulation results [16,17]. Although the experimental and simulation results cannot be directly compared because of their different conditions (e.g. whole blood viscosity, tube diameter, and velocity profile), the experimental and simulation results presented similar trends.

4.4. Mean velocity dependence on RBC aggregation

As shown in Fig. 7, the peak MAS decreased as mean velocity increased. As mean flow velocity increased, the mean shear rate increased, consequently reducing rouleaux formation. The smaller variation and peak of MAS with high mean velocity could be attributed to high mean shear rate, as also observed in previous in vitro experiments. An in vitro experiment showed that echogenicity and temporal variation decrease as mean flow increases during flow acceleration [17]. Notably, when we applied pure oscillatory flow in the experiment, the overall echogenicity weakened as mean steady flow increased from 10 cm/s to 20 cm/s at the center of the tube. Therefore, the mean shear rate became an important parameter that affects RBC aggregation whether in a pulsatile flow or steady flow system. The simulation results indicated that high mean shear rate should decrease RBC aggregation. Furthermore, the peak MAS shifted to earlier phase from the time of the lowest speed as the mean flow speed increased from 2 cm/s to 6 cm/s. The shift of the peak MAS with peak speed and mean flow speed resulted from the combined effects of shear rate and flow acceleration during a pulsatile cycle. Further systematic studies should be conducted to confirm this result.

4.5. Stroke rate dependence on RBC aggregation

The results of an in vitro porcine blood experiment showed that peak velocity increases as stroke rate increases because stroke rate and flow velocity are related [17]. Analyzing each parameter independently is difficult given their close interrelationships. However, numerical simulation can completely separate stroke rate from flow velocity. In the present study, we selected three stroke rates (60, 120, and 180 bpm) in consideration of the periodic boundary condition, wherein each particle that leaves the right boundary should be positioned at the same position with the same velocity and phase. The peak MAS was high at low stroke rate and decreased when stroke rate increased. Another previous experimental result showed that as stroke rate increases from 20 bpm to 60 bpm, the peak of the Doppler power at the tube center shifts toward the peak systole [13]. Nevertheless, the present simulation results did not provide any clue for the occurrence of phase shift with increased stroke rate. In vitro experiments on Doppler power measurements and porcine blood echogenicity revealed that the BRCR phenomenon weakens as stroke rate increases from 20 bpm to 60 bpm. The weakening of the BRCR phenomenon with limited variation in echogenicity at high stroke rates is in agreement with the present simulation results despite the different stroke rates used in the present numerical simulation and previous blood experiments. When mean flow velocity dominates and is sufficiently large, cyclic variation due to high mean shear rate is negligible. Even when mean flow velocity remains the same over one stroke period, high stroke rates resulted in low MAS variation in simulation and weak echogenicity in porcine blood.

4.6. Limitations and future

In nature, physiological pulsatile flow has a short duration and high acceleration during systolic phase and a long relaxation time during diastolic phase [29]. However, this study simulated the sinusoidal flow profiles of one frequency component with nonzero mean instead of physiological pulsatile flow with additional frequency components. Given that this velocity profile ignored volume conservation and incompressible fluid characteristics, future studies could focus on these limitations to provide rigorous analytical solutions. In this study, hematocrit was defined as a surface ratio of RBCs to whole blood within a rectangular area and not as a volumetric ratio given the two-dimensionality of the model. RBCs were designed as spheres instead of three-dimensional biconcave and flexible objects. Consequently, the simulated hematocrit did not precisely correspond to a true hematocrit of 30% in RBCs. Therefore, the parameters of this model need to be adjusted when translated to a three-dimensional model. Given the limitations of computation time in MATLAB, the assumed viscosity values in this study were lower than the real viscosity values. Specifically, the viscosity of human blood at 37^∘C is normally 0.003 Pa ⋅ s to 0.004 Pa ⋅ s, and in this simulation, the viscosity was 0.0027 Pa ⋅ s. Moreover, because hydrodynamic force is highly related with viscosity, RBCs can accelerate within a very short period under low viscosity. The step time was low to maintain RBCs to less than half of the tube length. For instance, in unit time, the program loop is more than 10,000 times. In future studies, high computation speed would solve the running time problem. At the same time, the ROI and the number of RBCs should be increased.

5. Conclusions

In this simulation study, we observed the kinetics, detailed processes, and cyclic and radial variation patterns of RBC aggregation under sinusoidal pulsatile flow. The cyclic and radial variation in RBC aggregation resulted from interaction forces and hydrodynamic force. We found that the MAS of RBCs and the cyclic variation in RBC aggregation increased with velocity amplitude but decreased with mean flow velocity and stroke rate. Furthermore, we simulated the dynamic process of RBC aggregation in a longitudinal tube. This process exhibited a parabolic shape. The simulation results for cyclic and radial variation patterns, including the parabolic shape of rouleaux, were in agreement with previous in vitro experimental results obtained with porcine blood. Our findings expand our understanding of the mechanisms of RBC aggregation in human blood vessels under physiological pulsatile flow.

Footnotes

Acknowledgement

This research was supported by the 2018 scientific promotion program funded by Jeju National University.

References

Baskurt

Neu

Meiselman

. Red Blood Cell Aggregation. CRC Press; 2011.

Lanotte

Mauer

Mendez

Fedosov

Fromental

Claveria

Red cells’ dynamic morphologies govern blood shear thinning under microcirculatory flow conditions. Proc Natl Acad Sci. 2016;113(47):13289–94.

Wells

Merrill

. Influence of flow properties of blood upon viscosity-hematocrit relationships. J Clin Invest. 1962;41(8):1591–8.

Vlachopoulos

O’Rourke

Nichols

. McDonald’s Blood Flow in Arteries: Theoretical, Experimental and Clinical Principles. CRC Press; 2011.

Paeng

Cao

Shung

. Doppler power variation from porcine blood under steady and pulsatile flow. Ultrasound Med Biol. 2001;27(9):1245–54.

Cloutier

Shung

. Study of red cell aggregation in pulsatile flow from ultrasonic Doppler power measurements. Biorheology. 1993;30(5–6):443–61.

Lin

Shung

. Ultrasonic backscattering from porcine whole blood of varying hematocrit and shear rate under pulsatile flow. Ultrasound Med Biol. 1999;25(7):1151–8.

Chang

Lee

Liao

Huang

. Study of the effects of flow acceleration on blood aggregation by high frequency duplex ultrasound: Under pulsatile flow. 2010;620–3.

Cao

Paeng

Shung

. The “black hole” phenomenon in ultrasonic backscattering measurement under pulsatile flow with porcine whole blood in a rigid tube. Biorheology. 2001;38(1):15–26.

10.

Huang

Liao

Lee

Shih

. The effect of flow acceleration on the cyclic variation of blood echogenicity under pulsatile flow. Ultrasound Med Biol. 2013;39(4):670–80.

11.

Huang

. High-frequency attenuation and backscatter measurements of rat blood between 30 and 60 MHz. Phys Med Biol. 2010;55(19):5801.

12.

Paeng

Nam

Shung

. Cyclic and radial variation of the echogenicity of blood in human carotid arteries observed by harmonic imaging. Ultrasound Med Biol. 2010;36(7):1118–24.

13.

Paeng

Shung

. Cyclic and radial variation of the Doppler power from porcine whole blood. IEEE Trans Ultrason Ferroelectr Freq Control. 2003;50(6):614–22.

14.

Thompson

Trudinger

Cook

. Doppler ultrasound waveforms in the fetal umbilical artery: quantitative analysis technique. Ultrasound Med Biol. 1985;11(5):707–18.

15.

Shung

. Cyclic variation of Doppler power from whole blood under pulsatile flow. Ultrasound Med Biol. 1996;22(7):883–94.

16.

Paeng

Chiao

Shung

. Echogenicity variations from porcine blood I: The “bright collapsing ring” under pulsatile flow. Ultrasound Med Biol. 2004;30(1):45–55.

17.

Paeng

Chiao

Shung

. Echogenicity variations from porcine blood II: The “bright ring” under oscillatory flow. Ultrasound Med Biol. 2004;30(6):815–25.

18.

Fenech

Garcia

Meiselman

Cloutier

. A particle dynamic model of red blood cell aggregation kinetics. Ann Biomed Eng. 2009;37(11):2299–309.

19.

Vincent

. The calculation of depletion layer thickness as a function of bulk polymer concentration. Colloids Surf. 1990;50:241–49.

20.

Joanny

Leibler

De Gennes

. Effects of polymer solutions on colloid stability. J Polym Sci. 1979;17(6):1073–84.

21.

Neu

Meiselman

. Depletion-mediated red blood cell aggregation in polymer solutions. Biophys J. 2002;83(5):2482–90.

22.

Liu

Zhang

Wang

Liu

. Coupling of Navier–Stokes equations with protein molecular dynamics and its application to hemodynamics. Int J Numer Methods Fluids. 2004;46(12):1237–52.

23.

Lovalenti

Brady

. The hydrodynamic force on a rigid particle undergoing arbitrary time-dependent motion at small Reynolds number. J Fluid Mech. 1993;256:561–605.

24.

Brenner

. Dynamics of neutrally buoyant particles in low Reynolds number flows. 1972;509–74.

25.

Paeng

Nam

Choi

Shung

. Three-dimensional reconstruction of the “bright ring” echogenicity from porcine blood upstream in a stenosed tube. IEEE Trans Ultrason Ferroelectr Freq Control. 2009;56(4):880–5.

26.

Nam

Paeng

Choi

Shung

. Ultrasonic observation of blood disturbance in a stenosed tube: Effects of flow acceleration and turbulence downstream. Ultrasound Med Biol. 2008;34(1):114–22.

27.

Qin

Durand

Cloutier

. Kinetics of the “black hole” phenomenon in ultrasound backscattering measurements with red blood cell aggregation. Ultrasound Med Biol. 1998 February 01;24(2):245–56.

28.

Allard

Cloutier

Durand

. Effect of the insonification angle on the Doppler backscattered power under red blood cell aggregation conditions. IEEE Trans Ultrason Ferroelectr Freq Control. 1996;43(2):211–19.

29.

Yellin

. Laminar-turbulent transition process in pulsatile flow. Circ Res. 1966 Oct;19(4):791–804.