Research on Human Erythrocyte's Threshold Free Energy for Hemolysis and Damage from Coupling Effect of Shear and Impact Based on Immersed Boundary-Lattice Boltzmann Method

Researches on the principle of human red blood cell's (RBC) injuring and judgment basis play an important role in decreasing the hemolysis in a blood pump. In the current study, the judgment of hemolysis in a blood pump study was through some experiment data and empirical formula. The paper forms a criterion of RBC's mechanical injury in the aspect of RBC's free energy. First, the paper introduces the nonlinear spring network model of RBC in the frame of immersed boundary-lattice Boltzmann method (IB-LBM). Then, the shape, free energy, and time needed for erythrocyte to be shorn in different shear flow and impacted in different impact flow are simulated. Combining existing research on RBC's threshold limit for hemolysis in shear and impact flow with this paper's, the RBC's free energy of the threshold limit for hemolysis is found to be 3.46 × 10−15 J. The threshold impact velocity of RBC for hemolysis is 8.68 m/s. The threshold value of RBC can be used for judgment of RBC's damage when the RBC is having a complicated flow of blood pumps such as coupling effect of shear and impact flow. According to the change law of RBC's free energy in the process of being shorn and impacted, this paper proposed a judging criterion for hemolysis when the RBC is under the coupling effect of shear and impact based on the increased free energy of RBC.


Introduction
A high-performance blood pump is of critical importance, because it can function as an artificial heart for a short or long time to maintain the body's blood circulation [1]. Normalized index of hemolysis (NIH), which is the height of free hemoglobin in the hematocrit of 100 L blood being pumped out by a blood pump in unit time, is the blood pump's most important performance index [2]. In natural blood, the hemoglobin is intracellular. However, the hemoglobin will be dissociative because the RBC will be damaged by shear, impact, turbulence, and other mechanical factors when it passes through the high-speed spiral flow field of the blood pump [3,4].
Experimental methods had been used in some papers for the research of RBC's damage in the blood pump; for example, a study of the RBC's threshold limit for hemolysis in shear flow by exposing the blood in a given shear rate of flow is proposed [5,6]. Based on these experimental data, some empirical formulas have been proposed for the calculation of the degree of RBC's injury when the RBC is exposed in different shear stresses and times [7][8][9]. The RBC's threshold limit for hemolysis in impact flow by dropping the blood from a given height to the glass plate was researched by Yun and Tan [10]. However, the experimental method cannot explain the criteria of RBC's injury.
The RBC was often been modeled as a capsule that consisted of a capsule membrane and cytoplasm. The Skalak model (SK), the zero-thickness (ZT) shell model, and the neo-Hookean (NH) model are widely used to describe the properties of the capsule membrane [11]. However, three models can only respond to elastically for all based on the elasticity. Moreover, the surface area and volume of capsule cannot be conserved when no restraints are imposed on in NH model and ZT model. In a recent research, the nonlinear spring model was used to simulate the capsule membrane's resisting area and volume change [12,13]. Considering the biological cell membrane's viscoelastic characteristics, the worm-like chain (WLC) model which is based on the spring model can perform better simulations of the membrane [14,15].
In the study of RBC's injury model, Yun [3] has built the model of RBC's mechanical damage and researched the principle of RBC's rupture because of being impacted, pressure, shear stress, and turbulence. There are also some studies on RBC's behavior in microchannel, tank-treading motion and large deformation, and so on. Omori et al. [16] numerically investigated the dynamics of an RBC in a thin microprobe by combining the immersed boundary method and the finite element method. Xu et al. [17] presented IB-LBM to research the RBC's deformation and aggregation in the blood flow. Li et al. [18] employed dissipative particle dynamics in the study on RBC's tank-treading motion and motion frequency.
The RBC will be damaged when the increased free energy of RBC can offset the energy needed for RBC to be broken. In order to form a criterion of RBC's mechanical injury from the aspect of RBC's free energy, this paper takes the following approaches. Firstly, a nonlinear WLC model of RBC was presented and so was the immersed boundary and lattice Boltzmann method (IB-LBM), based on which the motion of RBC in plasma was numerically simulated. Secondly, the free energy, shape, and other parameters of RBC in shear and impact flow are analyzed. Third, the free energy of RBC that changed in the coupling effect of different shear stresses and different impact velocities was researched. And then, by combining the existing research on the RBC's threshold limit for hemolysis in shear and impact flow with this paper's, we get the threshold of free energy for RBC's damage which can be used as criteria for RBC's damage. The method for judging RBC's damage in coupling effects of shear and impact was proposed too. The research of this paper can provide the study of hemolysis in the blood pump with basis and foundation.

