Validation of Numerical Simulations of Thoracic Aorta Hemodynamics: Comparison with In Vivo Measurements and Stochastic Sensitivity Analysis

Computational fluid dynamics (CFD) and 4D-flow magnetic resonance imaging (MRI) are synergically used for the simulation and the analysis of the flow in a patient-specific geometry of a healthy thoracic aorta. CFD simulations are carried out through the open-source code SimVascular. The MRI data are used, first, to provide patient-specific boundary conditions. In particular, the experimentally acquired flow rate waveform is imposed at the inlet, while at the outlets the RCR parameters of the Windkessel model are tuned in order to match the experimentally measured fractions of flow rate exiting each domain outlet during an entire cardiac cycle. Then, the MRI data are used to validate the results of the hemodynamic simulations. As expected, with a rigid-wall model the computed flow rate waveforms at the outlets do not show the time lag respect to the inlet waveform conversely found in MRI data. We therefore evaluate the effect of wall compliance by using a linear elastic model with homogeneous and isotropic properties and changing the value of the Young’s modulus. A stochastic analysis based on the polynomial chaos approach is adopted, which allows continuous response surfaces to be obtained in the parameter space starting from a few deterministic simulations. The flow rate waveform can be accurately reproduced by the compliant simulations in the ascending aorta; on the other hand, in the aortic arch and in the descending aorta, the experimental time delay can be matched with low values of the Young’s modulus, close to the average value estimated from experiments. However, by decreasing the Young’s modulus the underestimation of the peak flow rate becomes more significant. As for the velocity maps, we found a generally good qualitative agreement of simulations with MRI data. The main difference is that the simulations overestimate the extent of reverse flow regions or predict reverse flow when it is absent in the experimental data. Finally, a significant sensitivity to wall compliance of instantaneous shear stresses during large part of the cardiac cycle period is observed; the variability of the time-averaged wall shear stresses remains however very low. In summary, a successful integration of hemodynamic simulations and of MRI data for a patient-specific simulation has been shown. The wall compliance seems to have a significant impact on the numerical predictions; a larger wall elasticity generally improves the agreement with experimental data.


