Validation of 4D Flow based relative pressure maps in aortic flows

While the clinical gold standard for pressure difference measurements is invasive catheterization


Introduction
Measures of the blood pressure distribution are important diagnostic indicators for cardiovascular diseases and disturbed blood flow in large vessels. A prominent example, the common congenital heart disease coarctation of the aorta (CoA; a narrowing of the proximal descending aorta) is clinically evaluated by measurements of the pressure difference across the narrowing ( Warnes et al., 2008;Feltes et al., 2011 ). Catheterization is the current gold standard to assess pressure differences, but it is invasive, expen-ages of the heart and the great vessels for the assessment of the anatomy, function and flow. Three-dimensional and timedependent flow measurements, referred to as 3D cine PC-MRI or 4D Flow, consist in the acquisition of an anatomical image and velocity-encoded images in three orthogonal directions ( Markl et al., 2003;Uribe et al., 2009;Dyverfeldt et al., 2015 ). Its capacity of measuring complex flow patterns enables 4D Flow the quantification of different hemodynamic parameters Sotelo et al., 2021 ). 4D Flow allows to infer 3D and timedependent pressure maps using the NSE. The classical pressure recovery method is the Pressure Poisson Estimator (PPE), obtained by taking the divergence of the NSE and inserting the velocity measurements in the right-hand side ( Ebbers et al., 2002;Krittian et al., 2012 ). More recently, several additional methods have been introduced, see Bertoglio et al. (2018) for a comprehensive review. The Stokes Estimator (STE) ( Švihlová et al., 2016 ) computes 3D pressure maps using a Stokes equation based on the physical pressure and an auxiliary, non-physical velocity field. In addition to PPE and STE, less computationally expensive methods have been proposed, like the Work Energy-derived Relative Pressure (WERP) method ( Donati et al., 2015 ) based on an integral energy balance of the Navier-Stokes equation, the integral momentum relative pressure estimator (IMRP)  based on integral linear momentum conservation, or the virtual WERP ( v WERP) method ( Marlevi et al., 2019;, based on a different treatment of the convective term than in the IMRP. Studies comparing these methods are scarce, particularly using real data. Bertoglio et al. (2018) found using numerical data that the WERP method yielded more accurate results than PPE, but was less accurate than the STE and IMRP methods, which were shown to have similar accuracy. The v WERP method was reported to be more versatile, accurate and robust than the PPE and the WERP methods ( Marlevi et al., 2019 ), but has not been compared to STE or IMRP. A particular limitation of WERP is the assumption that the studied vessel segment does not include any bifurcations. As a consequence, the method cannot be used to estimate the pressure difference between the ascending aorta (AAo) and the descending aorta (DAo)-a common location of CoA-due to the presence of the supra-aortic branches. Moreover, instead of 3D relative pressure maps, WERP, IMRP and v WERP can only deliver mean pressure differences between two planes intersecting the vessel of interest. They are therefore difficult to apply in cases when pressure spatial variations are present as in post stenotic areas.
The PPE method has been analyzed in several validation studies ( Rengier et al., 2015;Riesenkampff et al., 2014;Goubergrits et al., 2019 ). An in vitro validation study showed a good correlation between PPE and catheter pressure differences ( r = 0 . 89 , p < . 001 ) in the simple setting of an elastic straight tube phantom ( Rengier et al., 2015 ). The PPE method was further assessed in 13 patients with moderate CoA in ( Riesenkampff et al., 2014 ), where instantaneous peak pressure differences from 4D Flow were found to be slightly underestimated on average in comparison to the catheterization data, with a bias of 1.5 mmHg and a variability of ±4.6 mmHg (two standard deviations). In Goubergrits et al. (2019) , PPE pressure differences showed a good agreement with catheter measurements in CoA patients in cases with sufficient spatial image resolution (at least 5 voxels/diameter). However, systematic underestimation of the pressure difference was found for lower resolutions ( 3 . 4 ± 0 . 64 voxels/diameter).
To the best of the authors' knowledge, no validation studies have been reported for the STE method using in vitro or in vivo 4D Flow despite its promising features. The aim of this work is to fill this gap by analyzing the PPE and STE methods for in vitro and in vivo 4D Flow data of cases of CoA. Invasive catheterization measurements constitute the reference against which the 4D Flow pressure estimates are validated. The in vitro study uses data from a realistic CoA phantom, designed and studied previously by the authors , and analyzes the effects of the MR image resolution and segmentation, noise, cardiac output and the severity of the CoA on the accuracy of the PPE and STE methods. Both methods are also compared for 4D Flow and catheterization data of two patients.
Given the limitations mentioned above of the WERP, v WERP and IMRP methods, they are not included in this study.

Problem statement
elative pressure can be computed directly from the velocity measurements by evaluating the linear momentum conservation equation of the incompressible Navier-Stokes model, i.e., where ρ is the density of the fluid and μ its dynamic viscosity,  ( Ebbers et al., 2002;Ebbers and Farnebäck, 2009;Krittian et al., 2012 ) and the Stokes Estimator method (STE) ( Švihlová et al., 2016;Cayco and Nicolaides, 1986 ), which will be described in this section.
Discretizing Eq. (1) in time, here with the first order backward difference formula, gives the following expression for the pressure gradient: The indices 1 ≤ k ≤ N denote the time snapshot of the measurements and t the temporal offset between two consecutive measurements or cardiac phases, with time stamps t k = k t. For the first step, k = 1 , a forward difference has to be used instead since no previous measurements are available. Evaluating the right hand side of Eq. (2) for spatially undersampled and noisy velocity measurements u m , yields a pressure estimate from its approximate gradient ∇ ˆ Higher order time schemes, while more accurate in theory for small time steps, are not beneficial in the present context due to the coarse time sampling of the measured velocities. Note that in previous works, for instance in Bertoglio et al. (2018) ; Marlevi et al. (2019) a second-order mid-point scheme was used. However, this leads to stronger underestimations of the pressure differences. This can be explained from the nature of time under-sampling in MRI, namely that u k m is reconstructed by assuming the flow velocity as constant within the interval [ t k − t / 2 , t k + t / 2] rather than being an instantaneous measurement at t k .
It is important to remark that in all methods derived from the Navier-Stokes equations, e.g., Bernoulli-based, PPE, STE, and in CFD simulations, at any instant of time, the pressure is uniquely defined up to a constant (with respect to the spatial coordinates). Therefore, only instantaneous pressure differences between different locations can be compared at different times. Catheterization or sphygmomanometer pressure measurements are taken relative to the atmospheric pressure. Hence, the pressures are calibrated with respect to a global reference and pressure values can be compared at different times and among patients. A common measure in the clinical practice are the so-called peak-to-peak pressure differences, which compares the largest pressure difference registered between two locations at any time during the cardiac cycle, thus taking into account time shifts due to the vessel elasticity. Peak-to-peak pressure differences can only be determined by means of catheterization or with the models described above when calibrated with catheterization data, which however violates the non-invasiveness of the estimation methods. For this reason, the present work focuses on instantaneous pressure differences instead of peak-to-peak values.

Pressure poisson estimator (PPE)
Assuming sufficient regularity (i.e., assuming that all required derivatives exist), a Poisson equation for the pressure estimation can be obtained by taking the divergence of (4) , Eq. (5) with BCs (6) can be discretized in space with the finite element method (FEM). In order to ensure that the resulting algebraic problem admits a unique solution, the pressure can be fixed arbitrarily at one point of the mesh via a Dirichlet boundary condition, without changing the pressure gradient.

Stokes estimator (STE)
The Stokes Estimator introduces a divergence-free auxiliary function w with w = 0 on ∂ . The Laplacian of w is subtracted from Eq. (4) as a regularization term (with unitary viscosity here for simplicity) and we obtain The auxiliary function w holds no physical interest, and it is expected to be negligible compared to the pressure term as long as the right-hand-side R k is the gradient of a scalar (irrotational). The advantages of the STE with respect to the PPE method are (1) that no artificial BCs for the pressure are necessary and (2) that it has lower regularity requirements, since no additional derivatives are applied on the measurements R k . In fact, in contrast to the PPE method, the STE method searches the pressure in the natural energy space of the pressure in the original Navier-Stokes equations ( Temam, 2001 ). As for the PPE method, the pressure constant has to be fixed for ensuring solvability of the algebraic problem.

Phantom setup
The experimental setup of the aortic phantom study is described in detail in Urbina et al. (2016) and Montalba et al. (2018) . The phantom represents the thoracic aortic circulation with a closed circuit, consisting in a MR-compatible pulsatile unit pump with a control unit (CardioFlow 50 0 0 MR, Shelley Medical Imaging Technologies, London, Canada) and a realistic aortic model built with flexible silicone (T-S-N 005, Elastrat, Geneva, Switzerland). The control unit was configured to simulate two physiological aortic flow conditions-rest conditions at 75bpm and stress conditions at 136bpm-, calibrated with average data of ten healthy volunteers. The resulting peak flow in the AAo was 4.5mL/min under rest conditions and 5.9mL/min under stress conditions. These flow conditions correspond to peak Reynolds numbers Re = ρUD /μ at the inlet of Re ∼ 3500 at rest and Re ∼ 50 0 0 under stress conditions, at the time corresponding to peak systole. Re is based on the diameter of the inlet tube D (as recovered from the 4D Flow images) and the average velocity U over a cross-section upstream of the aortic phantom.
Different degrees of CoA were placed in the DAo just after the left subclavian artery (at the isthmus level). CoA's were built with Technyl with an effective orifice of 13; 11; 9mm and a length of 10mm, leading to degrees of stenosis of 40;50;60% with respect to the native DAo distal to the CoA. The liquid used in the system consisted of a homemade volume-mixing blood mimicking fluid with 60% distilled water and 40% glycerol (Orica Chemicals, Watkins, CO), with a density of 1.119 g/cm 3 , viscosity of 4.83e-3Pa.s ( Bock et al., 2011;Holton et al., 2005 ), and a T1 value of 900ms, which are representative values for human blood. The density and viscosity values of the mixture were confirmed using an empirical formula reported by Cheng (2008) with an ambient temperature of 22 • .

Catheterization
The phantom was equipped with a catheterization unit to measure invasively and simultaneously the pressure gradient across the CoA. For this purpose, two catheters (5 French, Soft-Vu, AngioDynamics, Latham, NY) of side-hole type with transducers (AngioDynamics) were placed in the AAo and 2cm after the CoA and were connected to a patient monitor (Contec Medical Systems, Hebei, China). The pressure catheters were zeroed at the same height of the phantom.
Pressure information from the two catheters was recorded simultaneously during 1 minute in the CoA phantoms (40;50;60% degrees of stenosis) at rest and at stress conditions, using the commercial software Central Monitor System V3.0 (Contec Medical System). The pressure difference is obtained by subtracting the averages of both signals over the cardiac cycles. The average cycle (the so-called phase average) was obtained by first upsampling and filtering the signals, then determining the instantaneous phase angles by means of applying the Hilbert transform ( Luo et al., 2009 ) to the band-pass filtered signal around the cardiac rate. The phase angle interval ] − π , π ] was split into 52 (rest) or 28 (stress) segments, and the original signal was associated to these segments according to the instantaneous phase angles. The phase averages and phase standard deviations were found by averaging (or computing the standard deviation) within each phase segment. The phaseaveraged pressure signals and the corresponding phase variabilities are illustrated in Figs. A.10 . At the peak time, the value of two phase standard deviations, 2 σ (95.45% confidence interval), ranges between 10 mmHg (40% at rest) and 27 mmHg (40% at stress). For the pressure difference, 2 σ lies in the interval [0 . 8 , 11] mmHg.

