Simulating Subject-Specific Aortic Hemodynamic Effects of Valvular Lesions in Rheumatic Heart Disease Free

Acute rheumatic fever followed by rheumatic heart disease (RHD) are conditions that may develop from untreated Group A β-hemolytic streptococcal pharyngitis [1,2], occurring commonly in low- to middle-income countries. RHD results from an auto-immune reaction from epitopes of streptococcus pyogenes cross-reacting with proteins in the heart, in patients with genetic susceptibility [2]. RHD is characterized by endocardial inflammation and valvular lesions which can lead to eccentric hypertrophy of the left ventricle [3–5]. Fibrotic stiffening and thickening of the valve occurs when the interstitial cells of the valve are stimulated due to a T-cell predominant immune response [6]. RHD valve lesions can cause stenosis or regurgitation in the mitral and aortic valves [6–8], which may eventually lead to heart failure. Heart failure is estimated to be present in 33–40% of individuals living with RHD [9,10], while other complications include strokes or atrial fibrillation [11,12], all of which can be exacerbated during pregnancy [13–15]. The actual number of patients living with RHD might be much higher than the estimated 40 × 10⁶ [12] due to subclinical cases [16]. Therefore, there remains an urgent need to further characterize RHD, find preventive solutions, and improve current treatments.

Valvular lesions in RHD are typically diagnosed using M-mode and Doppler transthoracic echocardiography through the quantification of leaflet thickness, shape, and flow across the valves over the cardiac cycle. Computed tomography, X-ray fluoroscopy, and cardiovascular magnetic resonance (CMR) may also be used to further diagnose RHD severity [7,17]. While mitral valve regurgitation and stenosis are the most prevalent valve lesions in RHD [18], aortic regurgitation (AR) and stenosis (AS) are also common. Due to substantial differences in age, disease pathology, and socioeconomic circumstances of the patients, treatment guidelines established for valve replacements in older populations with degenerative aortic valve disease are not relevant in patients with RHD [19,20]. RHD may also worsen during pregnancy, further complicating treatment plans for pregnant patients [14,15].

Valve lesion severity determines the type and level of intervention and is classified by measuring various parameters including valve area, regurgitant fraction, leaflet thickness, and pressure gradient across the valve [7,21]. It is important to characterize the severity of RHD on a patient-specific basis. In patients with AR and AS, eccentric hypertrophy of the left ventricle often develops [7,22] and in severe cases, dilation of the ascending aorta may occur [22]. To our knowledge, the hemodynamic effects of AR and AS in patients with RHD on the ascending aorta have not been characterized, although Swanson et al. recently published a computational fluid dynamics (CFD) pipeline using echocardiography for limited-resource settings [23]. Using CFD to study hemodynamics is a well-established approach that has the potential to improve risk assessment and guide surgical intervention in patients with cardiovascular disease [24].

A common assumption in CFD simulations is a generic velocity inlet profile (e.g., plug, parabolic, Womersley) [25–27]. This assumption might be suitable for vessels with simple geometries, but complex flow is often observed in the thoracic aorta, particularly in the ascending region [28,29]. Morbiducci et al. quantified these flow differences between simulations using subject-specific MRI, plug profile, and an inlet flow extension, demonstrating inaccuracies from the idealized flow in shear and bulk flow structures [30]. Regarding the Womersley profiles, studies have shown that due to vessel wall deformation and aortic valve motion, these profiles may be poor estimates in this region of the aorta [31]. Recent studies investigating the effect of valve lesions on the ascending thoracic aorta have not specifically quantified the differences between subject-specific velocity profiles and generic profiles [32–34]. Rather, these studies focus on differences between diseased and nondiseased states. There are also currently no studies specifically modeling subjects with aortic lesions caused by RHD, a condition that presents differently when compared to other forms of aortic stenosis and regurgitation [6,22].