Model and Method
2.1. RBC Model. The membrane of RBC was modeled as a series of nonlinear wormlike chains as shown in Figure 1. The length change of the chain, the angle change of adjacent chains, and the change of area surrounded by all the chains can represent any deformation of RBC. When the shape of RBC was deformed, the energy which is defined as free energy of RBC is stored in chains. Based on the free energy change, the governing equations of the RBC deformation are [19,20] E L , E B , and E S are stretched/compressed energy, bending energy and area energy, respectively. E is the sum of E L , E B , and E S . k L , k B , and k S are the chain constant for the change in length, bending angle change, and area change, respectively. The elastic chain force acting on the ith membrane node is F i . The total number of the chain elements is n. x i = l i /l i max , set to 2.2 for the model [14], l i is the length of the ith chain, and l i max is the maximum extension of the ith chain. θ i and θ i0 are the changed angle and reference angle between two chain elements, respectively. S 0 is the reference area of the RBC, and S is the changed area. The second term of Equation (4) expresses the repulsive force between the neighboring nodes, and K pi is calculated by using [21] Vicoelasticity is an important property of RBC; the chains' viscosity in the WLC model can be expressed by [14] F D ij is the force acting on the particle. γ T and γ C are the coefficients and can be calculated from membrane viscosity η m = 0:022 Pa•s. Adding Equation (6) to Equation (11), the membrane viscosity can be analyzed using an algorithm.
The parameters of the RBC model are shown in Table 1. The model of RBC was verified by comparing the simulation result of RBC being stretched and experiment results of RBC being stretched by optical tweezers in our previous work [19,20].

IB-LBM Method. The deformation of RBC in different
flow relates to fluid-structure interaction problem. IB-LBM is applied to solve the interaction between the membrane and flow. The membrane was discretized to Lagrangian points, and the flow field was discretized to Eulerian points. An external forcing term is needed in order to adopt an immersed boundary method in the LBM; in the paper, we choose explicit-type split-forcing single-relaxation-time  Figure 1: The simulation model of RBC in plasma. Two coordinate systems were used: Eulerian grid was used for the plasma and RBC cytoplasm, and Lagrangian grid was used for the RBC membrane. The RBC membrane was modeled as a series of particles (the number is 40) connected by a chain network which can resist compression/stretch and bending. 2 Applied Bionics and Biomechanics lattice Boltzmann equation with a forcing term [22]: Δt is the streaming time step, t is the time, and τ is the single relaxation parameter. f α is the density distribution function of the Eulerian point x along the α direction and f ðeqÞ α is its corresponding equilibrium state: The discrete force distribution function F α ðx, tÞ is The membrane of RBC, which is a physical boundary, was discretized to Lagrangian points. And then, the boundary force is transformed into the forcing term f [23]: F is the force term at the Lagrangian points X, and δðr − X ðs, tÞÞ is a delta function. The forcing term in the LBM is distributed to the Eulerian points by using a Dirac delta function [24]: where F ij and F b are the force at the Eulerian point ij and The delta function is The velocity interpolation for the nearby Eulerian points is used for the position update of the Lagrangian points according to Here, u b and u ij represent the velocities of the Lagrangian point b and the Eulerian point ij, respectively.
2.3. The RBC in Shear and Impact Flow. Two parallel plates moving in the same velocity u but in the opposite direction generated the shear flow, and the RBC was placed in the center of the flow field. In the Couette flow, when the distance between the two plates is δ (20 μm in the paper), the kinetic viscosity of the fluid is μ, the shear stress τ 1 is 2μu/δ. The range of shear stress is 0 to 800 Pa, and the range of the corresponding Reynolds number is 0 to 226. In shear flow, Ladd's momentum correction term [25] was used to solve the velocity boundary of the moving plates at the top and bottom of a fluid field: where α is the opposite direction of α, u w is the velocity of the boundary, and c s is the Lattice velocity, which can be calculated from The inlet and outlet were periodic boundary condition.
The blood cells will be impacted when flowing through the blood pump because of a high-speed rotating impeller. In impact flow, the RBC was impacted on the wall at certain velocity with fluid. The range of impact velocity is 0 to 10 m/s. Half-way bounce-back boundary was applied to solve the rebound of nodes of the membrane on the wall. Similar to shear flow, the inlet and outlet were periodic boundary condition.
The parameters of IB-LBM which is employed for solving the RBC in shear and impact flow are shown in Table 2.  Figure 2. The higher the shear stress is, the longer and narrower the RBC is, and the shape of RBC varies from biconcave to fusiform. The diameter of RBC in the long axis direction is about 13 μm when the shear stress is 255 Pa. When the stress is 800 Pa, the diameter can be up to 17.5 μm.
The change of RBC's free energy under different shear stresses is demonstrated in Figure 3. When the stress is larger than 150 Pa, the free energy increases obviously and has approximate linear relationship with shear stress. This case indicates that the membrane of RBC cannot resist the shear stress higher than 145 Pa. The hemolytic threshold value of RBC in shear flow is different according to past research. The most commonly used is 400 Pa [14]; others are 255 [13], 600 [7], and 800 Pa [5]; and the corresponding free energy added is 3:45 × 10 −14 , 2:04 × 10 −15 , 5:18 × 10 −15 , and 7:81 × 10 −15 J, respectively, in this paper. The hemolytic threshold value of free energy is wide-ranging. In order to make sure of the hemolytic threshold value of free energy, other references and research are needed. The threshold free energy value of RBC's damage by impact was used in this paper. Figure 4 illustrates the time needed for RBC's free energy to increase in different shear flow. The higher the shear stress is, the shorter the time spent is. When the shear stresses are 255 Pa, 400 Pa, 600 Pa, and 800 Pa, the corresponding time spent is 7:18 × 10 −1 , 6:51 × 10 −1 , 4:73 × 10 −1 , and 4:04 × 10 −1 s, respectively.
As shown in Figure 5, the free energy increased over time when the RBC was shorn in 255, 400, 600, and 800 Pa shear flow. The tendency of free energy varying over time in different shear stresses was similar.
3.2. The RBC in Impact Flow. The biggest deformation of RBC in different impact flow is displayed in Figure 6. When the velocity of impact flow is 6, 8, and 10 m/s, the axial diameter of RBC in the long axial direction can be up to 14, 16, and 20 μm, respectively. Not like the shape of RBC in shear flow, the deformation of RBC in impact flow seems to be compressed along the vertical direction of impact velocity.
The change of RBC's free energy, which comes to the biggest deformation in different impact velocities, is presented in Figure 7. As the impact velocity increases, the flatter RBC is squeezed and the free energy change is increased nearly linearly. When the velocity of impact flow is 10 m/s, the change of RBC's free energy is 4:135 × 10 −15 J.
The processes of RBC's free energy changing in three kinds of impact velocity are demonstrated in Figure 8. The higher the impact velocity was, the bigger the increase of RBC's free energy was, and the time needed was shorter. Compared to Figure 5, the time needed in impact flow is much smaller than in shear flow, about an order of magnitude.

