Experimental and simulation study on shear stress-induced erythrocyte damage based on vortex oscillator

Background: Shear stress-induced erythrocyte damage, namely hemolysis, is an important problem in the development of blood-contacting medical devices such as mechanical circulatory support devices. Computational fluid dynamics simulation combined with hemolysis prediction models have been widely used to predict hemolysis. With the development of hemolysis prediction models, the new hemolysis prediction model requires more experimental data to verify. In addition, the difference of in vitro blood-shearing device also affect the accuracy of hemolysis prediction. Methods: To address these problems, a new in vitro blood-shearing device (vortex oscillator) was used to further verify the accuracy of the hemolysis prediction models, and to guide the optimal design of blood-contacting medical devices such as mechanical circulatory support devices. Firstly, the flow field information such as wall stress and velocity of the vortex oscillator under different speeds was analyzed. Secondly, different hemolysis prediction models were used to calculate hemolysis, and the predicted data was compared with the experimental data. Results and Conclusion: In this study, the flow field information inside the vortex oscillator at high rotational speeds was systematically investigated, and the prediction of hemolysis was carried out. The results showed that the predicted data of hemolysis was significantly different from the experimental data, which indicated that it was urgent to establish a standardized in vitro blood-shearing platform to provide a reference for accurate hemolysis prediction.


Background
Super-physiological shear stress causes an accumulation of damage to erythrocytes that can result in their destruction or lead to impaired deformability [1,2].Such mechanical stress frequently occurs inside mechanical circulatory support (MCS) devices, including blood pumps, percutaneous devices, and extracorporeal membrane oxygenators, which are used in patients with severe heart or lung failure.In addition, super-physiological shear stress causes complications such as impaired oxygen transport due to hemolysis [3][4][5][6].Therefore, to ensure the activity of erythrocytes, the influence of stress in the MCS devices must be considered.In the physiological flow field, the average stress on the blood vessel wall is about 1 Pa, but in the MCS devices, the shear stress caused by complex geometric structure and mechanical movement is often much higher than that in the physical flow field, and the stress in some areas even exceeds 1,000 Pa [7].If erythrocytes are subjected to a super-physiological flow field for a long time, severe hemolysis will occur.
To explore the law of hemolysis caused by the shear stress in the MCS devices, researchers have built several in vitro blood-shearing devices for experimental studies [8,9].Most of these researches have studied erythrocyte damage in the laminar flow state, while the flow within the MCS devices is a complex turbulent state.Therefore, to study the law of hemolysis under complex flow and provide more accurate and reliable information for hemolysis prediction.It is urgent to build a new in vitro blood-shearing device to better simulate the real conditions in the MCS devices.
In the development of the MCS devices, computational fluid dynamics (CFD) simulation combined with hemolysis prediction models have been widely used to predict hemolysis and assist the optimal design of the MCS devices [10,11].At present, there are three main types of hemolysis prediction models: Stress-based, Strain-based, and Energy-dissipation-based.The mainstream model is based on stress, whose prediction error is large.In recent years, researchers have found that the model based on energy dissipation can significantly improve hemolysis prediction in turbulent flow fields.To evaluate the accuracy of CFD flow field prediction and hemolysis prediction, Food and Drug Administration designed two benchmark models: "Medical Device Flow Model" and "Benchmark Nozzle Model", and carried out PIV flow field measurement and hemolysis experiment [12,13].52 teams submitted CFD calculation results, and there are great differences between the teams.Some hemolysis calculation results are even several orders of magnitude different from the experimental date [14].This is because the accuracy of flow field information (such as velocity, wall stress, etc.) obtained by CFD needs to be verified, and the errors during in vitro hemolysis experiments also directly affect the accuracy of hemolysis prediction models.
In this study, a new in vitro blood-shearing device (vortex oscillator) was adopted to further verify the accuracy of the hemolysis prediction model based on energy dissipation, which provides a new idea for the verification of hemolysis prediction models.It has been shown that the flow field of the vortex oscillator has an important effect on the biological properties.Salek et al. explored flow information in vortex oscillator, CFD was used to simulate flow and the distribution of wall shear stress in six-well culture plates [15].The flow pattern and wall shear stress of two rotational speeds (100 and 200 rpm) and two fluid volumes (2 and 4 mL) were investigated.The results showed that the distribution of wall shear stress was related to the topological structure of the fluid.The study is useful for interpreting the results of biological experiments.However, although the flow field information in this study is sufficient, the flow field at high rotational speed was not studied.Hemohydrodynamic factors are crucial to hemolysis, so the accurate acquisition of flow field information is a prerequisite to verify hemolysis prediction models.As a widely used instrument in the biomedical field, the vortex oscillator can explore its more complex flow field properties at high rotational speeds and be used to verify hemolysis prediction models.In this study, the flow field of the vortex oscillator at high speed is analyzed.The focus of this study is to analyze the quantitative relationship between hemolysis and non-physiological shear stress under different flow mechanisms to provide more accurate information for the development of hemolysis prediction models and helpful for the optimal design of the MCS devices.