INTRODUCTION
Techniques based on computational fluid dynamics (CFD) have been extensively used in the last few years to investigate hemodynamics inside arteries in both healthy and diseased subjects (see e.g. Refs. 8,9, and 13 and the references therein). CFD enables the investigation of pressure and flow field at a time and space resolution unachievable by any in vivo measurement. As a consequence, CFD permits to compute a variety of quantities and indicators that are difficult to be obtained from medical imaging, as e.g. wall shear stresses. In addition, since medical images are now capable to provide accurate morphological features and CFD techniques are mature to tackle complex 3D geometries, numerical simulations permit to investigate hemodynamics in patient-specific geometries. This is important since it has been shown in the literature that the actual geometric details may deeply affect the flow features, as e.g. the regions where the wall shear stress takes very-high or very-low values. 7 On the other hand, different sources of uncertainties are present in CFD models which can propagate and affect the output quantities of interest. Among them are boundary conditions assumed at the inflow and outflow of the computational domain and modeling of the vessel wall compliance properties. In previous works, we focused on inlet and outlet boundary conditions and we assessed their impact on the hemodynamics of ascending thoracic aortic aneurysms. [3][4][5] Once again, the integration of patient-specific medical imaging data could be useful to reduce the modeling assumptions. In particular, 4D-flow magnetic resonance imaging (MRI) 26 is a non-invasive technique that allows not only the geometry to be acquired, but also the velocity field to be measured at different points in space and at different time instants, although, as previously said, at a resolution significantly lower than that achievable in CFD simulations. MRI data can thus be used to provide patient-specific boundary conditions to numerical simulations (see e.g. Refs. 14,20, and 27-29).
In the present work, a framework integrating 4Dflow MRI data into the numerical simulation of a healthy thoracic aorta is presented. The numerical simulations were carried out with the open-source software SimVascular. First, the geometry was acquired by MRI; then, the experimental information extracted by 4D-flow MRI was used to provide/calibrate inflow and outflow boundary conditions. As for inflow conditions, a common practice is the imposition of idealized flow rate waveform. 20,25 However, studies in the literature have highlighted the impact of inlet boundary conditions (see e.g. Refs. 4,6,14, and 28). First, 4D-flow MRI data are used herein to obtain a patient-specific flow rate waveform to be imposed at the inlet. As for outflow boundaries, we adopt the widely used 3-element Windkessel model, in which the effect of downstream organs and vessels is described by means of an electric circuit scheme consisting of a proximal resistance in series with a parallel arrangement of a capacitance and of a distal resistance (see e.g. Refs. 2,25,40, and 42 and the reference therein). In particular, analogy is identified between the voltage difference and the drop in pressure, and between the current and the flow rate. This approach has the following positive features: (i) it can be naturally implemented in numerical methods as the one used in SimVascular through Neumann boundary conditions, (ii) it is simpler and less computational demanding than coupling 3D simulations with 1D models and (iii) it is more accurate than classical outflow boundary conditions, such as prescribed pressure or traction-free conditions, which only take into account the resistance to flow in the branches of the computational domain. The 4D-flow MRI data are integrated herein within the previously described approach by calibrating the Windkessel model parameters in order to match the experimental fractions of flow rate exiting each domain outlet during an entire cardiac cycle.
Clearly, 4D-flow MRI data can also be used for comparison against numerical results, to provide a cross validation (see e.g. Refs. 30 and 36). In the present work, comparisons are provided for flow rate waveforms and velocity distributions at selected sections in the ascending aorta, aortic arch and descending aorta.
This comparison was first carried out assuming the vessel walls to be rigid. From the comparison of flow rate waveforms, it appeared that the time lag present between each outlet flow rate waveform and the flow rate imposed at the inlet in the MRI data could not be captured by the rigid-wall model. This could have been anticipated since the lag is indeed an effect of wall compliance. We therefore evaluated the effect of wall compliance by using a linear elastic model with homogeneous and isotropic properties. This model is clearly simplified compared to the actual mechanical properties of the aortic wall. 21 Despite advanced imaging techniques enable both accurate morphological and functional information, however, the evaluation of the material properties of biological tissue still remain difficult to be estimated. It is well know that the material properties affect the results of a FSI simulation and that the biological tissue shows differences from healthy and pathological subjects as well as during aging. 44 However, since the mechanical and the fluid-dynamic problems must be coupled together, the use of more sophisticated mechanical models would lead to huge computational costs in patient-specific numerical simulations, which may become not affordable in the perspective of using CFD simulations in routine clinical follow-up. Thus, the question here is whether a linear elastic model may give accurate results and which is the impact of the value of the Young's modulus, E. In order to provide a systematic quantification of the effects of varying E on the output quantities of interest, namely flow rate, velocity distributions and wall shear stresses, a stochastic approach was adopted. In more details, E was considered as a random uncertain parameter and a continuous response surface in the parameter space was obtained through generalized polynomial chaos (gPC). 45 Stochastic sensitivity analysis plays a fundamental role in different fields of biomedical research, 6,10,11,12,37 in particular when patient data are used or when a comparison between experimental and simulated data is established. 1 Focusing on uncertainty quantification in blood flow simulations, different numerical approaches have been used starting from sampling techniques like Monte Carlo method 6,22,38 to more compact projection-based methods as the polynomial chaos expansions. 4,5,17,[31][32][33]35 As a recent example, Bozzi et al. 6 used the Monte Carlo method to propagate the uncertainty in measured PC-MRI velocity profiles applied as inflow BCs, in the CFD model of the aorta. A total number of 600 numerical simulations were carried out and the hypothesis of steady-state flow was made. The simulation number arises from the need to achieve a desired level of accuracy and to ensure the convergence of the probability density functions (PDFs) of the output variables; the simplified flow condition comes from the necessity to reduce the computational time. In our work the gPC approach was adopted because it converges much faster than Monte Carlo based methods, 16 allowing thus to obtain accurate response surfaces of the quantities of interest in the parameter space starting from a few deterministic simulations. In this way, computational costs can be reduced and more accurate and complex computational models can be used, as e.g. unsteady flow and compliant walls.
The remaining of the manuscript is organized as follows. Information about the acquisition of MRI data is provided in ''MRI Data Acquisition and Processing'' section, whereas the numerical methodology, the computational set-up, and the stochastic approach for sensitivity analysis are presented in ''Modeling and Numerical Methods''-''Methodology for Stochastic Sensitivity Analysis'' sections. ''Flow Rate Waveforms and Wall Motions'' section focuses on the comparison of the flow rate waveforms obtained in representative sections along the aorta and of the wall motions. Velocity maps are compared and analyzed in ''Velocity Maps'' section and in ''Wall Shear Stresses'' section the effects of E on the wall shear stress distributions are presented. Finally, ''Discussion and Concluding Remarks'' section provides some concluding remarks and highlights possible objects of future work.

