Investigation of hemodynamic bulk flow patterns caused by aortic stenosis using a combined 4D Flow MRI-CFD framework

TL;DR

This study uses advanced imaging and simulations to show how aortic stenosis changes blood flow patterns, increasing red blood cell damage.

Contribution

A validated 4D Flow MRI-based CFD model reveals altered turbulent and helical flow structures in aortic stenosis that may damage red blood cells.

Findings

Pathological flow in aortic stenosis increases shear stress on red blood cells by 125% during late systole.

A physiological bihelical flow structure is replaced by a left-handed helix in aortic stenosis, increasing turbulent kinetic energy.

The 4D Flow MRI-based CFD model shows excellent agreement with in-vivo data (R² = 0.9).

Abstract

Aortic stenosis (AS) leads to alterations of supra-valvular flow patterns which can cause increased damage of red blood cell (RBC) membranes. We investigated these patient specific patterns of a severe AS patient and their reversal in healthy flow through a 4D Flow MRI-based CFD methodology. Computational models of subject-specific aortic geometries were created using in-vivo medical imaging data. Temporally and spatially resolved boundary conditions derived from 4D Flow MRI were implemented for an AS patient and a healthy subject. After validation of the in-silico results with in-vivo data, a healthy inflow profile was set for the AS patient in the CFD model. Pathological versus healthy flow fields were compared regarding their blood flow characteristics, i.e., shear stresses on RBCs and helicity. The accuracy of the 4D Flow MRI-based CFD model was proven with excellent agreement between in-vivo and in-silico velocity fields and R² = 0.9. A pathological high shear stress region in the bulk flow was present during late systole with an increase of 125% compared to both healthy flow. The physiological bihelical structure with predominantly right-handed helices vanished for the pathological state. Instead, a left-handed helix appeared, accompanied by an overall increase in turbulent kinetic energy in areas of accumulated left-handed helicity. The validated 4D Flow MRI-based CFD model identified marked differences between AS and healthy flow. It suggests that altered turbulent and helical structures in the bulk flow are the cause for increased, potentially damaging forces acting upon RBCs in AS.

Aortic stenosis (AS) is a condition that alters the flow of blood through the aorta, potentially increasing damage to red blood cells (RBCs). It is therefore crucial to understand the precise nature of these changes, especially in the bulk region where the RBCs are mainly circulating. To investigate this, we developed a computational model based on patient-specific 4D Flow MRI data to simulate blood flow in a patient with severe AS and a healthy individual. Our high-fidelity models show high agreement with real-world data. We then imposed synthetically generated healthy flow conditions in the AS patient to eliminate age- and geometry-related uncertainties. Our results revealed significantly higher shear stresses and changes in helical flow patterns in the AS patient compared to healthy flow, suggesting that these alterations could increase the risk of RBC membrane damage. These findings underscore the importance of deepening our understanding about the interplay between bulk flow alterations and their impact on RBCs, as previous studies were mostly focused on near wall regions. This methodology could also be extended to other cardiovascular conditions, aiding in the development of patient-specific treatment strategies.

Figures

Location of ROIs in the aortic model.Shown for the AS case, the ROIs are used for validation and post-processing of the numerical results.

Validation results of the 4D Flow MRI-based CFD model.Each plot compares the 4D Flow MRI measured values with the corresponding numerically determined velocities. (a-b) Area-averaged temporal evolution of velocity for H-25 case on ROI B and C, respectively. (c) Linear regression plot for H-25 velocities on ROI A-C and SAG for area-averaged velocities throughout the cardiac cycle. (d-e) Area-averaged temporal evolution of velocity for AS-78 case on ROI B and C, respectively. (f) Linear regression plot for AS-78 velocities on ROI A-C and SAG for area-averaged velocities throughout the cardiac cycle. Dashed grey line represents the x = y line at 45°, revealing no obvious bias in the data. Abbreviations: H-25 – Healthy, 25 year old volunteer, AS-78 – Aortic Stenosis, 78-year old patient, ROI – region of interest.

Velocity contours at systolic peak.Shown on sagittal and cross sectional plane ROI. (a) Aortic stenosis case AS-78, (b) synthetic healthy case H-78, (c) healthy case H-25.

Temporal area-averaged shear stress evolutions and contours.Shown at their respective peak time frames on ROI A. Upper left: Shear stresses for the aortic stenosis case AS-78. Upper right: Shear stresses for the healthy case H-25. Bottom left: Shear stresses for the synthetic healthy case H-78.

Turbulent kinetic energy distribution contours.Displayed on sagittal plane SAG during deceleration for (a) aortic stenosis case AS-78, (b) synthetic healthy case H-78 and (c) healthy case H-25.