Test blood samples
Porcine whole blood was collected and mixed 1:20 with 3.24 wt% sodium citrate (Yuanye Bio-Technology, Shanghai, China) from a healthy pig.This anticoagulated whole blood was stored at a temperature of 0-10 °C and transported to the laboratory 2 h after blood collection (n = 3).

Experimental design
To avoid the undesired additional damage of blood, the inner surface of the vortex oscillator was wetted with saline before the blood-shearing experiment.Then, each vortex oscillator was filled with 4 mL of blood.The flow-dependent parameters shear stress and exposure time were controlled by changing the rotational speed and the running time.The rotational speeds in this study were set to 500 rpm, 1,000 rpm, and 3,000 rpm.The running time at each rotational speed was set to 60 min.A summary of the specific working condition information of the vortex oscillator is shown in Table 1.Sheared blood samples under each parameter condition were collected.Each sample was then centrifuged at 1,500 g for 10 min.After centrifugation, the supernatant was taken out and subjected to hemolysis measurement.

Measurement of shear-induced hemolysis
A Harboe assay was performed to measure pfHb levels [16].Briefly, 15 µL of plasma was diluted with 1,485 µL of 0.1% sodium carbonate as a measurement sample [17].Absorbances at 380, 415, and 450 nm were measured using a microplate reader (Molecular Devices, CA, USA).The amount of pfHb was calculated using the following formula and the acquired absorbance: Here A415, A450, and A380 represent the absorbance at each wavelength while Vtotal (µL) and Vplasma (µL) shows the total volume of the measurement sample and the plasma volume, respectively.

Hemolysis prediction models
At present, the mainstream hemolysis prediction model is based on stress, and it is believed that hemolysis has an exponential relationship with the equivalent shear stress and exposure time [18]: Here represents hemolysis index; represents the equivalent shear stress, which can be calculated from the stress tensor; Hb represents the total concentration of hemoglobin; hb indicates an increase in free hemoglobin in the plasma; A, and an empirical constant, which was generally obtained by fitting the hemolysis experimental data.The empirical constants of the three hemolysis prediction models used in this study are listed in Table 2 [19][20][21].
It has been shown that Reynolds stress has a large error in predicting the hemolysis in turbulent flow and overestimates the hemolysis.The most significant characteristics of turbulence are energy dissipation.Erythrocyte damage can essentially be regarded as a part of the energy dissipation in the flow field, and some of the energy dissipated in the flow field acts on the erythrocytes, causing erythrocyte damage.Therefore, it is more suitable to estimate hemolysis based on energy dissipation than stress.The form of the hemolysis prediction model based on the energy dissipation model is as follows: Here is the total energy dissipation, which can be calculated from the stress tensor, the dynamic viscosity of blood, and the density of blood.
The empirical constants of three models (Heuser model, Giersipen model, and Zhang model) were used in this study, each of which has two forms of strain-based and energy-dissipation-based models.