4D Flow data acquisition
Phantom data were acquired in a 1.5T MRI system (Achieva, Philips, The Netherlands) using a 4-channel body coil and retrospective cardiac gating. The control unit of the pulsatile pump generated a trigger signal to synchronize the MR data acquisition. In order to provide static tissue for phase correction algorithms used in PC-MRI, 6L of 1% agar were placed around the aortic phantom at least 6 hours before scanning.
4D Flow images were acquired with an isotropic voxel size of 0.9mm for all phantoms (CoA with 40;50;60% degrees of stenosis) under rest and stress conditions. The acquisition parameters are summarized in Table 1 (see column "in-vitro experiments"). A detailed discussion of the dominant flow features in the 4D Flow data is included in A.2 .

Patient study: 4D-Flow MRI and catheterization
In addition to the phantom study, we include two subjects with CoA ( Subject 1 : 12 years, weight 47 kg, height 151cm; Subject 2 : 35 years, weight 63 kg, height 205cm). Subject 1 presented a native CoA and mild aortic valve stenosis, mild left ventricular hypertrophy and systemic hypertension at rest. Subject 2 presented a repaired CoA using a subclavian flap. In addition, there was a very mild narrowing at the level of the transverse arch, close to the isthmus and a mild dilatation in the proximal DAo. Angiographic contrast agent enhanced images of the patients' aortas with anatomical information and 4D Flow streamlines are displayed in Fig. A.13 . The cardiac output obtained from the 4D Flow data for Subjects 1 and 2 were 3.27L/min and 6.49L/min, respectively.
The clinical patient data were acquired in a combined MRI/Catheter interventional suite (XMR, see Sotelo et al. (2015) for a more detailed description), equipped with a 1.5T Achieva MR scanner and a BT Pulsera cardiac radiography unit (Philips, Best, Netherlands) ( Razavi et al., 2005;Moore, 2005 ). Patients had general anesthesia according to institutional protocol. Two femoral artery vascular accesses by percutaneous puncture were performed. A heparin bolus of 50 IU/Kg was given with activated clotting time monitoring once vascular access was obtained. A MRI compatible multi-purpose catheter was advanced under fluoroscopic guidance from the right femoral artery to the AAo just above the aortic sinus for continuous hemodynamic pressure monitoring. A second multi-purpose catheter in the left femoral artery was advanced to the abdominal aorta at the level of the diaphragm. The catheter positions are visually determined during the fluoroscopy procedure by the interventionist. We used these positions in the 4D Flow based pressure computation as shown in Fig. A.13 . The floating table was then moved to transfer the patient to the MRI scan to acquire the 4D Flow MRI data. Following the acquisition of the MR data the patient was transferred back to the catheter table. Catheter pullbacks were performed using a biplane system (Siemens Axiom-Artis d-TA, Siemens, Germany) for evaluating the pressure distribution along the aorta. X-rays images were acquired with a frame rate of 15 images per second. Phaseaveraging of the catheter pressure signals was achieved by identifying and splitting the cardiac cycles according to the pressure peak locations and subsequent averaging over the resampled cycles.
The data was acquired at St. Thomas' Hospital, London, UK. The local research ethics committee approved this retrospective study and informed consent was obtained from all patients. Table 1 summarizes the acquisition parameters for the in-vivo study (see column "in-vivo experiments").

