Five patients with a normo-functional BAV were recruited prospectively at a single university hospital from cardiology clinics, cardiac MRI, and echocardiography sessions. Detailed demographics of the enrolled patients are reported in Table 1. 2D cine-MRI images were used to analyze valve function and measure aortic diameters on double oblique planes positioned along the aortic root. Aortic diameter was measured next to the right pulmonary artery. This measurement was repeated at 3-year follow-up to assess potential aortic remodeling. Ten age-matched healthy volunteers (HV, female males, age of 23 ± 7) were enrolled as control group and were free from any known significant medical problem. Both BAV patients and HVs underwent 4D Flow acquisitions. The outcome of follow-up measurements on BAV patients was not known when 4D Flow data were processed. The study complies with the declaration of Helsinki and was approved by the West Berkshire ethics committee, United Kingdom. All participants gave written informed consent.

4D flow acquisitions were performed at John Radcliffe Hospital (Oxford, United Kingdom) on a 3.0 T MR system (Trio, Siemens, Erlangen, Germany). Flow-sensitive gradient-echo pulse sequence MRI was used to characterize and quantify blood flow in the thoracic aorta. Datasets were acquired with prospective ECG-gating during free-breathing, using a respiratory navigator. The volume of acquisition was oriented along an oblique-sagittal plane encompassing the ascending aorta, the aortic arch, and the thoracic aorta. 4D Flow sequences were performed with the following parameters: echo time = 2.29–2.63 ms, repetition time = 38.4–41.6 ms, flip angle = 7°, voxel size = 1.46–1.67 × 1.46–1.67 × 1.8–2.2 mm³, temporal resolution = 50.72–76.91 ms. The velocity-encoding range (VENC) was properly determined using the lowest non-aliasing velocity on scout measurements, resulting in 1.3–1.5 m/s and 1.5–1.9 m/s for HVs and BAV patients, respectively.

Through ad-hoc in-house software implemented in MATLAB (The MathWorks Inc., Natick, MA, United States), 4D Flow data were analyzed through a semi-automated pipeline that consisted of five steps.

The first step consisted of 4D flow velocity data encoding and pre-filtering. Phase-contrast images of each encoded flow direction (Figure 1A) were scaled by the corresponding VENC value to obtain the Cartesian velocity components (Markl et al., 2011). Velocity components were automatically processed to correct from eddy current errors, velocity aliasing and to decrease MR noise through median filtering (Bock et al., 2007). The Euclidean norm of the 3D velocity vectors was then computed based on the Cartesian components for all acquired time frames to obtain the magnitude intensity of the velocity vector (Figure 1B).

In the second step, through visual inspection of the dataset, the time-frame with the highest velocity-to-noise ratio was identified as the most representative of peak systole. At the selected time-frame, four landmarks were manually positioned inside the acquired volume at locations corresponding to the left ventricle outflow tract (LVOT), right pulmonary artery (PA), brachiocephalic artery (BrA), and the most distal visible point of the descending aorta (DeA; Figure 1C). Moreover, the 3D geometry of the aortic lumen was defined. The boundaries of the lumen were manually traced to define the region of interest (ROI) on each para-sagittal slice of the 4D Flow volume (Figure 1D). Tracing resulted in a 2D binary mask associating a value equal to 1 to the voxels inside the traced boundary, and equal to 0 otherwise. A 3D binary mask was obtained merging the 2D masks (Figure 1E). Subsequently, the aortic wall and the centerline of the aorta were automatically defined. The wall was identified as the closed iso-surface encompassing the borders of the 3D binary ROI, and it was smoothed through Laplacian diffusion filtering (Desbrun et al., 1999) to obtain a regularized surface mesh (Figure 1F). The centerline of the vessel wall was obtained through skeletonization of the ROI (Shih, 2010); the points defining the skeleton were ordered end-to-end via Dijkstra algorithm (Dijkstra, 1959), specifying the projections of LVOT and DeA landmarks on the skeleton as starting and ending points, respectively (Figure 1G). The resulting ordered points were fitted with non-uniform rational bi-cubic spline (NURBS) curve interpolation that was smoothed via 5th order moving average filtering to obtain the centerline of the vessel wall. PA and BrA landmarks were finally projected to the closest points on the centerline (Figure 1H).