Many studies of RHD have examined the effects of aortic valve lesions on the left ventricle, but have not extensively studied the effects on the thoracic aorta or investigated any trends between aortic biomechanics and the aortic valve, despite evidence that some patients do have mild aortic dilation [22]. Therefore, we present a study simulating hemodynamic effects of RHD using CFD and subject-specific velocity profiles with previously collected CMR data from Groote Schuur Hospital in Cape Town, South Africa. We hypothesized that there would be measured differences in estimated WSS, oscillatory shear index (OSI), and relative residence time (RRT) between subjects with RHD and non-RHD in the thoracic aorta. Furthermore, we aimed to quantify the differences between using subject-specific velocity profiles and Womersley velocity profiles. We also included a brief morphological evaluation that estimated the angle between the left ventricle and the ascending aorta. This metric may provide additional information on the effect of the aorta in RHD. By using half-Fourier acquisition single-shot turbo spin echo (HASTE) and phase contrast magnetic resonance imaging (PC-MRI) data, we demonstrated the feasibility of using routinely collected image data to gain a more comprehensive understanding of the biomechanical effects of RHD (Fig. 1).

We used previously collected CMR data from 12 subjects (power analysis: $α=0.05, β=0.15$⁠) at the Cape Universities Body Imaging Center (CUBIC, Faculty of Health Sciences, University of Cape Town, Cape Town, South Africa) that we selected from 43 subjects with RHD and 28 control subjects without RHD. Written informed consent was obtained from all subjects, which was approved by the Human Research Ethics Committee in the Faculty of Health Sciences at the University of Cape Town. The study was carried out in accordance with the Declaration of Helsinki, International Code on Harmonization—Good Clinical Practice and South African Good Clinical Practice 2006 guidelines. The subjects were both male and female, ages 21 to 44 yrs, and have varying BMIs and ethnicities. Each control subject is age- (±3 yrs) and sex-matched. Subjects with RHD had mixed mitral and aortic valve disease with stenosis and regurgitation ranging in severity on a scale of one to three (one = mild, two = moderate, three = severe) as noted in Table 1. Severity was determined using the European Society of Cardiology and American Heart Association Guidelines [35], using parameters including valve area, aortic jet velocity, mean pressure gradient, and pulmonary artery pressure.

We used previously acquired CMR images from a 3T Siemens Skyra magnet at CUBIC and pulse sequences with parameters highlighted in Table 2. The PC-MRI datasets were used to obtain the subject-specific velocity profiles and were processed using open-source software (Medviso Segment,²) [36]. To obtain geometry for our simulations, we used transverse and coronal HASTE scans, which can be quickly obtained for each subject. Due to the large slice thickness of 8 mm, it was not possible to accurately segment the aorta using one acquisition type. We therefore aligned and overlayed the transverse and coronal HASTE images, using a custom matlab code (MATLAB 2018a, The MathWorks, Inc., Natick, MA), to create an image with a new voxel size of 1.5625 mm³. We then used ITK-SNAP³ [37] to manually segment the aorta beginning a centimeter above the aortic valve and following the vessel around the arch and into the descending thoracic aorta, while also including the brachiocephalic, left common carotid, and left subclavian arteries. The segmentations were inspected for errors and exported as steriolithography files for the computational analysis. We used three-chamber view cine CMR scans of the left ventricle to estimate the angles between the ascending thoracic aorta (ATA) and left ventricle in two-dimensional (2D).

To analyze the cine data, we used open source software, horos (Version 3.3.6) and the angle measurement tool. We based our method on a recent study by Roule et al. which looked at the angles between the ATA and the left ventricular outflow tract (LVOT) and left ventricle [38]. We estimated the planes from sagittal (three-chamber view) cine data. The ATA and LVOT planes were drawn parallel to the vessel walls, while the left ventricle plane was drawn from the apex to the center of the mitral valve. We performed this analysis on the 12 subjects in the study cohort as well as the entire CMR data collection (RHD: n = 43, control: n = 28).