MRI DATA ACQUISITION AND PROCESSING
Two different types of MRI acquisitions have been performed by means of a 3T MR-scanner on one healthy subject (28 years, male) with tricuspid aortic valve. Firstly a 2D balanced turbo field echo (BTFE) sequence has been acquired for an axial section including a projection of both ascending and descending aorta. This type of acquisition (Fig. 2a) enables an increased signal from fluid, allowing a more accurate wall motion and cross sectional area evaluation. Secondly, phase contrast sequences with threedirectional velocity encoding (4D-flow MRI) have been acquired. The 4D flow volume of acquisition was oriented along a sagittal plane encompassing the ascending aorta, the aortic arch and the thoracic aorta. An isotropic voxel resolution of 2 Â 2 Â 2 mm 3 was prescribed. The remaining parameters were: echo time ¼ 3 ms; repetition time ¼ 5:32 ms; flip angle ¼ 10 : Data were acquired with prospective ECGgating during free-breathing, using a respiratory navigator. The velocity encoding range (VENC) was properly set to 250 cm s À1 after scouting on cross-sections positioned in the ascending aorta. The MRI dataset volume was retrospectively reconstructed with phase contrast angiography (PC-MRA) technique. 15 Informed consent was obtained from patient. The study protocol conforms to the ethical guidelines of the 1975 Declaration of Helsinki as reflected in a priori approval by the institution's Human Research Committee. The patient-specific inner wall of the aorta was extracted from the obtained PC-MRA volume through a level set algorithm available in VMTKLab (Orobix Srl, Bergamo, Italy) for semi-automated segmentation. The aortic inner-wall surface was clipped proximally, at the sino-tubular junction, to define the inlet surface and distally to create the outlet surfaces for the three sovra-aortic branches and the descending aorta (Fig. 2b). The resulting volume was meshed with tetrahedral cells in VMTKLab and then a linear interpolation between voxel and cell was implemented in order to map the MRI data into the meshed domain. The centerline of the aorta was computed and equispaced planes were defined along the centerline (Fig. 2b). The processing of the 4D flow MRI data was performed with a custom Python-VMTK code. Three specific planes were used to compare experimental data with results from numerical simulations; the planes and the 3D model share the same coordinate system. This approach was used to guarantee the correspondence between the two platforms (4D flow processing and CFD simulations in SimVascular) without any additional registration step. Finally, the 4D flow data from MRI sequences have been processed to estimate the Young's Modulus (E) by using the Moens-Korteweg equation 24 : where h is the wall thickness, r the vessel radius, q the blood density and PWV the pulse wave velocity. In this work two different PWVs have been derived: one for the ascending and one for the descending portion of the aorta in order to evaluate the corresponding E asc and E disc modulus.

MODELING AND NUMERICAL METHODS
In the simulations, blood was considered as a Newtonian and incompressible fluid, as widely accepted for blood in large vessels, with the following values of density and kinematic viscosity: q ¼ 1:06 g cm À3 and m ¼ 3:77 Â 10 À2 cm 2 s À1 : The threedimensional Navier-Stokes equations for incompressible flows (omitted here for the sake of brevity) were thus considered as governing equations. No explicit turbulent eddy viscosity was introduced in the momentum equation. The viscosity is only the molecular one, together with the contribution added by the numerical method for stability reasons.
As for boundary conditions, at the inlet section of the computational domain we specified a Dirichlet boundary condition on the velocity, imposing a flow rate waveform (see e.g. Refs. 20 and 25). At the outflow boundaries, we used a 3-element Windkessel model, also called RCR model 2,25,40 and sketched in Fig. 1. In this scheme the proximal and distal resistances (R p and R d ) represent the viscous resistance the blood flow undergoes inside large and small vessels, respectively. The capacitance C quantifies the compliance, and thus the capability of storing fluid, of the major arteries (see also Ref. 5 for more details). On the arterial wall a no-slip condition was imposed between the fluid and the wall, i.e. u f ¼ u w ; being u f the velocity of the fluid at the wall and u w the velocity of the wall. Both rigid (u w ¼ 0) and deformable (u w 6 ¼ 0) simulations were performed.
The hemodynamic simulations were carried out with SimVascular, an open-source comprehensive package specific for cardiovascular problems. 39 The governing equations are discretized by means of a finite-element method, including SUPG/PSPG stabilizing terms. 43 The stabilized formulation makes it possible to choose linear shape functions for both velocity and pressure (so-called P1-P1 elements). Additional stabilizing terms are introduced to prevent also the numerical instabilities due to backflow at the outlets, which is typical of cardiovascular applications. 18 When taking into account the wall compliance, SimVascular uses a linearized kinematics approach, in which the fluid mesh is kept fixed but the nodes at the interface with the solid wall have in general non-zero velocities. The strong form of the mechanical problem for the vessel wall in the domain X w is the following: where s w is the displacement field, q and r w the vessel wall density and stress tensor, respectively, and b w is the prescribed volume force. The thin wall approximation allows b w to be expressed in terms of the wall thickness f and of the surface traction t w acting on the wall, and thus also in terms of the traction t f acting on the fluid: This permits to couple the solid and the fluid problems and to obtain a single variational formulation in which the effect of the vessel wall results in a Neumann-type boundary condition. 19 The stress tensor is expressed as r w ¼ D : e w ; where e w is the infinitesimal strain tensor, i.e. e w ¼ ru w ¼ 1=2ðru w þ ðru w Þ T Þ and the second-order tensor D is the following:  in which E is the Young's modulus, m w is the Poisson's ratio, and k w is a parameter accounting for a parabolic variation of the transverse shear stress through the membrane. 19 For advancing in time, the generalized a-method is used, which is an implicit technique that allows to achieve second order accuracy in time (see Ref. 23 for further information).

