Numerical study of shear-based hemolysis in aorta with left ventricular assistance device

Ventricular assistance devices (VADs) for heart failure treatment have been paid high attention among researchers for decades. However, the follow-up complications such as hemolysis and thrombosis require further optimization for this technique. Shear stress has been demonstrated to be significantly related to the hemolysis because of the rupture of red blood cells membrane with a leaking of hemoglobin in the plasma. This issue has already been investigated inside the pump of VAD, but estimations are still lacking regarding hemolysis generation in the aorta itself after VAD implantaion. Thus, the present study aims to evaluate the hemolysis in aorta through establishing the 3D numerical model of aorta with left ventricular assistance device (LVAD). Non-Newtonian Carreau model has been adopted. Comparisons of hemolysis evaluation have been made with two different mathematical models existing in literature. Moreover, the flow topology and hemodynamic variations have been studied. Different working conditions of LVAD have been considered corresponding to different heart failure severities. The results reveal a relatively low level of hemolysis risks in aorta. The thrombosis is more prone to occur in the case of severe heart failure condition. ⁎ Corresponding author: smaine.kouidri@ensam.eu


Introduction
Heart failure is a constantly increasing health problem to date affecting over 30 million people globally [1]. The heart transplantation as the best treatment is limited mainly by the lack of donor organs. Thus along with the advances in cardiac surgery, VADs have been developed and employed successfully for patients suffering from advanced heart failure since 1966 [2]. The theory of VADs is an implantable pump aiming to sufficiently restore the cardiac output through assisting a dysfunctional ventricle to pump blood to the distal organs as shown in Figure 1 [3]. Then, two sources of flow from ventricle and LVAD are generated to deliver blood into aorta. However, the follow-up adverse events such as thrombosis and hemolysis have been observed resulting in the VADs malfunction [4]. The relationship between local hemodynamics and VADs malfunction has caused significant interest among researchers as shear stress may lead to the rupture of red blood cells membrane, resulting in a leaking of hemoglobin in the plasma [5,6]. Numerical studies have been commonly adopted to investigate hemodynamic characteristics affected by VADs implantation. Different models considering this device have been established within literature in order to better understand the relevant hemodynamics in the specific physiological conditions [7][8][9]. The works carried out by A. Assmann et al. [10] have demonstrated the remarkable increase of aortic wall shear stress (WSS) during the pulsatile extracorporeal circulation compared to the nonpulsatile case. R. Mazzitelli et al. [11] found out that when LVAD is located in the ascending aorta, a better contribution can be obtained compared to the case where LVAD is located in the descending aorta resulting in more surfaces subjected to disturbed WSS. Y. Wang et al. [12] demonstrated that hemolysis index (HI) was higher in the cases of sinusoidal and pulsatile speed signal compared to constant speed signal. Moreover, among the existing literature within this topic, blood with Newtonian fluid model has been commonly considered regardless its shear thinning property [13,14]. In the present study, the 3D model of aorta with LVAD branch at the ascending location has been established referring to the patient-specific case. Blood flow with Carreau-model has been considered instead of Newtonian fluid model. As the previous investigations mainly focus on the blood cell damage inside pump of LVAD, the potential risks of hemolysis and intimal hyperplasia in aorta after LVAD implantation have been taken into account based on the hemodynamic evaluations. Different working conditions of LVAD have been considered corresponding to different heart failure severities. Aorta Fig. 1. Post-implantation of LVAD from ventricle to ascending aorta [3].

Geometrical model establishment
The simplified geometrical model of the aorta with LVAD branch at ascending location has been established regardless of the wall rugosity. Moreover, a second model without LVAD has been also built. Both geometrical models are visible in Figure 2. The constant mean diameters for aorta and its branches are adopted. The junction structure composed of LVAD and aorta can be observed along with three artery branches of brachiocephalic artery, left carotid artery and left subclavian artery at the ascending aorta. The corresponding diameters of these branches are 8 mm, 6 mm and 6mm respectively and the diameter of aorta is 30 mm. Regarding LVAD, the diameter is 1.6 cm with length of 7 cm so that fully developed flow can be achieved. The angle between centerlines of LVAD and aorta is 82.5 o to make the flow direction from LVAD more complying to the orientation of aorta geometry referring to literature [11].