We imported steriolithography files of the vessel geometries into simvascular (Version 2020.04.06) [39] and performed smoothing functions consistently across all 12 vessels to remove sharp edges and segmentation artifacts. After smoothing, we created a mesh using the meshing kernel, TetGen, which uses an adaptive meshing approach [39]. We conducted a mesh independence study to ensure that our velocity, wall shear stress, oscillatory shear index, and relative residence time values were not influenced by the mesh resolution. We used our case with the highest peak velocity to determine the optimal maximum mesh element size so that we could confidently apply that size across all subjects. Figure 2 below shows average velocity values of three probe locations in the ascending thoracic aorta over one cardiac cycle using six element sizes ranging from 0.07 to 0.27 cm. Based on these results, we concluded that a global maximum element size of 0.10 cm safely lies within a 5% margin of the average values computed using the finest mesh, without substantially increasing computational costs. We used meshing for the simulations with element sizes that do not exceed the chosen 0.10 cm. Furthermore, we incorporated boundary layer meshing along the wall with three layers to account for the higher velocity gradient in that region. The number of nodes (445,390 ± 141,977) and elements (2,663,721 ± 867,690) for each subject is shown in Supplemental Table 1, available in the Supplemental Materials on the ASME Digital Collection.

We found the mean Womersley numbers to be 16.5  ±  3.5 and $15.8±3$ at the inlet for RHD and control groups, respectively. We prescribed a flowrate over one cardiac cycle for each case using the mean values obtained from the PC-MRI analysis. This subject-specific flowrate waveform was then mapped to the inlet using the spatial velocity profile corresponding to the estimated Womersley number.

For the second set of simulations of each case, we applied a subject-specific velocity inlet profile obtained from our PC-MRI analysis. We kept the rest of the simulation parameters consistent, but instead of using mean flowrate data and Womersley spatial velocity distribution, we used subject-specific distribution. We aligned the inlet nodes to the inlets extracted from the PC-MRI data and assigned each node the corresponding PC-MRI velocity value at 50 time steps over one cardiac cycle. We took the Medviso Segment analysis mentioned previously and imported the results into a custom code in matlab (Version 2018a, The MathWorks, Inc., Natick, MA). This allowed us to replace the values in the Womersley profile distribution file with new PC-MRI values.

where $ΔT$ represents the time from peak systole to end diastole, Q_max and Q_min are maximum and minimum flow rates, and P_s and P_d are systolic and diastolic pressures. Mean arterial pressure, MAP, divided by the mean inlet flowrate, Q_in, was used to estimate resistance. simvascular was used to apply an automated flow split for the resistance and compliance values between all outlets. Mean outlet boundary condition values are listed in Table 3 for both groups.

As an initial pressure, we used estimated MAP values based on age and BMI from literature [40]. simvascular solves the Navier–Stokes equations and automatically applies the traction (no-slip) boundary condition on the surface of the vessel lumen. In addition, we assumed the blood to be a laminar Newtonian fluid with constant density (1.06 g/cm³) and viscosity (0.04 poise) using values from literature [39].

We ran each simulation for five cycles and confirmed that the simulation results converged by the fourth cycle. Simulations for each case were run on the simvascular Supercomputing Gateway, which uses XSEDE resources [41]. Each simulation took between 1.5 and 5.5 h to run. Results of the last cardiac cycle are presented in this section.

To measure flow asymmetry, Flow_asym, within the inlet, we subtracted the maximum velocity coordinates, V_m, from the inlet centroid coordinates, I_c, and divided that distance by the maximum effective radius, R, of the inlet. We compared the inlet flow profile in subject-specific simulations to Womersley by calculating the spatial velocity difference. Mean $ΔVD$ is calculated by the absolute velocity difference between matching node locations (a_i and b_i), then taking the sum of those values and dividing by the total number of nodes, n.