COMPUTATIONAL SET-UP
As already mentioned, boundary conditions based on 4D-flow MRI were specified. Experimental data showed some inconsistencies along the aorta, with non-balanced flow rates entering and exiting the domain during a cardiac cycle. To avoid these errors, probably related to the limited spatial resolution of the experimental technique, we decided to focus on a smaller segment where flow rate balance was found (Fig. 2c). Thus, the hemodynamic simulations were performed on the whole thoracic aorta shown in Fig. 2b, but only the experimental data in the reduced region were actually used for comparison. In particular, we extracted three representative aortic sections in the reference region above presented. These sections are shown in Fig. 5. Moreover, in order to perform comparison between patient data and numerical simulation in terms of aortic wall motions, an additional section has been considered among those analyzed during image acquisition through 2D BTFE MRI sequences (see Fig. 2b). In order to guarantee the best level of accuracy in the post processing of the experimental data, only the section in the ascending aorta has been used for the comparison because the acquisition has been set up on the parameters of this section, due to the better contrast of the image (see Fig. 2a).
In all the numerical simulations, at the inlet section we directly imposed the flow rate measured at the proximal section of the region of interest in Fig. 2c. As will be shown in the following, the imposition of the flow rate waveform measured not exactly in correspondence of the inlet section but a little downstream does not affect the overall results and, thus, it can be considered a good approximation. The flow rate waveform is shown in Fig. 3: the velocity distribution in the inlet section was assumed to be uniform. We decided not to use the MRI data to impose a patientspecific velocity distribution, as conversely done in Ref. 14, because this would have introduced additional sources of possible errors. Indeed, the direct implementation of the velocity map from 4D PC-MRI data requires a previous estimation of the error introduced by the movement of the inlet section during the cardiac cycle, which consists in both translation and rotation. In particular, the rotation is difficult to be tracked by MRI. Moreover, since the space and time resolution of MRI is significantly lower than that required for CFD simulations, interpolation is needed both in space and time and certainly, interpolation errors would be introduced. Finally, as explained above, some inconsistencies with non-balanced flow rates are present in the experimental data outside the area in which the data have been specifically used for comparison; therefore, the flow rate measured at the proximal section of the region of interest was used as inlet condition. Such an extrapolation cannot clearly be considered accurate for local values of velocity. Note also that, as previously said, the analysis is focused on a portion of the aorta significantly downstream the inlet, where the effects of the inlet velocity distribution may be considered small, especially if compared with those of the flow rate. The Reynolds number, based on the diameter of the inlet section and on the bulk inlet  velocity averaged during an entire cardiac cycle, is Re ' 1200: To specify the RCR parameters at the outlet sections, we first used the same simplified 0D procedure as described in Ref. 5. The resulting values are: R p ¼ 78:0 g cm À4 s À1 ; C ¼ 1:3 Â 10 À3 cm 4 s 2 g À1 and R d ¼ 1:36 Â 10 3 g cm À4 s: As pointed out in Ref. 20, the MRI data can be used to calibrate the flow rates at the different outlet sections. Then, the resistances of each outlet section were tuned so that the flow rate fractions during a cardiac cycle could match the experimental data. In other words, the previously computed equivalent R p and R d were distributed into the various outlets according to the experimental flow rate fraction exiting each section in a complete cardiac cycle, i.e. 74% for the descending aorta, and 16:0; 4.5 and 5:5% for the small branches, from proximal to distal respectively. Although the MRI flow rates are not exactly imposed at the outlets, this kind of approach allows experimental information to be incorporated and at the same time possible numerical/accuracy problems, as those occurring for instance when imposing the mass fractions at all the outlets, to be avoided. The effectiveness of this strategy will be analyzed in ''Flow Rate Waveforms and Wall Motions'' section.
When considering the wall compliance of the thoracic aorta, the arterial wall was assumed to be a linear elastic material, with homogeneous and isotropic properties, and a uniform wall thickness. The values of the mechanical properties are those reported in Table 1, where the Young's modulus is not present since it was selected as uncertain parameter in the stochastic analysis (see ''Methodology for Stochastic Sensitivity Analysis'' section).
The computational domain was discretized by means of a grid consisting of 4:0 Â 10 6 tetrahedral elements. The grid was refined at the wall: in detail, we set a maximum spacing of 0.04 cm at the wall and a growth rate of 1.1 until a maximum spacing of 0.1 cm was reached inside the aorta. Consequently, the grid counts 7:3 Â 10 5 nodes, 1:4 Â 10 5 of which on the surface. The number of the elements at the inlet and at the descending aorta outlet section is 5:2 Â 10 3 and 1:3 Â 10 3 ; respectively. The total number of elements in the outlet sections of the three small branches is approximately 1:6 Â 10 3 : To better appreciate the computational grid, Fig. 4 shows two sections: the inlet section (Fig. 4a), in which the above-discussed refinement is not present since it is a boundary, and a slice in the descending region (Fig. 4b).
The physical time step was set to 0.0005 s for both rigid and deformable cases, resulting in 2000 time steps every cardiac cycle (T ¼ 1 s).
In order to make sure that the flow reached periodicity after an initial transient, for each computation we evaluated the L 2 -norms of the normalized differences between two successive pressure waveforms and between two successive flow rate waveforms. These quantities were evaluated in correspondence of the descending aortic outlet section, and the simulations were run until the above L 2 -norms were smaller than The slice, which contains tetrahedral elements sectioned, is shown to appreciate the grid refinement near the arterial wall. 10 À3 : We usually needed to simulate 5-10 cardiac cycles in order to satisfy this criterion.