Temporal and spatial helicity distribution.Upper row: Area-averaged helicity evolution on ROI A and B for AS-78, H-78 and H-25. (a) LNH contour and tangentially projected streamline vectors on ROI A during deceleration peak for AS-78. (b) LNH contour and tangentially projected streamline vectors on ROI A during deceleration peak for H-78. (c) LNH contour and tangentially projected streamline vectors on ROI A during deceleration peak for H-25.

Equations

Full text

Introduction

The development of cardiovascular diseases is linked to altered hemodynamic environments of the volumetric blood flow and therefore, gaining a deeper insight into both pathological and physiological flow conditions is of undeniable importance [1, 2]. A pathology of high prevalence and clinical relevance is aortic valve stenosis (AS) and the American Heart Association (AHA) predicted an increase of more than 100% in patients suffering from this disease by 2050 [3]. A better understanding of AS pathophysiology can help to develop safer and more reliable therapies to tackle these upcoming challenges. On this front, aortic flow patterns during AS play a decisive role in the pathology [4]. Valvular pathologies lead to disturbed blood flow in the aorta which in turn has the potential to damage red blood cell (RBC) membranes. In a recent study by Quast et al., it has been shown that extracellular hemoglobin is elevated in AS patients which in turn promotes endothelial dysfunction [5]. In order to capture these in-vivo flow features non-invasively, 4D Flow Magnetic Resonance Imaging (4D Flow MRI) is commonly used. 4D Flow MRI alone, however, does not provide sufficient temporal and spatial resolution due to its space- and phase-averaging nature to reliably derive all velocity scales and gradient-based hemodynamic parameters, such as shear stress fields on circulating particles [6], which are essential quantities to investigate the cause of blood damage [2,7].

In order to overcome this limitation, patient-specific Computational Fluid Dynamics (CFD) has been gaining popularity in the field of medical engineering research. Here, the advantages of high spatial and temporal resolution of CFD simulations are combined with subject-specific information about geometry and boundary conditions. Different approaches have been developed over time to increase patient-specificity and therefore clinical credibility of these models, e.g., the use of realistic geometries segmented from medical images. More recent research is focusing on the definition of accurate boundary conditions at the inlets and outlets of vessel [2]. These boundary conditions can be defined in varying degrees of complexity and accuracy. Some studies used generic curves and uniform or parabolic waveforms at the inlet [8–11]. One-dimensional modeling strategies at the inlet, however, neglect the spatial in-plane distribution of velocity which is highly complex and particularly characteristic and therefore essential for non-physiological flow states, such as AS [12]. Pirola et al. [13] emphasize the importance of including secondary flows in the prescription of inlet profiles in order to achieve realistic velocity, stress and helicity fields, in particular in patients with aortic valve pathologies. Highest subject-specificity and agreement with in-vivo flow can be achieved by using flow measurement techniques, such as Doppler echocardiography, time-resolved 2D Flow MRI or 4D Flow MRI. While Doppler echocardiography can only provide one-directional flow information and the current gold standard for volumetric flow analysis of 2D Flow MRI is limited by a two-dimensional imaging plane, 4D Flow MRI is known to be the only existing method for non-invasive, time-resolved, and three-dimensional measurement of blood velocity in-vivo [2,8]. Similarly, the complexity of outlet boundary conditions can be increased from simple constant pressure or pressure waveforms by coupling 1-D models of the downstream vasculature and 0-D models of the capillary bed, as has been summarized in a review by Quarteroni et al. [8–10,14,15]. 1-D and 0-D models simplify arterial networks by neglecting spatial variations of single arterial segments, therefore enabling their mathematical description using algebraic and differential equations. Alastruey et al. have carried out a study on the impact of different modelling assumptions in 1D/0D and 3D/0D aortic flow models, with average waveform inconsistencies of < 7% between in-silico and in-vivo results [16].

When it comes to comparing patient-specific flow patterns, the definition of a control group poses challenges. Using age-matched individuals introduces anatomical differences as an unavoidable confounding factor. Furthermore, it remains challenging to recruit age-matched scans of healthy volunteers without any other diseases or pathologies [17]. According to the AHA, 48% of the American population aged ≥ 20 years suffer from a cardiovascular disease, with increasing prevalence with age in both sexes [18]. A potential approach to overcome this limitation could be the generation of a patient-specific synthetic healthy model as in-silico physiological reference for a patient, using pre-acquired knowledge about physiological flow profiles and patterns. According to Youssefi et al. [12], patient-specific profiles can be sufficiently substituted by idealized profiles when investigating healthy flows.

