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

Medical Image Analysis 74 (2021) 102195 Contents lists available at ScienceDirect Medical Image Analysis journal homepage: www.elsevier.com/locate/media Challenge Report Validation of 4D Flow based relative pressure maps in aortic flows David Nolte a , b , Jesús Urbina c , d , e , Julio Sotelo c , e , f , g , Leo Sok a , Cristian Montalba c , e , Israel Valverde h , Axel Osses b , e , Sergio Uribe c , d , e , Cristóbal Bertoglio a , b , ∗ a Bernoulli Institute, University of Groningen, Groningen, 9747AG, The Netherlands b Center for Mathematical Modeling, Universidad de Chile, Santiago, 8370456, Chile c Biomedical Imaging Center, Pontificia Universidad Catolica de Chile, Santiago, 7820436, Chile d Department of Radiology, School of Medicine, Pontificia Universidad Católica de Chile, Santiago, 833002, Chile e Millennium Nucleus for Cardiovascular Magnetic Resonance, Santiago, 7820436, Chile f School of Biomedical Engineering, Universidad de Valparaíso, Valparaíso, Chile g Department of Electrical Engineering, Pontificia Universidad Catolica de Chile, Santiago, 7820436, Chile h Hospital Universitario Virgen del Rocío, Sevilla, 41013, Spain a r t i c l e i n f o Article history: Received 29 October 2020 Revised 11 June 2021 Accepted 22 July 2021 Available online 15 August 2021 Keywords: 4D Flow Pressure difference Catheter Clinical and experimental validation a b s t r a c t While the clinical gold standard for pressure difference measurements is invasive catheterization, 4D Flow MRI is a promising tool for enabling a non-invasive quantification, by linking highly spatially resolved velocity measurements with pressure differences via the incompressible Navier–Stokes equations. In this work we provide a validation and comparison with phantom and clinical patient data of pressure difference maps estimators. We compare the classical Pressure Poisson Estimator (PPE) and the new Stokes Estimator (STE) against catheter pressure measurements under a variety of stenosis severities and flow intensities. Specifically, we use several 4D Flow data sets of realistic aortic phantoms with different anatomic and hemodynamic severities and two patients with aortic coarctation. The phantom data sets are enriched by subsampling to lower resolutions, modification of the segmentation and addition of synthetic noise, in order to study the sensitivity of the pressure difference estimators to these factors. Overall, the STE method yields more accurate results than the PPE method compared to catheterization data. The superiority of the STE becomes more evident at increasing Reynolds numbers with a better capacity of capturing pressure gradients in strongly convective flow regimes. The results indicate an improved robustness of the STE method with respect to variation in lumen segmentation. However, with heuristic removal of the wall-voxels, the PPE can reach a comparable accuracy for lower Reynolds’ numbers. ©2021 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license ( http://creativecommons.org/licenses/by/4.0/ ) 1. 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- ∗Corresponding author. E-mail address: [email protected] (C. Bertoglio). sive, non-exempt of risk, exposes the patient to X-rays and is difficult to reproduce. Non-invasive measurement techniques exist for the blood flow velocity, namely Doppler echocardiography and phase-contrast magnetic resonance imaging (PC-MRI). Since the velocity is related to the pressure via physical laws, i.e., the Navier–Stokes equations (NSE), estimates of the pressure can be computed from velocity data. The Simplified Bernoulli formula (semi-empiric, derived from the NSE) is commonly used with Doppler echocardiography data to estimate CoA pressure differences, but has severe limitations: it is operator dependent, has poor acoustic windows and spectral broadening, and in addition overestimates the peak velocities by up to 25% ( Hoskins, 1996 ). MRI is able to obtain imhttps://doi.org/10.1016/j.media.2021.102195 1361-8415/© 2021 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license ( http://creativecommons.org/licenses/by/4.0/ ) D. Nolte, J. Urbina, J. Sotelo et al. Medical Image Analysis 74 (2021) 102195 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 ( Markl et al., 2012; 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) ( Bertoglio et al., 2018 ) based on integral linear momentum conservation, or the virtual WERP ( v WERP) method ( Marlevi et al., 2019; 2021 ), 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 ( Urbina et al., 2016 ), 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. 2. Theory 2.1. 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., ρ∂  u ∂t + ρ(  u ·∇ )  u = −∇ p + μ u (1) where ρis the density of the fluid and μits dynamic viscosity,  u : ×[0 , T ] → R 3 denotes the velocity vector field (to be obtained by means of 4D Flow) and p : ×[0 , T ] → R is the pressure field.  represents the computational domain obtained from segmenting the 4D Flow images (see Section 3.3 and Fig. 1 ), its boundary— corresponding to the vessel wall—is denoted ∂. The considered time interval [0 , T ] represents one cardiac cycle. Relative pressure maps are reconstructed from 4D Flow data by means of the Pressure Poisson Estimator (PPE) approach ( 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: −∇p k = ρ  u k − u k −1 t + ρ(  u k ·∇)  u k −μ u k . (2) The indices 1 ≤k ≤Ndenote the time snapshot of the measurements and tthe 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 , R k := ρ  u k m − u k −1 m t + ρ(  u k m ·∇)  u k m −μ u k m . (3) yields a pressure estimate from its approximate gradient ∇ ˆ p k ≈ ∇p k , given by −∇ ˆ p k = R k . (4) 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 ( Markl et al., 2012 ). 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). 2 D. Nolte, J. Urbina, J. Sotelo et al. Medical Image Analysis 74 (2021) 102195 Fig. 1. Pressure recovery toolchain: (a) Segmentation of anatomic images (shown: 2D slice of 3D image) identifies the 3D lumen ( ) and the vessel walls ( ∂), averaged over the cardiac cycle. (b) A structured tetrahedral mesh is created from the segmentation. (c) The 4D Flow velocity data is represented on the mesh. (d) Relative pressure maps are computed from velocity data (shown: cuts though approx. centers of the vessels). Example data of the 60% CoA phantom at rest, at time of peak systole, for a voxel size of 2mm. 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. 2.2. 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) , −ˆ p k = ∇ ·R k . (5) Solving Eq. (5) requires boundary conditions (BCs) on the entire boundary surface ∂of the computational domain . A priori, no physical BCs for the pressure are known. An artificial Neumann BC can be obtained by projection of Eq. (4) , restricted to ∂, to the outward unit normal vector on the wall,  n ,  n ·∇ ˆ p k =  n ·R k . (6) 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. 2.3. 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 − w −∇ ˆ p k = R k in  ∇ · w = 0 in   w =  0 on ∂. (7) 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. 3. Materials and methods 3.1. Aortic phantom study 3.1.1. 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 Uover 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 ◦. 3 D. Nolte, J. Urbina, J. Sotelo et al. Medical Image Analysis 74 (2021) 102195 3.1.2. 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. 3.1.3. 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 . 3.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 Table 1 MR acquisition parameters. in-vitro experiments in-vivo experiments FOV (mm) 200 ×200 ×114 270 ×270 ×125 Matrix size 224 ×224 ×127 144 ×144 ×50 Recon. voxel (mm) 0 . 9 ×0 . 9 ×0 . 9 1 . 9 ×1 . 9 ×2 . 5 TFE factor rest: 2, stress: 1 rest: 2 Cardiac phases 25 28 Time resolution (ms) rest: 32, stress: 18 rest: 26 VENC (cm/s) rest: 160–400,stress: 250–500 rest: 300 TE/TR (ms) 3.7/6.4 2.4/3.8 Flip angle (deg) 6.5 5 Scan time (min) 18–22 19 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”). 3.3. Segmentation and mesh generation The 4D Flow data sets were processed using an in-house MATLAB library (The MathWorks, Natick, MA), similarly to previous studies ( Sotelo et al., 2016; 2018 ). 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 4 D. Nolte, J. Urbina, J. Sotelo et al. Medical Image Analysis 74 (2021) 102195 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). voxel size Mean pressure difference, PPE Mean pressure difference, STE h (mm) P 1 P 2 P 3 P 1 / P 1 stab P 2 / P 1 P 3 / P 2 2.0 1.5593 1.6714 ( + 7%) 1.6746 ( + 0.2%) 2.0755 2.0816 ( + 0.3%) 2.0790 (-0.1%) 1.4 1.7252 1.7986 ( + 4%) 1.8011 ( + 0.1%) —— — 0.9 1.8661 1.9108 ( + 2%) 1.9127 ( + 0.1%) —— — 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. 3.4. 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 ( Alnæs 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). 3.5. 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. 3.5.1. 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 velocity- to-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. Table 3 Numbers of degrees of freedom with different PPE and STE discretizations and voxel sizes. Exemplary data of the 50% CoA at rest, V +0 segmentation. voxel size PPE, degrees of freedom STE, degrees of freedom h (mm) P 1 P 2 P 3 P 1 / P 1 stab P 2 / P 1 P 3 / P 2 2.0 18,135 134,860 443,882 72,540 422,715 1,466,506 1.4 53,426 406,316 1,347,732 213,704 —— 0.9 203,345 1,575,443 5,259,104 813,380 —— 5 D. Nolte, J. Urbina, J. Sotelo et al. Medical Image Analysis 74 (2021) 102195 Table 4 Computation times with different PPE and STE discretizations and voxel sizes. LU was used unless noted otherwise. Exemplary data of the 50% CoA at rest, V +0 segmentation. voxel size PPE, CPU time (s) STE, CPU time (s) h (mm) P 1 P 2 P 3 P 1 / P 1 stab P 2 / P 1 P 3 / P 2 2.0 4.8 11.3 42.6 10.3 88.3 17627 b 1.4 7.5 39.1 465.8 33.4 —— 0.9 21.9 332.0 a 1720.4 a 245.7 —— a iterative solver CG/AMG. b iterative solver GMRES/AMG–Schur/pressure mass matrix. 4. Results 4.1. Phantom study 4.1.1. 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. 4.1.2. 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- Fig. 2. Pressure fields of the 40% CoA phantom obtained with the STE and PPE methods under (a) rest and (b) stress conditions, using the standard segmentation V +0 at the original voxel size of 0.9mm. Data is shown for the time of the observed peak pressure difference, t = 0 . 22 s at rest and t = 0 . 14 s at stress. The blue spheres show the locations between which the pressure difference is computed, averaged over the sphere volume. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Fig. 3. Pressure fields of the 50% CoA phantom obtained with the STE and PPE methods under (a) rest and (b) stress conditions, using the standard segmentation V +0 at the original voxel size of 0.9mm. Data is shown for the time of the observed peak pressure difference, t = 0 . 22 s at rest and t = 0 . 12 s at stress. The blue spheres show the locations between which the pressure difference is computed, averaged over the sphere volume. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Fig. 4. Pressure fields of the 60% CoA phantom obtained with the STE and PPE methods under (a) rest and (b) stress conditions, using the standard segmentation V +0 at the original voxel size of 0.9mm. Data is shown for the time of the observed peak pressure difference, t = 0 . 22 s at rest and t = 0 . 16 s at stress. The blue spheres show the locations between which the pressure difference is computed, averaged over the sphere volume. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) 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). 6 D. Nolte, J. Urbina, J. Sotelo et al. Medical Image Analysis 74 (2021) 102195 Fig. 5. Comparison of 4D Flow peak pressure differences, δp, obtained with the STE (left column) and the PPE (right column) methods with catheterization, for all studied in-vitro severities, segmentations and resolutions. Experiment indices: 1: 40% CoA, rest; 2: 40% CoA, stress; 3: 50% CoA, rest; 4: 50% CoA, stress; 5: 60% CoA, rest; 6: 60% CoA, stress. V +0 refers to the reference segmentation; the segmentation V +1 extends V +0 by 1 voxel in each direction, in V −1 the outer layer of voxels of V +0 is deleted. Top, middle, bottom row: results obtained for 0.9mm, 1.4mm, 2mm voxel size. The error bars indicate two standard deviations (95.45% confidence interval) of the catheterization pressure differences (cf. Fig. A.10 for the phase-averaged catheter data and associated uncertainties). 4.1.3. 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 δpestimates 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 δpis 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 δpestimates and therefore reduces the error in this case with respect to catheterization. In experiment 1 Fig. 6. 40% CoA phantom pressure differences obtained with STE and PPE compared to ground truth catheter data (line: phase average; gray shaded area: ±2 standard deviations), under (a) rest and (b) conditions, for voxel size 2.0mm and the V −1 segmentation. Left figures: original data. Right figures: pressure differences resulting under additional synthetic noise with indicated V ENC/σ(lines: mean over 30 realizations; shaded areas: ±2 standard deviations). (40% CoA at rest), the mean catheter δpis underestimated by 50%. However, the catheter is associated with a large variability, and the δpestimate 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 δpin 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 δpestimates 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. 4.1.4. 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 7 D. Nolte, J. Urbina, J. Sotelo et al. Medical Image Analysis 74 (2021) 102195 Fig. 7. 50% CoA phantom pressure differences obtained with STE and PPE compared to ground truth catheter data (line: phase average; gray shaded area: ±2 standard deviations), under (a) rest and (b) conditions, for voxel size 2.0mm and the V −1 segmentation. Left figures: original data. Right figures: pressure differences resulting under additional synthetic noise with indicated V ENC/σ(lines: mean over 30 realizations; shaded areas: ±2 standard deviations). Fig. 8. 60% CoA phantom pressure differences obtained with STE and PPE compared to ground truth catheter data (line: phase average; gray shaded area: ±2 standard deviations), under (a) rest and (b) conditions, for voxel size 2.0mm and the V −1 segmentation. Left figures: original data. Right figures: pressure differences resulting under additional synthetic noise with indicated V ENC/σ(lines: mean over 30 realizations; shaded areas: ±2 standard deviations). is generally underestimated by 3 to 7 mmHg with both methods under rest conditions, without added noise. The variability of the catheter is large under rest conditions, and in the positive lobes, the 4D Flow results lie within the ±2 σbounds. 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. 4.2. 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. 8 D. Nolte, J. Urbina, J. Sotelo et al. Medical Image Analysis 74 (2021) 102195 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. 5. 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 9