METHODOLOGY FOR STOCHASTIC SENSITIVITY ANALYSIS
For the stochastic analysis we used the gPC approach. In its non-intrusive form, it is based on a truncated projection of a given random process over an orthogonal known basis. 45 Let XðxÞ be a stochastic response and x a random event, the truncated gPC expansion can be written as follows: where nðxÞ is the vector consisting of the independent random variables in the event space (corresponding thus to the set of considered uncertain parameters), U k ðnÞ is the gPC polynomial of index k and a k is the corresponding Galerkin projection coefficient. Expansion (3) is truncated to a finite limit T. The truncation index T depends on the number of the considered uncertain parameters and on the maximum degree of the polynomials retained in the truncated expansion for each parameter (see e.g. Ref. 5). Using the maximum polynomial order for all one-dimensional polynomials (i.e., full tensor-product polynomial expansion), T is obtained as follows: where M is the number of the uncertain parameters and P i is the maximum polynomial order for the ith parameter. The coefficient a k can be computed as follows: where hÁ; Ái denotes the usual L 2 scalar product involving a weight function qðnÞ depending on the chosen polynomial family. The integrals in the scalar products are computed numerically, here we used Gaussian quadrature. The deterministic simulations have to be carried out for the values of the parameters corresponding to the quadrature points.
The polynomial family, U k ; must be a priori specified and its choice affects the speed of the convergence of the gPC expansion: a suitable polynomial family is able to approximate the stochastic response by means of fewer degrees of freedom. When using Gaussian quadrature, a optimal family has a weight function similar to the probability measure of the random variables. The choice of the polynomial family thus depends on the PDF shape of the uncertain parameters.
We aim to investigate here the sensitivity of the output hemodynamic quantities of interest to the Young's modulus E. For this purpose, we considered the following variation range: E 2 ½0:37; 2:28 MPa: The considered variation contains both the average value of the Young's modulus calculated from our experimental datasets for the healthy subject, which is the mean between E asc ¼ 0:375 MPa and E disc ¼ 0:625 MPa obtained from Eq. (1) and the E values of aneurysmatic subjects. We used a uniform PDF distribution of the random parameter: this seems a reasonable choice in that it is the least informative distribution with the highest variance in given intervals. Legendre polynomials thus represent the optimal polynomial family for the gPC basis, whose expansion was truncated to the third order. The four resulting quadrature points, corresponding to the values of E for which deterministic simulations must be carried out, are shown in Table 2. The lowest value corresponds to the average value estimated from experiments, while the highest values corresponds to those typical of pathological patients.

