Impact of heart failure severity on ventricular assist device haemodynamics: a computational study

This computational fluid dynamics study investigates the necessity of incorporating heart failure severity in the preoperative planning of left ventricular assist device (LVAD) configurations, as it is often omitted from studies on LVAD performance. A parametric study was conducted examining a common range of LVAD to aortic root flow ratios (LVAD/AR-FR). A normal aortic root waveform was scaled by 5–30% in increments of 5% to represent the common range of flow pumped by the left ventricle for different levels of heart failure. A constant flow rate from the cannula compensated for the severity of heart failure in order to maintain normal total aortic flow rate. The results show that LVAD/AR-FR can have a significant but irregular impact on the perfusion and shear stress-related haemodynamic parameters of the subclavian and carotid arteries. Furthermore, it is found that a larger portion of the flow is directed towards the thoracic aorta at the expense of the carotid and subclavian arteries, regardless of LVAD/AR-FR. The irregular behaviour found in the subclavian and carotid arteries highlights the necessity of including the LVAD/AR-FR in the preoperative planning of an LVAD configuration, in order to accurately improve the effects on the cardiovascular system post implantation.


Introduction
Left ventricular assist devices (LVADs) are an established form of treatment for patients with heart failure and result in clear survival benefit. Although LVADs demonstrate significant survival improvements, they are also associated with postoperative complications such as ischemic and haemorrhagic stroke due to abnormal haemodynamics (Boyle et al. 2014;Menon et al. 2012;Parikh et al. 2016). With a growing proportion of LVADs being implanted as a destination therapy, the associated risk of stroke (especially for women) remains too high (Kirklin et al. 2015(Kirklin et al. , 2014(Kirklin et al. , 2013Parikh et al. 2016). Although there is currently a wide range of research surrounding the optimization of the LVAD outflow cannula design to address this problem, there appears to be no emphasis on the haemodynamic effect of the flow ratio (FR) between the LVAD cannula outlet and aortic root (LVAD/AR-FR). Furthermore, this ratio is often omitted from studies on LVAD performance. The exact contribution of the remnant native cardiac output to the aortic flow is generally unknown at any given time. However, it can vary from patient to patient anywhere between approximately 5 and 30%, based on the average cardiac output of 5.08 L min − 1 (Brandfonbrener et al. 1955;Karmonik et al. 2014). Furthermore, LVAD/AR-FR varies from patient to patient depending on the severity of heart failure and amount of physical activity of the patient. Valuable insight may be gained by incorporating patient-specific LVAD/AR-FRs into preoperative LVAD configuration planning using simulation technologies.
There is a scope for work that defines the necessity of incorporating the LVAD/AR-FR into the preoperative planning of an LVAD configuration, in order to accurately improve the effects on the cardiovascular system post implantation. Therefore, the aim of this work is to quantify the cardiovascular effects associated with the common range of LVAD/ AR-FRs using computational fluid dynamics (CFD).
In the current study, a parametric analysis of a range of LVAD/AR-FRs is examined and compared against a reference solution. A three-element Windkessel model is used in order to accurately capture the physiological back-pressure response of the outlets. Metrics of the perfusion, as well as shear stress-based parameters are assessed and discussed, showing the importance of including LVAD/AR-FR in LVAD performance studies as well as preoperative LVAD planning.

Geometrical model
The most common and successful second-generation pump worldwide as bridge to transplant and destination therapy is the 'HeartMate II' (Thoratec Corp, Pleasanton, CA, USA) (Baldwin et al. 2015;Garbade et al. 2011;Smith et al. 2015). Therefore, in this study, a typical outflow cannula geometry (10-mm diameter) from a HeartMate II device configuration is digitally anastomosed to a patient-specific model of a 23-year-old female with normal aortic function and an AR diameter of 20 mm. The patient-specific geometry was produced using SimVascular software package (https://simtk.org) to process gadolinium-enhanced MRA image data (LaDisa et al. Jr 2011). The simplified cannula centreline path used by Osorio et al. (Osorio et al. 2013) with an inclination angle of Φ ≈ 30 degrees was reconstructed and implemented into this study and is anastomosed to the ascending aorta (AAo) approximately 10 mm (δ) from the innominate artery (see Fig.  1). An azimuthal angle of ≈ 17°(β) is implemented in order to minimize the jet effect from the cannula blocking flow from entering the aortic arch branches; this angle is similar to patient-specific azimuthal angles found in the literature (Karmonik et al. 2014(Karmonik et al. , 2012a.

Computational domain
The geometry was meshed using ANSYS-Meshing (Version 17.1). Following a series of mesh independency tests, the mesh with ≈ 3.3 million elements was used, which consisted of tetrahedral elements within the core region and prism layers (3 elements thick) near the wall to allow for large spatial velocity gradients.

Boundary conditions
As shown in Table 1, seven tests are studied in this work. 'Test 0' is used as a representative of normal aortic haemodynamics for comparison purposes, which does not include a cannula. 'Tests 1-6' represent a patient with an implanted LVAD suffering from various extents of heart failure. In all tests, a no-slip BC is applied to all walls and a rigid wall model is assumed (Brown et al. 2012;Moore et al. 1994). Details of the inlet and outlet BCs are shown in Table 1 and Fig. 2 and will be described next.
For all tests, a uniform velocity profile was imposed at the inlets. In vivo measurements in the literature have demonstrated that the velocity profile distal to the aortic valve is essentially flat in the AAo and only consists of a weak helical component (Caballero and Laín 2013). As shown in Fig. 2, for 'Test 0', the patient-specific PC-MRI waveform by LaDisa et al. (LaDisa et al. Jr 2011) was applied to the aortic root (AR). 'Tests 1-6' scaled the AR waveform by 5 to 30% in increments of 5% (corresponding to LVAD/AR-FRs between 70 and 95%). A constant flow rate from the cannula compensated for the unfulfilled total aortic flow to bring the total flow rate for all tests to 3.2 L min −1 . Each outlet is coupled to a three-element Windkessel model in this study, which has been implemented in the solver by the present authors. The threeelement Windkessel (Westerhof et al. 2009) has become the most widely used and accepted lumped parameter model of the circulatory system as it tends to produce realistic aortic pressures and flows (Capoccia 2015;Deyranlou et al. 2020;Kaufmann et al. 2012).

Numerical procedure
The finite-volume CFD solver ANSYS-CFX (Version 17.1) was used for the simulations and the commercial visualization tool, Ansys Ensight (Version 10.1.6) was used to post-process the results.
Blood was defined as a non-Newtonian, laminar, three-dimensional, incompressible, and isothermal fluid. Dynamic viscosity and density of the blood are set to 0.004 Pa s and 1060 kg m −1 , respectively (Osorio et al. 2013). Fig. 1 Schematic of the present configuration consisting of an aorta (that of a 23-year-old female) with normal aortic function with an LVAD outflow graft connected to the AAo. Inlets: AR and LVAD outflow graft. Outlets: RSA, RCCA, LCCA, LSA and TAo. Inclination angle, Φ ≈ 30°. Azimuthal angle, β ≈ 17°. Distance from innominate artery, δ ≈ 10 mm, similar to that used by Osorio et al. (2013)  Constant LVAD cannula flow rates were obtained by time-averaging the healthy pulsatile flow rate and scaling accordingly. Percentage values represent the boundary's contribution to the total aortic flow rate. 'Test 0' is a control model, in which heart function is normal and there is no LVAD outflow graft A fully implicit second-order backward Euler differencing scheme is used. Each simulation is run for 4 cardiac cycles (4 s) with a timestep size of 1.25 ms and 10 coefficient loops per timestep, which proved adequate in order to allow the boundary flow rate and pressure waveforms to converge (see Fig. 3), similar to previous studies in the literature (Satta et al. 2019;Wellnhofer et al. 2010). Results are recorded during the fourth cardiac cycle. For 'Test 0', the fourth cardiac cycle results correlate to producing a physiologically realistic blood pressure of 15131/8781 Pa (113/66 mmHg) (Stamler et al. 1993).

Results
In this section, the results from the simulations including timeaveraged velocity and 3 different haemodynamic parameters are presented. Note that the haemodynamics in 'Test 0' are assumed to be ideal as it is representative of the haemodynamics in a healthy patient without any invasive devices or cardiovascular dysfunction. Therefore, any comparisons will use 'Test 0' as a benchmark.

Time-averaged velocity
Through a range of LVAD/AR-FRs, the effects of the LVAD implementation on the perfusion of the aortic branches for different levels of heart failure are examined against 'Test 0'. In order to measure the aortic perfusion, TAFR and TAV are calculated over one cardiac cycle in order to examine the perfusion of vessels and general haemodynamics. Figure 4 provides a quantitative approach to assess the TAFR, whilst Figs. 5 and 6 show a more qualitative illustration of TAV using iso-surfaces and streamlines, respectively.
The results show that none of the geometries including an LVAD outflow cannula (i.e. 'Tests 1-6') experiences normal perfusion; this also supports common findings in the literature (Kaufmann et al. 2009;Laumen et al. 2010;Menon et al. 2013).

Haemodynamic parameters
The localization of various lesions including atherosclerosis and thrombosis has been extensively studied and have shown to be related to different local haemodynamic metrics including some of the work of the present authors Kabinejadian et al. 2016;Keshmiri et al. 2016;McElroy et al. 2016a;McElroy and Keshmiri 2018;Ruiz-Soler et al. 2017;Swanson et al. 2020). These haemodynamic parameters can be directly derived from the flow velocity fields obtained by CFD-based simulation tools. In the present study, these parameters are used as a tool to study whether the haemodynamics are favourable or detrimental to the health of the patients with an LVAD and to help identify some  to Table 1) characteristics of disturbed flow that are of pathobiological importance. The wall-based parameters chosen for the present analyses include WSS, time-averaged WSS (TAWSS) (He and Ku 1996), oscillatory shear index (OSI) (He and Ku 1996), and relative residence time (RRT) (Lee et al. 2009), calculated as: where τ ! W is the WSS vector and T is the time period of the flow cycle.
Time-averaged wall shear stress When comparing 'Test 0' with 'Tests 1-6' in Fig. 7, the latter all experience much less evenly distributed WSS, with generally higher WSS in aortic arch and DAo, whilst the WSS is lowered in all other locations. This is due to the cannula flow inducing nonphysiological haemodynamics and disturbing the WSS distribution (May-Newman et al. 2006). The most prominent feature to emerge as a result of increasing the LVAD/AR-FR is the effect of the impingement jet onto the DAo which can be clearly seen in the form of elevated TAWSS. To a lesser extent, proximal LSA and LCCA are also affected by the impingement jet exiting the LVAD cannula.
The flow in 'Test 0' may be a typical example of 'form follows function', as the human AAo is said to be insusceptible to atherosclerotic plaque formation (Liu et al. 2014). 'Test 0' does not produce any extreme distributions of WSS as the form has likely developed in such a way to maintain a relatively smooth distribution of WSS.
Similar to the TAV and TAFR results presented above, Fig. 8 quantitatively shows that different branches exhibit different degrees of sensitivity to increasing the flow from the LVAD. TAWSS within the subclavian and carotid arteries changed in an irregular manner as LVAD/AR-FR increased (locations 4-8 in Fig. 8). The TAWSS generally reduced by 22% within the innominate artery (location 4) and stayed relatively constant through 'Test 2 to 6', with the exception of 'Test 1', which increased slightly (+ 2%). Although the TAWSS within the subclavian and carotid arteries fluctuating as LVAD/AR-FR increases, an overall decrease in TAWSS of 40% was generally experienced.
Oscillatory shear index The subclavian and carotid arteries (as shown in Figs. 7 and 8) all experience much lower OSI compared with 'Test 0'. The innominate artery experiences a clear negative relationship between OSI and LVAD/AR-FR, due to its close proximity to the anastomosis, hence, experiencing a direct effect of the cannula flow rate. The other peripheral arteries however experience an almost zero value of OSI, most likely because (i) continuous flow from the LVAD dominates the pulsatile flow from the AR and (ii) the three-element Windkessel model has a high autoregulatory effect on the flow within the smaller vessels as previously discussed.
Unlike most locations, the AAo (location 3) experiences a significant increase in OSI from 'Test 1-3', peaking at 'Test 3' (+ 116%), and then significantly decreasing again until reaching almost healthy values of OSI at 'Test 6'. The elevated OSI magnitudes in location 3 are representative of a recirculating and unordered retrograde flow in the aortic arch near the cannula attachment point, where the LVAD and AR flows are mixing. Higher retrograde velocity represents an unnatural haemodynamic environment, hence a higher risk for the patient.
Relative residence time RRT in the subclavian and carotid arteries (locations 4-8) slightly increases when compared with 'Test 0' with no clear correlation to LVAD/AR-FR. Figure 7 also shows insignificant differences in the subclavian and carotid arteries for 'Tests 1-6'. The lack of any trend in the areaweighted averaged RRT is likely due to insignificant values in these locations. On the other hand, Figs. 7 and 8 show that the AR and AAo (locations 1 and 3) experience significantly elevated RRT with strong correlations to LVAD/AR-FR. The elevated RRT within the AR and AAo is again due to the cannula flow blocking and diverting the pulsatile flow from the AR and increasing the recirculating retrograde flow within the aortic arch. Furthermore, the DAo and TAo experience Fig. 4 Ratio (%) between the TAFR through the corresponding artery and the total TAFR through the system. Positive and negative ratios represent inlets and outlets, respectively. (a) Simulations that incorporated a three-element Windkessel model at the outlet boundaries. (b) Simulations that implemented a zero-pressure at the outlet boundaries regions of high RRT streaks which increase in size as LVAD/ AR-FR increases (this is shown most evidently in Fig. 7).

Subclavian and carotid arteries
The present results demonstrated that the main aortic outlets including the right subclavian artery (RSA), right common carotid artery (RCCA), left common carotid artery (LCCA) and left subclavian artery (LSA) for 'Tests 1-6' all experienced an insufficient amount of flow (i.e. malperfusion) compared with 'Test 0', with the exception of the RSA in 'Tests 1-2' (shown more clearly in an inset in Fig. 4a). The average perfusion of the subclavian and carotid arteries in most of the tests presented in the "Results" section was approximately 70% of what is deemed to be healthy. The general drop in flow rate to the subclavian and carotid arteries is mainly due to the jet from the cannula forming a blocking effect against the flow from the AR entering the subclavian and carotid arteries, as was shown in Figs. 5 and 6. A clear distinction was observed in Fig. 4a between the left and right carotid and subclavian arteries in relation to their sensitivity to different LVAD/AR-FRs. Specifically, both the LCCA and LSA maintained an almost constant outflow rate, showing minor changes in TAFR percentages for different LVAD/AR-FRs. On the contrary, the RCCA and RSA flow rates produce an irregular pattern as LVAD/AR-FR increases. The fluctuation of the RSA and RCCA flow rate is mainly due to their close proximity to the anastomosis, hence experiencing more disturbed flow (was shown more clearly in Fig. 6). Another reason behind the insensitivity of the LCCA and LSA to LVAD/AR-FR is likely to be due to the three-element Windkessel model used here creating a stronger autoregulatory effect on flow due to their smaller diameter (greater proximal and distal resistances).

Aorta
As was shown in Fig. 4a, adding the LVAD ('Tests 1-6') resulted in an average 75% increase in the thoracic aorta (TAo) flow rate. This is due to the underperfusion observed increasing from 'Test 1' to 'Test 4', peaking at 'Test 4', and then beginning to decrease.
It is worth noting that the perfusion experienced in this study should not be assumed to be experienced by all patients. Over a study of 956 patients with the 'HeartMate II' LVAD, Fig. 6 Time-averaged blood flow velocity streamlines coloured by TAV Fig. 7 Distributions of different haemodynamic parameters calculated for the aortic arch and descending TAo, viewed from outside the vessel. LVAD/ AR-FR increases from left to right (as shown by the percentages) the risk of ischemic stroke differed substantially by demographics; men and women over 65 years of age experienced an incidence rate of 0.028 and 0.067 events per patient-year, respectively (Boyle et al. 2014). The abnormal perfusion experienced in this study combined with the unclear pattern of ischemic events leads to the need of more thorough ways to predict the risk of such events on a patient-by-patient basis.

Shear stress-based haemodynamic parameters
The results of 3 haemodynamic parameters were presented in the "Haemodynamic parameters" section. These parameters have been chosen because (i) localized distribution of low-WSS and high-OSI strongly correlates with the locations of atheroma (He and Ku 1996), (ii) platelet activation may be induced by the combination of large exposure time and high shear stress (Hellums 1993;Ramstack et al. 1979;Wurzinger et al. 1986Wurzinger et al. , 1985, and (iii) stagnant and recirculating flow regions can cause platelet aggregation and thrombogenesis (Cito et al. 2013).
TAWSS The DAo (location 9) initially experienced a drop in TAWSS of 28% but then began experiencing a steady increase in TAWSS of approximately 13% per test, until equalling the healthy TAWSS at 'Test 3', and continuing to 'Test 6' with a relatively constant gradient. In agreement with the observations from Fig. 7, the more severe the heart failure, the larger and more intense the shear stress streak on the DAo due to the higher cannula jet velocity, which is required to maintain adequate overall flow rate. This phenomenon was also illustrated by the TAV and streamlines shown in Figs. 5 and 6, respectively.
The AAo (location 3) experienced a relatively linear negative relationship between TAWSS and LVAD/AR-FR. 'Test 1' experienced a sudden decrease in TAWSS by 62% and then steadily decreased by ≈ 3%. This shows that the AAo is less sensitive to changes in TAWSS by altering the LVAD/AR-FR but experiences a much lower TAWSS by simply introducing an LVAD into the model. TAWSS within the TAo (location 10) in 'Test 1' initially reduces by 60% and then remains relatively constant regardless of LVAD/AR-FR.
OSI The AR (location 1) experienced the greatest increase of OSI when compared against the reference case, with a strong positive correlation to LVAD/AR-FR. The main reason for such haemodynamic modification is that the high-velocity flow stream from the cannula tends to block and divert the pulsatile flow from the AR. The above observation, along with the findings of He and Ku (1996), suggests that the AAo may suffer from an increased probability of atheroma Fig. 8 Area-weighted average values of TAWSS, OSI and RRT at ten monitored locations within the geometry. LVAD/AR-FR increases from 'Test 0' to 'Test 6' development. It is also noted that the extent of this phenomenon is unlikely to be improved by the cannula tip positioning and configuration, based on positioning restrictions.
The DAo and TAo (locations 9 and 10) both experienced a negative relationship between LVAD/AR-FR, with the TAo being more sensitive; the OSI value within the TAo approximately halved each time the LVAD/AR-FR increased by 5%.
RRT The RRT distributions shown in the present study suggest that AR, AAo, DAo, and TAo are at serious risk for 'Tests 5-6' as these regions are experiencing significant stagnant and recirculating flow, which can cause platelet aggregation and thrombosis (Cito et al. 2013).

Effects of using conventional zero-pressure boundary condition
Despite more reliable methods existing, it is still common amongst the cardiovascular research community to use simpler outlet BCs such as zero-pressure (Cheng et al. 2014;Deyranlou et al. 2020;Janiga et al. 2015;Kabinejadian et al. 2016;Karmonik et al. 2014Karmonik et al. , 2012bMcElroy et al. 2016b;Moon et al. 2014;Ruiz-Soler et al. 2017;Shang et al. 2015). Running 'Tests 0-6' using a conventional zero-pressure BC for all outlets produces significantly different results and misleading conclusions when compared against using a threeelement Windkessel model. This is observed by comparing the TAFRs in Fig. 4a and b. Using a three-element Windkessel model concludes that the TAFR to the subclavian and carotid arteries is highly irregular and is overall reduced when compared against a healthy patient without an LVAD. In contrast, when using the zero-pressure BC, the subclavian and carotid arteries were less sensitive to the LVAD/AR-FR and behaved in a much more regular manner. Using the zeropressure BC also suggests that the LCCA and LSA are overperfused regardless of LVAD/AR-FR, which is in contradiction with the findings from the Windkessel model. Perfusion to the TAo using the zero-pressure BC was relatively unaffected by LVAD/AR-FR. This comparison highlights the importance of using a physiologically accurate BC when investigating haemodynamics within the aorta. Similar findings were found by McElroy and Keshmiri (2018) when studying the flow in an accurately represented rabbit aorta.

Limitations and future work
Due to the LVAD and aorta configuration in reality changing vastly from patient to patient, the results of this work only act as an indicator of the importance of incorporating the patientspecific LVAD/AR-FR when building simulations to model flow through the aorta during LVAD use or cardiopulmonary bypass surgery. Further studies are required in order to test the value of the LVAD/AR-FR on a wider range of LVAD configurations. Due to the sparsity of relevant data required to conduct a rigorous quantitative study of multiple patientspecific configurations, a more qualitative approach has been taken in the present study to assess the significance of the LVAD/AR-FR. Further value could be presented in this field by developing a methodology to accurately capture patientspecific LVAD/AR-FR. Clinical studies would need to be performed to validate hypothesized flow ratios. Then, with the knowledge that it is necessary to incorporate heart failure severity into CFD simulations, along with having a validated model to define the severity of heart failure on capable output flow rate, future models could be built to improve the accuracy of simulation-based preoperative planning techniques of the optimum cannula angle, design and location for LVAD and cardiopulmonary bypass surgery. Doing so will improve desired effects and reduce unfavourable risks on the circulatory system post-surgical intervention.

Conclusion
A patient-specific aorta model with normal aortic function was used in a CFD study to investigate the importance of incorporating the severity of heart failure in the computational preoperative planning of an LVAD configuration. Results were compared against a control model without an LVAD graft, representing healthy/normal aortic haemodynamics.
The results showed that in terms of overall perfusion, the smaller vessels (RSA, RCCA, LCCA and LSA) responded in an irregular manner as LVAD/AR-FR increased, with a general deficiency of perfusion by approximately 30%. On the other hand, the larger vessels' (aortic arch and TAo) perfusion response was generally more linear and predictable with a general increase of 75%. Moreover, similar patterns are experienced when comparing shear-based haemodynamic metrics including OSI, RRT and TAWSS. The fundamental reasons for these discrepancies in behaviour is mainly due to (i) the jet from the cannula blocking flow to the aortic arch branches (RSA, RCCA, LCCA and LSA), (ii) the aortic arch branches better autoregulate flow rate due to their smaller diameter, and (iii) the entrances of the aortic branches experiencing highly disturbed and retrograde flow due to their close proximity to the anastomosis.
The general irregular behaviour of the subclavian and carotid arteries to LVAD/AR-FR indicates that the perfusion and WSS-based haemodynamic metrics within these arteries cannot be accurately predicted unless the LVAD/AR-FR is incorporated into the computational preoperative planning of the optimal LVAD configuration. Engineering (MACE) at the University of Manchester, provided to the first author.

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.