In the third step, aorta cross-sections were defined along the vessel's centerline. Starting from LVOT, a set of equally spaced points P_i was automatically generated on the vessel centerline, prescribing a distance Δx between consecutive points P_i and P_i+1 equal to the average voxel dimension of the dataset (Figure 2A). Each point P_i and its corresponding unit vector x'₃ locally tangent to the centerline were used to automatically define a cross-section π_i passing through P_i and normal to x'₃ (Figure 2B). For each voxel whose center C_j was located within 0.5^*Δx from π_i, the projection (C'_j) of C_j was computed (Figure 2B), and the three velocity components (u_1j, u_2j, u_3j) measured at C_j were associated to C'_j. The discrete velocity values obtained over π_i were interpolated through cubic splines to increase the spatial resolution of data, and velocity components were transformed into the local reference system (x'₁, x'₂, x'₃; Figure 2C). In addition, for each point C'_B,l of π_i belonging to the boundary of the vessel wall Γ_π, a second local reference system was defined, which consisted of the local inward normal to the wall surface, assumed as the radial direction x'_r,l, the normal to the cross section x'₃, assumed as the longitudinal direction, and the cross product of x'_r,l by x'₃, assumed as the circumferential direction x'_t,l (Figure 2D).

In the fourth step, the aortic bulk flow was automatically quantified. Velocity data of each cross-section were processed to compute the variables of interest:

Vessel cross-sectional area A_i, as the surface of the piecewise linear polygon defined by the boundary points C'_B,l;

|V|_{i_max}, |V|_{i_mean}, α_πi and d_πi were computed at three consecutive time-frames centered in T_ps. The analysis was focused on two specific sectors of the ascending aorta: the proximal tract, running from LVOT to PA landmarks, and the distal tract, running from PA to BrA.

The values of the listed hemodynamic variables computed in BAV patients were compared off-line vs. those computed for HVs with non-parametric Mann–Whitney t-test, adopting a p = 0.05 as statistically significant, using GraphPad Prism 7 (GraphPad Software, Inc., La Jolla, CA, USA).

In the fifth step, the distribution of WSS vectors (WSS⃗) acting on the aortic endothelium was quantified. To this aim, the stresses tensor τ_ij was computed at the boundary of the fluid domain under the assumption of Newtonian blood rheological behavior (Baskurt and Meiselman, 2003):

where viscosity μ was set to 3.7cP (Baskurt and Meiselman, 2003), and spatial derivatives of velocities were computed numerically. In the local reference system x'_r,l, x'_t,l, x'₃, where x'_r,l is the local inward normal to the aortic wall, the 3D vector T was computed as:

The 3D shear stress vector WSS⃗ was then computed as:

The axial (WSS⃗Axial) and circumferential (WSS⃗Circ) components of WSS⃗ were finally obtained by projecting WSS⃗ on the local axes x'₃ and x'_t,l, respectively.

WSS⃗ was computed throughout the systolic phase, i.e., from the first time-frame to T_Es. As in a recent work of our group (Piatti et al., 2016), the numerical computation of the spatial derivatives of velocity was performed through two novel methods.