The RBC under Coupling Effect of Shear and Impact.
Shear stress and high velocity simultaneously exist in the blood flow of the blood pump. The coupling effect of shear and impact should be analyzed. It is hard to simulate the RBC being shorn and impacted simultaneously. The  Figure 2: The shapes of RBC in different shear stresses (255 Pa, 400 Pa, 600 Pa, and 800 Pa) and the natural shape of RBC. The lager the shear stress was, the longer and narrower the RBC stretch was, and the shape of RBC varies from biconcave to fusiform.
The grid in Figure 9 is the result of Equation (16). The coupling effect of shear and impact was smaller than the addition of two effects, but larger than anyone alone.

Discussion
The RBC will be broken when the free energy of RBC is bigger than the threshold free energy of RBC for hemolysis. There are different threshold shear stresses for hemolysis in past researches for the different accuracy of research methods and expose time. In impact flow, the expose time is meaningless for the process of RBC being impacted must be complete. According to Yun and Tan [10], the threshold limit for hemolysis in impact flow is between 8 m/s and 10 m/s. According to fitting and interpolation of the data in Figures 3 and 7 In the research of Lu et al. [5], we know that when the shear stress is 800 Pa, the exposure time is 1 ms, and the RBC will be damaged. As shown in Figure 5, when the shear stress is 800 Pa, the exposure time is 2 ms, and the free energy can reach to E 0 . Therefore, the research is reasonable for the result of this paper comes close to the experiment result [5].
Equation (16) and E 0 can be used for judgment of RBC being damaged in the coupling effect of shear and impact, and the RBC is damaged if the value of F e is larger than E 0 .

Conclusions
Conclusions should clearly explain the main findings and implications of the work, highlighting its importance and relevance. The work of this paper can be summarized as follows: (1) The free energy, shape, and time needed for RBC when being shorn in different shear stresses are researched in the paper by using IB-LBM based on the WLC model The research of RBC's damage judgment is very important for the design and optimization of the blood pump. It is hard to judge the damage of RBC (hemolysis) in the blood pump in theoretical research, especially when the factors, which play a role on RBC, is not single, for example, the coupling effect of shear stress and impact. The research in this paper provides a way of judgment for RBC being damaged. Comparing the free energy of RBC in certain effects to the threshold free energy for hemolysis, if the free energy of RBC is larger, the RBC is damaged. The three-dimensional model will be applied in future work for more reasonable research on the mechanism of RBC being damaged.

Data Availability
The simulation result of RBC being stretched for verification of the RBC model is available at 10.3938/jkps.74.607 and is cited at relevant places within the text as [20]. The threshold limit of shear stress for hemolysis is available at 10.1016/S0021-9290(01)00084-7, 10.1007/s10237-005-0005y, 10.1002/fld.3939, and 10.1016/j.mvr.2011.05.006 and is cited at relevant places within the text as [5,7,13,14], respectively. The free energy of RBC in different shear stress, impact velocity, and coupling effects of shear and impact which is used to support the findings of this study is available from the corresponding author upon request.