Mathematical model and boundary conditions
The blood flow in aorta is considered to be incompressible with constant density of 1060 kg/m 3 . The dynamic viscosity of blood follows non-Newtonian Carreau-model as described below [15]: Where  is the shear rate, 1/s, The blood in aorta model is assumed to be laminar as illustrated in [16]. The corresponding fluid behavior is governed by the Navier-stokes equations: Where V is the flow velocity, m/s. p is the pressure, Pa.  is the flow density, kg/m3.  is the dynamic viscosity, Pa.s. The aortic inlet flow in the reference case without LVAD is shown in Figure 3-(a) representing the normal state of one person without need of artificial assistance [17]. For the cases with LVAD, a pulsatile flow rate and a continuous flow rate are employed separately at the inlets of aorta and LVAD. Different inlet flow rates of LVAD and aorta are considered corresponding to different heart failure conditions as shown in Figure 3-(b) and (c) [11]. The aortic flow rate in Figure 3-(b) indicates a relatively moderate condition of heart failure which requires intervention of LVAD with continuous flow rate of 1 l/min. Comparatively, the aortic flow shown in Figure 3-(c) represents the more severe condition of heart failure with increased requirement of flow rate from LVAD of 2 l/min. The inlet velocity profile has been set to be uniform which has been demonstrated in-vivo measurements for the ascending aorta distal to the ventrical valve [18]. The reconstruction of the flow profile is achieved based on the nineteen-coefficient Fourier equations described in equation (4).
 is the angular velocity, radian/s. At the outlets for brachiocephalic artery, left carotid artery, left subclavian artery and descending aorta outlet, the physiological pressure waveforms have been imposed termed as P1, P2, P3 and P4 respectively. The relationship between pressure and flow rate follows equation (5) referring to the literature [11].
Where r C is the resistance at the specific branch. Q is the flow rate in that branch, m 3 /s. Referring to the literature [19], the flow to branches of brachiocephalic artery, left carotid artery and left subclavian artery are all estimated to be 5%. Thus, the flow to the aortic outlet is 85% of the volumetric flow from inlets. P0 is the physiological pressure level (11,332 Pa). The constants for calculating P1, P2, P3 and P4 are 1.8e8, 1.44e8, 1.44e8 and 9e6 respectively. Therefore, the outlets pressure waveforms are displayed in Figure 4 for different cases. Moreover, the artery wall is considered to be rigid.

Computational mesh establishment
In order to achieve the mesh independence, the effects of different mesh elements on local WSS have been investigated as shown in Figure 5. The mesh sizes adopted are 0.47 million, 0.7 million, 1 million and 1.4 million respectively. As observed, the accuracy can be controlled below 2% when the mesh elements is over 1 million. Thus, the final mesh size is determined to be 1 million with mesh distribution around the junction structure between aorta and LVAD displayed in Figure 6.

Numerical methodology
Comsol 5.1 has been adopted for modeling, meshing and computing with finite element method. The generalized minimal residual method (GMRES) is used for iteration. Regarding the spatial and temporal discretization of governing equations, the streamline upwinding Petrov-Galerkin (SUPG) scheme and the backward differential formula (BDF) scheme with variable order are employed. The residual value of 10 -5 has been set for continuity and momentum equation. The initial time step of 0.001 s is set which can be adaptively adjusted by Comsol during computing process depending on the physics and scheme used. The simulations are performed on a desk computer GEDEON with 3.00 GHz Intel Xeon 64, 2 Processor and 128 GB of RAM. One simulation takes approximately 15h to compute two cardiac cycles (around 2 seconds of physical time for each case). At least two cycles have been carried out for each simulation in order to make sure that the final results will not be affected by the initial boundary conditions.