The first method, named Global Volumetric (GV), is based on the use of 3D Sobel filters to compute voxelwise the components of the strain rate tensor and yielded the distribution of WSS⃗GV for each point of the triangular elements of the surface representing the aortic wall. We showed previously that this method correctly captures the spatial trends in WSS⃗, despite underestimating pointwise values (Piatti et al., 2016). Additionally, the oscillations in time of WSS⃗GV were computed by means of the Oscillatory Shear Index (OSI), of the Relative Residence Time (RRT), and of the Transversal WSS (TransWSS), as defined in Table 2 (Soulis et al., 2011; Jahangiri et al., 2015; Mohamied et al., 2017). Specifically, the OSI quantifies changes in WSS⃗ direction during time, and ranges from 0.0 to 0.5 as WSS⃗ rotates by 0° to 180°. The OSI was computed also for the single components WSS⃗Axial and WSS⃗Circ, obtaining the two indexes OSI_Ax and OSI_Circ (Supplementary Material). RRT combines the effect of time-averaged WSS vector magnitude at a given region of the endothelium with its corresponding OSI value: RRT gets lower as WSS⃗ increases and OSI decreases, i.e., when locally the wall is loaded by fast and constant flow jets, and it gets greater as WSS⃗ decreases and OSI increases, i.e., when the wall is loaded by slow and highly time-dependent flow jets. TransWSS measures the component of the WSS⃗ that is perpendicular to the time-averaged WSS⃗, representative of the average direction of the flow. This metric is hence defined in a reference frame based on the flow direction and not on the geometry of the vessel wall.

The second method, named Local Planar (LP), was adopted to obtain the distribution of WSS⃗LP for specific locations of the vessel wall based on the velocity values over a single vessel cross-section. As demonstrated in Piatti et al. (2016), this approach allows for improving the pointwise quantification of WSS⃗. It is based on the application of 2D derivative filters; in Piatti et al. (2016), 2D Sobel filters were exploited to compute the components of the strain rate tensor, while in the present study B-Spline based filters (Hast, 2014) were used, owing to the improved results verified through preliminary benchmarking (unpublished data). WSS⃗LP was characterized in terms of magnitude and angle that quantifies the deviation of WSS⃗LP from the longitudinal, i.e., axial, direction of the vessel (for a 0° angle, the WSS vector is purely axial and points downstream from the considered cross-section; for 90° angle, the WSS vector is purely circumferential; for a 180° angle, the WSS vector is again purely axial, but points upstream of the considered cross-section).

The analysis of WSS⃗GV, OSI and RRT data was focused on the ascending aorta, i.e., on the tract of the vessel running from LVOT to BrA. WSS⃗GV values were averaged over three time-frames centered in T_ps. To rule out the effects of inter-subject anatomical variability when comparing different aortas, WSS⃗GV, OSI and RRT data were represented in terms of heat maps over a normalized 2D template (Figure 3A). The position of points in the relevant tract of the aortic wall was expressed in terms of: (i) normalized longitudinal position s along the centerline with s = 0 at LVOT and s = 1 at BrA; (ii) given a cross section at s = ŝ, normalized circumferential position r(ŝ) along the boundary of the cross-section (r = 0 at the most anterior point, and increases when moving counterclockwise over the boundary). The 2D template was sub-divided into 8 sectors with respect to LVOT, PA and BrA landmarks and anatomical short axis orientations (Antero: A; Posterior: P; Right: R; Left: L; Figure 3B).

No clinically relevant increase in diameter was measured at 3-year follow-up. On patients BAV02-BAV05, an increase in diameter by <1.5 mm was measured, which is below the threshold indicated by guidelines as symptom of disease progression (0.5 mm/year, Della Corte et al., 2013). On patient BAV01, an increase in diameter by 1.8 mm was measured, which was equivalent to a yearly growth rate of 0.6 mm/year, only slightly above the threshold. However, 1.8 mm is comparable to the in-plane resolution of the cine-MRI images used to perform the measurements, and indeed no macroscopic remodeling of the aorta could be observed.

For all the analyzed datasets, streamlines over the whole acquired domain within the ROI were computed and visualized using ParaView (Sandia Corporation, Kitware Inc., Albuquerque, NM, USA) at t = T_ps and t = T_dec (Figure 4).