FLOW RATE WAVEFORMS AND WALL MOTIONS
This section focuses on the analysis and comparison of the flow rates at the chosen reference section of the descending aorta (see Fig. 5) and at the outlets of the side branches. Figure 6 shows the comparison between experimental and simulated flow-rate waveforms at the descending aorta reference section and at the three outlet sections of the branches. The tuning of the resistances appears to be effective both for the descending aorta reference section and for the supraaortic branches; indeed rather good agreement with MRI data is obtained. The difference between the peak value given by the simulation and that measured by MRI is equal to À8 cm 3 s À1 for the descending aorta reference section (minus indicates an underestimation in numerical simulations) while the differences are equal to + 12, + 5 and À7 cm 3 s À1 for Branches 1, 2 and 3 respectively. The main difference with the experimental data is a time lag, approximately equal to 0.04 s, between the experimental and numerical peaks of the flow rate, observable in the descending aorta section (Fig. 6a). The measurement error on the experimental waveforms can be estimated to 5% of the average value over the cardiac cycle 41 and thus negligible compared to the observed differences.
Since the time delay of flow rate waveforms along the aorta is a peculiar characteristic of a compliant arterial wall, we performed two deformable simulations with elastic modulus respectively equal to E ¼ 0:50 MPa; which is the value of the Young's modulus estimated from our experimental dataset of the healthy subject, and E ¼ 1:00 MPa: The tuned RCR parameters were imposed at the outlets for both values of the elastic modulus. Figure 7 indicates that, besides matching the measured flow rate fractions, the simulated waveforms better reproduce the time delay found in the experimental data when the Young's modulus is reduced (see in particular in the simulation with E ¼ 0:50 MPa in Fig. 7a). At the descending aorta reference section, the time lag between the flow rate peaks obtained from MRI data and the simulation carried out by using E ¼ 0:50 MPa is reduced to 0.02 s, practically half the one found for the rigid simulation. On the other hand, smaller values of the Young's modulus lead to a reduction of the peak values of the flow rate waveform for all the considered sections (see Fig. 7). This improves the agreement with MRI data for Branches 1 and 2, in which the flow rate peak was slightly overestimated in rigid simulations, while the underestimation for the descending aorta reference section and for Branch 3, already present in rigid simulations, slightly increases. In particular, the difference between the peak value of the flow rate waveform for the simulation and the MRI results is equal to À15 cm 3 s À1 for the descending aorta reference section and À2; 0 and À11 cm 3 s À1 for the three branches.    Based on the previous observations, it is interesting to carry out a systematic analysis of the effect of the wall compliance (quantified by the Young's modulus) on the output quantities of interest. As previously described, a gPC stochastic approach was used which allows continuous response surfaces to be obtained starting from a few deterministic simulations. In the present case, four hemodynamic simulations were therefore performed, in correspondence of the values of the Young's modulus reported in Table 2.
The flow rate waveforms at the three sections introduced in Fig. 5 are presented in Fig. 8 in terms of their PDF and compared with the corresponding experimental data. Figure 8a shows that, in the section at the ascending aorta, the PDF distribution contains the experimental data for almost the whole time period. This indicates that the experimental waveform can be obtained by modifying the parameter E in its variation range. In addition, this confirms that the imposition at the inlet section of a flow rate waveform measured a little downstream does not affect the results, as discussed in ''Computational Set-Up'' section. Figures 8b and 8c show that the experimental time delay at the aortic arch and at the descending aorta is underestimated, at least in the most probable part of the distribution (represented in blue), i.e. for most values of E. The time delay can be reproduced by the tail of the PDF, but at the same time the peak value reduces; this happens for the lowest considered values of E, as already noticed in Fig. 7 for E ¼ 0:50 MPa: For the rest of the cardiac cycle the agreement is good in the descending aorta and fair in the aortic arch section.
Furthermore, a comparison between aortic wall motions from MRI and simulations has been made for the section W in the ascending portion shown in Fig. 2c. The time behavior of the cross-section area and of the mean radius of the selected section during a cardiac cycle for MRI data is presented in Fig. 9a, whereas the comparison between the obtained displacements of mean radius in the cardiac cycle obtained from MRI data and simulations is shown in Fig. 9b. The qualitative behavior is well captured but all the simulations underestimate the mean wall displacement; the best agreement is found, once again, for the lowest value of the considered Young's modulus, which is also the closest one to the patient-specific value estimated from the experimental dataset.

VELOCITY MAPS
Now let us focus on the velocity maps in the three reference sections along the aorta. For each time instant shown with markers in Fig. 3 the comparison is presented among the MRI data, the rigid simulation and the deformable simulation with E ¼ 1 MPa: In addition, the stochastic standard deviation is also reported to give information about the variability of the results due to variations of E. Figure 10 refers to peak systole and indicates that the rigid and deformable results are close. Some small differences can be detected in the aortic arch and in the descending region (compare Figs. 10f and 10g, and Figs. 10j and 10k), which can be explained by noting that the more one moves downstream, the larger is the time delay obtained in case of deformable simulations. This is consistent with the stochastic standard deviation (Figs. 10d, 10h and 10l), which in the ascending aortic section is very low except for a limited region   Fig. 3) in the ascending aorta (a-c), aortic arch (e-g) and descending aorta (i-k). Data from: MRI (a, e, i), rigid simulation (b, f, j) and deformable simulation with E = 1 MPa (c, g, k). Maps of the stochastic standard deviation of the velocity component normal to the cross-section at peak systole for the deformable simulations in the ascending aorta (d), aortic arch (h) and descending aorta (l). near the wall, while in the aortic arch and in the descending aorta is slightly larger. The simulation results are in good qualitative agreement with the experiments, even if quantitative differences are present especially in the shape and extension of the highvelocity region. Since, as previously said, the results do not change significantly with E, the same agreement is found for rigid and deformable simulations.
The time instants of maximum deceleration is the object of Fig. 11. In this case, the computed velocity is underestimated in the aortic arch and in the descending aorta, which is consistent with the flow rate waveforms (and, in particular, the smaller time delay of the simulations) presented in Figs. 8b and 8c. The other main difference is that the simulations predict flow reversal, whereas in the experimental data these regions are generally characterized by low velocity, but not by reverse flow. The lack of reverse flow in experimental data may be explained by the coarser resolutions in MRI, both in space and time, which tend to smooth the velocity gradients (see e.g. the discussion in Ref. 34). Similar discrepancies between numerical results and MRI data can also be found in the literature concerning flow reversal occurring in both ascending and descending regions during the deceleration phase (see, e.g., Ref. 25). The stochastic variability with the elastic properties in this case is larger than in the previous one, and this explains also why a larger difference can be found between rigid and deformable simulations. For all the three sections, the region of largest variability is in correspondence of the region of predicted reverse flow, whose extent tends to decrease as the wall elasticity increases, improving thus the agreement with MRI data.
At early diastole (Fig. 12), the velocities are significantly lower than in the previous instants and small reverse flow regions are present also in MRI data. As previously, both rigid and deformable simulations overestimate the regions of reverse flow compared to MRI. The variability of the numerical results with E is significant especially in the zones of negative velocity, whose extent, as previously observed, is reduced as E decreases, leading to a slightly better agreement with the experiments.
Finally, the time instant of mid diastole is characterized by low velocities and by an increasing variability of deformable results by moving from the ascending to the descending aorta. This is also responsible of the differences between the rigid and the deformable simulations, which is significant in the aortic arch (compare Figs. 13f and 13g) and in the descending aorta (compare Figs. 13j and 13k). The qualitative agreement with the experimental results is quite good in the ascending aorta, whereas in the aortic arch more significant differences can be found, even if the wall compliance seems once again to improve the comparison with experiments (compare Figs. 13f and 13g with Fig. 13e).