Segmentation and mesh generation
The 4D Flow data sets were processed using an in-house MAT-LAB library (The MathWorks, Natick, MA), similarly to previous studies . The library contains a homemade segmentation toolbox, which consists in (a) a contrast adjustment of the images to increase the intensities in the vessel of interest, prior to (b) performing a 3D threshold and 3D labeling of the different regions of interest in the images. If the user agrees with the selected thresholding level, the next step (c) is the iterative manual disconnection of different objects inside the image. 3D labeling at this step is used to verify that the disconnection was performed properly.
The segmentation procedure of the patient's data used the angiographic image, i.e., the time average of the anatomic images multiplied by the magnitude of the velocities images, as described in Bock et al. (2010) . For the phantom data, the good contrast between the lumen of the vessel and the agar-agar used in the phantom reservoir enabled the segmentation based on the time average of the anatomic images.
The toolchain proceeds as illustrated in Fig. 1 . Structured tetrahedral meshes ( Fig. 1 (b)) were created from the segmented images, such that the mesh vertices matched the centers of the image voxels. The velocity vectors attributed to each voxel of the 4D Table 2 Mean pressure difference over the cardiac cycle for finite elements with varying order. Exemplary data of the 50% CoA at rest, V +0 segmentation. In parentheses: percent difference of the mean pressure difference with respect to lower order results (previous column). Flow images were transferred to the corresponding mesh vertices ( Fig. 1 (c)). The last step in the toolchain consisted in the pressure map reconstruction ( Fig. 1 (d)), described in the next section.