Qualitatively, HVs were systematically characterized by a centered and nearly laminar flow in the entire ROI, as exemplified by dataset HV02 (Figure 4). At T_ps, macroscopic secondary flows were absent, but were observed at T_dec close to the intrados of the ascending aorta. All normo-BAV patients presented altered secondary flows at T_ps, owing to the BAV deflected jet (Figure 4, white asterisks). The jet was directed anteriorly, toward the extrados of the aorta, in patients characterized by RL leaflet fusion (BAV01-BAV04). The jet was directed posteriorly, toward the intrados, in case of RN fusion (BAV05). Also, aberrant flows were observed in the distal curvature of the aortic arch in BAV01 and BAV05 patients. At T_dec, the same pattern of secondary flows observed in HVs was observed also in BAV01, BAV03, and BAV04. The pattern was somewhat complementary, with disturbances displaced toward the extrados of the ascending aorta, in the only case of RN fusion (BAV05). More evident flow disturbances represented by helicoidal structures were observed all over the ascending aorta and arch in BAV02.

In BAV patients QmaxLVOT and SV values ranged from 12.6 to 26.3 l/min and from 80 to 156 ml/beat, respectively. These values were slightly higher as compared to those obtained for HVs, which ranged from 18.4 to 30.3 l/min and from 62 to 135 ml/beat, respectively. However, these differences were not statistically significant (p = 0.44 and p = 0.68, respectively).

For HVs and BAV patients, the peak value of |V| _{i_max} and D_i and the average value of |V|_{i_mean}, α_πi and d_πi over the centerline of the ascending aorta were considered for the two tracts LVOT-PA and PA-BrA (Table 3). No statistically significant differences were observed between the two groups in terms of peak D_i in the LVOT–PA tract (p = 0.31) nor in the PA–BrA tract (p = 0.25). Statistically significant differences between HVs and BAVs (p = 0.04) were found only for the peak values of |V|_{i_max} in the LVOT–PA tract, and for the average value of α_πi in the PA–BrA tract. Non-negligible though not statistically significant differences were observed for the average value of |V|_{i_mean} in both tracts (p = 0.05 in the LVOT-PA tract, p = 0.07 in the PA-BrA tract), and for the peak value of |V|_{i_max} (p = 0.05) in the PA–BrA tract.

For HVs, the heat maps of the 10th, 50^th, and 90th percentiles of the |WSS⃗GV|, OSI and RRT distribution in this group (Figures 5A, 6A, 7A) were used as term of comparison for the heat maps obtained for each BAV patient (Figures 5B, 6B, 7B).