Computational fluid dynamic simulations
The vortex oscillator used in this study rotates around an axis 4 mm parallel to the central axis, where the blood volume accounts for 50%, and the sustained non-physiological shear stress generated by the rotation causes erythrocyte damage.The dimensions of the vortex oscillator are shown in Figure 1A.To visualize flow patterns and shear stress distribution, a CFD solver, Fluent v6.3.26 (Ansys Fluent, Inc., Lebanon, NH, USA) was used to calculate the flow field information.The grid was meshed using Ansys ICEM.The mesh quality was above 0.3, and y+ was below 1.5.The mesh generation effect diagram is shown in Figure 1B.Blood was specified to be incompressible Newtonian fluid with a dynamic viscosity of 3.5 mPa•s and density of 1,060 kg/m 3 in simulation.Because the model has a gas-liquid two-phase flow, a mixture model was adopted.The viscous model adopted the shear-stress transport model.Three different mesh sizes were fabricated by the tetrahedral mesh method to study the mesh independence.

Mesh independence analysis
Numerical calculation requires grid division.Obviously, the denser the number of mesh, the more accurate the calculation results.But what the crypto grid brings is a greater amount of computation.In practical applications, it is necessary to consider the relationship between computational accuracy and computational resources, and choose an appropriate mesh number to achieve a balance between them.Taking the condition of 3,000 rpm as an example, the relationship between the wall shear stress of the model and the  3.

Table 2 Empirical constants of the hemolysis prediction models
The instantaneous distribution of the wall shear stress of the vortex oscillator under different mesh sizes is shown in Figure 2. As shown in these figures, due to the eccentric rotation of the vortex oscillator, the wall shear stress on one side of the tube will be significantly higher than the other areas, that is, the red area in each figure.The wall shear stress cloud map of 180 thousand meshes was significantly different from that of the other two mesh sizes, which can be inferred that using a mesh size of 180 thousand for numerical calculation, the analysis of the flow field information may not be accurate enough.In addition, from the average value of wall shear stress on the entire wall surface shown in Table 3, it can be seen that the error of wall shear stress obtained by 2.4 million meshes and 1 million meshes was within 5%, which met the engineering error requirements, and it can be considered that the accuracy of calculation results of 2.4 million meshes and 1 million meshes is comparable.
Combining the above qualitative and quantitative analysis results, it can be concluded that to ensure that the number of mesh will not affect the flow field information, minimize the errors caused by spatial dispersion, and take into account the factors of computing resources, 2.4 million meshes were determined as the last meshes to be used.

Vortex oscillator flow field information
It can be seen from the Reynolds number in Table 1 that the change in rotational speed will bring about the change in Reynolds number, and the difference in Reynolds number between different rotational speeds was very large.These inevitably lead to differences in the flow of field information.The analysis results of flow field information at different rotational speeds are shown in Table 4.
 is the wall stress of 2.4 million mesh)  The instantaneous distribution of the wall shear stress of the vortex oscillator under different rotational speeds is shown in Figure 3.The difference in rotational speeds makes the difference between the wall shear stress greatly.With the increase in rotational speeds, the value of wall shear stress also increases.
The instantaneous cloud image of blood velocity in the vortex oscillator is shown in Figure 4.As the rotational speed increases, the vector arrows become more dense, demonstrating that the magnitude of the blood velocity increases with the vortex oscillator rotational speed.The three working conditions all show that the blood velocity speed was fast and the middle blood velocity was low.