Pressure maps computation
Pressure maps were computed from all 4D Flow data sets with the PPE and STE methods. The pressure differences, to be compared with the corresponding catheter values, were defined as differences of the pressure averages over two spheres with a radius of 4mm at locations proximally and distally to the CoA. Averaging over several voxels rendered the pressure estimate more robust to local perturbations, e.g., induced by noise in the 4D Flow data.
The finite element method (FEM) was used to discretize the partial differential equations of the PPE and STE methods. For the PPE, continuous elements P k of first, second and third order, k = 1 , 2 , 3 , were considered. For the STE, we consider linear, continuous polynomials for both the auxiliary function and the pressure, combined with standard pressure stabilization (PSPG) in order to ensure the solvability of the problem ( Elman et al., 2014 ), and denote it P 1 / P 1 stab . In addition, inf-sup stable so-called Taylor-Hood elements are considered, consisting in continuous polynomials of order k for the auxiliary function and order k − 1 for the pressure. We use k = 2 and k = 3 , and denote the discretizations P 2 / P 1 and P 3 / P 2 , respectively. Direct solvers are ideal for solving the linear systems, taking advantage of re-using the LU factorization of the system matrix at the first time-step for all subsequent steps. Memory limits required iterative solvers to be used for higher-order simulations, which have a significantly smaller memory footprint, but lead to increased computation times. In practice, on a standard desktop workstation (8 CPU cores, 32 GB RAM memory), we could employ direct solvers for system with < 1 e 6 unknowns and used iterative solvers for larger problems (see Table 3 in Section 4.1.1 ). In the latter case, for the PPE system we used the conjugate gradient (CG) method preconditioned with algebraic multigrid (AMG). The Stokes saddlepoint system of the STE method was solved with a flexible GMRES method ( Saad, 2003 ), preconditioned with a Schur decomposition approach using AMG for the velocity block and performing 10 inner GMRES iterations on the Schur complement, preconditioned with the diagonal of a pressure mass matrix ( Elman et al., 2014;Benzi et al., 2005 ). These methods were implemented using the FEM library FEniCS ( Alnaes et al., 2015 ).
In a preliminary step, the effect of the element order is analyzed and optimal discretizations are selected for the subsequent analyses. We focus on the (arbitrary) example of the 50% CoA at rest (standard segmentation).

Test cases
The 4D Flow phantom dataset was augmented synthetically in order to study the effects of image resolution, segmentation and noise, in addition to the influence of severity and cardiac load, on the pressure difference estimation.

Phantom study
Resolution study. Lower resolution data with isotropic voxel sizes h = 1 . 4 mm and h = 2 . 0 mm were generated from the original highresolution image (0.9mm isotropic voxel) using linear interpolation.
Segmentation study. Additional segmentations were created for the phantom velocity images by adding or subtracting one voxel at the boundary of a reference segmentation, thus extending or decreasing the lumen cross-section. The reference segmentation will be referred to as V +0 , the segmentations with one voxel added or subtracted as V +1 and V −1 .
Noise study. The phantom 4D Flow data had a very high velocityto-noise ratio (VNR) in comparison to in vivo data. In order to study the effect of realistic noise intensities on the numerical methods, synthetic noise was added to the phantom velocity measurements (voxel size 2.0mm, matching typical in vivo data resolutions, and the V −1 segmentation) and enters the computations in nonlinear fashion, due to the nonlinearity of Eq. (3) . The noise was defined as additive white Gaussian noise with zero mean, independently sampled for every time step, with a standard deviation σ proportional to the VENC value ( Dyverfeldt et al., 2015 ) with V ENC/σ = 10 (we use this measure in analogy to the VNR). In comparison, the noise intensity estimated from the patient 4D Flow data, by computing the standard deviation of the velocity components perpendicular to the flow direction within the DAo, during diastole at negligible mean flow, was approximately V ENC/σ ≈ 20 . Hence, the synthetic noise applied to the phantom data is "twice as bad". The original 4D Flow phantom data at voxel size 0.9mm had on average V ENC/σ ≈ 50 , which is reduced by the subsampling process to ≈ 150 on the synthetic datasets with h = 2 . 0 mm .
Statistics of the pressure drop were analyzed using 30 independent noise realizations.