In summary, studies on ascending aortic blood flow exist using in-vitro and in-vivo measurements and imaging for validation [8,16,19–21]. In-vitro techniques for validation include, for example, particle imaging velocimetry, and in-vivo measurements involve for instance echocardiography, cardiac catheter measurements or 4D Flow MRI. Further, these existing studies often investigate flow behavior in idealized phantoms or for physiological cases. Thorough investigations of pathological flows with modeling parameters derived from in-vivo imaging data pose as a more recent advance in in-silico analyses. However, investigations of pathological aortic bulk flow structures during AS, and their relationship to potential RBC damage are lacking. In the past, the focus has mostly lied on the near-wall hemodynamics, e.g., wall shear stresses, residence times and pressure distributions, and only little knowledge is available on volumetric bulk properties within the vessel lumen. However, bridging this gap is indispensable to continuously deepen the clinical knowledge and understanding of in-vivo aortic hemodynamics and their interplay with cardiovascular diseases.

Therefore, the goal of the present work is to develop a 4D Flow MRI-based CFD methodology to generate patient-specific models of healthy and pathological aortic flow and to show its ability to accurately characterize blood bulk flow structures and alterations caused by AS in the aortic arch. Applying this methodology to large patient cohorts and by comparing healthy and pathological flow structures, hemodynamic differences between these flow fields will aid the identification of those flow features with potentially damaging effect on circulating erythrocytes during AS.

Materials and methods

Ethics statement

All acquisitions in this study were conducted for research purposes according to the Declaration of Helsinki following local ethical approval by the respective Institutional Review Boards of the University Hospital Düsseldorf (HHU, No. 2018-86, No. 5761R, No. 4080, R5761R) and the University Hospital Aachen (EK 038/22), Germany. Formal written consent of all study participants was obtained before the study.

The complex aortic flow fields of a healthy adult and an adult suffering from severe AS were determined using a methodology of CFD analysis based on 4D Flow MRI data and are described in detail in the following. The two datasets served as input data for the generation of a synthetic healthy model of the AS patient. The derived hemodynamics were compared using qualitative and quantitative methods to identify the largest differences in the aortic flow structures between the physiological and pathological conditions measured in this study.

Study design

4D Flow MRI and multisliced computed tomography (CT) measurements of a 78 year-old, female patient diagnosed with severe AS (AS-78), originally acquired at the University Hospital Düsseldorf (UKD) were anonymized and made available for retrospective analysis. In addition, a 4D Flow MRI scan of a healthy (H-25) volunteer (female, age = 25 years) was acquired at the University Hospital Aachen (UKA) as physiological reference. Approvals were granted for all acquisitions by the respective Institutional Review Boards. 4D Flow data were acquired on 1.5 T MRI systems (Achieva at UKD and Ingenia Ambition at UKA, Philips Healthcare, The Netherlands) using the torso and posterior coil. Relevant scan parameters are listed in S1 Table. Single source CT data were acquired using a SOMATOM Definition AS+ (Siemens Healthcare, Forchheim, Germany), with a temporal resolution of 150 ms and collimation of 128 x 0.6 mm for the AS-78 patient.

In-vivo data processing

The handling of in-vivo data followed general guidelines and recommendations presented in the cardiovascular 4D Flow MRI consensus statement by Dyverfeldt, Bissell et al. [2].

The phase contrast DICOM images were processed using in-house developed MATLAB R2021b (MathWorks, Natick, MA) scripts. The data were first sorted according to timeframe and slice number. Furthermore, the stored greyscale values were transformed to their flow values in cm/s by multiplication with the velocity encoding value (V_enc_). The magnitude images were used to segment the aortic geometry (aortic arch and supraaortic vessels) during systolic peak ejection with the open-source software ITK Snap (www.itksnap.org) [22]. The segmentation was carried out by a used-guided automatic segmentation process based on active contours and region competition, followed by manual postprocessing to account for errors in the captured contours. Notably, the AS case required more manual refinement due to retrospective use of the underlying MRI scans from a previous study by Quast et al. [5]. Thus, the supraaortic vessels in particular were not the main region of interest and therefore did not provide sufficient quality for the active contour method to act successfully. Finally, the magnitude-based binary segmentations were superposed with the volumetric velocity data for each timeframe from the 4D Flow MRI scans using an in-house MATLAB script. This yielded a segmented velocity volume for the patient-specific lumen of the thoracic aorta. Resulting flow fields were visualized in ParaView 5.10.1 (KitWare, Clifton Park, NY).

Computational model

For subsequent generation of the numerical models, the commercially available CFD software package ANSYS CFX 2021 R2 was used, including the solid modeling Computer-Aided Design (CAD) software ANSYS SpaceClaim 2021 R2 and the meshing tool ANSYS ICEM CFD 2021 R2, all developed by ANSYS, Inc. (Canonsburg, PA).