Hemodynamic parameters
Wall shear stress (WSS): Where T is the total time for one flow cycle, s. TAWSS represents the shear load over time acting on the arterial wall and the threshold value is 0.4 Pa [20].
Normalized artery surface of TAWSS below 0.4 Pa ( TA S ): Where T S is the total artery wall surface. The "if" function is used for integrating the artery wall surfaces where TAWSS<0.4 Pa (=1 if so, =0 if not) [21]. Thus TA S represents the artery surface percentage of potential disease sites.
Normalized index of hemolysis (NIH) [22]: As introduced above, hemolysis occurs because of the rupture of red blood cells membrane, followed by a leaking of hemoglobin in the plasma. The mechanical stress imposed to the cell's membrane can result in a rupture. Hence, a direct relation between stress and the amount of plasma free haemoglobin exists. The hemolysis index is a way to evaluate the amount of mechanical blood damage. Indeed, this index calculates the variation of hemoglobin contained in blood plasma. In this study, the NIH has been adopted to evaluate the hemolysis risks with critical value of 0.01 as defined in equation (9).  (9) Where Hct corresponds to the patients' hematocrit (percentage of volume occupation of red blood cells in the blood). We used the average value of hematocrit which is 40% and the k constant is the normal concentration of hemoglobin in the plasma, 1.5 g/100L. Therefore, NIH in g/100L allows us to conclude on the amount of red blood cells damage. Hb Hb  represents the change of hemoglobin.
Two lagrangian models have been taken into account to evaluate the amount of blood damage with the definition shown in equation (10) and (11) respectively. Model 1 is the most commonly used model because of its wider range of applicability and lower complexity of use developed by M. Grigioni et al. [23]. In addition, model 2 has been proved to be more accurate among the existing models developed by A. Garon and M. I. Farinas [24,25].
Model 2: Hb (11) Where  is the local shear stress. t  is the exposure time for that shear stress value which is calculated by the ratio between cubic root of mesh cell volume and flow velocity within the cell. C ,  and  are constants which are determined by the regression of experimental data. Different values of these constants can be found in literature corresponding to different ranges of shear stress and exposure time. In the present study, two sets of constants are considered shown in Table 1 termed as GW constants and HO constants respectively. The applicability ranges of these two sets of constants are the most consistent with our case. To calculate NIH, 200 path lines are generated in flow domain. Then, the total hemolysis along path lines from the inlet to the outlet can be calculated. The final NIH value is the averaged one across these 200 path lines.