Discretization order
First, the convergence of the pressure difference under prefinement (increasing the polynomial order of the finite element space) was studied for the example of the 50% CoA at rest (standard segmentation) to establish the optimal order for the subsequent studies. Table 2 shows the temporal mean of the pressure difference over the cardiac cycle, obtained with the PPE method at the three voxel sizes h = 0 . 9 , 1 . 4 , 2 . 0 mm and using first, second and third-order finite element discretizations.
Corresponding results are shown for the STE method, using P 1 / P 1 stab , P 2 / P 1 , P 3 / P 2 elements, for the coarsest grid only. The effect of higher-order approaches on the STE mean pressure difference is negligible at h = 2 mm . In contrast, the PPE pressure differences change significantly for all voxel sizes between P 1 and P 2 , whereas P 3 does not lead to a significant change with respect to P 2 . The number of unknowns (degrees of freedom) are summarized in Table 3 , for the examples reported above. Iterative solvers, here required for the cases with 1 e 6 degrees of freedom, are more costly and more complicated for STE than for the PPE system. Table 4 lists the computation times required for solving the PPE and STE problems for the complete cycle of 25 measurements. Most configurations were solved within seconds up to several minutes, depending on the problem size. The largest PPE scenario terminated after approximately 30min, with CG/AMG converging after 26 iterations. In comparison, the STE case requiring iterative solution, P 3 / P 2 elements at 2mm voxel size, took approximately 5h on a single CPU core ( ∼ 270 iterations per time step).
The results presented in the following sections were computed using P 1 / P 1 stab elements for STE and P 2 elements for the PPE method.

Pressure maps
Maps of the relative pressure, computed by means of the STE and the PPE methods, are displayed in Figs. 2-4 , for the times of the maximum pressure differences (see below). The figures show slices through the approximate centers of the AAo, DAo, the aortic arch and the CoA. Since not all parts are aligned in a 2D plane, two slices are joined, hence the kinks visible in the 3D figures. For the purpose of comparing STE with PPE pressure maps, the arbitrary pressure constant was set such that the pressure is zero in the center of the outflow section. In addition, the figures include the spheres over which the pressure is averaged for the evaluation of the pressure differences.
All cases exhibit similar qualitative features. Along the length the inlet hose, the pressure slightly increases and reaches a maximum downstream of the diameter expansion which remains approximately constant through the aortic arch, until the constriction. At the junction of the inlet hose with the AAo, local pres- sure minima occur related to flow recirculation. Significant pressure differences are present in all scenarios across the CoA. Within the DAo the pressure shows little variation. However, local pressure minima and maxima are visible at the outer phantom wall downstream of the coarctation, related to the impingement of the jet (see Fig. A.12 ). In the 60% CoA, pressure fluctuations surround the jet. The most striking feature in all figures (albeit to a lesser extend in Fig. 2 (a)) is that the PPE method consistently computes significantly reduced CoA-related pressure gradients than the STE method, resulting in a lower AAo pressure (because the pressure was 'zeroed' in the DAo).

Catheter vs. 4D Flow pressure differences -effects of segmentation and resolution
A quantitative comparison of the PPE and STE results is undertaken with respect to ground truth catheterization data in terms of pressure differences between locations in the AAo and DAo (spheres in Figs. 2-4 ). In particular, we consider the effects of different segmentations and image resolution as described in Section 3.5.1 . Fig. 5 compares the maximum instantaneous pressure differences, denoted δ p, obtained with the PPE and the STE methods with catheterization data for all investigated phantoms, segmentations and resolutions. Focusing the analysis first on the results obtained on the original resolution of 0.9mm, STE δ p estimates show a very good agreement with the catheterization data ( Fig. 5 (a)). For indices 1-4, i.e., the 40;50% CoA cases, the influence of the segmentation on δ p is weak. A larger variability due to segmentation is noted for indices 5 and 6 (60% CoA at rest/stress) where δ p ≥ 30 mmHg . Under rest conditions (experiment 5), the V +0 result matches perfectly with the mean catheter value, however, in this case, the phase standard deviation of the catheter is very large, as indicated by the error bars. At stress (experiment 6), the STE method clearly underestimates the catheter result. Using the reduced V −1 segmentation consistently yields larger δ p estimates and therefore reduces the error in this case with respect to catheterization. In experiment 1 (40% CoA at rest), the mean catheter δ p is underestimated by 50%.
However, the catheter is associated with a large variability, and the δ p estimate lies at the boundary of the ±2 σ interval.
In comparison, the PPE method at voxel size 0.9mm ( Fig. 5 (b)) leads to strongly reduced δ p in the 60% CoA and features a larger variation due to the segmentation, in particular in the 50% CoA. On V +0 , the catheter mean value is underestimated for experiments 4-6 and V −1 gives a closer match. For the 40% CoA, the differences between STE and PPE are negligible.
Under lower image resolutions 1.4;2mm, the δ p estimates of both methods for the standard segmentation are consistently lower than with h = 0 . 9 mm ( Figs. 5 (c)-(f)). Most affected are the cases 5 and 6 (60% CoA). As the voxel size increases, the variability of δ p due to the segmentation changes increases, in particular for the PPE method. V −1 significantly improves the PPE results in experiments 4-6 ( δ p ≥ 20 mmHg ), and the STE results in experiments 5 and 6 ( δ p ≥ 30 mmHg ). In the other cases, the segmentation effect is weak and no clear improvement can be noted. The dilated segmentation V +1 leads to strong underestimation, in particular for the PPE method and at low resolution.