The vascular geometry for the AS-78 case was semi-automatically segmented from CT images of the same patient. For the healthy case, the segmentation obtained from the 4D Flow MRI scans was used. Subsequently, the segmentation surfaces were manually smoothed using ANSYS SpaceClaim and all inlet and outlet surfaces were cut perpendicularly. Finally, the patient-specific geometries were meshed with suitable grid settings ensuring mesh independency of the numerical results. Meshing was performed with ANSYS ICEM CFD using unstructured tetrahedral meshes and prism elements for the near wall regions. A mesh independence study (see S1 Method for further information) yielded a number of 4.65 million elements and 12 prism layers with a refinement region in the ascending aorta.

Subsequently, the meshes were imported to ANSYS CFX. Here, the numerical model was built using patient-specific boundary conditions derived from the in-vivo measurements.

The Carreau-Yasuda model was implemented for modeling of the Non-Newtonian blood behavior in terms of a shear rate $[eqn]$ dependent viscosity η function, utilizing parameter values defined by Abraham et al. [23] and a constant density ρ of 1060 kg/m^3^.

$$[eqn]$$

Here, $[eqn]$ is the viscosity when the shear rate tends to zero, $[eqn]$ is the viscosity for high shear rates, λ is a time constant, a and n are dimensionless parameters. Turbulent flow has been observed in AS in numerous studies, e.g., in a turbulence analysis by Manchester et al. [24]. According to Benim et al. [25], the SST model can efficiently adapt to transitional effects and is widely used in industrial and academic context. Therefore, the Reynolds-averaged SST model was chosen to describe the expected formation of turbulent eddies.

The inlet of the numerical model was located shortly above the aortic valve and the 4-dimensional velocity profile serving as input data was extracted from the corresponding plane in the 4D Flow MRI velocity volume, see S1 Fig. The velocity values were interpolated linearly in both space and time to adjust to the high-resolution CFD grid and time using an in-house MATLAB script. A preliminary analysis of cycle-dependency showed, that the results converged for the third heart cycle. Thus, three cardiac cycles were simulated and the third one was post-processed, with a numerical time step size of 1/40 of the corresponding MRI timeframe (0.00111563s for H-25, 0.000255s for AS-78)

The supra-aortic vessels and the descending aorta were treated as openings defined by loss coefficients, which control the resulting mass flow in the respective branch i and serve as a substitute of the downstream vasculature and its impact on the flow structure. These loss coefficients can only be determined iteratively, following equation (2) from Benim et al. [25].

$$[eqn]$$

Here, ρ quantifies to fluid density, $[eqn]$ the mean pressure at outlet branch i, and $[eqn]$ the mean velocity at the respective outlet. The pressure $[eqn]$ denotes the pressure present at the farthest end of the downstream vasculature of each considered branch. Benim et al. [25] postulated for this pressure to be the same for all branches. According to numerous studies, the percentage distribution of mass flow among the four branches of the aortic arch (brachiocephalic artery, left common carotid artery, left subclavian artery, descending aorta), remains constant over a cardiac cycle. Therefore, for each numerical time step, an inner feedback loop was implemented with 100 iterations, which compared the percentage mass flow distribution through one opening for each iteration step with the prescribed percentage flow. According to the deviation, the loss coefficient was adapted for the next iteration until the two values were synchronized. The required pressure $[eqn]$ , resulting for this iterative process, was then implemented as new outlet boundary condition. This process was repeated for each numerical time step and a detailed implementation workflow can be found in S2 Method. The percentage flow distributions were calculated from values provided in Benim et al. [25] and are listed in Table 1. Finally, the vessel wall was modeled as a rigid, no-slip wall. Convergence control was set for a root mean square residual level of 10e-5 for the mass and momentum conservation equations.

Table 1: Mass flow of aortic outlets, defined as percentages of inlet flow.

Synthetic healthy model

It is known, that aortic hemodynamics change with increasing age of the individual [26, 27]. Therefore, in order to identify the hemodynamic changes caused by the presence of aortic stenosis, an age-matched healthy reference model of the aortic flow is required. However, since the risk for developing cardiovascular diseases increases with age, it poses a large challenge to recruit age-matched healthy volunteers that do not show any signs of cardiovascular pathologies or other age-related comorbidities. Moreover, by using a second individual as healthy reference, inter-individual variabilities, e.g., in shape, are introduced which can be misleading. This is the case for the H-25 and the AS-78 subjects, showing differences in age and geometry. Thus, a synthetic healthy model of the AS patient was developed and checked for credibility by comparing the resulting flow patterns with the healthy model. For the development of the synthetic healthy model of the AS patient (H-78), a parabolic profile was assumed, which, according to Youseffi et al. [12], is sufficiently accurate as inlet conditions to be a simplified substitute for medical data-derived profiles in case of physiological flow. A similar approach was also implemented by Zolfaghari et al. [28]. The parabolic profile was adapted to match the 4D Flow MRI-derived volume flux of the AS patient at each MRI acquisition timestep. The resulting geometry-specific, parabolic velocity profiles were interpolated between these time points by fitting a Fourier series. According to Morbiducci et al. [29], in-plane velocity components are non-negligible for bulk flow structure investigations, adding up to approximately 32% of axial velocity during systole and 111% during diastole. Therefore, secondary flows were prescribed as percentage amounts of the current axial flow magnitude. The percentage values were derived from the healthy case H-25 for each MRI timestep and used to determine the Fourier series coefficients to fit these datapoints as well. This resulting patient-specific, three-dimensional velocity waveform was prescribed as inlet boundary condition to mimic synthetic healthy flow of the AS patient. The detailed calculations and Fourier series fittings can be found in S3 Method.

