Abstract
Most of the existing hemolysis mechanism studies are carried out on the macro flow scale. They assume that the erythrocyte membranes with different loads will suffer the same damage, which obviously has limitations. Thus, exploring the hemolysis mechanism through the macroscopic flow field information is a tough challenge. In order to further understand the non-physiological shear hemolysis phenomenon at the cell scale, this study used the coarse-grained erythrocytes damage model at the mesoscopic scale based on the transport dissipative particle dynamics (tDPD) method. Combined with computational fluid dynamics the hemolysis of scalarized shear stress () in the clearance of “Impella 5.0” was evaluated under the Lagrange perspective and Euler perspective. The results from the Lagrange perspective showed that the change rate of scaled shear stress () was the most critical factor in damaging RBCs in the rotor region of “Impella 5.0”and other transvalvular micro-axial blood pumps. Then, we propose a dimensionless number with time integration based on to evaluate hemolysis. The Dissipative particle dynamics simulation results are consistent with the evaluation results, so may be an important factor in the hemolysis of VADs. Finally, we tested the hemolysis of 30% hematocrit whole blood in the “Impella 5.0” shroud clearance from the Euler perspective. Relevant results indicate that because of the wall effect, the RBCs near the impeller side are more prone to damage, and most of the cytoplasm is also gathered at the rotor side.
1 Introduction
Heart failure (HF) is a complex clinical syndrome caused by structural and functional damage to the filling or ejection of blood from the ventricles. According to statistics, the prevalence of HF with clinical symptoms is 1.3%–1.8%, the prevalence of asymptomatic HF is 1.5%–2%, and the prevalence of HF in people over 65 years old is as high as 10% (Basil et al., 2017). Cardiogenic shock (CS), which occurs in the late stages of HF, is very lethal, and patients must receive timely cardiac perfusion therapy to save their lives. Thoracotomy implantation of an artificial heart can be very dangerous for patients with end-stage HF. The “Impella” series transvalvular axial blood pump can be implanted from the patient’s femoral artery (similar to the implantation of a vascular stent) to provide additional flow to the heart ejection, thereby extending the patient’s life. Currently, transvalvular micro-axial blood pumps can provide cardiac perfusion at various flow rates in a minimally invasive manner. The “Impella” pump series (Abiomed, Danvers Massachusetts, US), which represents this category, is exclusive to “Abiomed” company in the United States and is the world’s largest supplier of such left ventricular assist devices (LVADs). Available clinical data showed that the “Impella” series could achieve complete left ventricular support in a short period, providing up to 5.5 L/min of flow with a good prognosis (Mehra et al., 2019). In FDA testing, the mean time to implantation of Impella 5.0 was 16.3 ± 4.7 days, 9/22 patients had a ventricular recovery, 6/22 patients were free from the risk of acute heart failure, and the 30-day survival rate was 62% (Pahuja et al., 2021). However, various blood-contacting medical devices (BCMDs), including “Impella 5.0,” produce non-physiological shear stress that can damage red blood cells (RBCs) and the freeing of large amounts of hemoglobin (Hb) into the plasma. It can cause damage to the liver and kidneys, leading to organ failure (Olsen, 2000). Several studies have shown that most of the dangerous effects caused by BCMDs are associated with blood damage (Ellis et al., 1998; Maraj et al., 1998). During the design phase of VADs, the hydraulic properties must be carefully weighed against the blood compatibility. The empirical power-law index of hemolysis (IH) formulation is the most common method when using computational simulation techniques to assess hemolysis in VADs. The researchers first use computational fluid dynamics (CFD) to evaluate various flow field information in the vicinity of VADs and then process the flow field results to predict the IH. Eq. 1 is an empirical IH formula obtained through in vitro experiments summarized by Giersiepen et al. (Giersiepen et al., 1990)where Hb is the total amount hemoglobin, ∆Hb is the amount of variation of hemoglobin in plasma. are empirical constants, is the exposure time of blood in a non-physiological flow field, and is the equivalent shear stress of the flow field. It can be found that the empirical power-law model greatly simplifies the hemolysis process, and the factors affecting hemolysis are exposure time and shear stress. Song et al. improved the IH power-law model by integrating along the blood flow trajectory, then summarized several cumulative exponential formulas for blood damage (Song et al., 2003; Yano et al., 2003). However, there are limitations to the use of macroscopic scale methods to assess hemolysis: 1. The intrinsic nature of hemolysis is erythrocyte damage, so assessing hemolysis requires consideration of the mechanical properties of the erythrocyte membrane and erythrocytes of different loads should not be subjected to the same incremental damage values (Grigioni et al., 2004; Grigioni et al., 2005). The clearance between the rotor and shroud is usually the location where VADs’ hemolysis is most severe. The link between erythrocyte kinetics and blood rheology in small-size flow becomes important (Lanotte et al., 2016).
Dissipative particle dynamics (DPD) has shown great potential for development in various simulations at the cellular scale. Fedosov et al. have developed a spring network erythrocyte model and “erythrocyte-erythrocyte” and “erythrocyte-platelet” aggregation force models based on DPD. Their study explains the Rouleaux stacking of erythrocytes and the Fahraeus-Lindquist effect of blood (Fedosov et al., 2010; Fedosov et al., 2011). Ye et al. have developed a cellular-scale model of thrombogenesis. They explored the adsorption properties of platelets (PLTs) and the process of thrombus formation in stenotic vessels using smooth dissipative particle dynamics (SDPD) (Ye et al., 2020). We have established an erythrocyte damage model and a transport dissipative particle dynamics (tDPD) based hemoglobin spillover model based on a coarse-grained erythrocyte model (Xu et al., 2022; Xu et al., 2023).
This study combined CFD and tDPD to explore the sensitive factors that lead to blood damage in non-physiological shear stress flow around “Impella 5.0". It is expected that, the cell-scale hemolysis analysis of “Impella 5.0” may hopefully guide the optimization design of VADs in the future.
The “Method” section introduces the transport dissipative particle dynamics (tDPD) and the RBC damage model in DPD systems. This section also provides the aggregation model of the RBCs population and the mapping relationship between DPD units.
In the “Result” section, the scalarized shear stress in Impella 5.0 trace lines was obtained from CFD as the boundary condition for DPD simulation to evaluate the sensitive factors of RBC damage at the cellular scale.
We propose a new hemolysis index based on cell-scale, which can calculate the damage rate of RBC membrane in a trace line and reflect the hemolysis situation. Then, this study collected a large amount of trace data and reflected it in the flow field of “Impella,” in order to observe the high hemolytic regions in the flow field. Finally, we used the traditional hemolysis index to observe the hemolysis details in the “Impella” clearance region.
2 Methods
2.1 Transport dissipative particle dynamics
Dissipative particle dynamics (DPD) was developed on the basis of molecular dynamics (MD) and combined with statistical mechanics theory (Hoogerbrugge and Koelman, 1992; Groot and Rabone, 2001). In the DPD method, the flow field is described as a series of interacting particles. The DPD uses a number of discrete particles to represent the state of the flow, representing the evolution of the system by recording the position and velocity of the particles at each time step. Eq. 2 is the basic equation of DPD.where the interaction force between two adjacent particles consists of three components, where is time and are the velocity of the particle and the external force on it, respectively. The conservative force is a soft repulsive force along the relative direction between the particles, and the magnitude of the conservative force affects the intensity of the momentum exchange within the system. The dissipative force is used to hinder the relative motion between the particles and characterizes the viscous dissipation within the system. Its magnitude is related to the particles’ relative velocity. The random force characterizes the random momentum exchange due to the thermal fluctuation of the system (Espanol and Warren, 1995). The expressions for the three basic forces of DPD are shown in Eq. 3.
is the conservative force coefficient, which needs to meet to ensure that the compressibility coefficient of water is close to the real physical scale (Groot and Warren, 1997). is the particle distance between particle and particle, is the unit direction vector. and are parameters of dissipative force and random force, and are the distribution function. To satisfy the Boltzmann distribution and diffusion-fluctuation theorem theory of the microcanonical ensemble system, the dissipative force and random force coefficients have the following relationship: is the Boltzmann constant, is the temperature of the system, is the particles cutoff radius, is a correction factor proposed by Fan et al. (Fan et al., 2006). When using the DPD method to solve fluid dynamics problems, the relationship between particle diffusion and viscous diffusion is not physical. Reducing the numerical value of can improve this nonphysical relationship, our paper took .
Li et al. (2015) added the discrete diffusion equation to the basic DPD equation. They proposed the transport dissipative particle dynamics model (tDPD), which can accomplish the advection-diffusion-reaction process (ADR) in flow The DPD-based discrete diffusion equation is derived from the continuous medium diffusion equation, and the equation framework is essentially the same as the DPD basic equation, as shown in Eq. 5.
In the diffusion equation for a continuous medium, is the medium concentration, is the diffusion coefficient, and is the source term. In the diffusion equation for the discrete phase of DPD, is the component concentration of the particle, is the concentration source term, and is the concentration flux of the corresponding particle exchange, where the Fickian flux represents the diffusion of the medium during convection, and the random flux is the medium diffusion during random momentum exchange. Eq. 6 is the basic equation of concentration diffusion of tDPD, is the component concentration of the particle, is the concentration source term, is the concentration flux of the corresponding particle exchange, which can be divided into Fickian flux and Random flux :where refers to Fickian diffusion flux intensity and random diffusion flux intensity. These two coefficients also have a conversion relationship, as shown in Eq. 7.where is the molecular weight of the solute, is the number density of solvent particles. Not only can tDPD represent stochastic diffusion due to thermal fluctuation phenomena, but the particle-based computational method can also save significant computational resources when dealing with ADR problems. This study will use tDPD to simulate erythrocyte dynamics in the flow around “Impella 5.0” and to evaluate the processes of erythrocyte damage and Hb leakage to plasma. We compared the simulation results with other scholars’ hemolysis experiments in vitro many times and empirically set the Fickian diffusion coefficient in the DPD system (Zhang et al., 2011).
2.2 Coarse-grained erythrocyte damaged model
Prior to the twenty-first century, researchers exploring hemodynamics would typically treat blood as a Newtonian fluid. With the advancement of microfluidic techniques and rheological theories, blood-related research is set to enter a new stage. Incorporating erythrocyte modeling into blood flow has become a topic of great interest in recent years. Many optical tweezers experiments have been conducted to test the deformation of erythrocyte membranes under various states of stress and to derive the material parameters of erythrocyte membranes (Zhu et al., 2020). Dao et al. developed a viscoelastic erythrocyte membrane constitutive model. They completed the erythrocyte’s center stretching simulation using the three-dimensional finite element method, and their simulation results were also in full agreement with the optical tweezer stretching experiments (Dao et al., 2003). There are hundreds of millions of RBCs in the blood, and the RBCs have large deformation and viscoelastic characteristics. Therefore, it is very difficult to construct the fluid-structure coupling problem of RBCs-blood flow based on the traditional mesh method. Therefore, in recent years, researchers have proposed various erythrocyte models to accomplish blood flow simulation at the cellular scale with fewer computational resources. A low-dimensional model of erythrocytes consisting of ten colloidal particles was proposed by Pan et al. (Pan et al., 2010). However, the erythrocyte shape differs too much from the actual one to reproduce the deformed state of erythrocytes in flow. It is exciting that the spring network-based model of coarse-grained erythrocytes proposed by Fedosov et al. is able to achieve various deformation features of erythrocytes with a fewer number of particles (Fedosov et al., 2010). In their research tests, the error between the 500-particle constructed erythrocyte model and the real optical tweezer stretching experiment is less than 1%. The energy of the spring network-based erythrocyte model can be expressed as:where are the surface elastic energy, bending energy, surface area conservation penalty function and volume conservation penalty function of RBCs, respectively. The erythrocyte spring network model uses the wormlike chain nonlinear spring model, and the surface elastic energy can be expressed as:where is the number of the bonds of the erythrocyte system, is the persistence length, kp is the spring constant, is the energy unit, is the length of the spring , is the maximum spring extension, and . The bending strength of erythrocyte membrane can be expressed as:where is the bending constant, is the instantaneous angle between two adjacent triangles having the common edge , and is the spontaneous angle. Finally, and are penalty functions that constrain the area and volume of RBCs, which can prevent non-physical deformation:where is the number of triangular mesh in the erythrocyte model, is the equilibrium value of a triangle mesh area, and , and are the local area, global area and volume constraint coefficients, respectively. and are area and volume of RBCs at the present time step, and are area and volume in equilibrium. The RBC particle forces corresponding to the above energies are derived from the following formula:
The DPD parameters of RBCs and PLTs used in this study were referenced from Yazdani et al. (2021), as shown in Table 1.
TABLE 1
| 500 | 100 | 132.87 | 92.45 | 6.30 | 5000 | 5000 | |
| 48 | 10000 | 19.63 | 6.02 | 100 | 5000 | 10000 |
Model parameters of RBCs and PLTs. is the number of particles, is the shear modulus, the RBC’s initial area and initial volume are , , the bending constant is , the constraint constant of the surface area penalty function is , and the constant of the volume constraint penalty function is .
Erythrocyte destruction is based on the ultimate force between the cytoskeleton and the cell membrane tetramer. When the force between the cell membrane tetramer is greater than 24.9 pN, the erythrocyte membrane at this location is considered to have reached a state of local maximum strain (
Koshiyama and Wada, 2011). In addition, when the local maximum strain region of the erythrocyte membrane exceeds 6%, the erythrocyte membrane will produce irreversible damage and generate large-scale pores (
Leverett et al., 1972). Then cytoplasm such as Hb will diffuse from the pores into the plasma. After erythrocyte damage, cytoplasmic begins to spill out of the cell membrane, and the volume of erythrocytes is no longer conserved (
). Li et al. studied that the binding energy between the spectrin and cytoskeleton would decrease by 9 orders of magnitude when the intramembrane protein “Band3” embedded in the phospholipid bilayer and the anchor protein “Band4.1″connected to the spectrin and cytoskeleton were separated (
Li et al., 2007). Therefore, when the cell membrane produces pores, the cytoskeleton at the corresponding location would also be fractured and relaxed, the surface elastic energy of the cell membrane becomes lower (
). Therefore, we made the following assumptions in the RBC damage model:
1. If the spring force in the RBC system is higher than , the corresponding particles will be marked as “local maximum strain region particles”.
2. When the local maximum strain region of the erythrocyte membrane exceeds 6%, the corresponding area will produce large-scale holes, the cytoplasm will overflow the cell membrane, and the erythrocyte volume will no longer be conserved .
3. When the erythrocyte membrane is destabilized and damaged, the spectrin skeleton with immense stress begins to disconnect, manifested as connecting bonds relaxation in the coarse-grained model. Once the spring force is greater than 24.9pN, the maximum length of the spring will increase with the increase of the spring force.
The surface elastic energy provided by the spring will be changed :
The parameters are set to , is the spring force at the current moment, and is the critical force for bond relaxation. The energy composition of the damaged erythrocyte membrane can be expressed as (Xu et al., 2022):
In the coarse-grained erythrocyte damage simulation, if a certain erythrocyte membrane particle and other adjacent particles reach the local maximum strain, this particle will be labeled as the “Hemoglobin Diffusion Pore Particle”, implying that the cytoskeleton within this region will undergo relaxation. Meanwhile, the pore particle will generate a source term to diffuse Hb into the plasma, as shown in Figure 1B.
FIGURE 1
2.3 Aggregation force model and units mapping relationships
There is a significant aggregation effect of RBCs in healthy blood. The RBC population exhibits a rouleaux coin-like stack during low-speed movement. Therefore, modelling the aggregation force of red blood cells is necessary. Besides, in a non-physiological shear stress flow, the collision among erythrocytes is frequent, the aggregation force model can provide a great repulsive force over short distances (similar to LJ potential), preventing the overlap of RBC membranes. Fedosov et al. used Lennard-Jones (LJ) potential to represent the short-range repulsive force and Morse potential to represent the long-range attraction between RBCs (Fedosov et al., 2011). However, the time step of the DPD simulation is larger than MD, and the LJ potential may generate an extremely large force between two close particles during the simulation, causing the computing system to crash. Therefore, Yazdani et al. used only the Morse potential in DPD to build cell aggregation model (Yazdani and Karniadakis, 2016). Our study also uses the Morse potential to model the aggregation force between RBCs. The basic form of the Morse potential is shown in Eq. 8.
The equilibrium distance is set . are Morse potential parameters. We tested the parameters using the aggregation force test proposed by Fedosov et al. When , a uniform stretch of or a unilateral stretch of will depolymerize two Rouleax-like stacked erythrocytes. The parameters of the aggregation force model between RBCs and PLTs are shown in Table 2.
TABLE 2
| 1.6 | 1.5 | 0.3 | |
| 2.0 | 2.0 | 1.0 | |
| 2.0 | 2.0 | 1.0 |
Model parameters for aggregation forces between RBCs and PLTs.
The units used in the DPD method are not based on real scale units, and the results obtained from DPD simulation need to be mapped to real units through dimensional analysis. The study by Groot et al. shows that, when using the Verlet-Velocity integration algorithm, the DPD time steps less than 0.04 can stabilize the Boltzmann temperature (Groot and Warren, 1997). This study takes a time step size of 0.001, with each step corresponding to s of the real-time unit.
The focus of this study is on the motion and deformation state of RBCs, so the length , cell membrane shear modulus , and characteristic viscosity are set as the basic units of the DPD calculation system (Yazdani et al., 2021). The corresponding relationships are shown in Table 3.
TABLE 3
Correspondence between DPD basic units and real physical units.
The scaling relationship between the basic physical quantities and other physical quantities can be obtained according to the dimensional principle, as shown in Eq. 16.where are the force, velocity, time and energy respectively, the subscript represents the real physical unit, and the subscript represents the DPD unit. The DPD particle mass is set to , the energy unit is , and the cutoff radius is set to . The calculation process involves a total of four types of particles, RBC particles, plasma particles, cytoplasmic particles and PLT particles. Table 4 shows the conserved force coefficients and dissipative force coefficients between different types of particles.
TABLE 4
DPD coefficients between different types of particles, where is red blood cell particle, is plasma particle, is cytoplasmic particle, is platelet particle, is the conserved force coefficient, is the dissipative force coefficient.
3 Result
3.1 Non-physiological shear stress flow field analysis for “Impella 5.0”
The two most commonly used computational methods for CFD evaluation of the hydraulic performance of VADs are the multiple reference coordinate method (MRF) and the slip mesh method (SMM). The MRF method “freezes” the rotor in a phase. At the same time, the boundary is given a rotational speed. The VAD’s rotation is reflected by the relative rotational speed of the boundary and the rotor, thus calculating the quasi-constant flow field of the rotor at a certain phase. In the SMM method, the fluid mesh in the rotor region will rotate with time, more closely resembling the true rotation state. Wu et al. found that the time step significantly impacts the SMM method’s computational accuracy, the calculation time step must be less than for the calculation results to converge (Wu et al., 2019). The “Impella 5.0” rotor in this study is constant at 30000 rpm, which can be considered as a steady-state calculation. Therefore, our study adopted the MRF method that consumes less computational resources.
As shown in Figure 2, the “Impella 5.0” can be divided into four components. The “Shroud” is able to hold the heart valve opened and provide space for the rotor to work. The “Impeller” is the device’s rotor, which at 30,000 rpm is able to provide an additional of ejection volume to the patient’s heart. The “Diffuser” is the device’s rear impeller, which serves as a guidance and buffering, can reduce energy loss due to vortex generation, and the “Motor” provide energy for device. This section evaluates the flow field and hydraulic performance of “Impella 5.0” using FLUENT 21.0 (ANSYS, Inc., United States). We divided the simulation area into three parts. The “Rectification Zone” enables the boundary flow to develop fully before entering the rotor area, which helps to improve the numerical simulation stability; The “Rotor Zone” is an important area for our analysis in this study. The high-speed rotor is able to boost the flow rate of 4L/min, meanwhile generating a large amount of non-physiological shear stressful flow; The “Tailrace Zone” can observe the development of tail flow and the formation of vortices. The inlet boundary condition uses a volume flow inlet boundary condition of , the outlet boundary condition uses a pressure outlet boundary, the impeller rotation computational method uses the MRF, and the other areas are set to a no-slip wall boundary. Our CFD model uses the SST for the turbulence model, the “SIMPLE” algorithm for the mesh difference method, and no-slip wall boundary conditions. The CFD numerical calculation can be considered to have converged, when the residual of the continuous equation is less than the residual of the equation is less than . The average results of the last 1,000 steps are calculated as the final result.
FIGURE 2
The domain segmentation is shown as Figure 3A. The whole domain has a total length of 15.5D, which is 93 mm. The Mesh independence was carried out by rescaling the basic mesh in Figure 3. It can be observed that when the number of grids exceeds 5 million. The pressure change is very small, and the calculation results are independent of the number of meshes. This study used 8.29 million meshes for CFD calculations. In addition, the quality of the final meshing is shown in Figure 3D, which was judged by mesh skewness:where , are the largest and smallest angle in the face or cell, is the angle for an equiangular face/cell (such as 60 for a triangle, 90 for a square). Thus, the smaller the skewness is, the better the mesh is.
FIGURE 3
Scalarized the shear stress tensor not only makes it easier to quantify hemolysis, but also greatly simplifies the analysis process. Thus, tensor scalarization is usually required for the VADs hemolysis evaluation, and the most popular method of scalarization method was proposed by Bludszuweit et al. (Bludszuweit, 1995):where is the scalarized shear stress and is the shear stress tensor in the Cartesian coordinate system . The blood shear thinning model adopts the Carreau-Yasuda non-newtonian fluid model (Gijsen et al., 1999), which can be expressed as: is the dynamic viscosity of blood, is the shear rate of the flow domain, other parameters are constants: . Arvand et al. proved that Reynolds stress would produce a non-physical shear stress value during numerical calculations, which is usually not considered in the scalarization of shear stress (Arvand et al., 2005; Ge et al., 2008). In addition, Tullio et al. explored the effects of viscous shear stress and turbulent Reynolds stress on hemolysis of the BCMDs. They found that the scale of Reynolds stress is very large, and in micro-scale blood flow, viscous shear stress is the hemolysis dominant factor (De Tullio et al., 2009). There is no unified conclusion on how to incorporate turbulence into the process of hemolysis evaluation. Most of the serious hemolysis regions of the transvalvular micro-axial blood pumps are in the impeller clearance. The viscous stress effect on hemolysis is far stronger than Reynolds stress, so our study does not consider Reynolds stress in the VADs shroud clearance.
Figure 4A extracts the velocity contour and two-dimensional streamlines of the section. The flow characteristic dimension from the rotor zone to the tailrace zone is relatively large so that two vortices are formed at the corner, and these two vortices are not symmetrical. After careful observation, we can find that the tail flow field on the right side is interfered by the rear impeller, thus suppressing the formation of vortices. Figure 4B shows the three-dimensional streamline. The gap at the tip of the rotor is the position with the maximum flow velocity, the blood flow velocity can reach . Figure 4C shows the contour of section. Areas with high shear stress are concentrated on the gap, and the of most rotor areas is .
FIGURE 4
The of the blood flow is highly positively correlated with IH. We intercepted the data of three sections in the Y direction to facilitate a clearer analysis of the hemolytic region, as shown in Figure 5. Section A-A is the top section, and the high shear stress area is concentrated at the impeller edge. In section B-B, the high area of the clearance is significantly increased, and is concentrated in the clearance between the impeller and shroud. When the flow developed to the C-C position, the high stress area has already spread to most of the domain.
FIGURE 5
3.2 Evaluation of erythrocyte damage based on the Lagrange analysis
Section 3.1 carried out a CFD-based simulation on the rotor region of “Imeplla 5.0.” We analyzed the distribution of in the flow field at the macro-scale to evaluate the region prone to hemolysis. In order to observe the hemolysis details of blood in the non-physiological shear stress flow from the cell-scale, we will track the movement details of RBCs in the flow field. This section will extract four characteristic trace lines from the rotor zone of “Impella 5.0,” then use the coarse-grained RBC model and DPD method to simulate the RBC motion state on the trace lines, as shown in Figure 6A. It is noteworthy that, the time step in the DPD simulation is far less than the shear stress sampling frequency on the trace line. Suppose the CFD trace data is directly imported into the DPD program. In that case, it will result in the Dirichlet velocity boundary condition being unable to derive from time, which brings more difficulties for further analysis. Therefore, we should fit a complete function curve based on the CFD data points, even if the fitting operation loses some accuracy.
FIGURE 6
The relationship curve between and time on the four trace lines of red, blue, purple and green is also displayed in corresponding colours in Figure 6B. The red and purple traces contact the impeller edge more frequently. The change of is intense, and the curve has 2-3 peaks, which can change from 10Pa to 150Pa in a short time. The degree of RBC damage may also be relatively serious. The green curve represents the flow state of most of the blood. The RBCs will enter the rear impeller along the flat position of the vane without passing through the areas prone to hemolysis, and the is between 20 Pa and 40 Pa.
This study uses the “Large-Scale Atomic/Molecular Massively Parallel Simulator” (LAMMPS) as the DPD particle simulation platform and uses code C++ to write specific RBC damage model functions. In DPD simulation, a 500-particle RBC model is used for simulation, and the number density of liquid particles is 3 (Feng et al., 2012). Figure 7A shows the schematic diagram of the shear test calculational settings, and the RBC’s deformation state under stable shear stress. In reality, RBCs will be subjected to multi-directional stress, and the flow characteristics are more complex. Scalarization is the most commonly used method in hemolysis analysis, so we adopted scalarized shear stress () to reflect hemolysis and RBC’s damage. This section used the plane Couette flow formed by the movement of the boundary to characterize the on RBCs. To reflect the stress of micro-scale vortices on RBCs, this section evaluated the Kolmogorov scale of “Impella 5.0”, and the DPD calculation scale is also set near the Kolmogorov scale:
FIGURE 7
Is the Kolmogorov scale, the characteristic scale is set as the pipe diameter, is the kinematic viscosity of blood. When the flow rate , is the characteristic velocity. From Eq. 20, we can calculate .
Figures 7B–E shows the deformation state of RBCs under different . The characteristics of the red trace (Figure 7B) and the purple trace (Figure 7D) are similar, both have a great change rate of . The increases from 10 Pa to 150 Pa, taking only . Although the peak of these two traces can reach 160 Pa, the exposure time of RBCs in high shear stress is also relatively short. Therefore, RBC has not been fully sheared under the action of the first peak . With the increase of exposure time, RBC after the second peak will change into an elongated ellipsoid-shaped under the action of shear stress flow. Although the peak in the blue trace (Figure 7C) is only 120 Pa, the time of RBCs exposed to shear stress flow above 60 Pa is much longer than that in the red and purple trace lines. Therefore, under the 120Pa peak, the RBCs also change into tank-treading motion form. The peak of the green trace is relatively low. In the low shear stress, the RBC motion form presents a rolling, folding and multilobe shape, and the changing shape is relatively rich. According to the simulation results based on Lagrange method. In addition to and exposure time, the shear stress change rate may also be one of the sensitive factors of hemolysis.
To further verify the effect of shear stress change rate () on erythrocyte damage. In the DPD simulation, we recorded the ratio of the local maximum stress area (RBC’s damage area) to the total area of the cell membrane () every 500 computing steps. So as to evaluate the damage degree of RBCs quantitatively, as shown in Figure 8A–C. The RBCs on the green trace line are not damaged, so this paragraph does not discuss the green trace. The RBC damage curve simulated by DPD has a stepwise growth trend. When the trace line passes through the high shear stress flow in the impeller clearance every time, the cytoskeleton breaks suddenly after exceeding the maximum stress, and the damage to the erythrocyte membrane will increase sharply. Interestingly, although the peak of the blue trace can reach 120 Pa and the exposure time is relatively long, the RBCs have been fully stretched into an ellipsoid-shaped, and the damage rate is the lowest, less than 10%. The red and purple traces have two peak regions, and the corresponding is also much higher than the blue trace, and the RBC damage rate is as high as 50%. So far, in “Imeplla 5.0” and other transvalvular micro-axial blood pumps, we may judge that the is one of the most important factors affecting hemolysis. Moreover, there are more regions of high in the trace passing through the impeller clearance frequently, resulting in more serious hemolysis. Further, we take as the main variable to construct a dimensionless number () related to RBC damage, and its expression is as follows:
FIGURE 8
The dimensionless number of RBC damage should be expressed as an integral relationship with time, is the average energy of blood flow in the rotor zone. Our previous studies have shown that RBCs can be damaged when the shear stress is above 80 Pa. Thus, in order to correctly express the accumulation of on blood damage, we add two limiters to the dimensionless integration formula. indicates that only flow acceleration can injure RBCs, while flow deceleration does not affect RBC damage; represents that the shear stress flow above 80 Pa can affect the damage amount of RBCs (Xu et al., 2022).
Figures 8D–F shows the relationship between the dimensionless number and time, which can directly evaluate the blood injury status through flow field information. Comparing the curve obtained from the flow field information with the curve obtained from the DPD simulation, although there are some errors in the damage peak data, the damage time point can coincide. This phenomenon also explains that the may be an important factor leading to hemolysis in the high-speed rotating impeller zone. Thus, in designing VADs, it may be necessary to consider adding a flow buffer layer to reduce the shear stress change rate to inhibit hemolysis.
3.3 Hemolysis evaluation in position of rotor clearance
To analyze erythrocyte damage from Euler perspective, this section will monitor the of a fixed area in the gap. The monitoring track will also change into a complete circular edge when the impeller rotates for one circle. Figure 9A shows the intercepted data area and is marked as a magenta curve. At the same time, we extracted the distribution contour in the rotor clearance (Figure 9B). The at the contact position with the impeller is higher than 100 Pa, and the average in the clearance is between 60 Pa and 70 Pa. Figure 9C shows the -t curve of the monitoring area in three rotating cycles, lasting for .
FIGURE 9
Next, we will use tDPD to analyze the RBCs’ population morphology and 30% hematocrit whole blood hemolysis in the area. There are 270 RBCs and 50 PLTs in the calculation area. The simulation involves 406200 particles, 410280 bonds, 273520 angles and 410280 dihedrals. It takes about 48 h to compute every 100000 calculation steps. The figure hides plasma particles and cytoplasmic particles to show the simulation results more clearly. Figure 10 shows the process of tDPD simulating whole blood hemolysis under Euler’s perspective, and the simulation are divided into two steps. When the blood has not entered the rotor area, it will be Poiseuille flow driven by blood pressure. First, we should simulate the RBC movement and random distribution in Poiseuille flow. We use a pressure gradient to simulate the whole blood flow in the rectification zone, as shown in Figure 10A. This simulation step lasted 400000 steps, and the movement of RBCs has become stable. It can be found that the RBCs in the centre still maintain the “Rouleux” coin-like stack, which is consistent with the findings of Fedosov et al. and Yazdani et al. (Fedosov et al., 2014; Yazdani and Karniadakis, 2016). Figure 10B shows the second step of the simulation. When blood enters the rotor area, the flow in the impeller clearance can be simplified as Couette flow. We regard the tDPD simulation area as the monitoring region in the clearance between the impeller and the shroud. Because of the tiny monitoring area, the influence of curvature on flow is ignored in the tDPD simulation process. The moving wall at the lower side of the simulation area represents the rotor side, which can generate shear stress and injure RBCs.
FIGURE 10
We show the tDPD simulation results of the first, third and fifth peaks at the monitoring region, and the corresponding time points are marked in Figure 9C. IH calculation formula is as follows (Sohrabi and Liu, 2017):where is hematocrit, is the total hemoglobin, is the hemoglobin in plasma, and is the hemoglobin in healthy whole blood plasma. In healthy blood, no hemoglobin is free in plasma . Thus, this hemolysis index formula can be simplified. As shown in Figure 11A, most RBCs are still in the aggregation state in the first peak because of the short exposure time. Due to the interface non-slip effect, the RBCs near the wall have large deformation, while a small number of RBCs on the side of the rotor have been damaged. With the increased exposure time, most RBCs have shown a tank-treading motion state and were no longer clustered. The Hb overflowing from the damaged RBCs on the rotor side and the IH value near the wall increased significantly (Figures 11B, C). When the whole blood enters the rotor zone, RBCs are more likely to be damaged on the side of the rotor. Other cytoplasms, such as Hb, are also concentrated on the side of the rotor. It is worth noting that the more Hb and damaged RBCs, PLTs, and coagulation factors would be gathered on the side of the rotor, which may also be an indirect reason for the frequent attachment of thrombosis and gel on the VADs impeller.
FIGURE 11
4 Discussion
4.1 Conclusion
Combining computational fluid dynamics (CFD) and transport dissipative particle dynamics (tDPD), this paper evaluates the hemolysis in the rotor zone of “Impella 5.0” at the macro-scale and cell-scale. Then, we analyze the damage of RBCs and whole blood from the “Lagrange perspective” and the “Euler perspective” and summarize a dimensionless number used to evaluate the blood compatibility of “Impella 5.0” and other transvalvular micro-axial blood pumps.
We first used CFD to focus on the velocity and scalarized shear stress () distribution in the “Impella 5.0” rotor region. The velocity gradient and at the clearance between the impeller and the shroud are the highest, which are the areas prone to hemolysis in the rotor zone.
Then, we extracted four trace lines of the rotor region and analyzed the movement and deformation of RBCs on the four traces using the DPD method. Under low shear stress, RBCs can show various deformation states, such as rolling, folding, and multilobe. However, RBCs exposed to high shear stress for a long time will present ellipsoid-shaped, and the cell membrane will also be irreversibly damaged.
Notably, in the rotor zone of “Impella 5.0,” the change rate of scalarized shear stress may be one of the most important factors for RBC damage and hemolysis. Therefore, we proposed a dimensionless number of hemolysis based on in the Lagrange perspective. After the comparison, it was found that the result of RBC damage simulation was consistent with the result of dimensionless number evaluation, which can preliminarily confirm our conclusion. Therefore, in designing VADs, it may be necessary to add buffer segments to reduce the RBC’s damage of to reduce hemolysis symptoms effectively.
Finally, we continuously observed the of the three rotation cycles in the “Imeplla 5.0”clearance region from the Euler perspective. We used tDPD to simulate the RBC motion form and IH distribution with a hematocrit of 30% of the whole blood. Due to the boundary effect, the RBCs at the rotor side boundary are easily damaged, and the Hb and other cytoplasm are concentrated on the rotor side. In the same way, damaged RBCs, PLTs, and coagulation factors would also gather on the side of the rotor, which may also be the reason why the VADs impeller often produces thrombus and gel.
4.2 Limitations
Nevertheless, there are still many unsolvable limitations in this study. 1. Considering that the viscous shear stress dominates the non-physiological shear stress in the clearance, this study simplifies the model and does not consider the role of turbulent Reynolds stress. 2. DPD method is a mesoscopic scale calculation method which cannot produce complex flow. Thus, we simplify the blood flow in the clearance and use the plane Couette flow to characterize and quantify the damage of shear stress to RBCs. 3. The hemolysis index proposed in this study is an integral formula of Lagrange’s perspective. In addition, is also an index that is difficult to measure in the process of in vitro tests. The construction of the hemolysis formula based on the hemolysis principle from Euler’s perspective is still a huge challenge.
In addition, it is necessary to use tensors to characterize the damage of red blood cells. Currently, many scholars have obtained the relationship between red blood cell deformation and flow tensor (Faghih and Sharp, 2020), but no researchers have considered the mapping relationship between flow tensor and red blood cell damage. In future work, we need to construct a flow stress tensor based on complex flow in DPD method, then further improve mesoscale hemolysis research by combining flow tensor with blood damage models.
Multi-scale blood injury research is a huge challenge, as it is difficult to couple macroscopic hemolysis phenomena with microscopic damage to red blood cells. The RANS model used in this study is obviously not enough to characterize the Kolmogorov scale vortex that causes damage to red blood cells (Antiga and Steinman, 2009). Therefore, in future, we plan on combining the DNS method and the DPD method which may make the blood damage research more accurate.
Statements
Data availability statement
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.
Author contributions
All authors listed have made a substantial, direct, and intellectual contribution to the work and approved it for publication.
Funding
We gratefully acknowledge the Funds of the National Key R&D Program of China (Grant No. 2022YFC2402600) and National Natural Science Foundation of China (Grant No. 11972215, 12072174).
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
References
1
AntigaL.SteinmanD. A. (2009). Rethinking turbulence in blood. Biorheology46 (2), 77–81. 10.3233/bir-2009-0538
2
ArvandA.HormesM.ReulH. (2005). A validated computational fluid dynamics model to estimate hemolysis in a rotary blood pump. Artif. Organs29 (7), 531–540. 10.1111/j.1525-1594.2005.29089.x
3
BasilS.BryanS.ClydeW. (2017). Heart failure and sudden cardiac death. Card. Electrophysiol. Clin.9, 709–723. 10.1016/j.ccep.2017.07.010
4
BludszuweitC. (1995). Three-dimensional numerical prediction of stress loading of blood particles in a centrifugal pump. Artif. Organs19 (7), 590–596. 10.1111/j.1525-1594.1995.tb02386.x
5
DaoM.LimC. T.SureshS. (2003). Mechanics of the human red blood cell deformed by optical tweezers. J. Mech. Phys. Solids51 (11-12), 2259–2280. 10.1016/j.jmps.2003.09.019
6
De TullioM. D.CristalloA.BalarasE.VerziccoR. (2009). Direct numerical simulation of the pulsatile flow through an aortic bileaflet mechanical heart valve. J. Fluid Mech.622, 259–290. 10.1017/s0022112008005156
7
EllisJ. T.WickT. M.YoganathanA. P. (1998). Prosthesis-induced hemolysis: Mechanisms and quantification of shear stress. J. Heart Valve Dis.7 (4), 376–386.
8
EspanolP.WarrenP. (1995). Statistical mechanics of dissipative particle dynamics. Europhys. Lett.30, 191–196. 10.1209/0295-5075/30/4/001
9
FaghihM. M.SharpM. K. (2020). Deformation of human red blood cells in extensional flow through a hyperbolic contraction. Biomechanics Model. Mechanobiol.19 (1), 251–261. 10.1007/s10237-019-01208-3
10
FanX.Phan-ThienN.ChenS.WuX.Yong NgT. (2006). Simulating flow of DNA suspension using dissipative particle dynamics. Phys. Fluids18 (6), 063102. 10.1063/1.2206595
11
FedosovD. A.CaswellB.KarniadakisG. E. (2010). Systematic coarse-graining of spectrin-level red blood cell models. Comput. Methods Appl. Mech. Eng.199 (29-32), 1937–1948. 10.1016/j.cma.2010.02.001
12
FedosovD. A.DaoM.KarniadakisG. E.SureshS. (2014). Computational biorheology of human blood flow in health and disease. Ann. Biomed. Eng.42 (2), 368–387. 10.1007/s10439-013-0922-3
13
FedosovD. A.PanW.CaswellB.GompperG.KarniadakisG. E. (2011). Predicting human blood viscosity in silico. Proc. Natl. Acad. Sci.108 (29), 11772–11777. 10.1073/pnas.1101210108
14
FengR.XenosM.GirdharG.KangW.DavenportJ. W.DengY.et al (2012). Viscous flow simulation in a stenosis model using discrete particle dynamics: A comparison between dpd and cfd. Biomechanics Model. Mechanobiol.11 (1-2), 119–129. 10.1007/s10237-011-0297-z
15
GeL.DasiL. P.SotiropoulosF.YoganathanA. P. (2008). Characterization of hemodynamic forces induced by mechanical heart valves: Reynolds vs. Viscous stresses. Ann. Biomed. Eng.36 (2), 276–297. 10.1007/s10439-007-9411-x
16
GiersiepenM.WurzingerL. J.OpitzR.ReulH. (1990). Estimation of shear stress-related blood damage in heart-valve prostheses - invitro comparison of 25 aortic valves. Int. J. Artif. Organs13 (5), 300–306. 10.1177/039139889001300507
17
GijsenF. J. H.van de VosseF. N.JanssenJ. D. (1999). The influence of the non-Newtonian properties of blood on the flow in large arteries: Steady flow in a carotid bifurcation model. J. Biomechanics32 (6), 601–608. 10.1016/s0021-9290(99)00015-9
18
GrigioniM.DanieleC.MorbiducciU.D'AvenioG.Di BenedettoG.BarbaroV. (2004). The power-law mathematical model for blood damage prediction: Analytical developments and physical inconsistencies. Artif. Organs28 (5), 467–475. 10.1111/j.1525-1594.2004.00015.x
19
GrigioniM.MorbiducciU.D’AvenioG.BenedettoG. D.GaudioC. D. (2005). A novel formulation for blood trauma prediction by a modified power-law mathematical model. Biomechanics Model. Mechanobiol.4 (4), 249–260. 10.1007/s10237-005-0005-y
20
GrootR. D.RaboneK. L. (2001). Mesoscopic simulation of cell membrane damage, morphology change and rupture by nonionic surfactants. Biophysical J.81 (2), 725–736. 10.1016/s0006-3495(01)75737-2
21
GrootR. D.WarrenP. B. (1997). Dissipative particle dynamics: Bridging the gap between atomistic and mesoscopic simulation. J. Chem. Phys.107 (11), 4423–4435. 10.1063/1.474784
22
HoogerbruggeP. J.KoelmanJ. (1992). Simulating microscopic hydrodynamic phenomena with dissipative particle dynamics. Europhys. Lett.19, 155–160. 10.1209/0295-5075/19/3/001
23
KoshiyamaK.WadaS. (2011). Molecular dynamics simulations of pore formation dynamics during the rupture process of a phospholipid bilayer caused by high-speed equibiaxial stretching. J. Biomechanics44 (11), 2053–2058. 10.1016/j.jbiomech.2011.05.014
24
LanotteL.MauerJ.MendezS.FedosovD. A.FromentalJ-M.ClaveriaV.et al (2016). Red cells’ dynamic morphologies govern blood shear thinning under microcirculatory flow conditions. Proc. Natl. Acad. Sci.113 (47), 13289–13294. 10.1073/pnas.1608074113
25
LeverettL. B.HellumsJ. D.AlfreyC. P.LynchE. C. (1972). Red blood cell damage by shear stress. Biophysical J.12 (3), 257–273. 10.1016/s0006-3495(72)86085-5
26
LiJ.LykotrafitisG.DaoM.SureshS. (2007). Cytoskeletal dynamics of human erythrocyte. Proc. Natl. Acad. Sci.104 (12), 4937–4942. 10.1073/pnas.0700257104
27
LiZ.YazdaniA.TartakovskyA.KarniadakisG. E. (2015). Transport dissipative particle dynamics model for mesoscopic advection-diffusion-reaction problems. J. Chem. Phys.143 (1), 014101. 10.1063/1.4923254
28
MarajR.JacobsL. E.IoliA.KotlerM. N. (1998). Evaluation of hemolysis in patients with prosthetic heart valves. Clin. Cardiol.21 (6), 387–392. 10.1002/clc.4960210604
29
MehraM. R.UrielN.NakaY.ClevelandJ. C.YuzefpolskayaM.SalernoC. T.et al (2019). A fully magnetically levitated left ventricular assist device — final report. N. Engl. J. Med.380 (17), 1618–1627. 10.1056/nejmoa1900486
30
OlsenD. B. (2000). The history of continuous-flow blood pumps. Artif. Organs24 (6), 401–404. 10.1046/j.1525-1594.2000.06652.x
31
PahujaM.Hernandez-MontfortJ.WhiteheadE. H.KawaboriM.KapurN. K. (2021). Device profile of the Impella 5.0 and 5.5 system for mechanical circulatory support for patients with cardiogenic shock: Overview of its safety and efficacy. Expert Rev. Med. Devices19, 1–10. 10.1080/17434440.2022.2015323
32
PanW.CaswellB.KarniadakisG. E. (2010). A low-dimensional model for the red blood cell. Soft Matter6 (18), 4366. 10.1039/c0sm00183j
33
SohrabiS.LiuY. (2017). A cellular model of shear-induced hemolysis. Artif. Organs41 (9), E80–E91. 10.1111/aor.12832
34
SongX. W.ThrockmortonA. L.WoodH. G.AntakiJ. F.OlsenD. B. (2003). Computational fluid dynamics prediction of blood damage in a centrifugal pump. Artif. Organs27 (10), 938–941. 10.1046/j.1525-1594.2003.00026.x
35
WuX.ZhangX. W.HaoP. F.HeF. (2019). Comparison of three control strategies for axial blood pump. J. Mech. Med. Biol.19 (6), 1950058. 10.1142/s0219519419500581
36
XuZ.WangC.HeF.HaoP.ZhangX. (2023). Coarse-grained model of whole blood hemolysisand morphological analysis of erythrocytepopulation under non-physiological shear stressflow environment. Phys. Fluids35 (3), 031901. 10.1063/5.0137517
37
XuZ.WangC.XueS.HeF.HaoP.ZhangX. (2022). The erythrocyte destruction mechanism in non-physiological shear mechanical hemolysis. Phys. Fluids34 (11), 111901. 10.1063/5.0112967
38
YanoT.SekineK.MitohA.MitamuraY.OkamotoE.KimD. W.et al (2003). An estimation method of hemolysis within an axial flow blood pump by computational fluid dynamics analysis. Artif. Organs27 (10), 920–925. 10.1046/j.1525-1594.2003.00034.x
39
YazdaniA.DengY.LiH.JavadiE.LiZ.JamaliS.et al (2021). Integrating blood cell mechanics, platelet adhesive dynamics and coagulation cascade for modelling thrombus formation in normal and diabetic blood. J. R. Soc. Interface18 (175), 20200834. 10.1098/rsif.2020.0834
40
YazdaniA.KarniadakisG. E. (2016). Sub-cellular modeling of platelet transport in blood flow through microchannels with constriction. Soft Matter12 (19), 4339–4351. 10.1039/c6sm00154h
41
YeT.ShiH.Phan-ThienN.LimC. T. (2020). The Key events of thrombus formation: Platelet adhesion and aggregation. Biomechanics Model. Mechanobiol.19 (3), 943–955. 10.1007/s10237-019-01262-x
42
ZhangT.TaskinM. E.FangH-B.PamporiA.JarvikR.GriffithB. P.et al (2011). Study of flow-induced hemolysis using novel Couette-type blood-shearing devices. Artif. Organs35 (12), 1180–1186. 10.1111/j.1525-1594.2011.01243.x
43
ZhuR.AvsievichT.PopovA.MeglinskiI. (2020). Optical tweezers in studies of red blood cells. Cells9 (3), 545. 10.3390/cells9030545
Summary
Keywords
ventricular assist devices, cell-scale hemolysis simulation, dissipative particle dynamics, erythrocyte damage model, erythrocyte motion form
Citation
Xu Z, Chen C, Hao P, He F and Zhang X (2023) Cell-scale hemolysis evaluation of intervenient ventricular assist device based on dissipative particle dynamics. Front. Physiol. 14:1181423. doi: 10.3389/fphys.2023.1181423
Received
07 March 2023
Accepted
26 June 2023
Published
05 July 2023
Volume
14 - 2023
Edited by
Brent A. Craven, United States Food and Drug Administration, United States
Reviewed by
Nicolas Tobin, The Pennsylvania State University (PSU), United States
Katharine Fraser, University of Bath, United Kingdom
Updates
Copyright
© 2023 Xu, Chen, Hao, He and Zhang.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Xiwen Zhang, zhangxiw@tsinghua.edu.cn
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.