Sensitivity to noise
Building on the results of the previous section, the discussion of the noise study is limited to the V −1 segmentation, and to voxel sizes of 2mm as found in clinical practice. Figs. 6-8 show the time profiles of the pressure differences obtained with the STE and PPE methods in comparison with averaged catheter measurements for all phantoms at rest and stress, with no additional synthetic noise (original data) and with additional noise as described in Section 3.5.1 .
The gray shaded areas indicate the band of two standard deviations σ of catheterization data over all measured cardiac cycles. For the additional noise case, ±2 σ bands are also shown for the PPE and STE results, colored accordingly, and the solid lines represent the mean of 30 realizations of noise. For the 40% CoA phantom ( Fig. 6 ), while showing a qualitatively correct behavior, the averaged catheter amplitude of the oscillation  Under stress conditions, the peak pressure difference is recovered with a good accuracy. There is a negative pressure difference lobe that was not correctly recovered by any of both methods. It must be noted that in the particular case of the 40% CoA phantom, artifacts appeared in the velocity measurements, most likely connected to issues with the experimental setup, like air bubbles.
The 40% CoA phantom configuration was repeatedly scanned at different resolutions (results not reported) and the corresponding estimated pressure differences showed similar characteristics in all cases. The difference in the width of the systolic peak interval between the 4D Flow data and the catheter data is also likely to be caused by such issues with the experimental setup.
The effect of noise at rest conditions (lowest VENC value, hence weakest noise standard deviation) is negligible. Under stress, the variability of the noisy results is slightly increased. The mean value accurately matches the time profile obtained for the original data. The STE and PPE results are indiscernible.
In the case of the 50% CoA phantom, Fig. 7 , the PPE and STE results are near identical, and slightly overestimate the catheter results under rest conditions during systole.
Note that in this case, the standard segmentation V +0 would lead to less overestimation (cf. Fig. 5 ). Under stress conditions, similarly to the 40% stress results, the 4D Flow-derived pressure differences show a double peak. The temporal duration of the peak appears to be overestimated, but the systolic amplitude and slope leading to the peak is recovered accurately. The effect of noise results in an large variability of the results, in particular at stress, but the mean pressure difference matches precisely with the original data case.
Results from the 60% CoA phantom under rest conditions show an excellent agreement between the catheter data and the pressure difference computed with the STE method ( Fig. 8 (a)).
The PPE method underestimates the peak pressure difference. The discrepancy between the pressure difference reconstruction and catheter measurements increases for stress conditions. Both methods significantly underestimate the catheter pressure differences and yield flattened profiles around the peak location. Particularly the estimation of the PPE method severely deteriorates under the given flow conditions. The addition of noise results in a large spread of the results during systole, but again, the mean pressure difference matches well with the results obtained for the original data.

Patient data
Pressure differences obtained by catheterization and from 4D Flow are shown in Fig. 9 for both patients. The locations where the pressure difference is evaluated is indicated by the green spheres in Fig. A.13 (right column). An excellent agreement of the STE pressure difference with catheter data was found for Subject 1 during systole. The local extrema after t = 0 . 4 s are underestimated. While similar qualitative agreement was found with the PPE method, it significantly underestimates the pressure difference during systole.
Subject 2 exhibits good qualitative and quantitative agreement between catheter data and numerical pressure difference reconstruction. However, the pressure difference peak observed by catheterization is too steep to be captured by the time resolu- Fig. 9. 4D Flow pressure differences computed with STE and PPE compared to catheter data (line: mean profile over all measured cycles, shaded area: ±2 standard deviation bands) for Subject 1 and Subject 2. tion of the 4D Flow protocol. The resulting maximum value lies below the catheter value, possibly because no velocity image was recorded matching exactly with the maximum pressure difference.
Note that the lengths of the cardiac cycles differ significantly between catheterization and MRI scans of both patients. This indicates a change in the heart rate which is likely to contribute to the differences between 4D flow and catheter pressure differences.
Reducing the diameter of the segmentation by one voxel was not possible in one of the patients due to the large voxel size with respect to the diameter. The original segmentations were thus not modified.