In total, three CFD models were created and the setup is summarized in Table 2.

Table 2: CFD models overview.

Numerical evaluation

Finally, several parameters of interest were extracted from the CFD models for the validation of the model as well as for the subsequent investigation of the aortic flow behavior at characteristic timeframes. The following cross-sectional and sagittal planes served as regions of interest (ROI) for investigation: ROI A, ROI B, ROI C and SAG, as shown in Fig 1.

Location of ROIs in the aortic model.Shown for the AS case, the ROIs are used for validation and post-processing of the numerical results.

The temporal and spatial velocity fields were subjects of investigation for validating the generated CFD model. Therefore, the area-averaged velocity values on ROI A-C were compared between the CFD and MRI data for all MRI timeframes by means of linear regression analysis. Further, spatial velocity distribution within these ROIs was analyzed qualitatively for the systolic timeframe.

In order to locate the region of highest forces in the bulk flow acting upon erythrocytes potentially causing damage for the AS case, acting stresses and turbulence parameters were processed in both their temporal and spatial properties. The following quantities were of interest: Total shear stresses (TSS), together with Reynolds shear stresses (RSS) and viscous shear stresses (VSS). Turbulent kinetic energy (TKE), helicity and local normalized helicity (LNH). Their individual derivation and implementation is explained in S4 Method.

Results

Validation of the in-silico model

For a quantitative evaluation of the temporal evolution of velocity predicted by the CFD simulation, the area-averaged velocity at ROI A, B and C are plotted over the cardiac cycle for MRI and CFD for H-25 and AS-78, respectively. The temporal evolution for ROI B and C are shown in Fig 2A, 2B, 2D and 2E, the corresponding evolution for ROI A can be found in S2 Fig. The noise of each acquisition was calculated by determining the standard deviation of the measured velocity in the surrounding non-vessel regions during diastole, where no signal should exist. The resulting noise values of approximately 57 mm/s for AS-78 and 46 mm/s for H-25 correspond to approximately 5% of the V_enc_ value, which Papathanasopoulou et al [30] determined to be an approximation of the noise in MRI measurement. The noise is represented by the shaded areas in Fig 2. Hence, the numerical results show very good agreement concerning the shape and magnitude of the area-averaged velocity. Additionally, a linear regression analysis for ROI A-C and SAG deliver an R^2^-value of 0.9 (p < 0.001) for both the physiological and the pathological flow, indicating excellent agreement and statistically significant correlation between 4D Flow MRI-based CFD results and the measured in-vivo 4D Flow MRI velocities. The corresponding regression plots are presented in Fig 2C and 2F. The x=y curve reveals no obvious biases.

Validation results of the 4D Flow MRI-based CFD model.Each plot compares the 4D Flow MRI measured values with the corresponding numerically determined velocities. (a-b) Area-averaged temporal evolution of velocity for H-25 case on ROI B and C, respectively. (c) Linear regression plot for H-25 velocities on ROI A-C and SAG for area-averaged velocities throughout the cardiac cycle. (d-e) Area-averaged temporal evolution of velocity for AS-78 case on ROI B and C, respectively. (f) Linear regression plot for AS-78 velocities on ROI A-C and SAG for area-averaged velocities throughout the cardiac cycle. Dashed grey line represents the x = y line at 45°, revealing no obvious bias in the data. Abbreviations: H-25 – Healthy, 25 year old volunteer, AS-78 – Aortic Stenosis, 78-year old patient, ROI – region of interest.

Further, the numerically calculated and measured velocity fields also show high agreement spatially. Examples of spatial flow patterns and depictions of momentary regions of higher and lower velocities on ROI A-C are appended in the Supporting Information (S1 Result) for the pathological flow field during early systole. Here, the difference between the high resolution CFD grid and the corresponding low-resolution MRI grid also becomes obvious.

Physiological and pathological in-silico velocity fields

With excellent agreement between the in-vivo and in-silico results, the in-silico velocity fields are investigated in more detail for the three numerical models.