WALL SHEAR STRESSES
Let us now consider the TAWSS distributions, again in terms of results of the rigid simulation (Figs. 14a and 14b), results of the deformable simulation with E ¼ 1 MPa (Figs. 14c and 14d) and stochastic standard deviation (Figs. 14e and 14f). The differences between the deterministic rigid and deformable simulations is practically negligible; consistently, the stochastic variability of TAWSS is very low, as can be appreciated in Figs. 14e and 14f.
The situation is instead slightly different when considering the instantaneous WSS distributions, shown in Fig. 15 for the first three relevant time instants presented in Fig. 3; only the results of the deterministic simulation with E ¼ 1 MPa are shown together with the stochastic standard deviation.
At peak systole, the deterministic values of WSS are characterized, in both ascending and descending regions, by large spatial variations from internal and external parts of the vessel: in the ascending aorta, the internal region is subjected to larger wall shear stresses, whereas the opposite occurs in the descending aorta. The stochastic variability of WSS with E takes lower values in the ascending region and in the aortic arch, except for the posterior-lateral part of the ascending aorta, in which variability is the largest of all the domain (Fig. 15d). This behavior can be related to the large velocity variability with E observed in that region in Fig. 10d. On the other hand, the descending aorta is characterized by a stochastic standard deviation larger than average, also consistently with previous observations on the velocity maps.
At maximum deceleration time, the deterministic values of WSS are smaller than those at the systolic peak, as expected, but the variability with E in the descending aorta seems to reach even larger values (Fig. 15h). This is in accordance to the high velocity variability detected near the internal wall (blue region in Fig. 11l).
At early diastole, the deterministic values are further reduced. The maximum values of the deterministic simulation and of the standard deviation are in the internal part of the descending aorta, which is consistent with the large nucleus of reverse flow, and the large variability with E of velocity in that region, found in Fig. 12l. Significant values of standard deviation are also found in the posterior part of the aortic arch (Fig. 15l).