To estimate the effects of aortic valve lesions on the ascending thoracic aorta, we quantified time-averaged wall shear stress (TAWSS), OSI, and RRT for RHD and control subjects, including simulations using either subject-specific or Womersley flow profiles. We analyzed a region of the ascending aorta of consistent size across vessels of each subject that incorporated where we found maximum hemodynamic parameter values due to the influence of the aortic valve. This region was placed approximately in the same location for all subjects, a few millimeters above the inlet.

Statistical analyses were conducted on each group. Mean and standard deviations are represented in the tables and figures below. We used QQ plots to test for normality, then compared group means using paired and unpaired t-tests. A p-value of <0.05 was considered significant. All statistical analyses were performed with commercial software, prism (version 9.2; GraphPad Software, Inc., Boston, MA).

We first examined the flow asymmetry of our subject-specific inlet velocity profile simulations. We did this by measuring the difference between the inlet centroid coordinates, I_c, and maximum velocity coordinates, V_m, at peak systole, as shown in Fig. 3(a). Comparing mean flow asymmetry values, we found no significant difference between RHD and control groups (Fig. 3(b)) although the RHD group mean is higher.

In Fig. 4, each subject-specific simulation inlet is visualized at peak systole in two different views. Observing individual subject results, we see that RHD subjects N2, N3, and N4 have the highest flow asymmetry among both groups. In addition, if we examine other metrics listed in Table 2, we see that, again, RHD subjects N2, N3, and N4 all have left ventricular ejection fractions (LVEFs) outside of the normal range of 50–75% [43]. These results confirm the idea that aortic valve lesions can affect both left ventricular function and aortic blood flow.

Observing the control group flow asymmetry, we qualitatively can conclude that velocity distribution is much more equal; however, we see that maximum velocity points are still not overlapping with the inlet centroid as they would in generic velocity profiles.

To quantify the spatial velocity profiles in a different way, we did a point-wise comparison by comparing each inlet node's velocity value at peak systole between Womersley and the subject-specific simulation parameters. Representative inlets in Fig. 5(a) show how the spatial velocity difference was calculated. In Fig. 5(b), we found that the mean $ΔSV$ of the RHD group was significantly higher (p < 0.05) than that of the control group. Qualitatively, we can see more variation between cases within the RHD group, in Fig. 5(c), with N3 having the highest difference.

We followed our inlet flow assessment by quantifying the effects of flow on the lumen wall. The objective for this analysis was to measure differences between the RHD and control groups, but also to measure differences between subject-specific and Womersley inlet velocity profile simulations. We summarized these values in Fig. 6 and Table 4, while values for each simulation case can be found in Supplemental Table 2 available in the Supplemental Materials on the ASME Digital Collection.

As mentioned previously, time-averaged wall shear stress provides an estimation of drag caused by blood flow on the wall of the vessel lumen during one cardiac cycle. This metric can be useful for understanding the effect of the inlet flow on the wall in the ascending region of the thoracic aorta. We plotted the mean TAWSS among both groups and both simulation types (i.e., subject-specific and Womersley). We found that mean TAWSS of the subject-specific inlet simulations was significantly higher than the subject-specific inlet control group (p < 0.05) as well as the Womersley RHD group (p < 0.05, Fig. 3). Due to the aortic stenosis and regurgitation present among the RHD subjects, we can assume that the resulting flow asymmetry would also increase TAWSS values as higher velocities are observed near the wall. Therefore, we can anticipate significantly higher TAWSS in RHD subjects with aortic valve lesions. However, the Womersley simulations fail to reveal these differences and even show significantly lower values when comparing to subject-specific simulations within the RHD group (See Supplemental Fig. 1 available in the Supplemental Materials). Therefore, these findings highlight the importance of imposing a subject-specific velocity distribution when analyzing the ascending region of the thoracic aorta.