Numerical calculation and analysis of Hemolysis
The setting of hemolysis numerical calculation is as follows： a)Mesh: According to the mesh independence analysis above, the number of meshes selected was 2.4 million.
b) Solver setup: The internal iteration time of 500 rpm, 1,500 rpm, and 3,000 rpm was 0.0012 s, 0.0004 s, and 0.0002 s, respectively, that is, the vortex oscillator rotated once every 100 steps calculated.c) Convergence standard: The convergence standard of the flow field in the three working conditions was unified to 1E-04.The convergence standard of the hemolysis calculation is that the error of the rate of change of the hemolysis date was within 1% every 100-time steps.When the flow field converges stably, the hemolysis numerical solution is performed.The user defined functions were used to calculate the hemolysis based on energy dissipation and strain models, respectively.After ensuring the convergence of hemolysis values, the data were collated for comparative analysis.Figure 5, Figure 6, and Figure 7 show the comparison between the CFD hemolysis predicted date and experimental date using the Heuser model, Giersipen model, and Zhang model at 500 rpm, 1,500 rpm, and 3,000 rpm, respectively.Judging from the results of the predicted date of hemolysis, both the predicted date of hemolysis based on energy dissipation method and those based on strain model are far from the experimental date, so it is impossible to strictly judge the merits of the predicted models of hemolysis based on energy dissipation method and those based on strain model.This may be because the rotation condition of the vortex oscillator is not stable, and the vibration in the rotation process causes the change of shear stress.However, the predicted date of the hemolysis prediction model based on the energy dissipation method is closer to the experimental data.In addition, the data of the same hemolysis prediction model based on energy dissipation and based on strain are also quite different.Among the three models, the predicted date of the Zhang model is the closest to the experimental date.These results indicated that it was urgent to establish a standardized in vitro blood-shearing platform to provide a reference for accurate hemolysis prediction.

Limitations and future directions
There were also some study limitations and future directions.In the flow field simulation part, there is still room for further encryption of the grid number, which can continue to reduce the error.There are some improvements in the setting of Fluent, such as the selection of the turbulence model and the improvement of the accuracy of the momentum equation, which can more accurately analyze the flow field.If possible, PIV can also be used to test the flow field and compare the simulated flow field.In the part of numerical calculation of hemolysis, the obtained result differs greatly from the experimental date.The rotation condition of the vortex oscillator is not stable, which may cause a change of shear stress due to vibration during rotation.The accuracy of the flow field is a very important factor which directly affects the numerical calculation of blood damage.Therefore, further simulation of the accurate flow field is the key to future research.

Conclusion
Three typical models of hemolysis prediction were verified based on the results of the hemolysis experiment with the vortex oscillator.The flow field of the vortex oscillator is studied by numerical simulation, and the flow field information, such as wall stress and velocity, is investigated and summarized.At the same time, it is used to verify the quality of the hemolysis prediction model based on the energy dissipation method and the stress-based hemolysis prediction model, as well as the comparison of different models and experimental data, to judge the accuracy of different hemolysis prediction models.Through numerical calculation, it is found that the predicted date of hemolysis is very different from the experimental date, and it is impossible to strictly judge the quality of the hemolysis prediction model based on the energy dissipation method and the strain model.It was urgent to establish a standardized in vitro blood-shearing platform to provide a reference for accurate hemolysis prediction.

Figure 1
Figure 1 The computational Domain.(A)The dimensions of the vortex oscillator.(B) The effect diagram of mesh.

Figure 2
Figure 2 Shear stress distribution of the vortex oscillator under different mesh sizes at 3,000 rpm.(A) 180 thousand meshes, (B) 1 million meshes, (C) 2.4 million meshes.

Figure 4
Figure 4 Cloud image of blood velocity in the vortex oscillator under different rotational speeds.(A) 500 rpm, (B) 1,000 rpm, (C) 3,000 rpm.

Figure 5 Figure 6 Figure 7
Figure 5 Comparison between the predicted date of different models and the experimental data at 500 rpm

Table 1 A summary of the specific working condition information of the vortex oscillator
where ρis the blood density (1,060 kg/m 3 ); µ is the dynamic viscosity (3.5 mPa•s); n is the rotational speed of the vortex oscillator; D is the model diameter (27 mm).Submit a manuscript: https://www.tmrjournals.com/bmec