First, a credibility comparison between H-25 and H-78 shows that both models reproduce an overall homogeneous, parabolic shaped and rotationally symmetric velocity distribution in the ascending aorta, as shown in Fig 3. In the descending aorta, however, large differences can be observed with high velocities for H-25, while the velocity magnitudes remain relatively constant along the aortic arch for H-78.

Velocity contours at systolic peak.Shown on sagittal and cross sectional plane ROI. (a) Aortic stenosis case AS-78, (b) synthetic healthy case H-78, (c) healthy case H-25.

In contrary, the numerically derived systolic aortic stenotic flow displays an eccentric high velocity jet flow along the outer wall of the ascending aorta. When comparing the maximum and averaged velocities of cross-sectional planes in the ascending aorta over a cardiac cycle, the area-averaged systolic value is higher in the healthy cases, whereas the peak velocity is significantly higher for the aortic stenosis case, due to the stenotic jet. In the descending aorta, both average and maximum velocity are higher in physiological flow. Thus, the highest velocities are located in the ascending aorta for the AS patient and further downstream for both H-25 and H-78. The volume-averaged (maximum) Reynolds number (Re) at systolic peak in the ascending aorta is Re = 2350 (5126.5) for AS-78, Re = 2027 (7887.19) for H-78 and Re = 2447 (6950.8) for H-25.

Shear stresses acting onto erythrocytes

First, the total shear stress (TSS), viscous shear stress (VSS) and Reynolds shear stress (RSS) are compared between the healthy and the synthetic healthy flow. The synthetically generated shear stress field on ROI A displays similar temporal and spatial distribution as the reference field in H-25. In both cases, RSS remains negligible and the TSS and VSS peak coincides with the moment of systolic velocity peak, as shown in Fig 4. The peak TSS contour on ROI A is rather homogeneous in both cases with the highest values occurring at the vessel wall. Further, Fig 4 depicts results of the pathological flow on ROI A. The area-averaged, temporal evolution shows, that TSS for the AS-78 case is significantly higher than for the H-25 and H-78 case, with a second peak occurring during deceleration and being the global stress peak. At this point, it becomes 125% higher than the TSS value in the healthy case. The first shear stress peak also shifts into the deceleration phase whereas it occurred at the systolic peak for both healthy cases. Further, RSS as an indicator for turbulent flow is now clearly contributing to the overall acting stresses for AS. The AS-78 TSS contour in Fig 4 shows the eccentric region of high TSS in the pathological bulk flow at its peak, whereas the physiological stress fields do not show any obvious characteristics at its respective peak.

Temporal area-averaged shear stress evolutions and contours.Shown at their respective peak time frames on ROI A. Upper left: Shear stresses for the aortic stenosis case AS-78. Upper right: Shear stresses for the healthy case H-25. Bottom left: Shear stresses for the synthetic healthy case H-78.

With the significant rise in bulk RSS for AS-78, bulk flow behavior is investigated in detail. First, turbulent kinetic energy (TKE) as a quantity of turbulence is compared qualitatively on the sagittal plane SAG. Both H-25 and H-78 display negligible elevations of TKE values in the entire domain, with a small region of higher TKE in the central aortic arch, see Fig 5. The analysis shows strong alteration for AS flow with swirling structures in the ascending aorta.

Turbulent kinetic energy distribution contours.Displayed on sagittal plane SAG during deceleration for (a) aortic stenosis case AS-78, (b) synthetic healthy case H-78 and (c) healthy case H-25.

Second, helicity and its locally normalized value, namely Local Normalized Helicity (LNH), were analyzed in the ascending aorta over the heart cycle and during late systole, respectively, as shown in Fig 6. For both healthy models, the average helicity values remain positive over the entire cardiac cycle with their peak occurring shortly after the systolic velocity peak in ROI A and B. The area-averaged values are lower for H-78 and the cross-sectional contour of the synthetic healthy case displays less chaotic patterns and a smooth transition between a smaller dominantly left-handed to a larger dominantly right handed region, whereas the H-25 case shows more mixed patterns. Overall, if putting a convex hull around the predominantly right- and left-handed regions, both cases lead to similar distributions of positive and negative helicity areas. In AS, however, the mean helicity turns from right-handed to left-handed characteristics during late systole, and the magnitude reaches its peak later than in the healthy case. Fig 6C and 6D present the LNH contours on ROI A at their peak during deceleration. The healthy flow displays two counterrotating helices (b, c), while one nearly central, left-handed vortex can be identified for the AS flow (a). When comparing the LNH contour for AS-78 with the corresponding TKE contour in Fig 5, the region of highest TKE coincides with the region of left-handed helicity (LNH < 0).

