Cell-scale hemolysis evaluation of intervenient ventricular assist device based on dissipative particle dynamics

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 ($\mathrm{\tau }$) 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 ($\dot{\mathrm{\tau }}$) 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 ${\mathrm{D}}_{\mathrm{k}}$ with time integration based on $\dot{\mathrm{\tau }}$ to evaluate hemolysis. The Dissipative particle dynamics simulation results are consistent with the ${\mathrm{D}}_{\mathrm{k}}$ evaluation results, so $\dot{\mathrm{\tau }}$ 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)

$\begin{array}{c}IH=\frac{∆Hb}{Hb}=C{t}^{\alpha }{\tau }^{\beta }\end{array}(1)$

where Hb is the total amount hemoglobin, ∆Hb is the amount of variation of hemoglobin in plasma. $C,\alpha ,\beta$ are empirical constants, $t$ is the exposure time of blood in a non-physiological flow field, and $\tau$ 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.

$\begin{array}{c}\frac{d{v}_i}{dt}=\sum \frac{1}{m}\left(F_{ij}^C+F_{ij}^D+F_{ij}^R\right)+F_i^{ext}\end{array}(2)$

where the interaction force $F_{ij}$ between two adjacent particles $i,j$ consists of three components, where $t$ is time and $v_i,F_i^{ext}$ are the velocity of the $i-th$ particle and the external force on it, respectively. The conservative force $F_{ij}^C$ 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 $F_{ij}^D$ 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 $F_{ij}^R$ 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.