Results and discussion
The flow topology has been investigated in Figure 7 at systolic peak (t=0.14 s for reference case, t=0.3 s for other two cases) under different conditions. The streamlines tend to straightly go downward the artery in each case. Subtle swirls can be observed within the junction region for case 1 and case 2 which is affected by the coming flow conflict from LVAD and aorta. Higher velocity along aorta in case 1 is observed compared to the reference case and case 2 which is positively related to the systolic peak flow rate. The velocity vectors have been displayed as well and it is found that the flow direction is mainly dominated by the aortic inlet flow because of stronger inertial force . The peak velocity can be always found at upstream inner wall along ascending aorta in each case because of centrifugal forces directed from the center of curvature to the outer wall. Thus the pressure at outer wall and inner wall will be increased and decreased respectively along the aortic curvature. The NIH has been evaluated at different time instants in one cycle as shown in Figure 8. Comparisons of NIH value have been made between two different calculation models considering GW constants under different LVAD working conditions. As observed, a peak value of NIH tends to be obtained around the systolic peak for both calculation models because of the higher flow rate at this time instant. Larger NIH is observed with model 1 compared to model 2 in each case. Moreover, the peak values of NIH with model 1 and model 2 using GW constants are around the order of 10 -3 and 10 -4 respectively in each case. Compared to the NIH threshold value of 0.01, they all conform to the acceptable design to avoid hemolysis. This observation is also consistent with the works of P. A. Segalova [22]. Similarly, the comparisons of NIH between model 1 and model 2 with HO constants have been investigated in Figure 9 for different time instants and LVAD working conditions. Higher NIH can be observed with model 1 compared to model 2 which is consistent with the observations in Figure 10. The peak value of NIH tends to be generated around the systolic end for model 1 which is mainly affected by the stronger backward flow at this time instant. However, peak NIH tends to be obtained at systolic peak for model 2 corresponding to the higher flow rate. Compared to the NIH values calculated with GW constants shown in Figure 10, lower value of NIH can be observed with HO constants when focusing on the same calculation model. The peak values of NIH for model 1 and model 2 with HO constants are in the order of 10 -4 and 10 -5 separately in each case. Thus, based on the observations both in Figure 8 and 9, it can be concluded that hemolysis is not expected to occur in the aorta and mainly happens inside the pump.
(a) (b) (c) Fig. 9. Comparisons of NIH evaluation at different time instants between these two calculation models with HO constants for different LVAD working conditions.
As discussed above, the exposure time is calculated as the ratio of cubic root of mesh cell volume and local velocity within the cell which has a direct influence in the NIH evaluation. Thus, the distribution of calculated exposure time on the aortic slice (shown in Figure 10) has been explored at different time instants as displayed in Figure 11 with case 1 as the example. As observed, when the time goes from 0.3 s to 0.5 s (systolic phase), the exposure time along the aorta is increased corresponding to the decreased flow rate. Decreased flow velocity contributes to the longer exposure time. Consistently, higher level of exposure time at time instants 0.7 s and 0.9 s (diastolic phase) has been maintained compared to the systolic period. Moreover, compared to the ascending aorta, higher exposure time has been mainly kept at the descending aorta which is negatively related to the velocity distribution. Similar phenomena can be found with other two cases which are not shown here. TAWSS distribution has been investigated in Figure 12 for different working conditions of LVAD. High value of TAWSS in case 1 is found compared to reference case and case 2 which is mainly affected by the larger aortic inlet flow rate. Larger TAWSS exists along the junction part of aorta and LVAD with case 1 and case 2 compared to the downstream part of aorta. The presence of LVAD tends to increase the TAWSS value compared to the reference case without LVAD because of the decreased blood cycle period with reduced volumetric flow within normal conditions of aorta. As the potential risks of thrombosis and intimal hyperplasia are significantly related to the TAWSS with limit value of 0.4 Pa, the corresponding quantitative analysis of TA S has been made in Figure 13. The normalized artery surface where TAWSS below 0.4 Pa is integrated at these different cases with/without LVAD. As observed in Figure 13, the value of STA in case 1 keeps lower level compared to other cases. It is also noticed that TA S in reference case is still up to 18% which is influenced mainly by the aorta curvature and the lower aortic flow rate in one cycle.

Conclusions
In the present study, numerical investigations have been carried out regarding blood flow in aorta with LVAD of different working conditions. 3D simplified geometrical model of aorta with LVAD has been established. Possible blood cells damage (hemolysis) due to shear stress are investigated as well as the flow topology and hemodynamic variations. The results reveal that the existence of LVAD contributes to swirls within junction part and the blood flow is more disturbed during diastolic period compared to the systolic period. According to NIH value based on the evaluations of shear stress, there is no risk of hemolysis along the aorta for these different cases. Thus, it is concluded that the hemolysis mainly happens inside the pump of LVAD instead of the aorta. With the intervention of LVAD, the shear stress distribution at connecting region between LVAD and aorta is decreased with lower flow rate from aorta. Moreover, lower level of TA S can be found with larger aortic flow rate indicating that severe hear failure condition can increase the risks of thrombosis.