Discussion
Main findings. This study compared two relative pressure reconstruction methods, STE and PPE, in terms of accuracy and sensitivity with respect to image resolution, segmentation, noise, CoA severity and cardiac load (rest and stress). The main finding of this study is that the STE method applied to 4D Flow data provides a significantly closer agreement with catheter measurements in terms of instantaneous pressure differences than the PPE method at large pressure differences ( ≥ 20 mmHg ). In the clinical practice, the typical criterion for the necessity of surgical intervention is pressure differences ≥ 20 mmHg , hence the advantage of the STE method is of practical importance. Heuristic modification of the segmentation by deleting one layer of wall voxels significantly improves the results of the PPE method to match the accuracy of the STE method in most cases ( ≤ 30 mmHg ). In contrast, the STE method is considerably less sensitive to the segmentation. This may provide an advantage for clinical use in terms of robustness to inter-operator variability.

Effect of image resolution & segmentation.
Both methods were robust with respect to the image resolution for pressure differences below 30 mmHg, where the effect of resolution was minor. Estimates of larger pressure differences were slightly reduced under lower resolution.
The PPE results could be improved significantly in most cases by removing the outer layer of voxels at the vessel wall, and proved highly sensitive to the studied manipulation of the segmentation (i.e., comparing a high-fidelity segmentation with extensions or reductions by one layer of voxels at the walls). This is in line with Goubergrits et al. (2019) who also observed an improved accuracy of the PPE method by eliminating the outer layer of voxels. However, to the best of our knowledge there is no (theoretical) rationale supporting that the accuracy may be improved in this way, so these observations are purely empirical. In particular in the absence of catheter ground truth data, there is no clear rule for deciding when eliminating the boundary voxels is beneficial. In contrast, the STE results appeared to be less sensitive with respect to the segmentation, and are not always improved by choosing a heuristically manipulated segmentation. Indeed, at high resolution, the standard segmentation seems to be generally the best choice if extreme pressure differences ∼ 45 mmHg are not expected.
The theory offers an explanation for the different sensitivities of the methods: the PPE method explicitly uses the measured data in the boundary conditions. These boundary values determine the solution in the interior, and errors due to partial-volume effects or low velocity-to-noise ratio (VNR) can contaminate the solution. The STE method completely avoids such pressure boundary conditions, which is likely the reason for the improved robustness. Dilated segmentations (V +1 ) add no-flow voxels with insignificant VNR, which introduce spurious information into the estimation problem. The pressure gradient computation is required to accommodate to such unphysical conditions, hence the deterioration of the results with both methods.
Effects of noise. Both the PPE and the STE methods proved robust with respect to elevated noise in the sense that the mean pressure difference appeared to be unbiased and matched precisely with the low-noise original phantom data. This is consistent with the theoretical bias analysis performed in Bertoglio et al. (2018) . At high VENC values, the synthetic noise (its intensity proportional to the VENC) led to significant variability of the results around the mean.
The issue of high VENC values causing low VNR in regions of low flow velocities much smaller than the VENC could be alleviated by Dual-VENC techniques ( Nett et al., 2012;Ha et al., 2016;Callaghan et al., 2016;Carrillo et al., 2018 ), allowing for lower VENC-values (hence lower noise), but at increased scan times.
Regularity & turbulence. A further advantage of the STE method, guaranteed by the theory, is the fact that it searches the pressure in its natural space , while the PPE method imposes artificial higher regularity requirements, implying that its solutions are generally smoother than naturally required. Instead, the STE method allows for stronger spatial variations in the pressure and is capable of computing more realistic pressure fields. This fact is likely to cause the better agreement of the STE method with the catheter data, in relative and in absolute terms, in the regime of large pressure differences produced by the most severe 60% CoA. Under such conditions, the flow is strongly convection-dominated, resulting in sharp gradients and large changes in the flow. These are precisely the flow features under which the PPE method must be expected to deteriorate due to the mentioned model assumptions. The strong spatial oscillations visible in the pressure fields computed for the 60% CoA under stress, located at the edges of the jet ejected from the CoA, are presumably caused by the very large, underresolved velocity gradients. Furthermore, it seems likely that in this case the flow periodically transitions to turbulence around systole and subsequently relaminarizes. No measurements of turbulent quantities have been made in this study to probe the presence of turbulent flow, and due to their inherently phase-averaged character, the 4D Flow data ( Fig. A.12 ) do not allow deducing the presence of turbulences. Future work should both study the possible occurrence of turbulence in the present CoA phantom, i.e., by imaging the Reynolds stress tensor ( Haraldsson et al., 2018;Walheim et al., 2019 ), and address the extension of the STE method to account for turbulent effects as was done in other studies for other pressure estimation methods, see e.g. Marlevi et al. (2020) ; Haraldsson et al. (2018) ; Ha et al. (2019) .
In-vivo results. The STE and the PPE pressure reconstruction methods were also applied to real patient data. For one of the patients, the STE method showed a great improvement over the PPE method. Both methods showed satisfactory results for the second patient. From the findings in the phantom experiments, the differences between PPE and STE in patient one is most likely due to strong convective effects.