$\begin{array}{c}\left\{\begin{array}{l}{F}_{ij}^C=a_{ij}\left(1-r_{ij}\right)\stackrel{\sim }{{\mathit{r}}_{\mathit{i}\mathit{j}}}\\ F_{ij}^D=-\gamma {\omega }_D\left(r_{ij}\right)\left(\stackrel{\sim }{{\mathit{r}}_{\mathit{i}\mathit{j}}}\bullet {\mathit{v}}_{\mathit{i}\mathit{j}}\right)\stackrel{\sim }{{\mathit{r}}_{\mathit{i}\mathit{j}}}\\ F_{ij}^R=\sigma {\omega }_R\left(r_{ij}\right){\zeta }_{ij}\stackrel{\sim }{{\mathit{r}}_{\mathit{i}\mathit{j}}}d{t}^{-\frac{1}{2}}\end{array}\right.\end{array}(3)$

$a_{ij}$ is the conservative force coefficient, which needs to meet $a_{ij}=75k_{B}T/\rho$ to ensure that the compressibility coefficient of water is close to the real physical scale (Groot and Warren, 1997). $r_{ij}$ is the particle distance between $i-th$ particle and $j-th$ particle, $\stackrel{\sim }{{\mathit{r}}_{\mathit{i}\mathit{j}}}$ is the unit direction vector. $\gamma$ and $\sigma$ are parameters of dissipative force and random force, ${\omega }_D\left(r_{ij}\right)$ and ${\omega }_R\left(r_{ij}\right)$ 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:

$\begin{array}{c}\left\{\begin{array}{l}{\omega }_D\left(r_{ij}\right)={\left[{\omega }_R\left(r_{ij}\right)\right]}^2={\left(1-\frac{r_{ij}}{r_c}\right)}^s\\ \gamma =\frac{{\sigma }^2}{2k_{B}T}\end{array}\right.\end{array}(4)$

$k_B$ is the Boltzmann constant, $T$ is the temperature of the system, $r_C$ is the particles cutoff radius, $s$ 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 $s$ can improve this nonphysical relationship, our paper took $s=0.2$.

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.

$\begin{array}{c}\frac{dC}{dt}=D{\nabla }^{2}C+Q^S\Rightarrow \frac{d{C}_i}{dt}=\sum \left(Q_{ij}^D+Q_{ij}^R\right)+Q_i^S\end{array}(5)$

In the diffusion equation for a continuous medium, $C$ is the medium concentration, $D$ is the diffusion coefficient, and $Q^S$ is the source term. In the diffusion equation for the discrete phase of DPD, $C_i$ is the component concentration of the particle, $Q_i^S$ is the concentration source term, and $Q_{ij}$ is the concentration flux of the corresponding particle exchange, where the Fickian flux $Q_{ij}^D$ represents the diffusion of the medium during convection, and the random flux $Q_{ij}^R$ is the medium diffusion during random momentum exchange. Eq. 6 is the basic equation of concentration diffusion of tDPD, $C_i$ is the component concentration of the particle, $Q_i^S$ is the concentration source term, $Q_{ij}$ is the concentration flux of the corresponding particle exchange, which can be divided into Fickian flux $Q_{ij}^D$ and Random flux $Q_{ij}^R$:

$\begin{array}{c}\left\{\begin{array}{l}{Q}_{ij}^D=-{\kappa }_{ij}{\omega }_D{r}_{ij}\left(C_i-C_j\right)\\ Q_{ij}^R={ϵ}_{ij}{\omega }_R{r}_{ij}{\zeta }_{ij}∆t^{-1/2}\end{array}\right.\end{array}(6)$

where ${\kappa }_{ij},{ϵ}_{ij}$ refers to Fickian diffusion flux intensity and random diffusion flux intensity. These two coefficients also have a conversion relationship, as shown in Eq. 7.

$\begin{array}{c}{ϵ}_{ij}=m_s^2{\kappa }_{ij}\rho \left(C_i+C_j\right)\end{array}(7)$

where $m_s$ is the molecular weight of the solute, $\rho$ 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 ${\kappa }_{ij}={4*10}^{-8}$ 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:

$\begin{array}{c}{V}_{rbc}=V_s+V_b+V_a+V_v\end{array}(8)$

where $V_s,V_b,V_a,V_v$ 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 $V_s$ can be expressed as:

$\begin{array}{c}{V}_s={\displaystyle \sum _{j=\mathrm{1,2}\cdots N_s}}\left(\frac{k_{B}T\bullet l_m\left(3x_j^2-2x_j^3\right)}{4p\left(1-x_j\right)}+\frac{k_p}{l_j}\right)\end{array}(9)$

where $N_s$ is the number of the bonds of the erythrocyte system, $p$ is the persistence length, kp is the spring constant, $k_{B}T$ is the energy unit, $l_j$ is the length of the spring $j$, $l_m$ is the maximum spring extension, and $x_j=l_j/l_m$. The bending strength of erythrocyte membrane can be expressed as:

$\begin{array}{c}{V}_b={\displaystyle \sum _{j=\mathrm{1,2}\cdots N_s}}{k}_b\left(1-\cos \left({\theta }_j-{\theta }_0\right)\right)\end{array}(10)$

where $k_b$ is the bending constant, ${\theta }_j$ is the instantaneous angle between two adjacent triangles having the common edge $j$, and ${\theta }_0$ is the spontaneous angle. Finally, $V_a$ and $V_v$ are penalty functions that constrain the area and volume of RBCs, which can prevent non-physical deformation:

$\begin{array}{l}\left\{\begin{array}{c}{V}_a={\displaystyle \sum _{j=\mathrm{1,2}\cdots N_t}}\left(k_d\frac{{\left(A_j-A_0\right)}^2}{2A_0}\right)+k_a\frac{{\left(A^{rbc}-A_0^{tot}\right)}^2}{2A_0^{tot}}\\ V_v=k_v\frac{{\left(V^{rbc}-V_0^{tot}\right)}^2}{2V_0^{tot}}\end{array}\right.\end{array}(11)$

where $N_t$ is the number of triangular mesh in the erythrocyte model, $A_0$ is the equilibrium value of a triangle mesh area, and $k_d$, $k_a$ and $k_v$ are the local area, global area and volume constraint coefficients, respectively. $A^{rbc}$ and $V^{rbc}$ are area and volume of RBCs at the present time step, $A_0^{tot}$ and $V_0^{tot}$ are area and volume in equilibrium. The RBC particle forces corresponding to the above energies are derived from the following formula:

$\begin{array}{c}{F}_i=-\frac{\partial V\left(x_i\right)}{\partial r_i}\end{array}(12)$

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

TABLE 1. Model parameters of RBCs and PLTs. $N_v$ is the number of particles, ${\mu }_s$ is the shear modulus, the RBC’s initial area and initial volume are $A_0^{tot}$, $V_0^{tot}$, the bending constant is $k_b$, the constraint constant of the surface area penalty function is $k_d+k_a$, and the constant of the volume constraint penalty function is $k_v$.

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 ($V_v=0$). 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 ($V_s\to V_s^{\prime }$). Therefore, we made the following assumptions in the RBC damage model:

1. If the spring force in the RBC system is higher than $24.9\mathrm{p}\mathrm{N}$, 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 $\left(V_v=0\right)$.

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 $\left(V_s\to V_s^{\prime }\right)$:

$\begin{array}{l}\left\{\begin{array}{l}{V}_s^{\prime }={\displaystyle \sum _{j=\mathrm{1,2}\cdots N_s}}\left(\frac{k_{B}T\bullet l_m^{\prime }\left(3x_j^{{\prime }^2}-2x_j^{{\prime }^3}\right)}{4p\left(1-x_j^{\prime }\right)}+\frac{k_p}{l_j}\right)\\ l_m^{\prime }=l_m\left(r_0-D_e\left(e^{\beta \left(1-\alpha \left(\frac{F_b}{F_c}-1\right)\right)}-e^{\frac{\beta }{2}\left(1-\alpha \left(\frac{F_b}{F_c}-1\right)\right)}\right)\right)\end{array}\right.\end{array}(13)$

The parameters are set to $r_0=1.5,D_e=0.03,\beta =3,\alpha =20$, $F_b$ is the spring force at the current moment, and $F_c=24.9pN$ is the critical force for bond relaxation. The energy composition of the damaged erythrocyte membrane $V_{rbc}^{\prime }$ can be expressed as (Xu et al., 2022):

$\begin{array}{c}{V}_{rbc}^{\prime }=V_s^{\prime }+V_b+V_a\end{array}(14)$

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 $Q_i^S$ to diffuse Hb into the plasma, as shown in Figure 1B.

FIGURE 1

FIGURE 1. (A) Schematic diagram of the particles’ fundamental forces in DPD. The conservative force $F_C$ is a soft repulsive force indicating the intensity of the momentum exchange of the system, the dissipative force $F_D$ can reflect the viscosity of the system, and the random force $F_R$ can satisfy the dissipative rise and thermal fluctuation at the molecular scale. (B) Schematic diagram of the coarse-grained RBC damaged model. If all the springs around a particle reach the maximum value, this particle will be marked as a “Hemoglobin Diffusion Pore Particle” and diffuse cytoplasm into the adjacent plasma particles.

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.

$\begin{array}{c}{V}_{morse}=D_e\left(e^{2\beta \left(r_0-r\right)}-2\beta e^{\left(r-r_0\right)}\right)\end{array}(15)$

The equilibrium distance is set $r_0=0.3$. $D_e,\beta$ are Morse potential parameters. We tested the parameters using the aggregation force test proposed by Fedosov et al. When $D_e=1.6,\beta =1.5$, a uniform stretch of $7pN$ or a unilateral stretch of $3pN$ 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

TABLE 2. 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 $2.14\ast {10}^{-7}$ s of the real-time unit.

The focus of this study is on the motion and deformation state of RBCs, so the length $l$, cell membrane shear modulus $\mu$, and characteristic viscosity $\eta$ 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

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.

$\begin{array}{c}\left\{\begin{array}{c}\begin{array}{l}\frac{F_p}{F_m}=\frac{{\mu }_p}{{\mu }_m}\frac{l_p}{l_m}=4.73\hspace{0.17em}\ast \hspace{0.17em}{10}^{-2}\left(pN\right)\\ \frac{v_p}{v_m}=\frac{{\mu }_p}{{\mu }_m}\frac{{\eta }_m}{{\eta }_p}=4.64{\hspace{0.17em}\ast \hspace{0.17em}10}^{-3}\left(m/s\right)\end{array}\\ \frac{t_p}{t_m}=\frac{{\eta }_p}{{\eta }_m}\frac{{\mu }_m}{{\mu }_p}\frac{l_p}{l_m}=2.16{\hspace{0.17em}\ast \hspace{0.17em}10}^{-4}\left(s\right)\end{array}\right.\end{array}(16)$

where $F,v,t,E$ are the force, velocity, time and energy respectively, the subscript $p$ represents the real physical unit, and the subscript $m$ represents the DPD unit. The DPD particle mass is set to $m=1$, the energy unit is ${\left(k_{B}T\right)}_m=0.1$, and the cutoff radius is set to $r_c=1.58$. 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

TABLE 4. DPD coefficients between different types of particles, where $R$ is red blood cell particle, $S$ is plasma particle, $C$ is cytoplasmic particle, $P$ is platelet particle, $a_{ij}$ is the conserved force coefficient, ${\gamma }_{ij}$ 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 ${10}^{-5}s$ 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 $4L/\mathit{min}$ 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 $4L/\mathit{min}$, 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 $k-\omega$ 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 ${10}^{-2}$ the residual of the $k-\omega$ equation is less than ${10}^{-3}$. The average results of the last 1,000 steps are calculated as the final result.

FIGURE 2

FIGURE 2. Structure and flow area of “Impella 5.0”. The “Shroud” can protect the device from working properly, the high-speed rotor “Impeller” that provides additional energy for the heart ejection. The “Diffuser” has the function of buffering and rectifying and the micro “Motor” provide energy for device. In addition, the computational area is divided into rectification zone, rotor zone and tailrace zone.

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:

$\begin{array}{c}\mathrm{s}\mathrm{k}\mathrm{e}\mathrm{w}\mathrm{n}\mathrm{e}\mathrm{s}\mathrm{s}=\max \left[\frac{{\theta }_{\mathit{max}}-{\theta }_e}{180-{\theta }_e},\frac{{\theta }_e-{\theta }_{\mathit{min}}}{{\theta }_e}\right]\end{array}(17)$

where ${\theta }_{\mathit{max}}$, ${\theta }_{\mathit{min}}$ are the largest and smallest angle in the face or cell, ${\theta }_e$ 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

FIGURE 3. CFD mesh quality and mesh independence assessment. (A) Schematic diagram of CFD calculation area size $\left(D=6mm\right)$; (B) A total of 8.29 million meshes were used in the calculation domain. (C) Verification of mesh independence using three conditions and four mesh numbers. The horizontal axis represents the number of meshes, and the vertical axis represents the pressure difference between the inlet and outlet. (D) Mesh quality assessment of the 8.29 million meshes. The horizontal axis represents the skewness, and the vertical axis represents the volume percentage of different mesh quality.

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):

$\begin{array}{c}\tau ={\left(\frac{1}{6}\sum {\left({\tau }_{ii}-{\tau }_{jj}\right)}^2+\sum {\tau }_{ij}^2\right)}^{\frac{1}{2}}\end{array}(18)$

where $\tau$ is the scalarized shear stress and ${\tau }_{ij}$ is the shear stress tensor in the Cartesian coordinate system $\left(i,j=\mathrm{1,2,3}\right)$. The blood shear thinning model adopts the Carreau-Yasuda non-newtonian fluid model (Gijsen et al., 1999), which can be expressed as:

$\begin{array}{c}\frac{\eta -{\eta }_{\infty }}{{\eta }_0-{\eta }_{\infty }}={\left(1+{\left(\lambda \dot{\gamma }\right)}^a\right)}^{\frac{n-1}{a}}\end{array}(19)$

$\eta$ is the dynamic viscosity of blood, $\dot{\gamma }$ is the shear rate of the flow domain, other parameters are constants: ${\eta }_0=0.022Pa\bullet s,{\eta }_{\infty }=0.0022Pa\bullet s,\lambda =0.11s,a=0.644,n=0.392$. 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 $Y-Z$ 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 $10m/s$. Figure 4C shows the $\mathrm{\tau }$ contour of $Y-Z$ section. Areas with high shear stress $\left(\tau >100Pa\right)$ are concentrated on the gap, and the $\mathrm{\tau }$ of most rotor areas is $30Pa\sim 50Pa$.

FIGURE 4

FIGURE 4. (A) Velocity contour and two-dimensional streamline of $Y-Z$ section. There are two obvious eddies in the section; (B) Three-dimensional streamline distribution in rotor area and tailrace area. The maximum speed is at the impeller tip clearance position; (C) The $\mathrm{\tau }$ contour of $Y-Z$ section shows that the shear stress of most rotor areas is $30Pa\sim 50Pa$.

The $\mathrm{\tau }$ of the blood flow is highly positively correlated with IH. We intercepted the $\mathrm{\tau }$ 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 $\mathrm{\tau }$ 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

FIGURE 5. $\mathrm{\tau }$ contour of three sections in Y direction. The velocity gradient in the clearance area between the impeller and shroud is very large and prone to hemolysis. With the development of blood flow, the area of high shear stress gradually diffused.

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 $\mathrm{\tau }$ 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

FIGURE 6. (A) Distribution of four trace lines in the rotor zone; (B) The relationship curve between $\mathrm{\tau }$ and time on the four trace lines. The red and purple traces pass through the clearance, and the $\mathrm{\tau }$ changes sharply. The green trace does not pass through the area prone to hemolysis, the peak $\mathrm{\tau }$ is small.

The relationship curve between $\mathrm{\tau }$ 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 $\mathrm{\tau }$ 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 $\mathrm{\tau }$ 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 ($\mathrm{\tau }$) to reflect hemolysis and RBC’s damage. This section used the plane Couette flow formed by the movement of the boundary to characterize the $\mathrm{\tau }$ 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:

$\begin{array}{c}\frac{\sigma }{d}={\left(\frac{\upsilon }{d\ast u}\right)}^{\frac{3}{4}}=R{e}^{-\frac{3}{4}}\end{array}(20)$

FIGURE 7

FIGURE 7. (A) DPD calculation setting and RBC shear flow test; (B–E) The movement morphology of RBCs on the four traces. Under low shear stress, the RBCs showed rolling and folding shapes. In high shear stress, almost all RBCs will evolve into tank-treading motion form and present an ellipsoidal shape.

$\sigma$ Is the Kolmogorov scale, the characteristic scale $d=6mm$ is set as the pipe diameter, $\upsilon =3.5\ast {10}^{-6}{m}^2/s$ is the kinematic viscosity of blood. When the flow rate $Q=4L/\mathit{min}$, $u=2.2m/s$ is the characteristic velocity. From Eq. 20, we can calculate $\sigma \approx 12.5\mu m$.

Figures 7B–E shows the deformation state of RBCs under different $\mathrm{\tau }$. The characteristics of the red trace (Figure 7B) and the purple trace (Figure 7D) are similar, both have a great change rate of $\mathrm{\tau }$. The $\mathrm{\tau }$ increases from 10 Pa to 150 Pa, taking only ${10}^{-4}s$. Although the peak $\mathrm{\tau }$ 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 $\mathrm{\tau }$. 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 $\mathrm{\tau }$ and exposure time, the shear stress change rate $\dot{\mathrm{\tau }}$ may also be one of the sensitive factors of hemolysis.

To further verify the effect of shear stress change rate ($\dot{\mathrm{\tau }}$) 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 ($D_R$) 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 $\mathrm{\tau }$ 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 $\dot{\mathrm{\tau }}$ 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 $\dot{\mathrm{\tau }}$ is one of the most important factors affecting hemolysis. Moreover, there are more regions of high $\dot{\mathrm{\tau }}$ in the trace passing through the impeller clearance frequently, resulting in more serious hemolysis. Further, we take $\dot{\mathrm{\tau }}$ as the main variable to construct a dimensionless number ($D_k$) related to RBC damage, and its expression is as follows:

$\begin{array}{c}{D}_k=\int \left(\frac{\dot{\tau }}{\rho v^2}\bullet {\varphi }_1\bullet {\varphi }_2\right)dt\end{array}(21)$

FIGURE 8

FIGURE 8. (A–C) Use DPD to evaluate the RBC membrane surface damage rate $D_R$ on red (597 data points), purple (680 data points) and blue traces (847 data points). These three curves are characterized by stepwise growth. (D–F) Use the time integral dimensionless number $D_k$ to evaluate the RBC damage amount on the three trace lines. After comparison, it can be found that, although there are some errors in the peak value, the time point of RBC injury is coincident.

The dimensionless number of RBC damage $D_k$ should be expressed as an integral relationship with time, $\rho v^2$ 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 $\dot{\mathrm{\tau }}$ on blood damage, we add two limiters ${\varphi }_1,{\varphi }_2$ to the dimensionless integration formula. ${\varphi }_1$ indicates that only flow acceleration can injure RBCs, while flow deceleration does not affect RBC damage; ${\varphi }_2$ represents that the shear stress flow above 80 Pa can affect the damage amount of RBCs (Xu et al., 2022).

$\begin{array}{c}{\varphi }_1=\left\{\begin{array}{c}1\left(\dot{\tau }>0\right)\\ 0\left(\dot{\tau }<0\right)\end{array}\right.\end{array}(22)$

$\begin{array}{c}{\varphi }_2=\left\{\begin{array}{c}1\left(\tau >80Pa\right)\\ 0\left(\tau <80Pa\right)\end{array}\right.\end{array}(23)$

Figures 8D–F shows the relationship between the dimensionless number $D_k$ and time, which can directly evaluate the blood injury status through flow field information. Comparing the $D_k-t$ curve obtained from the flow field information with the $D_R-t$ 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 $\dot{\mathrm{\tau }}$ 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 $\mathrm{\tau }$ 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 $\mathrm{\tau }$ distribution contour in the rotor clearance (Figure 9B). The $\mathrm{\tau }$ at the contact position with the impeller is higher than 100 Pa, and the average $\mathrm{\tau }$ in the clearance is between 60 Pa and 70 Pa. Figure 9C shows the $\mathrm{\tau }$-t curve of the monitoring area in three rotating cycles, lasting for $6*{10}^{-3}s$.

FIGURE 9

FIGURE 9. (A) Magenta curve represents the $\mathrm{\tau }$ monitoring area during rotor one cycle rotation; (B) $\mathrm{\tau }$ distribution contour of “Impella 5.0″clearance. The $\mathrm{\tau }$ of the flow field in contact with the impeller is as high as 100 Pa. (C) Monitor the $\mathrm{\tau }$-t curve of the specific location within three rotation cycles of the rotor.

Next, we will use tDPD to analyze the RBCs’ population morphology and 30% hematocrit whole blood hemolysis in the $100\mu m\times 100\mu m\times 5\mu m$ 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 $∆p=80\mathrm{m}\mathrm{m}\mathrm{H}\mathrm{g}$ 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

FIGURE 10. Schematic diagram of tDPD whole blood hemolysis simulation. (A) The whole blood flow in the rectification area is regarded as the Poiseuille flow driven by blood pressure. (B) The flow in the rotor area is regarded as Couette flow. Then the velocity gradient generated by the movable wall can simulate the shear stress generated by the impeller in the clearance.

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):

$\begin{array}{c}IH=\left(1-H_{ct}\right)\frac{H{b}_P-H{b}_p^0}{H{b}_B}\approx \left(1-H_{ct}\right)\frac{H{b}_P}{H{b}_B}\end{array}(24)$

where $H_{ct}$ is hematocrit, $H{b}_B$ is the total hemoglobin, $H{b}_P$ is the hemoglobin in plasma, and $H{b}_p^0$ is the hemoglobin in healthy whole blood plasma. In healthy blood, no hemoglobin is free in plasma $\left(H{b}_p^0=0\right)$. 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

FIGURE 11. Simulation results of three shear stress peaks in the shroud clearance flow domain. (A) RBC’s deformation and hemoglobin distribution at 0.00082s; (B) RBC’s deformation and hemoglobin distribution at 0.00282s; (C) RBC’s deformation and hemoglobin distribution at 0.00482s.

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 $D_k$ 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 ($\mathrm{\tau }$) distribution in the “Impella 5.0” rotor region. The velocity gradient and $\mathrm{\tau }$ 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 $\dot{\mathrm{\tau }}$ may be one of the most important factors for RBC damage and hemolysis. Therefore, we proposed a dimensionless number $D_k$ of hemolysis based on $\dot{\mathrm{\tau }}$ in the Lagrange perspective. After the comparison, it was found that the result $D_R$ of RBC damage simulation was consistent with the result $D_k$ 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 $\dot{\mathrm{\tau }}$ to reduce hemolysis symptoms effectively.

Finally, we continuously observed the $\mathrm{\tau }$ 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, $\dot{\mathrm{\tau }}$ 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.

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

Antiga, L., and Steinman, D. A. (2009). Rethinking turbulence in blood. Biorheology 46 (2), 77–81. doi:10.3233/bir-2009-0538

PubMed Abstract | CrossRef Full Text | Google Scholar

Arvand, A., Hormes, M., and Reul, H. (2005). A validated computational fluid dynamics model to estimate hemolysis in a rotary blood pump. Artif. Organs 29 (7), 531–540. doi:10.1111/j.1525-1594.2005.29089.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Basil, S., Bryan, S., and Clyde, W. (2017). Heart failure and sudden cardiac death. Card. Electrophysiol. Clin. 9, 709–723. doi:10.1016/j.ccep.2017.07.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Bludszuweit, C. (1995). Three-dimensional numerical prediction of stress loading of blood particles in a centrifugal pump. Artif. Organs 19 (7), 590–596. doi:10.1111/j.1525-1594.1995.tb02386.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Dao, M., Lim, C. T., and Suresh, S. (2003). Mechanics of the human red blood cell deformed by optical tweezers. J. Mech. Phys. Solids 51 (11-12), 2259–2280. doi:10.1016/j.jmps.2003.09.019

CrossRef Full Text | Google Scholar

De Tullio, M. D., Cristallo, A., Balaras, E., and Verzicco, R. (2009). Direct numerical simulation of the pulsatile flow through an aortic bileaflet mechanical heart valve. J. Fluid Mech. 622, 259–290. doi:10.1017/s0022112008005156

CrossRef Full Text | Google Scholar

Ellis, J. T., Wick, T. M., and Yoganathan, A. P. (1998). Prosthesis-induced hemolysis: Mechanisms and quantification of shear stress. J. Heart Valve Dis. 7 (4), 376–386.

PubMed Abstract | Google Scholar

Espanol, P., and Warren, P. (1995). Statistical mechanics of dissipative particle dynamics. Europhys. Lett. 30, 191–196. doi:10.1209/0295-5075/30/4/001

CrossRef Full Text | Google Scholar

Faghih, M. M., and Sharp, M. K. (2020). Deformation of human red blood cells in extensional flow through a hyperbolic contraction. Biomechanics Model. Mechanobiol. 19 (1), 251–261. doi:10.1007/s10237-019-01208-3

CrossRef Full Text | Google Scholar

Fan, X., Phan-Thien, N., Chen, S., Wu, X., and Yong Ng, T. (2006). Simulating flow of DNA suspension using dissipative particle dynamics. Phys. Fluids 18 (6), 063102. doi:10.1063/1.2206595

CrossRef Full Text | Google Scholar

Fedosov, D. A., Caswell, B., and Karniadakis, G. E. (2010). Systematic coarse-graining of spectrin-level red blood cell models. Comput. Methods Appl. Mech. Eng. 199 (29-32), 1937–1948. doi:10.1016/j.cma.2010.02.001

CrossRef Full Text | Google Scholar

Fedosov, D. A., Dao, M., Karniadakis, G. E., and Suresh, S. (2014). Computational biorheology of human blood flow in health and disease. Ann. Biomed. Eng. 42 (2), 368–387. doi:10.1007/s10439-013-0922-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Fedosov, D. A., Pan, W., Caswell, B., Gompper, G., and Karniadakis, G. E. (2011). Predicting human blood viscosity in silico. Proc. Natl. Acad. Sci. 108 (29), 11772–11777. doi:10.1073/pnas.1101210108

PubMed Abstract | CrossRef Full Text | Google Scholar

Feng, R., Xenos, M., Girdhar, G., Kang, W., Davenport, J. W., Deng, Y., 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. doi:10.1007/s10237-011-0297-z

CrossRef Full Text | Google Scholar

Ge, L., Dasi, L. P., Sotiropoulos, F., and Yoganathan, A. P. (2008). Characterization of hemodynamic forces induced by mechanical heart valves: Reynolds vs. Viscous stresses. Ann. Biomed. Eng. 36 (2), 276–297. doi:10.1007/s10439-007-9411-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Giersiepen, M., Wurzinger, L. J., Opitz, R., and Reul, H. (1990). Estimation of shear stress-related blood damage in heart-valve prostheses - invitro comparison of 25 aortic valves. Int. J. Artif. Organs 13 (5), 300–306. doi:10.1177/039139889001300507

PubMed Abstract | CrossRef Full Text | Google Scholar

Gijsen, F. J. H., van de Vosse, F. N., and Janssen, J. 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. Biomechanics 32 (6), 601–608. doi:10.1016/s0021-9290(99)00015-9

CrossRef Full Text | Google Scholar

Grigioni, M., Daniele, C., Morbiducci, U., D'Avenio, G., Di Benedetto, G., and Barbaro, V. (2004). The power-law mathematical model for blood damage prediction: Analytical developments and physical inconsistencies. Artif. Organs 28 (5), 467–475. doi:10.1111/j.1525-1594.2004.00015.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Grigioni, M., Morbiducci, U., D’Avenio, G., Benedetto, G. D., and Gaudio, C. D. (2005). A novel formulation for blood trauma prediction by a modified power-law mathematical model. Biomechanics Model. Mechanobiol. 4 (4), 249–260. doi:10.1007/s10237-005-0005-y

CrossRef Full Text | Google Scholar

Groot, R. D., and Rabone, K. L. (2001). Mesoscopic simulation of cell membrane damage, morphology change and rupture by nonionic surfactants. Biophysical J. 81 (2), 725–736. doi:10.1016/s0006-3495(01)75737-2

CrossRef Full Text | Google Scholar

Groot, R. D., and Warren, P. B. (1997). Dissipative particle dynamics: Bridging the gap between atomistic and mesoscopic simulation. J. Chem. Phys. 107 (11), 4423–4435. doi:10.1063/1.474784

CrossRef Full Text | Google Scholar

Hoogerbrugge, P. J., and Koelman, J. (1992). Simulating microscopic hydrodynamic phenomena with dissipative particle dynamics. Europhys. Lett. 19, 155–160. doi:10.1209/0295-5075/19/3/001

CrossRef Full Text | Google Scholar

Koshiyama, K., and Wada, S. (2011). Molecular dynamics simulations of pore formation dynamics during the rupture process of a phospholipid bilayer caused by high-speed equibiaxial stretching. J. Biomechanics 44 (11), 2053–2058. doi:10.1016/j.jbiomech.2011.05.014

CrossRef Full Text | Google Scholar

Lanotte, L., Mauer, J., Mendez, S., Fedosov, D. A., Fromental, J-M., Claveria, V., et al. (2016). Red cells’ dynamic morphologies govern blood shear thinning under microcirculatory flow conditions. Proc. Natl. Acad. Sci. 113 (47), 13289–13294. doi:10.1073/pnas.1608074113

PubMed Abstract | CrossRef Full Text | Google Scholar

Leverett, L. B., Hellums, J. D., Alfrey, C. P., and Lynch, E. C. (1972). Red blood cell damage by shear stress. Biophysical J. 12 (3), 257–273. doi:10.1016/s0006-3495(72)86085-5

CrossRef Full Text | Google Scholar

Li, J., Lykotrafitis, G., Dao, M., and Suresh, S. (2007). Cytoskeletal dynamics of human erythrocyte. Proc. Natl. Acad. Sci. 104 (12), 4937–4942. doi:10.1073/pnas.0700257104

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Z., Yazdani, A., Tartakovsky, A., and Karniadakis, G. E. (2015). Transport dissipative particle dynamics model for mesoscopic advection-diffusion-reaction problems. J. Chem. Phys. 143 (1), 014101. doi:10.1063/1.4923254

PubMed Abstract | CrossRef Full Text | Google Scholar

Maraj, R., Jacobs, L. E., Ioli, A., and Kotler, M. N. (1998). Evaluation of hemolysis in patients with prosthetic heart valves. Clin. Cardiol. 21 (6), 387–392. doi:10.1002/clc.4960210604

PubMed Abstract | CrossRef Full Text | Google Scholar

Mehra, M. R., Uriel, N., Naka, Y., Cleveland, J. C., Yuzefpolskaya, M., Salerno, C. T., et al. (2019). A fully magnetically levitated left ventricular assist device — final report. N. Engl. J. Med. 380 (17), 1618–1627. doi:10.1056/nejmoa1900486

PubMed Abstract | CrossRef Full Text | Google Scholar

Olsen, D. B. (2000). The history of continuous-flow blood pumps. Artif. Organs 24 (6), 401–404. doi:10.1046/j.1525-1594.2000.06652.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Pahuja, M., Hernandez-Montfort, J., Whitehead, E. H., Kawabori, M., and Kapur, N. 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. Devices 19, 1–10. doi:10.1080/17434440.2022.2015323

PubMed Abstract | CrossRef Full Text | Google Scholar

Pan, W., Caswell, B., and Karniadakis, G. E. (2010). A low-dimensional model for the red blood cell. Soft Matter 6 (18), 4366. doi:10.1039/c0sm00183j

CrossRef Full Text | Google Scholar

Sohrabi, S., and Liu, Y. (2017). A cellular model of shear-induced hemolysis. Artif. Organs 41 (9), E80–E91. doi:10.1111/aor.12832

PubMed Abstract | CrossRef Full Text | Google Scholar

Song, X. W., Throckmorton, A. L., Wood, H. G., Antaki, J. F., and Olsen, D. B. (2003). Computational fluid dynamics prediction of blood damage in a centrifugal pump. Artif. Organs 27 (10), 938–941. doi:10.1046/j.1525-1594.2003.00026.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, X., Zhang, X. W., Hao, P. F., and He, F. (2019). Comparison of three control strategies for axial blood pump. J. Mech. Med. Biol. 19 (6), 1950058. doi:10.1142/s0219519419500581

CrossRef Full Text | Google Scholar

Xu, Z., Wang, C., He, F., Hao, P., and Zhang, X. (2023). Coarse-grained model of whole blood hemolysisand morphological analysis of erythrocytepopulation under non-physiological shear stressflow environment. Phys. Fluids 35 (3), 031901. doi:10.1063/5.0137517

CrossRef Full Text | Google Scholar

Xu, Z., Wang, C., Xue, S., He, F., Hao, P., and Zhang, X. (2022). The erythrocyte destruction mechanism in non-physiological shear mechanical hemolysis. Phys. Fluids 34 (11), 111901. doi:10.1063/5.0112967

CrossRef Full Text | Google Scholar

Yano, T., Sekine, K., Mitoh, A., Mitamura, Y., Okamoto, E., Kim, D. W., et al. (2003). An estimation method of hemolysis within an axial flow blood pump by computational fluid dynamics analysis. Artif. Organs 27 (10), 920–925. doi:10.1046/j.1525-1594.2003.00034.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Yazdani, A., Deng, Y., Li, H., Javadi, E., Li, Z., Jamali, S., 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. Interface 18 (175), 20200834. doi:10.1098/rsif.2020.0834

PubMed Abstract | CrossRef Full Text | Google Scholar

Yazdani, A., and Karniadakis, G. E. (2016). Sub-cellular modeling of platelet transport in blood flow through microchannels with constriction. Soft Matter 12 (19), 4339–4351. doi:10.1039/c6sm00154h

PubMed Abstract | CrossRef Full Text | Google Scholar

Ye, T., Shi, H., Phan-Thien, N., and Lim, C. T. (2020). The Key events of thrombus formation: Platelet adhesion and aggregation. Biomechanics Model. Mechanobiol. 19 (3), 943–955. doi:10.1007/s10237-019-01262-x

CrossRef Full Text | Google Scholar

Zhang, T., Taskin, M. E., Fang, H-B., Pampori, A., Jarvik, R., Griffith, B. P., et al. (2011). Study of flow-induced hemolysis using novel Couette-type blood-shearing devices. Artif. Organs 35 (12), 1180–1186. doi:10.1111/j.1525-1594.2011.01243.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, R., Avsievich, T., Popov, A., and Meglinski, I. (2020). Optical tweezers in studies of red blood cells. Cells 9 (3), 545. doi:10.3390/cells9030545

PubMed Abstract | CrossRef Full Text | Google Scholar