We also analyzed oscillatory shear index and relative residence time index to further our understanding of the hemodynamic effects from the valvular lesions. We found no significant difference in mean OSI values between simulation type or between RHD subjects and controls. Although we initially expected that the mean OSI of the RHD group would be higher than that of the control group, the lack of difference may be explained by averaging regions of high and low OSI together [33,44]. We also did not find any significant differences in RRT values between either RHD and control groups. RRT is a less important metric when it comes to RHD as it has been used as an indicator for atherosclerosis susceptibility and thrombosis formation [27,45]. However, we felt it was an important metric to compare simulation results between subject-specific inlet profiles and Womersley profiles for potential application to other aortopathies. We anticipated that RRT would be significantly lower in the RHD group because of the increased velocity near the walls. However, we did find a significant difference in mean RRT between subject-specific inlet simulations and Womersley inlet simulations (p < 0.01) among the RHD group. This finding demonstrates the substantial influence inlet boundary conditions can have on hemodynamic metrics. These differences are also demonstrated when we examine individual subjects. We see that from the subject-specific simulations, RHD subjects N2 and N3 have the lowest RRT, whereas the Womersley simulations do not show that N2 and N3 have the lowest values. Taken together, TAWSS, OSI, and RRT can all provide insights into the effect of valve lesions on the ascending aortic wall when using subject-specific velocity inlets.

To supplement our CFD analysis and characterize the relationship between the ATA, aortic valve, and left ventricle, we quantified the angles between them. Figures 7(a) and 7(c) show representative images of how the angles were calculated with the summarized data in Figs. 7(b) and 7(d). We found that the angle between the left ventricle and LVOT is significantly larger in RHD patients for both the total subjects and the study cohort (p < 0.01; p < 0.05). These findings are similar to Roule et al. which found that patients with AR had significantly greater ATA-LV angles and noted that patients with moderate to severe AR following transcatheter aortic valve replacement (TAVR) have increased one-year mortality and heart failure [38]. Despite their significant findings, we found it more relevant to measure the LVOT-LV angle since flow from the LV is first affected by the LVOT. To investigate further within the RHD group (n = 43), we divided the subjects into groups with moderate to severe mitral regurgitation (MR), mitral stenosis (MS), AR, and AS (⁠$≥2$⁠) and the resulting groups of mild (<2) or none present. No significant differences were found between MS severity nor AS severity groups. We also did not find a significant difference between AR severity (p = 0.21), but did find that patients with moderate to severe MR have significantly larger LVOT-LV angles than those with mild or no MR (p < 0.05). These angle differences may be important when considering TAVR for RHD patients and the performance of the valves. One recent study also considered the effect of aortic arch geometry on prosthetic valves [46]. They evaluated aortic pressure differences post-TAVR with varying aortic geometries and found elevated pressure differences in those with greater curvatures, which may compromise valve performance. All of these findings suggest that rheumatic heart disease has a complex effect on the cardiovascular system, with valve lesions that affect both left ventricular function and aortic hemodynamics.

There are some limitations to this study that should be noted. One limitation is the rigid-wall assumption in the simulations. Because we aimed to create a feasible CFD method through quick and affordable image acquisitions, we were not able to obtain wall deformation information. Furthermore, incorporating wall deformation would increase computational costs. However, modeling the vessel wall using a hyperelastic model, such as the one derived by Raghavan and Vorp [47], could allow the wall to expand and contract with pressure changes improving our estimation of hemodynamic metrics near the wall including TAWSS and OSI. These metrics may be overestimated when using rigid walls [48]. Since the effect of rigid walls is consistent among simulations and we were focusing mainly on comparisons between groups, the actual values are of less importance. Another limitation of this study is that the custom code to impose the subject-specific inlet boundary conditions does not account for vessel translation and is mapped to a stationary inlet. Although we ensured that the PC-MRI data were adequately overlayed at maximum systole, accounting for vessel translation throughout the cardiac cycle would further improve simulation accuracy, particularly in the diastolic phase [31]. One study even incorporated valve geometry into the model and compared the plug velocity profile to subject-specific profiles in the LVOT. The authors did not observe significant differences in WSS, arguing that subject-specific profiles are not needed [49]. However, without access to subject-specific valve geometry and software to model it, their presented method is not feasible. Moreover, our WSS findings do not agree with theirs as we did find significant differences between methods, which suggests that subject-specific is important.