DISCUSSION AND CONCLUDING REMARKS
The present study is an example of integration between 4D MRI flow data, hemodynamic numerical simulations and stochastic sensitivity analysis. Integration of in vivo measurements with CFD tools is FIGURE 13. Maps of the velocity component normal to the cross-section at mid diastole (point D in Fig. 3) in the ascending aorta (a-c), aortic arch (e-g) and descending aorta (i-k). Data from: MRI (a, e, i), rigid simulation (b, f, j) and deformable simulation with E = 1 MPa (c,g,k). Maps of the stochastic standard deviation of the velocity component normal to the cross-section at peak systole for the deformable simulations in the ascending aorta (d), aortic arch (h) and descending aorta (l).
clearly important to obtain patient-specific information at a level of detail which is hardly reachable in experiments only. On the other hand in vivo measurements may provide information to be used in the computational set-up in order to be able to obtain accurate patient-specific results. Different examples can be found in the literature pointing out the importance of using patient-specific data for boundary conditions 4,6,14,28 ; however these studies usually focus on one particular boundary condition. Clearly, experimental data can also be used to validate simulation results (see e.g. Refs. 30 and 36). In our study we tried to reach a higher level of integration by exploiting the information available from patient-specific experiments for different aspects of the computational set up as well as for validation. In particular, the 4D MRI flow data were used to: (i) provide inlet conditions, (ii) calibrate outflow conditions within the RCR Windkessel model, (iii) estimate the Young's modulus of the aorta and (iv) validate the simulations results. Another question at issue was whether a model of wall compliance assuming a linear elastic behavior together with homogeneous and isotropic properties, which is a rather rough simplification of the real wall behavior, is able to give reliable results for suitable values of the Young's modulus, and, more in general, to investigate the effects of E on the simulation results. To this aim a stochastic approach was used, based on gPC, which allows a systematic exploration of the parameter space starting from a few deterministic simulations, limiting in this way the computational costs.
The results of the stochastic analysis indicate that the experimental behavior of the flow rate waveform can be accurately reproduced in the ascending aorta; on the other hand, in the aortic arch and in the descending aorta, the experimental time delay can be matched with low values of the Young elastic modulus, close to the average value estimated from MRI data. However, by decreasing E, the peak flow rate tends to become underestimated. As for the velocity maps, we found a generally good qualitative agreement of simulations with MRI data. The main difference is that the simulations overestimate the extent of reverse flow regions or predict reverse flow when it is absent in the experimental data, as e.g. at the maximum deceleration phase of the cardiac cycle. On the other hand, the stochastic sensitivity analysis showed that the region of reversal flows is significantly affected by the elastic properties of the vessel walls; moreover, the extent of these zones decreases with E, leading to a better agreement with the experimental data for the lowest values of the Young's modulus. Note also that, the lower extent or the absence of regions of reverse flow may be, at least partly, due to the low resolution of MRI which tends to smooth the velocity fields. A comparison of wall displacements during the cardiac cycle has also been provided. As previously described, PC-MRI suffers from a low spatial resolution, consequently this technique is not able to accurately measure the variation during the cardiac cycle such as the radius and area of the aorta cross sections. As far as the wall motion is concerned, therefore comparison has been performed on one specific section (W) where a BTFE sequence was performed, which is able to pro- duce images with increased signal from fluid allowing an increase of accuracy in the area and motion evaluation. The qualitative behavior of the cross-section mean radius displacement during the cardiac cycle is well captured but all the simulations underestimate it; the best agreement is found for the lowest value of the considered Young's modulus, which is, as previously said, the closest one to the patient-specific value estimated from the experimental dataset. Finally, the sensitivity of the wall shear stress to wall compliance has been quantified. The wall shear stress maps indicate a significant sensitivity to wall compliance during large FIGURE 15. Instantaneous WSS distributions: deterministic results for the deformable simulation with E ¼ 1 MPa and tuned RCR parameters (a, b, e, f, i, j) and stochastic standard deviation (c, d, g, h, k, l). Time instants: peak systole (a-d), max deceleration (e-h) and early diastole (i-l).
part of the cardiac cycle period but, interestingly, a very low variability of the time-averaged wall shear stresses.
In summary, a successful integration of hemodynamic simulations and of MRI data for a patientspecific simulation has been shown. The wall compliance seems to have a significant impact on the accuracy of numerical predictions; a larger wall elasticity generally improves the agreement with experimental data. However, the simple model adopted here of a constantthickness elastic vessel wall does not seem to be able to quantitatively reproduce the MRI data.
In order to obtain a better quantitative agreement, different aspects of both numerical model and experimental set up could be improved. As previously described, the first limitation of the work is the assumption of a uniform material behavior along the aorta. The E values estimated in in vivo experiments have shown significant variations along the aorta with a more compliant behavior in the ascending portion. These results are in agreement with previous studies and stress out the importance to implement in the numerical simulations a model with varying material properties. The implementation and application of more sophisticated and accurate models of the vessel wall mechanical properties could thus be the object of future work. Another possible source of inaccuracy in numerical simulations is the assumption of a plug flow at the inlet. The use of a velocity distribution obtained from MRI would certainly be interesting. However, as previously discussed, it is not trivial to obtain accurate patient-specific velocity distributions from MRI due to the difficulty in characterizing in experiments rotation movements of the inlet section. Moreover, possible interpolation errors both in space and time may be introduced when passing from the MRI resolution to the one needed for numerical simulations. Further work is thus definitely needed to reach this higher level of integration between in vivo experiments and numerical simulations. Finally, it would be interesting to investigate the possible effects of the introduction of an explicit turbulence model in the simulations.
As for the experiments, besides the limited resolution of MRI, the second main limitation of the methodology described herein, comes from the use of a single VENC value during MRI acquisition. When the VENC value is set too high, the velocity-to-noise ratio, itself inversely proportional to the VENC setting, is affected by errors. The increased noise will lead to erroneous velocity readings, particularly in regions of slow flow. Finally, the error on MRI data could be better quantified by investigating, for instance, the data repeatability.