In the HVs population, |WSS⃗GV| values were higher in the posterior and anterior sides of the wall of the distal ascending aorta, i.e., from PA up to BrA (posterior: 1.08 ÷ 1.13 Pa, anterior 0.90 ÷ 1.22 Pa), while lower values were detected in the lateral sides (right: 0.58 ÷ 0.80 Pa, left: 0.40 ÷ 0.56 Pa (Figure 5A). In the proximal ascending aorta (i.e., sectors immediately downstream of LVOT), |WSS⃗GV| was uniformly redistributed over the circumferential extent of the vessel wall, but with a high inter-subject variability (range: 0.20 ÷ 1.19 Pa), as underscored by the differences between 10th and 90th percentile heat maps. In BAV patients, a high inter-subject variability of |WSS⃗GV| distribution was evident all over the wall of the ascending aorta, with different patient-specific localizations and extent of high and low |WSS⃗GV| regions (Figure 5B).

For both HVs and BAV patients, the distribution of OSI was consistent with the one of |WSS⃗GV|. For instance, in HVs the 50th percentile heat map showed regions of high OSI at the right and left side of the distal ascending aorta (right: 0.24 ÷ 0.38, left: 0.21 ÷ 0.28; Figure 6A), which corresponded to low WSS regions. A similar correspondence between low |WSS⃗GV| regions and high OSI values could be observed for most locations of the heat maps obtained for BAV patients, thus preserving a high inter-subject variability in the distribution of OSI values (Figure 6B). The decomposition of the OSI along the circumferential and axial directions of the vessel showed that in BAV patients extended regions of the aortic wall were characterized by OSI_Ax values exceeding the 90th percentile of the distribution observed in HVs (Supplementary Figure 1). These regions were consistent with those characterized by extremely high value of the OSI computed for the full WSS vector (Figure 6B). By contrast, in BAV patients wide regions of the aortic wall were characterized by OSI_Circ values that were below the 10th percentile of the distribution observed in HVs (Supplementary Figure 2).

RRT heat maps consistently combined the features observed for OSI and |WSS⃗GV| (Figure 7).

Differently from the previously reported metrics, the TransWSS was characterized by a lower inter-subject variability over the HVs: the only exception consisted in the posterior-left region, where TransWSS ranged from 0.06 to 0.39 Pa (Figure 8). Also, all BAV patients were characterized by spread regions with abnormally high or low TransWSSs. These differences were observed particularly in the distal ascending aorta, i.e., from PA up to BrA, with values ranging from 0.03 Pa to 0.51 Pa (Figure 8).

To better highlight the patient-specific deviations of BAV patients from the HVs distribution, HV-relative heat maps were also obtained and represented on the same 2D template for each wall-related variable, with an approach resembling the one used in Van Ooij et al. (2014b). Specifically, each point on the wall of the ascending aorta was classified based on whether the corresponding value of the relevant variable was (i) abnormally high, i.e., higher than HVs 90th percentile, (ii) abnormally low, i.e., lower than HVs 10th percentile, (iii) between HVs 50th and 90th percentiles, (iv) between HVs 10th and 50th percentiles (Figure 9).

When considering BAV patients with RL fusion (BAV01-BAV04), localized regions with abnormally high |WSS⃗GV| values were detected for patient BAV02. Wide regions with abnormally high |WSS⃗GV| were detected on the distal ascending aorta of BAV03 and BAV04, and on the proximal ascending aorta of BAV01. Detected values exceeded the HV 90th percentile by up to 27.3, 27.3, 30.5, 45.8% for BAV01-BAV04, respectively. Regions of abnormally low |WSS⃗GV| were localized distally in BAV03 or proximally in BAV04, as well as in large portions of the anterior and posterior ascending aorta of BAV02. For patient BAV05, characterized by RN fusion, a pervasive abnormal increase in |WSS⃗GV| was evident: around 72% of the wall of the ascending aorta was subject to abnormally high |WSS⃗GV| values, which exceeded the HV 90th percentile by up to 131.7% (Figure 9A).

Regions characterized by abnormally high OSI values were found in the anterior sectors of the ascending aorta of patients from the RL fusion cohort (BAV02, BAV03) and for the patient with RN fusion (BAV05). In BAV04, few regions of abnormally high OSI values were found proximally, while abnormally low OSI values were computed in the distal left-extrados of the aortic wall. Of note, BAV01 was free from abnormally high OSI values, but was characterized by abnormally low OSI values over 22% of the ascending aorta wall (Figure 9B). Overall, RRT HV-relative heat maps showed a complementary pattern as compared to the obtained for |WSS⃗GV| values, thus confirming the prevalent contribution of |WSS⃗GV| over OSI in determining the value of RRT. This predominant behavior of |WSS⃗GV| was particularly evident in the anterior regions of the ascending aorta of BAV05 characterized by high values of both |WSS⃗GV| and OSI (Figure 9C). Considering the HV-relative maps of TransWSS computed at peak systole, BAV patients were characterized by ubiquitous differences with respect to the HVs population, both in regions above HV 90th percentile and below HV 10th percentile. These differences appeared more widely spread on the aortic wall for all BAV patients as compared to the patterns obtained for the other computed indexes (Figure 9D).

Furthermore, six different hot-spots were selected (Figure 9A) on HV-relative |WSS⃗GV| heat maps of different BAV patients. These hot spots were selected because they were characterized by different combinations of |WSS⃗GV|, OSI and RRT HV-relative classifications. At each selected hot spot, the time-course of magnitude and orientation of WSS⃗LP was computed. For both features, the time-course was compared to the time dependent range (10th–90th percentile) obtained in HVs at the same normalized location (Figure 10). Of note, computed |WSS⃗LP| values were approximately one order of magnitude higher than the |WSS⃗GV| computed over the aortic wall. This result is consistent with the much better capability of the LP approach, as compared to the GV one, to realistically capture |WSS⃗| values, as demonstrated in Piatti et al. (2016). The analysis of the different time-courses of |WSS⃗LP| showed that regions of abnormally high values can be subject to overloading throughout the whole systolic phase (hot spots 4 and 5) or only for a part of it (hot spots 1 and 6). Considering the time-dependent courses of WSS⃗LP angle, similar ranges were found in BAV patients with respect to the corresponding hot-spots in HVs, especially considering the first part of the systole. However, the dynamics of angle variations showed altered patterns when associated either with high (hot spots 2, 3, and 6) or low (hot spot 1) OSI values. Also, hot spots characterized by abnormally high OSI values were systematically characterized also by time-oscillations of |WSS⃗LP| only during the systolic deceleration phase (hot spots 2, 3, and 6).

In this work, we exploited the in vivo measurements yielded by 4D Flow sequences to analyze blood fluid dynamics and WSS⃗ -related indexes in young normo-functional BAV patients (age 25 ± 10 years) without aortic stenosis or aortic dilation and to compare them vs. age-matched healthy controls. To the best of our knowledge, this is the first study of this kind specifically focused on a cohort of BAV-affected patients without any anatomical derangement of the aorta at the organ length-scale. From this standpoint, the present study differs from previous studies based on 4D Flow imaging, which were run on very heterogeneous cohorts of BAV patients that included elder patients and subjects already affected by aortic dilation or aortic stenosis (Barker et al., 2012a; Hope et al., 2014). Consequently, those studies typically aimed to testing possible correlations between the altered fluid-dynamic stimuli on the aortic wall and its microstructural/anatomical derangements (Barker et al., 2012a; Bissell et al., 2013; Hope et al., 2014; Mahadevia et al., 2014) following the onset of aortopathy: differently, the present study aims at gaining insight into the possible role of such altered fluid-dynamic stimuli as trigger of the response pathways that ultimately lead to aortopathy.

A second novel aspect of the study resides in the computation of WSS⃗ and of the related metrics. Regarding the computation of WSS⃗, we obtained the 3D vector based on the calculation of the strain rate tensor through derivative filters. To the best of our knowledge, this method is not used by other research groups in the context of 4D Flow data processing. As expected from our previous benchmarking (Piatti et al., 2016), this method, and namely its implementation on 2D cross-sections, yielded values of WSS⃗ magnitude that are greater than those reported in the literature by approximately one order of magnitude. For instance, in Barker et al. (2012a) systolic values of 0.3 ± 0.1 Pa are reported for the ascending aorta of healthy subjects with a cardiac output of 5.2 ± 1.2 l/min, and values ranging from 0.4 ± 0.2 Pa to 0.8 ± 0.3 Pa, depending on the aortic side, are reported for BAV patients with non-stenotic valve, moderate aortic dilation, and cardiac output of 6.3 ± 1.7 l/min. Similarly, in Hope et al. (2011), even when considering the wall region impacted by highly eccentric systolic jets, peak WSS⃗ magnitude values of 1.67 Pa were reported in BAV patients. Regarding the WSS⃗ -related metrics, we provided for the first time such a broad spectrum of indexes, in the attempt to test whether some indexes may be more relevant in the context of BAV-related aortopathies.

In the present study, we could find remarkable abnormalities in the space distribution, in the peak magnitude values, and in the time-course of WSS⃗ acting on the ascending aorta of the considered cohort of BAV patients characterized by normo-functional valves and non-dilated aortas. These abnormalities were even more evident when considering WSS⃗ -derived indexes; in particular, every considered BAV patient was characterized by pervasive alterations of the TransWSS and by broad regions with extremely low values of OSI_Circ, indicating that alterations affect not only peak |WSS⃗|, but also the time-course and the directionality of WSS⃗. Of note, the computed WSS⃗ -related alterations vs. healthy controls were remarkable and more evident than those appearing from the qualitative and quantitative analysis of the bulk flow. Also, within the sub-group of RL fusion patients, the flow field pattern appeared repeatable, in particular at peak systole. Instead, |WSS⃗|, OSI and RRT heatmaps highlighted a non-negligible inter-subject variability. These evidences should be cautiously evaluated owing to the small number of considered BAV patients, however they suggest that the quantitative analysis of WSS⃗ -derived indexes may have a greater potential to discriminate between healthy and aberrant fluid dynamics, and possibly to grade the severity of the BAV-related fluid-dynamic alterations on a patient-specific basis as compared to the analysis of the bulk flow field.

Even though we computed far greater values of |WSS⃗| as compared to previous studies focused on the ascending aorta and based on the processing of 4D Flow data, the alterations we found in our BAV patients vs. HVs were consistent, in terms of percentage variations, with those obtained on more heterogenous cohorts of BAV patients (Hope et al., 2010; Bissell et al., 2013), as well as on BAV patients characterized by severe aortic dilation and selected for immediate aortic resection (Van Ooij et al., 2014a). Yet, measurements of ascending aorta diameter at 3-year follow-up did not reveal any clinically relevant aortic anatomical remodeling at the organ length scale, despite the anomalies in WSS⃗ detected at baseline on the endothelium of the aortic wall. This evidence was true not only when abnormally high |WSS⃗| values affected limited regions and exceeded HVs values by up to 46% (BAV01-BAV04), but also when abnormally high |WSS⃗| values affected the entire ascending aorta and were characterized by a striking 132% increase as compared to HVs values (BAV05).

Hence, in the context of BAV-related aortopathy our results, if analyzed per-se, may suggest two mutually exclusive interpretations. On the one hand, one may conclude that there is no cause-effect relationship between altered fluid-dynamic stimuli and onset of aortopathy. On the other hand, one may conclude that the anatomical remodeling of the aortic wall due to altered fluid-dynamic stimuli is a slow process that develops over decades and may not be detected in young patients.

However, the analysis of our results in light of the evidences of previous studies from the literature suggests that the first interpretation is extremely unlikely, whereas the second one might be realistic. Indeed, recent studies investigated the effects of BAV-related fluid-dynamic alterations in non-dilated and initially healthy aortas, isolating these alterations from possible concomitant pathophysiological factors (Atkins et al., 2014; Atkins, 2016; Cao et al., 2017; McNally et al., 2017). Morphotype-dependent flow alterations in the aorta proved to generate different mechanical stimuli on the aortic wall in terms of WSS⃗ magnitude and directionality (Cao et al., 2017). These BAV-specific mechanical stimuli were different on the convexity and on the concavity of the ascending aorta. Interestingly, when applied in vitro to fresh specimens of healthy porcine aortas, these two region-dependent stimuli induced different responses: only when reproducing the WSS⃗ pattern computed at the convexity significant alterations in the expression of MMP-2 and MMP-9, which are involved in medial degradation, were detected already after 48 h (Atkins et al., 2014; Atkins, 2016). Similar evidences of dysregulations in gene expression and metalloproteinase concentration were also reported at the sites of highest |WSS⃗| in BAV patients undergoing aortic resection, thus in presence of clinically-relevant aortopathy at the organ-scale length (Guzzardi et al., 2015).

Combining these previous evidences with the ones we obtained in the study herein presented we may speculate that WSS⃗ alterations precede the onset of aortopathy, that their short-term effect is detectable only at the cell level or at the microscale, and that the subsequent anatomical remodeling occurs only on the very long term. Should this speculation be correct, its practical implication would be that any follow-up study aimed at assessing the prognostic potential of fluid-dynamic alterations should be planned over a time-frame of decades.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.