Based on the results of our study and the limitations, there are key future works to be mentioned. First, accounting for vessel translation within the code may improve the accuracy of simulations using subject-specific inlets. In addition, the method would be more accessible and clinically translatable if subject-specific velocity imposition could be directly incorporated within the simvascular software. To further validate our CFD approach using 2D PC-MRI, it could be useful to simultaneously collect 4D flow MRI. In addition to having lower temporal resolution [50], it is less clinically translatable due to system requirements, scan times, and advanced techniques, but could be used in a validation study to compare aortic flow downstream to CFD results from our current method.

Outside of these small improvements, extending this method to other diseases may provide valuable insights. Specifically, there is a need to determine whether patients with bicuspid aortic valve (BAV) develop dilated ascending thoracic aortas because of genetic causes or abnormal hemodynamic causes [51]. Patients with BAV often have regions of high WSS due to the asymmetric blood flow caused by the aortic valve [52] and may even develop ascending thoracic aortic aneurysms. Nightingale et al. recently showed strong evidence that aortic biomechanics are influenced by hemodynamics and that regions of elevated WSS have low stretch ranges in the vessel wall tissue [53], yet we do not see dilation of the ascending aorta in RHD patients as we do in BAV, despite having elevated WSS regions. Comparing RHD patients with similar aortic stenoses and WSS patterns may help elucidate the level of biomechanical effect the hemodynamics are having on aortic wall changes. A second direction could be to study the progression of RHD, which may help reveal any geometrical changes occurring in the LV or aorta, which may adversely effect treatment outcomes. The method presented here could be applied to a small animal study using a mouse model such as the one established by Galvin et al. [54], acquiring CMR data at multiple time points.

The hemodynamics of patients with aortic valve lesions due to RHD differ depending on the severity of stenosis and regurgitation, function and angle of left ventricle, and potentially other biomechanical factors. It is imperative to use subject-specific inlet velocities to model flow in the ascending thoracic aorta in patients with altered aortic valve function. Imposing subject-specific inlets revealed a significant increase in TAWSS and a decrease in RRT, while imposing a Womersley profile did not reveal these differences. We supplemented this hemodynamic analysis by providing simple 2D angle estimations of the ATA, LVOT, and left ventricle that revealed significant differences between RHD and non-RHD subjects, as well as differences between MR severity groups. Exposing differences in these metrics may help to stratify risk and personalize pre-operative information for RHD cases. Furthermore, we demonstrated how to accomplish this CFD approach through routinely collected CMR scans at a relatively low cost and faster acquisition, expanding the availability of this type of analysis for more patients.

• National Science Foundation (NSF) (Grant No. 1829436; Funder ID: 10.13039/100000001).

• South African Medical Research Council (Funder ID: 10.13039/501100001322).

• Medical Research Council UK (Funder ID: 10.13039/501100000265).

• Lily and Ernst Hausmann Trust.

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

The authors thank Alejandro Quintero Rendon from Universidad Santiago de-Cali for assisting with image segmentation and Joseph Muskat from Purdue University for computational guidance. We also acknowledge the funding support for HLC provided by IIE-GIRE. We also thank the patients with RHD and healthy control subjects who volunteered to participate in this study.

AR =

aortic regurgitation

AS =

aortic stenosis

ATA =

ascending thoracic aorta

CFD =

computational fluid dynamics

CMR =

cardiovascular magnetic resonance

HASTE =

half-Fourier acquisition single-shot turbo spin echo

LVEF =

left ventricular ejection fraction

LVOT =

left ventricular outflow tract

MR =

mitral regurgitation

MS =

mitral stenosis

OSI =

oscillatory shear index

PC-MRI =

phase contrast magnetic resonance imaging

RHD =

rheumatic heart disease

TAVR =

transcatheter aortic valve replacement

TAWSS =

time-averaged wall shear stress

ΔSV =

spatial velocity difference