Limitations.
A limitation of the study was the lack of availability of real low resolution MRI data for all scenarios, hence requiring synthetic subsampling of the high resolution data. Furthermore, only one segmentation was used for all cardiac phases. Another aspect to consider is the approach used for assessing the sensitivity to the segmentations, i.e., automatic manipulation instead of inter-operator variability. The latter was expected to be insignificant due to the excellent contrast and high SNR of the anatomic phantom images, hence an automatized voxel increment or decrement was adopted. Regarding the patient data, significant effects of inter-operator variability in the segmentation process are likely.
In addition, the comparison of catheter data with MRI scans is limited by the following observations. The locations where the catheter recorded the pressure during catheterization are only known approximately. A mismatch of the catheter positions with the locations selected for evaluating the computed pressure gradient can introduce significant errors. The pressure fields can be seen in Figs. 2-4 to vary near the extraction locations, both in the AAo and the DAo. The PPE and STE methods have the advantage, that the pressure difference can be evaluated in various locations, after inspection of the pressure field. Furthermore, the catheter measurements showed important variabilities among the cardiac cycles and their sampling rate and accuracy ( ∼ 1 mmHg ) were limited, which restricts the authority of catheterization as ground truth in our study. As an invasive technique, it is possible that the presence of the catheter in the vessel disturbs and alters the aortic flow during catheterization, while the 4D Flow data was acquired immediately after without the catheters.
In the patient study it was seen that the heart rate changed significantly between catheterization and the MRI scan, hence possibly also affecting the outcome of this comparison.

Conclusion
In conclusion, in our phantom study, the STE method delivered results that were more accurate and robust with respect to the segmentation than the PPE method, in particular in severe cases of CoA and under more challenging flow regimes. Very large pressure gradients were underestimated. Low image resolution increased this effect and otherwise was of minor importance.
By heuristically eliminating the outermost layer of voxels of the segmentation, the PPE method could be significantly improved and was on a par with the STE method for low to moderately large pressure differences, but remained more sensitive to the segmentation. The STE method still provided a superior agreement with catheterization, without the need of artificially altering geometrically accurate segmentations, which proved beneficial only in the most severe case with very large pressure gradients. Both methods proved robust under noise. The mean pressure differences were unbiased, considerable variability was found for very high VENC values.
The data will be made available upon request of the interested parties.

Declaration of Competing Interest
The authors declare no conflict of interest.

A1. Catheterization data
The phase-averaged catheterization data obtained in the AAo and DAo locations is displayed in Fig. A.10 for the phantoms and in Fig. A.11 for the patients. The graphs include the catheter uncertainty across the measured cycles in terms of ±2 standard deviation bands (shaded areas).
The AAo pressure profiles are similar for all CoA phantoms under the same cardiac load conditions. At rest, the peak amplitude slightly increases with increasing CoA severity. Under stress conditions, it stands out that the 40% CoA data has a higher minimum pressure than all other cases. The same case features a slight phase shift between the AAo and DAo signals. The DAo pressure visibly decreases in all cases with increasing severity. The catheter pressure variability across cycles also increases with CoA severity. The spread is smaller under stress conditions.
The catheter pressures measured in the patients are shown in Fig. A.11 . In contrast to the phantom data, marked phase shifts appear between the AAo and DAo pressures. This results in very large discrepancy between the instantaneous pressure differences as discussed here, and clinical peak-to-peak pressure differences, which only compare the difference between the overall maximum pressures. In Subject 2, the instantaneous pressure difference reaches 40 mmHg due to the phase shift, while the peak-to-peak pressure difference is close to 0 mmHg.

A2. 4D Flow data
Peak velocity streamlines in the vicinity of the CoA are shown in Fig. A.12 for all phantoms under rest and stress conditions. The displayed 4D Flow data correspond to the instant of maximum catheter pressure differences, see Figs. 6 -8 . Also included in the figures are isosurfaces of the Q-criterion, commonly employed for identifying vortices ( Haller, 2005;Dubief and Delcayre, 2000 ). The levels of the isosurfaces were adapted manually, such that the largest, dominant vortices were identifiable and fractional, cluttered regions minimized.  In all cases, the CoA produces a strong jet impinging on the outer wall of the DAo. Strong recirculation zones, driven by the deflected flow, develop below the jets. Large vortex rings envelop the jets. For each severity, under stress conditions the maximum jet velocity is increased by 23% to 30% (at the time of the maximum catheter pressure difference).
In the 40% CoA phantom, the vortex ring associated with the large recirculation bubble is comparatively smooth and coherent. The flow quickly realigns and remains approximately parallel a short distance downstream of the recirculation zone. In comparison, the 50% CoA exhibits larger recirculation zones, both under the rest and stress conditions, and increasingly complex circulating flow patterns. The streamlines and Q-criterion isosurfaces appear less smooth and coherent than in the 40% case. For the 60% CoA, the recirculation zones are significantly larger with more chaotic and less coherently oriented streamlines. At rest, a large vortex, associated with the large contiguous Q-criterion ring structure, can be appreciated, superimposed by less orderly features. Under stress conditions, in Fig. A.12 (f), a single large vortex cannot be identified, in contrast to all other scenarios. The Q-criterion isosurface is significantly fractured, forming a number of tight vortex rings around the jet core. The streamlines in the recirculation zone are highly disorganized. Fig. A.13 shows angiographic images of the patients paired with 4D Flow measurements corresponding to the peak systole. The flow appears to be highly regular, with high velocities in the narrow sections (positions (4)) and some recirculation in Subject 2 at position (5).