Temporal and spatial helicity distribution.Upper row: Area-averaged helicity evolution on ROI A and B for AS-78, H-78 and H-25. (a) LNH contour and tangentially projected streamline vectors on ROI A during deceleration peak for AS-78. (b) LNH contour and tangentially projected streamline vectors on ROI A during deceleration peak for H-78. (c) LNH contour and tangentially projected streamline vectors on ROI A during deceleration peak for H-25.

Discussion

The present study aimed to identify the alterations in aortic blood flow patterns caused by aortic valve stenosis which may be a leading trigger for erythrocyte membrane damage. This was achieved by developing a subject-specific and realistic numerical flow model using in-vivo 4D Flow MRI. After validating the CFD model results with in-vivo data of two subjects, differences between AS and healthy aortic bulk flow were analyzed. In order to omit confounding factors such as anatomy and age, a synthetic healthy profile was set for the AS patient geometry. Especially for healthy flows, such approach is justifiable, as outlined by Youssefi et al. [12]. In order to create more realistic flow conditions, we enhanced the synthetic flow profile by the same secondary flow components as seen in the healthy subject of 25 years (H-25). This aspect is crucial when one aims to realistically model vorticity and helicity of aortic flow [29,31].

Our comparison of the synthetic healthy (H-78) and realistic healthy (H-25) yielded overall good agreement in ascending aorta flow features such as homogeneous velocity distributions, very low turbulent kinetic energies and right-handed helices (see Figs 3-6). With the latter flow feature being underestimated by the synthetic flow case. Differences in the descending aorta with H-25 flow velocities being higher than H-78 flow velocities can be attributed to age-dependent changes. Interestingly, the in-vivo study of Garcia et al. [27] has identified similarly increased flow velocities in the descending aorta when comparing female subjects of 21-39 years with subjects of > 60 years. Overall, the synthetic healthy flow model H-78 can be considered sufficiently credible to serve as patient-specific, healthy reference for the AS patient in our study and was used for the subsequent comparison with the AS flow. Thus, in the remaining analysis pathological (AS-78) and healthy (H-78) were compared.

Flow alterations during aortic stenosis and erythrocyte damage

Differences in the cross-sectional velocity distribution were caused by the presence of a high-velocity jet due to the reduced orifice of the aortic valve in AS. In order to pinpoint the impact of aortic flow alterations on shear stress on circulating cells (e.g., erythrocytes), shear stress distributions, turbulence kinetic energy and helicity were investigated in detail. The TSS evolution displays a key novel finding: it develops a second peak for the pathological flow during early diastole, which shows an increase by 125% compared to the corresponding maximum physiological values. Erythrocyte damage is known to be caused by a combination of elevated shear stresses and exposure time. Therefore, the increased shear stresses acting in the bulk flow, which are indicative of the pathological state during AS, are likely to impact and apply damaging forces on circulating RBCs. Nonetheless, only the spatial characteristic of forces acting onto the erythrocytes is investigated in this present work. Flow-induced damage of erythrocytes, however, is an integrated effect of both spatial and temporal factors and Lagrangian-based analyses are required, since the instantaneous streamlines and particle pathlines do not align in complex aortic flow. Further, the interaction between the AS jet and the aortic wall can lead to tissue remodeling such as dilation of the ascending aorta and aneurysm formation, which has been shown in a previous clinical study [32].

Next, an additional investigation of the bulk flow structures during the TSS peak was carried out for deeper understanding. Here, an increase in turbulence characteristics in the bulk flow of the ascending aorta is present in AS, whereas TKE was negligible for the healthy case. Further, significant differences in helicity were identified for the ascending aortic bulk hemodynamics, which is a known parameter to evaluate blood flow alterations, particularly caused by valvular diseases [31,33]. However, only little is known about the helical properties during AS [34]. According to Kilner et al. [35] and Markl et al. [36], two consistent features present in physiological flow of the human aorta are: helicity and retrograde flow. The anatomical curvature of the aortic arch leads to right-handed, well-aligned helical outflow during systolic peak. During deceleration, the ascending aorta is still predominantly filled with right-handed helices with a small portion of retrograde flow developing over a short distance during end-systole, which Morbiducci et al. [31] describe as a bihelical pattern with two counterrotating, Dean-like vortices, which is also characteristic for fully developed flow in a bended pipe. This can also be detected in the healthy flow structure of H-25 in this study, and to some extent in H-78, given its nature to underestimate helicity as described earlier. The numerical results of the pathological flow in this study, however, show a decrease in right-handed helicity and an increase of the region of left-handed vortices. These areas of accumulated left-handed helices coincide with the areas of high TKE (see Figs 5 and 6), thus pointing towards a correlation between left-handed helices and high presence of RBC-damaging turbulence potentially impacting erythrocyte function.

Overall, the results allow for the formulation of two main findings that – to the authors’ best knowledge – have not been identified in existing literature so far:

