Patient-specific aortic geometries derived from MRI angiograms were obtained from the Vascular Model Repository (VMR), an open source database developed for clinically-motivated CFD investigations (Pfaller et al., 2022). We selected 12 patient aortas for simulation: six from patients after CoA repair REEA (“Repair” cohort, Figure 1, top), and six from age- and sex-matched control patients with healthy, anatomically normal aortas (“Control” cohort, Figure 1, middle). Each geometry consisted of the ascending aorta, the aortic arch/associated branch vessels (brachiocephalic, left common carotid and left subclavian arteries) and the descending aorta, terminating proximal to the iliac bifurcation. In order to only compare patients with similar aortic anatomies, we excluded patients who either 1) underwent CoA repairs other than REEA or 2) had collateralized arterial networks between the arch vessels and descending aorta, (e.g., Patient #3 in this sampling of different CoA variants (LaDisa et al., 2011b)). Since all data were anonymous and open source, this study did not require IRB approval or informed consent.

The aortic geometries of all 12 study patients were downloaded from the VMR database as stereolithography (STL) files. We will refer to these unmodified, anatomically-accurate, patient-specific geometries as the “baseline” geometries/models for the remainder of this study. Since we also wanted to model the interaction between aortic endothelial shear stress and the restenosis process, we then modified each baseline geometry to have progressively severe concentric stenoses of 10%, 50%, and 80% (Figure 1, bottom) at the CoA repair site (for the Repair cohort) or at the anatomically corresponding location immediately distal to the left subclavian artery where re-stenosis is most often seen (for the Control cohort). This resulted in a total of 48 aortic geometries for this study.

Geometry modifications were performed with the open source packages morphMan (version 1.4, Simula Research Lab, Norway) (Bergersen et al., 2020) and Blender (version 2.7.9, Amsterdam, Netherlands) where Voronoi diagrams derived from the discretized surface mesh were created to represent the geometry in terms of centerlines, maximum inscribed spheres, and their radius. Then, the inscribed spheres and the corresponding radius were scaled relative to the centerlines in the regions of interest to circumferentially “pinch” the aortic wall, thereby recreating a restenosis in the descending thoracic aorta distal to the left subclavian artery (Figure 1). Consistent with previous definitions (Ou et al., 2004; Donazzan et al., 2014; Ou et al., 2008), the severity of the restenosis lesion σ was defined as: σ = 1 − R stenosis R baseline (1) where R stenosis is the aortic diameter at the most severe point of stenosis, and R baseline is the corresponding baseline aortic diameter at that location. Modified geometries (in both cohorts) will be referred to by the degree of stenosis (sigma) introduced, as defined above in Equation 1.

Fluid simulations were performed with HARVEY (Peters Randles et al., 2013), a massively parallel CFD code that uses the lattice Boltzmann method (LBM) to solve the Navier-Stokes equations of fluid flow. A comprehensive description of the LBM methodology applied to fluid dynamics can be found elsewhere (Succi, 2003), and HARVEY’s development and validation for simulating cardiovascular flow have been previously described (Randles et al., 2015a; Randles et al., 2015b; Feiger et al., 2019; Vardhan et al., 2024). Briefly, instead of directly solving for fluid velocity and pressure, LBM-based solvers such as HARVEY take an alternative approach of modeling fluid as a probability distribution f ( x , t ) of particles within a discrete 3D Cartesian lattice. The temporospatial evolution of this particle distribution over a specified time interval Δ t is governed by collision interactions described by the lattice Boltzmann equation (Equation 2), where Ω represents the collision operator and f i e q ( x , t ) represents the particle distribution at equilibrium: f i x + c i Δ t , t + Δ t = f i x , t − Ω f i x , t − f i e q x , t (2)

HARVEY uses a standard D3Q19 velocity discretization pattern and a single-relaxation Bhatnager-Gross-Krook (BGK) collision operator. Due to the high shear rate ( > 100 s − 1 ) present in large arteries (Alexy et al., 2022), the blood is modeled as an incompressible Newtonian fluid with constant density ρ = 1,060 k g / m 3 and dynamic viscosity μ = 0.004 P a ⋅ s , with a no-slip boundary condition enforced at the fluid-wall interface using a halfway bounce-back method. All vessel walls are assumed to be rigid, which is commonly done in fluid simulations of large vessel blood flow (Goubergrits et al., 2015; LaDisa et al., 2011a; LaDisa et al., 2011b).

We wanted our simulations to best isolate the effect of vessel anatomy and CoA lesion severity σ on aortic hemodynamics. Consequently, to minimize the risk of confounding the results by implementing different boundary conditions for each simulation to achieve a physiological state selected a priori, we chose to apply simple and consistent boundary conditions to all simulations. To capture time-averaged hemodynamic metrics such as TAWSS and OSI, pulsatile blood flow at the aortic inlet was modeled as a time-varying velocity waveform with a period T of 0.75 s (corresponding to a heart rate of 80 beats per minute), a maximum systolic velocity of 0.45 m/s, and a parabolic wave profile consistent with previous studies using the same repository (LaDisa et al., 2011a; LaDisa et al., 2011b). The outlets were modeled with 0-pressure boundary conditions. Spatial convergence studies were performed at grid spacings of 15, 20, 25, 50, and 75 microns, with convergence results seen for all simulations at a resolution of 25 microns. Specifically, maximum values of TAWSS and OSI were monitored at the CoA region for incrementally decreasing grid spacings until ≤3% change was observed which was the case for all the simulations done at grid spacings <50 microns. Temporal convergence studies were performed for a duration of seven cardiac cycles, with simulation stability observed after the second cardiac cycle.

The primary outcomes of interest were endothelial WSS τ ( x , t ) , TAWSS (Equation 3), and OSI (Equation 4): T A W S S = 1 T ∫ 0 T | τ t | d t (3) O S I = 1 2 1 − ∫ 0 T τ t d t ∫ 0 T | τ t | d t (4)

These metrics were chosen based upon previous literature linking TAWSS/OSI to abnormal aortic remodeling and to allow direct comparison with the results of previous CFD investigations of CoA.

To facilitate comparing simulation results from aortas of different lengths, diameters, and CoA severities σ , all data were spatially normalized prior to analysis. First, every point on the wall of the vessel of each aorta was re-parameterized from Cartesian to cylindrical coordinates that were defined relative to the aortic centerline (Figure 2). For a given vessel wall point P ( r , θ , l ) , r signifies the shortest Cartesian distance to the vessel centerline, θ represents the circumferential position around the vessel wall (with θ = 0 and θ = ± π corresponding to the inner and outer curves of the aorta, respectively), and l denotes the normalized position along the vessel centerline, with the aortic inlet set to be l = 1 and the distal aortic outlet set to be l = 0 . This coordinate transformation, which has been previously employed by other researchers (Perinajová et al., 2021), allows each aortic wall to be unwrapped and assigned to a normalized rectangular two-dimensional domain in θ and l . In turn, this allows us to make quantitative comparisons of the magnitude and spatial distribution of shear metrics within/between cohorts and coarctation severities.

The initial analysis consisted of calculating the circumferentially averaged (i.e., across 0 ≤ θ ≤ 2 π ) TAWSS/OSI for each patient. Baseline hemodynamics for the entire cohort were calculated as median [interquartile range] and plotted against the normalized position of the aortic centerline l. The aorta was subdivided into five anatomic regions of interest: the ascending aorta (AscAo), the aortic arch, the coarctation repair site (CoA), the proximal descending aorta (pDA) and the distal descending aorta (dDA). Statistical comparisons in each region were performed with the Kolmogorov-Smirnov test for non-normal continuous distributions (Sheskin, 2003). Inter-patient variation in the spatial distribution of TAWSS/OSI was quantified by calculating the L2 norm of the difference in the magnitude of the shear metric, which was then plotted on heat maps.

Since no two patients have anatomically identical aortas, we expected some level of inter-patient variation in shear metrics, even between baseline patient models in the same cohort. Furthermore, these differences would carry over into the σ = 10%, 50%, and 80% variants generated from each baseline geometry. However, we needed to determine that any observed differences in TAWSS/OSI magnitude and spatial distribution between cohorts or different levels of σ were not due to these chance variations in aortic anatomy. Therefore, we needed to account for the variability introduced by these shared features and the interactions between anatomy and hemodynamics.

We did this by using a linear mixed effect (LME) model that incorporates random effects to account for the inherent variation between different patient anatomies (Faraway, 2006). The LME model was employed with the lme4 package in R (Bates et al., 2015). The fixed effect coefficient, as determined from the TAWSS and OSI results, captured the effect of group (control vs repaired) and severity (degree of re-coarctation). The baseline variability was then taken into account by a random intercept. Additionally, the random-effects intercept captures correlations among all measurements derived from a common geometry and can thus be interpreted as a geometry-specific intercept adjustment in the fixed-effects parameter.

We evaluated these factors at each of the five previously defined anatomic regions of the aorta (AscAo, arch, CoA, pDA, and dDA). In addition, we considered angular locations, including the inner, left, outer, and right sides of the aorta. The significant coefficients of the model for the group, severity and their interaction were visualized in Figure 10, with the regions where any of the coefficients were significant highlighted in dark gray.

Simulations showed clear differences in aortic hemodynamics between healthy and repaired patients. The median flow rates in the descending aorta of each patient were 70% of the inlet flow rate, representing a physiologic flow split in all investigated geometries. Healthy patients exhibited a uniform distribution of TAWSS in θ and l, which smoothly increased from a median of 1.07 Pa in the AscAo to 2.09 Pa in the dDA (Figure 3, top), due to the slightly narrower diameter of the distal aorta. For repaired patients, the median TAWSS was similar in the AscAo (1.20 Pa) but significantly higher in both the aortic arch (3.46 vs 1.24 Pa, p < 0.05) and the CoA repair site (4.34 vs 1.56 Pa, p < 0.05, Figure 3, top). This finding is reinforced by heat maps showing that the magnitude of TAWSS is largely isotropic between control patients (Figure 3, bottom left), but in repaired patients it is clearly elevated around the aortic arch and the outer curvature of the coarctation repair site (Figure 3, bottom middle). The elevation of TAWSS along the outer curvature of the aortic and pDA indicates that most repaired aortas exhibit flow separation, with the primary blood flow jet impinging on the distal descending aortic wall, leaving a region of recirculating flow along the inner aortic curvature. This finding is further demonstrated by calculating and plotting the ratio of the magnitude of TAWSS between repaired and control patients (Figure 3, bottom right).

OSI distributions are consistent with these findings, with control patients showing a smooth, uniform OSI distribution centered around θ = 0 , again indicating that normal aortic flow hugs the inner vessel curve from the inlet to the distal aortic outlet (Figure 4, left). In contrast, repaired patients have regions of high OSI centered around θ = 0 along the arch and CoA region (Figure 4, bottom center). Together with the TAWSS findings, this pattern of OSI distribution indicates regions of recirculation and low flow along the inner aortic curvature of repaired CoA patients, not strongly pulsatile flow.

We observed a clear, non-linear relationship between CoA severity σ and shear stress magnitude at the CoA site for all investigated geometries (Figure 5). Although TAWSS changed little when increasing σ from 0% to 10%, there was a significant increase from σ = 10 % to σ = 50 % (control, 2.17 → 12.19 Pa; repaired, 7.64 → 24.28 Pa), and even more dramatic from σ = 50 % to σ = 80 % (control, 12.19 → 27.99 Pa; repaired, 24.28 → 42.70 Pa). The nonlinear increase for peak systolic WSS was even greater (control: 8.33 → 45.81 → 95.89 Pa; repaired: 27.83 → 87.25 → 154.26 Pa).

Compared to TAWSS, the relationship between σ and OSI magnitude is more complex (Figure 6). Since OSI can only be calculated over the course of the cardiac cycle, we report the maximum and median OSI within the CoA region. For control patients, increasing σ led to a region of progressively higher OSI in the inner curvature of the AscAo and the outer curvature of the pDA (Figure 7). However, there is no correspondingly clear trend for repaired patients, who demonstrate a weaker trend towards increased aortic arch OSI.

TAWSS distributions in the repaired patient cohort showed considerably more variation than in control patients. For example, the TAWSS IQR was nearly an order of magnitude higher in repaired patients than in control patients (2.92 vs 0.32 Pa, aortic arch, and 1.82 vs 0.20 Pa, repair site). This indicates that a wide range of outcomes are possible after CoA repair, while TAWSS distributions in healthy patients remain fairly uniform.

To further examine this, we quantified intra-cohort variance by calculating the L2 norm of the difference in TAWSS and OSI magnitude between each pair of geometries within each cohort (Figure 8). From this, it is immediately apparent that the hemodynamics between the control geometries are much more similar than between the repaired patients. This technique could be used to potentially identify different “hemodynamic phenotypes” following coarctation repair. For example, in the repaired cohort, the greatest pairwise differences in the TAWSS distribution are between geometries R0, R2, and R3. Plotting these geometries shows clear differences in the TAWSS distribution pattern: R0 shows a high TAWSS mainly around the repair site, while R2 is centered around the narrowed aortic arch, and R3 is qualitatively uniform, although with a slight residual stenosis at the repair site (Figure 9).

The implemented LME model distinguished the relative importance of severity or residual stenosis in each studied cohort confirming that aortic coarctation, even after surgical repair, significantly disturbs the spatial distribution of TAWSS and OSI compared to healthy controls. Notably, these differences persist across various regions of the aorta and are not merely a function of the severity of re-coarctation, which is a separate, significant disruptor of TAWSS/OSI patterns in its own right.

The LME model revealed that the group effect (repaired versus control) remained significant in specific regions (18%–24% of the luminal surface of the aortic, regardless of the severity of the re-coarctation. This consistent pattern further highlights that repaired patients consistently experience abnormal hemodynamics. Furthermore, the interaction between group and severity was also significant in most regions, suggesting that the degree of re-coarctation exacerbates these differences but does not alter the fundamental disparity between the groups. Specifically, TAWSS was significantly affected in the AscAo and arch, particularly towards the outer and left regions (Figure 10). Additionally, significant TAWSS patterns were observed on the inner side of the pDA, extending to the inner and right side of the dDA. These findings indicate regions where repaired CoA patients are likely to experience increased shear stress, which can contribute to adverse remodeling and complications. In contrast, OSI exhibited a different pattern of significance. It was primarily affected in the AscAo, particularly on the right outer side. Sparse significant patterns of OSI were also observed in the arch and CoA region, mainly towards the outer and left side of the aorta (Figure 10). These results suggest that the oscillatory nature of shear stress in these regions may play a role in pathological remodeling processes distinct from those driven by TAWSS. In summary, our findings indicate that repaired CoA patients continue to experience adverse hemodynamics compared to healthy controls. The patterns of TAWSS and OSI distribution are distinct between the two groups and are significantly influenced by both the presence of re-coarctation and its severity. These results underscore the importance of considering both the spatial distribution of hemodynamic metrics and the impact of residual or recurrent stenosis when assessing long-term outcomes in CoA patients. The differences in TAWSS and OSI patterns across the aorta may provide insights into the mechanisms driving pathological remodeling and restenosis, especially if correlated with surveillance imaging findings and long-term clinical outcomes. This highlights the potential for improving patient outcomes by integrating CFD simulations and metrics into the routine clinical management of complex, chronic cardiovascular disease such as CoA.

The present study should be interpreted in relation to several limitations. This study provides valuable information by using open source patient geometries from the VMR, enhancing reproducibility and accessibility for future research. However, the VMR dataset does not include follow-up imaging or long-term patient outcomes, such as the development of restenosis over time. Future investigations could build on this work by using clinical data sets that include long-term follow-up data and multiple imaging time points, enabling simulations of evolving patient-specific hemodynamics.

Our use of 0-pressure outlet boundary conditions was specifically chosen to isolate the effects of vascular morphology on flow patterns and reduce confounding variables. Although this approach ensures consistency and comparability between patient cohorts, future applications in clinical settings would benefit from tunable boundary conditions to account for downstream vascular adaptation and variability, particularly in cases of severe CoA stenosis, where patient-specific boundary conditions have a more significant effect (Mariotti et al., 2023). For example, incorporating lumped parameter models, such as Windkessel networks, could better capture the dynamics of the vascular response, although these adjustments may introduce additional variability in flow predictions - reported to differ by up to 96% compared to clinical MRI data (Nicole Antonuccio et al., 2021). To balance these considerations, our use of standardized boundary conditions aligns with established practices in computational modeling (LaDisa et al., 2011a), providing a robust foundation for comparative analyses while minimizing potential confounding factors.

Recent comparisons of rigid vs deformable-wall CFD simulations in the ascending thoracic (Calò et al., 2023) and abdominal aorta (Peng et al., 2023) have found that rigid-wall models tend to overestimate aortic TAWSS and underestimate OSI, particularly in the compliant regions of the aorta. However, the site of CoA repair following REEA is known to be significantly less elastic and more rigid than normal aortic tissue (Verhaaren et al., 2001). This suggests that a deformable-wall model would yield even stronger differences, further supporting the robustness of our conclusions. Additionally, deformable-wall models of CoA after EEA also conclude that the presence of residual stenosis is the most significant factor in causing residual hemodynamic abnormalities (Taelman et al., 2016).

Finally, the pronounced variation in TAWSS and OSI patterns observed among repaired patients, compared to controls, appears to stem not only from residual anomalies at the repair site but also from anatomical variability. Factors such as arch angulation and other morphological anomalies, which are absent in the control cohort, significantly influence postoperative shear stress outcomes (Oliver et al., 2009; Goubergrits et al., 2015).