1)The presence of a second TSS peak in the bulk flow during early diastole was observed, which correlates with the onset of strong turbulent structures in AS flow.2)AS further causes a decrease in physiological right-handed helical structures. This also leads to the destruction of the bihelical structure, which is found to be characteristic for physiological aortic flow. Instead, a main left-handed helix is present, coinciding with the peak region of TKE.

These two hemodynamic alterations are hypothesized to cause damage on circulating cells and cellular dysfunction. However, our work only included a size of N = 2. Therefore, these findings need to be further validated with a larger cohort of patients with the developed framework.

4D Flow MRI-based CFD methodology

In order to validate the developed framework and the resulting findings, the contours of the numerical velocities profiles were analyzed on a set of planar, cross-sectional ROIs along the aortic centerline. Their temporal and spatial agreement with the corresponding velocity profiles obtained in the MRI process were compared qualitatively and quantitatively. A good qualitative agreement of patterns and contours is found on the ROIs. The linear regression analysis used to assess the overall correlation between the measured MRI velocity and the predicted CFD results yields excellent results with R^2^ = 0.9 for both the physiological and the pathological case.

To the best of the authors’ knowledge, this is the first study to compare patient-specific bulk flow patterns in healthy flow and aortic stenosis by combining spatially and temporally resolved 4D Flow MRI measurements and CFD modeling in order to identify the pathological impact on erythrocytes due to flow alterations during AS. It implements realistic, four-dimensional in-vivo boundary conditions and the numerical modeling of non-Newtonian and turbulent blood behavior. However, the 4D Flow MRI acquisition and post-processing suffered from limitations characteristic for this technique, e.g., the averaging nature of the scans, systematic errors [2], or segmentation uncertainties. Further, simplifying the vessel wall with rigid properties can lead to overestimations of some quantities such as the wall shear stress. According to Torii et al. [37], however, this overestimation becomes negligible (< 5%) when analyzing time-averaged values and was therefore accepted in this study to keep computational costs low.

Finally, the SST turbulence model was used in this study, which is a widely used Reynolds-averaging approach to model turbulent structures. Especially describing the transition from laminar to turbulent flow, however, is challenging. While the employed SST model is expected to possess abilities in capturing transition to turbulence, deeper investigations of more complex turbulence models, especially scale resolving ones, should be investigated. Examples exist in the literature investigating the overall impact of turbulence model choice, for instance [24,38].

Conclusion

Large differences between the given models of healthy flow and the altered hemodynamics present during aortic stenosis were identified: The occurrence of a second TSS peak in the bulk flow of the ascending aorta during early diastole and the loss of the physiological, predominantly right-handed bihelical structure. The left-handed helices coincided with areas of high TKE. These results were obtained by the development of a 4D Flow MRI-based CFD framework of high fidelity. These patient-specific models promise patient-individual treatment planning in the future, e.g., for aortic valve replacement, and improved long-term outcomes of patients undergoing such procedures.

Supporting information

S1 TableScanning parameters of the 4D Flow MRI measurements.(PDF)

S1 MethodsMesh independency study and mesh parameters.(PDF)

S2 MethodsImplementation of Loss Coefficient Model in Theory and ANSYS CFX.(PDF)

S3 MethodsGeneration of synthetic healthy inflow profile based on pathological in-vivo 4D Flow MRI measurements.(PDF)

S4 MethodsMathematical derivation and explanation of all hemodynamic parameters of interest quantified in the study.(PDF)

S1 Fig4D Flow MRI-derived velocity profiles prescribed at the aortic inlet for H-25 and AS-78.Prescribed inlet velocity profiles at systolic peak for (a) H-25 and (b) AS-78. The inlet boundary conditions of the H-25 and AS-78 case were defined by 4D Flow MRI-derived velocity profiles. The respective profiles at systolic peak are shown in S3 Fig. The profile for the healthy case shows the typical nearly parabolic profile, whereas the AS profile is characterized by a strong excentric velocity peak, even leading to a small portion of backflow.(TIF)

S2 FigValidation of 4D Flow MRI-based CFD models for H-25 and AS-78 on ROI A with corresponding 4D Flow MRI measurements.Validation results of the 4D Flow MRI-based model, each plot compares the 4D Flow MRI measured values with the corresponding numerically determined velocities. (a) Area-averaged temporal evolution of velocity for H-25 case on ROI A. (b) temporal evolution of velocity for AS-78 case on ROI A.(TIF)

S1 Result4D Flow MRI and CFD velocity contour maps on different cross-sectional planes for H-25 and AS-78.The velocity contour maps at systolic peak are shown for ROIs A, B and C for both subject cases. The large difference between low-resolution 4D Flow MRI scans and high-resolution CFD simulation can be observed.(PDF)