Estimation of Cardiovascular Relative Pressure Using Virtual Work-Energy

Many cardiovascular diseases lead to local increases in relative pressure, reflecting the higher costs of driving blood flow. The utility of this biomarker for stratifying the severity of disease has thus driven the development of methods to measure these relative pressures. While intravascular catheterisation remains the most direct measure, its invasiveness limits clinical application in many instances. Non-invasive Doppler ultrasound estimates have partially addressed this gap; however only provide relative pressure estimates for a range of constricted cardiovascular conditions. Here we introduce a non-invasive method that enables arbitrary interrogation of relative pressures throughout an imaged vascular structure, leveraging modern phase contrast magnetic resonance imaging, the virtual work-energy equations, and a virtual field to provide robust and accurate estimates. The versatility and accuracy of the method is verified in a set of complex patient-specific cardiovascular models, where relative pressures into previously inaccessible flow regions are assessed. The method is further validated within a cohort of congenital heart disease patients, providing a novel tool for probing relative pressures in-vivo.


Results
Virtually probing hemodynamic flow to compute relative pressure. The complete theory and process for computing relative pressure using the virtual Work-Energy Relative Pressure (νWERP) formulation is described in detail in the Online Methods, and in Fig. 1. Briefly, from an acquired 4D flow field, the method requires segmentation of the main vascular flow compartment. Inlet and outlet planes are then selected between which relative pressures are to be estimated. A virtual field to interrogate the anatomical region of interest is computed based on the solution to a Stokes flow problem, isolating the relative pressure of interest. Combining the acquired flow and virtual fields, νWERP evaluates virtual work-energy equations to identify relative pressure over time. While νWERP is not inherently connected to any specific imaging modality, here we focus on integration of MRI-based 4D flow. The method could however easily be adapted to other full-field imaging modalities such as refined 3D ultrasound 20 or contrast-enhanced computed tomography 33 . In-silico evaluation in cases of aortic coarctation and aortic dissection. To evaluate the νWERP method in a set of challenging realistic flow scenarios, a set of personalised in-silico tests were performed using Computational Fluid Dynamics (CFD) (for a full description of the models used, see Online Methods). Analysis using in-silico image data provides a powerful test-bed for new techniques, as velocity and pressure information can be defined at high spatiotemporal resolution throughout the entire vasculature, giving a clear comparative ground truth. Two pathological conditions are considered in this work: a model of a coarctation of the aorta (CoA) 34 and a model of an acute type B aortic dissection (AAD) 35 . Both reflect pathological cases where relative pressure is of direct diagnostic importance 6,36 . In-silico results were sampled using spatiotemporal resolution typical of 4D flow data (2 mm 3 voxel size images at 22 and 32 phases over the cardiac cycle for CoA and AAD, respectively) as illustrated in Fig. 2. Effects of realistic signal-to-noise (SNR) levels were evaluated by amending the sampled in-silico results with 100 different noise fields at low-(SNR = 30) and high-noise (SNR = 10), respectively. With the in-silico image data generated, the flow fields were fed into the νWERP work-flow ( Fig. 1) and relative pressures were estimated between sets of vascular planes. Figure 3 shows the CoA results, illustrating planes between which relative pressures were estimated and providing graphs comparing νWERP estimates to the simulated (CFD) ground-truth. In general, a distinct correspondence can be seen between νWERP and true (CFD-based) relative pressure, with errors consistently below 20% for all evaluated planes and SNR, (peak systolic error < 1.0 mmHg), with a mean error of 10% at SNR = ∞ (averaged over the entire cardiac cycle). In particular, the method shows the ability to accurately probe the relative pressure between multiple sites, with νWERP estimates showing no evidence of bias due to vascular branching. A slight increase in error (mean error average 13/18/18% at SNR = ∞/30/10) can be seen into the left common carotid artery (plane C in Fig. 3), potentially due to the smaller mean diameter of 6 cm (~3 voxels).
νWERP estimates demonstrate robustness to noise -important for clinical use -showing only minor increases in mean error (13/15% at SNR = 30/10). For all noise-levels the νWERP estimates consistently follow the general behaviour of the simulated pressure through the cardiac cycle. Focusing on the clinically important coarctation (D → E, Fig. 3), νWERP estimates give mean errors of below 16/18% for SNR = 30/10, indicating the potential value of this technique to compute quantities of important clinical value. Figure 4 similarly presents results for a selection of anatomical planes in the AAD model. In this instance, the flow is significantly impacted by the septum separating true and false lumen as well as the numerous connecting tears linking flow between the lumina. However, with the use of a virtual field, νWERP was able to isolate the relative pressure of interest (a visualisation of the virtual field into true and false lumina is given in Supplementary  Fig. 3). As with the CoA, νWERP illustrates the capacity to probe relative pressure between any of the selected points in the dissected aorta, with a mean error of below 14/16/21% for all evaluated regions at SNR = ∞/30/10. With noise, narrower regions show a slight deterioration in accuracy such as seen for the right subclavian artery (plane C, lumen diameter ~4 voxels) and the descending aortic part of the true lumen (plane E, lumen diameter ~3 voxels), showing mean errors of 23/24/31% and 16/20/32% at SNR = ∞/30/10, respectively. However, νWERP estimates still follow the temporal behaviour of the true (CFD-based) relative pressure.
A particularly complex case is provided for plane F: here, the right renal artery is fed by the true lumen, before continuing into the false lumen and the ascending thoracic aorta through connecting tears. Given sufficient spatial sampling to capture the narrow anatomy (in this case 1 mm 3 ), νWERP shows the ability to estimate relative pressure with a mean error of 12%. This again exemplifies the strength of the proposed approach, and underlines the versatility of using virtual fields to probe relative pressures even in highly complex cardiovascular morphologies. For the above evaluated in-silico patient-specific cases, the application of conventional flow based relative pressure estimates become problematic, specifically since the test cases considered expand beyond the usual intended use of these methods. The routinely used Bernoulli estimate originates from a singular measure of peak velocity, and similarly the extension into the unsteady Bernoulli 37 neglects a large portion of the flow field and is not optimised for traces into vasculatures of turbulent or low flow. Using work-energy estimates without the use of a virtual field 32 also deteriorate with vascular branching or low flow regimes. For these reasons, the above methods perform poorly for both CoA and AAD, with none of the evaluated relative pressures successfully assessed. Comparing with νWERP estimates into the true lumen of the AAD (A → E, Fig. 4) render mean errors of 100%, 88%, and > 200% for simplified Bernoulli, unsteady Bernoulli and conventional work-energy estimates, respectively (a complete numerical comparison is given in Supplementary Table 3). With more refined approaches such as the PPE showing reduced accuracy even in simplified analytical setups 29 , the added value of the proposed νWERP approach is highlighted, with a gained ability to probe relative pressure in regions inaccessible by other proposed methods.
In-silico evaluation of image resolution and SNR. A core practical concern for image-driven techniques is their dependence on spatiotemporal resolution and image quality. Using the cardiovascular simulations, where image spatial and temporal resolution can be freely adapted, it is feasible to explore the impact of resolution and noise on νWERP performance. For this analysis, two regions spanning from ascending aorta to descending true and false lumen of the AAD (A → D and A → E, Fig. 4) were selected, with temporal and spatial sampling within clinical acquired image ranges (for details of the model generation, see Online Methods).
In general, νWERP performance improves with increasing spatiotemporal sampling, as seen in Fig. 4e. For the mean error, this behaviour is seen for both the wide false and the narrow true lumen, respectively, despite the significant flow differences. Considering image noise, the results indicate that νWERP can assess relative pressure in all evaluated noise scenarios with regards to both peak and mean relative pressure, even if image noise exhibits a more pronounced effect on the estimates into the narrower true lumen (at 2 mm 3 and 32 time frames, mean error is 16/20/32% in the true lumen, compared to 10/10/12% in the wider false lumen at SNR = ∞/30/10). With a reduction in the number of voxels, increasing boundary layer effects as well as singular noise will influence the physics of the investigated flow and impact νWERP estimation accuracy. The sought estimate accuracy will thus require consideration of both the image acquisitions settings, as well as the vasculature of interest. Regardless, νWERP still shows reasonable behaviour in the narrow true lumen when using routine imaging settings (2 mm 3 voxel size and 32 temporal phases), where the peak relative pressure was estimated with an accuracy of 2 ± 2/6 ± 5% at SNR = 30/10. Analyses at the noise-free state (SNR = ∞) indicate a slightly higher dependence on temporal sampling, where a twofold increase in temporal sampling renders an average decrease in mean error of 32%, whereas a similar increase in spatial sampling gives an average decrease in error of 19%. However, the effect seems slightly reversed when image noise is added, where a twofold increase in spatial and temporal sampling renders a decrease in mean error of 34/39% and 31/31% at SNR = 30/10, respectively. For the peak error the effect of increased spatial does not appear as clear, even though temporal refinements seems to have a defined effect at all noise levels (twofold temporal refinement rendering a decrease in peak error of 26/30/32% and 39/40/43% at SNR = ∞/30/10 for false and true lumen, respectively). Consequently, the appropriate choice of spatiotemporal image sampling will be dependent on the physical nature of the investigated vasculature and demand for accuracy; narrow vessels require increased spatial sampling, whereas transient cardiac events (such as those during ventricular ejection in early systole) require increased temporal sampling.
An additional analytical spatiotemporal convergence analysis is provided in Supplementary Material.

In-vivo validation with invasive cardiac catheterisation.
To establish its clinical applicability, we compared νWERP estimates against gold-standard invasive pressure measurements acquired in-vivo during combined cardiac catheterisation and cardiac MRI (XMR). 4D flow MRI (voxel size ~2 mm 3 , 22 time samples) was collected in a cohort of five anaesthetised adolescent patients with complex congenital heart disease (for patient characteristics, see Supplementary Table 1). During the same XMR session, absolute invasive pressures were recorded through catheterisation at ascending and mid-descending or diaphragm level of the aorta, respectively (positions registered using fluoroscopy). For a complete description of data preparation, see Online Methods. All patients included in this study participated under informed consent, with data collection and study approved and true lumen (E). Note that the analysis is only performed when the selected spatial image sampling is able to resolve the anatomy. Data for SNR = ∞/30/10 is given in red/blue/green, respectively. Each table entry presents the mean error (lower left) and the error at peak relative pressure with standard deviation (upper right).
Scientific RepoRts | (2019) 9:1375 | https://doi.org/10.1038/s41598-018-37714-0 by the Regional Ethics Committee, South East London, United Kingdom (REC, 10/H0802/65). Additionally, all methods and experiments were performed in accordance with relevant guidelines and regulations. Results for all six subjects are presented in Fig. 5. On average, νWERP could capture the temporal shift in relative pressure, with estimates corresponding well to gold-standard catheter measurements (for a graphical comparison of all investigated patients, see Supplementary Fig. 4). A mean error of 16 ± 13% is reported for all subjects, corresponding to an absolute mean error < 1.0 mmHg. The error varied slightly over time, with the initial systolic pressure peak captured with defined accuracy in all patients (mean error < 18%, corresponding to an accuracy within 1.2 mmHg). The return from the systolic pressure peak is also well captured by νWERP as seen in both Fig. 5 and Supplementary Fig. 4. With the relative pressure in general, and the peak relative pressure specifically acting as an important clinical biomarker 1,3 , νWERP shows direct utility for clinical diagnostics. Additionally, with this accuracy reported over the aortic arch where several vessel branches are present, as well  Table summarising the results for the entire cohort, indicating mean error, root-mean square difference (RMSD), peak relative pressure by catheterisation and νWERP, respectively, as well as the corresponding absolute error. as into the narrow adolescent descending aortas (in average ~5 voxels over the vessel cross-section), νWERP enables recognised biomarkers to be applied in vascular regions where it has previously not been used, following limitations in the routinely and previously proposed methods.
Comparing alternative methods, these again perform with varying accuracy, following the unconventional application field for which none of the methods are optimised. Using simplified Bernoulli errors in maximum relative pressure of >100% were found, with the accuracy significantly deteriorating due to the simplifications imposed in the hemodynamic description of the flow. For the extension to Unsteady Bernoulli, the accuracy increases (error at peak systole of 67% or 6.2 mmHg), although it produces errors of above 100% in singular subjects due to the assumptions made on the assessed flow. With regards to the work-energy based approach without the virtual field mean errors in relative pressure of >300% were found, again following the branched anatomy as well as low diastolic flows. Lastly, the unsuitability of the PPE method to deal with the narrow adolescent lumina further highlights the novelty of the proposed νWERP method, which represents a definite improvement for hemodynamic pressure estimations and enables the probing of relative pressure in a range of vascular segments accessible by image-based 4D flow acquisition.
Cumulative analysis of νWERP estimates of relative pressure.  Relative pressure estimations for patient-specific coarcted and dissected aorta, respectively, comparing νWERP (y-axis) to simulated data (x-axis). (b) Relative error of νWERP estimates at peak relative pressure for all evaluated inlet/ outlet planes in patient-specific coarcted and dissected aortas. Data for SNR = ∞/30/10 is given by black circles, red crosses, and blue triangles, respectively. (c) Relative pressure estimation for the validation study, comparing νWERP estimates (y-axis) to catheter measurements (x-axis). The distribution of the catheter data is given by gray bars, with νWERP estimates superimposed by black crosses. (d) Absolute error of νWERP estimates at peak relative pressure for the entire validation study cohort. Data for estimates against the mean of all catheter measurements is given by the black circles, and relative errors of ±10% are given by the dashed black lines. Examining the in-vivo data, the deviation between νWERP and the catheter measurements shows more variability likely due to catheter acquisition issues 38 , inherent patient variability, and imaging artefacts. However, the correlation between νWERP and the catheter measures remains strong, with an average slope of m = 0.87 (Fig. 6c), representing a slight underestimation in relative pressure. Isolating the peak relative pressure (Fig. 6d), the absolute and mean error does seem to comply with a slight underestimation (0.6 ± 1.5 mmHg or 8 ± 26%). Importantly, considering the improvement in performance compared to other methods, and judging from the novelty in accessible regions of which the above νWERP validation is performed, the proposed approach represents a novel method for improved hemodynamic diagnostics in cardiovascular applications.

Discussion
In this work, we have presented a work-energy based method (νWERP) to derive relative pressure between any given vascular segment within an acquired flow field. This versatility is provided by the use of a virtual field to define the vascular region of interest, achieving accurate assessments of relative pressure. νWERP also provides a theoretical foundation for evaluating cardiovascular conditions with fewer assumptions than other methods, and can be further extended into anatomically complex scenarios currently inaccessible by today's techniques. Moreover, the technique works directly from images with minimal user input required. For clinical use, the validation study revealed that νWERP can capture important clinical biomarkers such as peak relative pressure drop at high accuracy, as well as the transient behaviour of the flow.
Relative pressure estimations are part of routine clinical practice when assessing scenarios such as the degree of valvular stenosis or the severity of aortic coarctation. However, current clinical techniques have shown significant limitations when applied to complex hemodynamic conditions. Doppler ultrasound systematically overestimates valvular pressure drops when severe stenosis skews the flow profile 16 . Conversely, invasive catheterisation of CoA patients is restricted to severe symptomatic cases due to the inherent risks of the procedure. Refined techniques such as the unsteady Bernoulli, iterative methods, PPE or WERP have fundamental limitations illustrated in our results and reported previously 29 . Comparatively, νWERP shows improved versatility, accuracy and robustness in-silico and in-vivo, extending the feasibility of an accurate and non-invasive approach into challenging anatomies such as the narrow subclavian and renal arteries of the in-silico evaluations; situations where implementation of previously proposed techniques would violate fundamental assumptions inherent to these methods. νWERP sensitivity also proved to be stable under realistic image noise scenarios (Figs 3 and 4).
In-silico evaluation also enabled the systematic analysis of the required spatiotemporal image resolution of νWERP. As noted previously, performance deteriorates slightly with small lumen diameters where spatial resolution can limit the capacity for images to represent vascular morphology. Moreover, limitations on data such as noise and boundary wall effects will affect νWERP performance. Optimal acquisition settings will also depend on the flow physics: when assessing domains of low flow with significant anatomical constriction, higher spatial resolution is needed to capture morphology and flow. In contrast, higher temporal resolution may be required if quick transient cardiac events occur. The versatility of νWERP allows for adaption to specific cases, providing advantages compared to other routine methods. This approach also allows for more complete hemodynamic analysis to be performed for improved disease risk stratification.
Our comparisons to invasive catheterisation revealed that νWERP can capture the temporal shifts of in-vivo aortic relative pressure. Results are similar to the in-silico evaluations, with the minor error increases potentially spurring from data acquisition and analysis (image segmentation or co-registration between imaging and catheterisation). Regardless, the clinically important biomarker peak pressure drop is still captured with high accuracy and comparing to other proposed methods, νWERP shows a defined improved performance.
The in-vivo cohort represents a particularly challenging study case, since the aortic anatomy of the adolescent subjects present narrow vascular segments (on average ~4 voxels through the diameter). Additionally, the Fontan procedure performed in 3 of the 5 subjects typically result in abnormal ascending aortic vessel morphology and flow profiles 39 , complicating the assessed hemodynamics. Based on this, the νWERP estimates are underlined, with potential applicability spanning over multiple clinical applications. Alternative methods performed poorly in this cohort, largely due to their inherent assumptions. Comparing with situations where simplified Bernoulli estimates are optimised 40 νWERP show similar or improved levels of accuracy. With that, νWERP represents a viable alternative to current techniques when applied in the same application field, and further enables investigations into more complex, previously inaccessible cardiovascular conditions. Even though not necessarily restricted to data acquired by 4D flow MRI, the presented νWERP formulation relies on full-field measures of blood flow in the vascular region of interest. An extension of similar work-energy based methods into planar 2D flow acquisitions has been proposed 16 , opening the possibility for similar strategies in the application of νWERP. However, coupling the refinement of νWERP with the rapid development of refined imaging techniques promises the potential of using 3D-based flow methods as standardised diagnostic tools.
To date, the diagnostic usage of relative pressure assessment has been applied for stenotic flow regimes in which the development of turbulent flow plays a role 41 . In severe cases of stenosis, methods incorporating turbulent flow behaviour have shown promising potential to non-invasively assess relative pressure 31,[41][42][43] . In its current form, νWERP does not include analysis of turbulent fluctuations, however the theoretical derivation does allow for the inclusion of a turbulence term originating from the incorporation of incoherent flow variations in the velocity field. Again using the concept of a virtual field, νWERP could thus enable even the assessment of turbulence-driven relative pressures into complex cardiovascular regimes. However, a full derivation, analysis, and method validation remains to be performed, and is beyond the scope of this current work. Instead νWERP is here presented as a method to expand on previous techniques by enabling hemodynamic probing into previously inaccessible vasculatures where the relative pressure drop is a result of combined kinetic, advective, and viscous flow behaviour. Improvements in medical imaging enable more complete and accurate investigation of hemodynamic biomarkers. Utilising these developments, νWERP represents a novel tool for assessing relative pressure between virtually any connected locations in cardiovascular flow images. The technique demonstrates accuracy and robustness for realistic image quality and resolution over a range of complex cardiovascular conditions in-silico and in-vivo. Providing improved flexibility and versatility, νWERP constitutes a promising tool opening the possibilities for expansion of hemodynamic biomarkers, both in research scenarios and standardised clinical care.

Methods
theoretical derivation of νWERP. The derivation of νWERP originates from the conservation of mass and momentum for an isothermal viscous Newtonian fluid given by the Navier-Stokes equations, i.e.
where v is the velocity of the fluid, p the fluid pressure, and μ and ρ are the dynamic viscosity and density, respectively. Derivation of the virtual work-energy equations for such a fluid can be achieved by multiplying (1) by an arbitrary velocity field -or virtual field -w and integrating over the entire fluid domain Ω, yielding: In this case, equation (3) describes the work-energy within the fluid in the direction of the virtual field w. In the case where w = v, equation (3) reduces to the classic work-energy equation and provides insight into the actual work and energy held within the fluid at a given time (see Donati et al. 32 ). Importantly, however, is that the sought relative pressure drop originates from the acquired v, and is not connected to the virtual field, w.
Using integration by parts (converting volume integrals over Ω to surface integrals over the domain boundary Γ with outward normal n), the solenoidal nature of v, and identifying separate virtual energy terms, equation (3) can be concisely written as Here, the intuitive understanding of the different virtual energy terms is most clearly understood in the case where w = v. In this case, K e denotes the current kinetic energy held within the fluid, A e indicates the rate at which kinetic energy enters, exits or grows in the region of interest Ω, S e indicates the power added or removed from the system due to shear, V e details the rate of viscous energy dissipation in the fluid and H(p) defines the hydraulic power.
The calculation in equation (4-9) varies from our previous IMRP estimatior 29 . The key difference resides in equation (6), where we previously applied integration by parts to reduce the numerical derivatives required on imaging data. This, however, lead to a more significant dependence on surface integrals, deteriorating method performance.
To isolate the relative pressure between two boundaries of the region of interest, we seek to assign attributes (or restrictions) on w. Using the desired regions to isolate the pressure drop, the boundary can be split into inlet Γ i , outlet Γ o , and wall Γ w regions (with Γ = Γ i ∪ Γ o ∪ Γ w ). Focusing on H(p), we observe that if the virtual field is a solenoidal field with zero velocity on Γ w , i o Under these assumptions, applying the Divergence theorem we observe the total inflow, Q, given by w is identical to its outflow, e.g.
is the w⋅n weighted mean pressure on the k th boundary. Further, the conditions on the virtual field reduce S e to more simplified forms In this work, we sought to minimise the added power due to S e . By ensuring that the virtual field acts in the surface normal direction, we focus on the added power due to flow gradients along the length of the artery of interest. As these variations are small, the added power becomes negligible.
Combining equations (5,6,8,(12)(13), we can express the relative pressure drop as e e e meaning that the relative pressure, Δp, is now expressed as a function of the virtual work-energy performed by w. Apart from acquiring the velocity field v, all that remains for deriving Δp is then to find a virtual field, w, abiding to the assumptions outlined above. In theory, any virtual field fulfilling these criteria can be selected and, in the absence of noise and sampling errors, will result in the same relative pressure estimate. In this study, we elect to construct w by solving Stokes problem, which ensures the existence of a unique solution satisfying the physical constraints outlined. Specifically, we solve the boundary value problem i w where λ is the virtual pressure field (Lagrange multiplier) corresponding to the virtual field w used to enforce equation (15), and n is the normal vector on the inlet plane Γ i . Information on the practical implementation of the above is given below under Data pre-processing and practical implementation of νWERP.
Data pre-processing and practical implementation of νWERP. Below we outline the pre-processing steps and implementation of νWERP for both in silico and in-vivo flow fields.
Flow field corrections. Using 4D flow MRI (as in the current study), eddy currents, field inhomogeneities, and potential regions of aliasing have to be corrected using customised pre-processing tools. For this specific study, 4D flow fields were reconstructed and corrected using GTFlow (GyroTools LLC).
Segmentation of fluid domain. The vascular region of interest required segmentation, separating it from the surrounding tissues and airways. A range of segmentation techniques exist in literature 44 . For this study, segmentation was performed by the following protocol: 1. An initial binary mask was created by collapsing the temporal dimension of the 4D flow velocity magnitude, leaving a spatial image of the average velocity magnitude of each pixel. A moving threshold was applied, removing voxels with velocity magnitude <25-40% of the image velocity encoding (or VENC), with the exact threshold defined by the user. 2. Isolated regions consisting of 10 voxels or less were removed, and remaining open enclosed regions within the remaining binary masks were filled. 3. Distinction of the vessel of interest was achieved by manual selection, separating the vascular segment from other disconnected regions. The vascular segment was isolated using a region growing approach within the binary mask generated from the previous step. 4. The remaining segmented velocity field was evaluated to identify potential erroneous voxels based on examination of: a. The temporal derivative of the velocity field. For this, temporal derivatives were evaluated throughout the entire mask, and ordered with respect to magnitude creating a velocity-rate histogram. Border voxels where subsequently removed from the mask when the velocity temporal derivative magnitudes were outside two standard deviations of the distribution mean.
Scientific RepoRts | (2019) 9:1375 | https://doi.org/10.1038/s41598-018-37714-0 b. The predominant flow direction in a given voxel. For this, the mean flow direction in a given voxel at peak systole was compared to the average flow direction in the neighbouring voxel area (within a spherical area of radius <5 voxels). If a voxel's mean flow direction deviated more than 120° from that of the neighbourhood, the voxel was deemed erroneous and removed. Subsequently, a second check was performed by evaluating the remaining masked voxels against a solved Stokes flow from inlet to outlet. Again, if a singular voxel's mean flow direction deviated more than 120° from that of the Stokes flow within a <5 voxel radius, it was deemed erroneous and removed. With this, a final binary mask was created, isolating the fluid domain of interest from static tissue or spurious noise in the acquired flow field.
Selection of inlet and outlet planes. Inlet and outlet planes within the binary mask were chosen by manually selecting an input/output point on the surface of the binary mask. The centre-line of the vessel was then automatically identified in the vicinity of the selected points. Using this, a spatial plane was then automatically identified, cutting the vascular cross-section perpendicular to the centre-line direction of the vessel, and including the selected input/output point on the mask surface.
Labelling the fluid domain. To facilitate for the definition of boundary conditions within the computation of the virtual field (see e) Computation of virtual fields below) a labelling of the vascular domain is performed, making use of the aforementioned domain, inlet, and outlet mask, respectively. Specifically, region growing technique was applied to create a mask nested within the original segmentation and bridging between the inlet and outlet planes. With this, labels are assigned to individual voxels to separate: interior voxels (entirely within the fluid), exterior voxels (entirely outside the fluid), inlet/outlet voxels, and wall voxels (separating interior and exterior domain). With this, appropriate boundary conditions can be set within e) Computation of virtual fields.
Computation of virtual fields. As described in Theoretical derivation of νWERP, the virtual velocity field w is chosen to satisfy a Stokes problem (given by equations (15)(16)(17)). In this work, the problem is set up and solved numerically using the Finite Difference Method (FDM) (see Supplement Material B). For this, an initial base-sampling of nodal points is generated at each voxel inside the segmented acquired velocity field. To improve on the accuracy of the solved virtual field, a subsampling of the data is performed, subdividing base-voxels by a chosen integer value, making a single image voxel consist of several uniformly distributed nodal points (within this work, all virtual fields were subsampled to an image resolution of 0.5 mm 3 ). Note that this subsampling only effects the resolution of w and all analysis of v is kept in the given image resolution of the performed acquisition.
With nodal points defined, boundary conditions of the FDM were set such that nodal points on the binary inlet and wall masks were set based on the velocity constraint in equation (17). The virtual field w can be computed by solving a linear system of equations generated by the FDM. Further details of the method and solver can be found in Supplementary Material. An overview of the pre-processing from segmentation to computed virtual field is also provided in Supplementary Fig. 1.
Noise filtering and data smoothing. To reduce the effects of spurious noise in the acquired flow field, spatial noise filtering is applied to the data prior to νWERP estimation. Specifically, a k-order 3D polynomial fitting is used over local kernels to approximate a noise-free signal (see Supplementary Fig. 2). The technique applied uses a standard Savitzky-Golay filtering 45 commonly applied in image processing. For a given flow field, optimal interpolation order and kernel size was determined by evaluating the filtering on a 1D sinusoidal signal, with discretisation and maximum spatial gradient identical to that of the input data. Truncated Gaussian noise resembling the SNR of the input data is added to the 1D signal, and noise filtering using polynomial order from 0-3, and kernel sizes from 1 3 to 9 3 were evaluated. The combination of polynomial order and kernel size minimising the error to the known 1D sinusoidal function is then chosen as the optimal filter parameters for the given flow field. This choice of optimal filtering parameters was performed once for an entire 4D dataset.
νWERP discretisation. The quantification of relative pressure from equation (14) for a spatiotemporally discrete flow field is defined as: being the flow and virtual field at a time half-way between t and t + 1. Note that computations are performed at the midpoint between time steps to improve temporal accuracy.
Surface and volume integral estimation. Energy components of νWERP in equation (18)  Generation of patient-specific in-silico image data. For evaluating vWERP, patient-specific in-silico data was generated from computational fluid dynamics simulations of a coarcted and an acute type B dissected aorta, respectively. Both models have been presented in detail elsewhere 34,35 , but for completeness a summary of the data generation process is outlined below.
Coarcted aorta. With anatomical data originating from 3D contrast enhanced MR angiography, a model of the aortic arch, including connecting carotid bifurcations and details of the pathological coarctation, was created using a combination of SimVascular (simtk.org) and MeshSim (Simmetrix, Clifton Park, NY). Using acquired 2D PC-MRI data, appropriate boundary conditions were set to simulate flow through the coarcted aorta model. Simulations were run in SimVascular using a Coupled-Momentum Method 46 developed for solving the incompressible Navier-Stokes equations in deformable arteries. For the νWERP analysis, the simulated dataset was sampled onto a uniform grid of 2 mm 3 voxels, with 22 time slices representing one cardiac cycle.
Dissected aorta. Using a combination of thoracic computed tomography angiography (CTA) and 2D PC-MRI, a patient-specific simulation model of an acute type B dissected aorta was segmented, with the model including vascular anatomy spanning from the carotid to the femoral arteries. The model was built using CRIMSON 47 by defining well-spaced transmural contours along each vessel in the acquired image data. A volume segmentation was produced by lofting the contour sets using non-uniform rational B-splines, and individual vessels were joined to create a full CAD segmentation, including connecting the identified segments 48 . Each lumen was created separately and the final model was visually checked to verify consistency to the image data. 17 connecting tears were introduced between true and false lumen. With that, the incompressible Navier-Stokes equations were solved using a stabilised finite-element formulation. The acquired 2D PC-MRI at the aortic root was used as inlet boundary condition. The distal vascular bed was modelled using 0D Windkessel RCR elements coupled to the 3D domain using a multi-domain method 49 and attached at all 17 outflow arteries. Parameters for each RCR were tuned to match 2D PC-MRI acquired in a number of positions in the descending aorta. The meshes were created with tetrahedral elements (global element size: 1 mm, local refinements at vascular wall between 0.1-0.05 mm), with iterative mesh adaptation performed at regions experiencing high velocity Hessians. Simulations were run for multiple cycles until a periodic solution was found for both flow and pressure. For the νWERP analysis, and particularly the spatiotemporal convergence analysis, the simulated dataset was sampled onto a uniform grid of varying spatiotemporal resolution. In silico images were generated for 1, 2, 3, and 4 mm 3 voxel size, with temporal sampling of between 8 and 64 time slices representing one cardiac cycle. Image noise. For both in-silico sets, image noise was applied frame-by-frame using the maximum velocity, v max , of the entire dataset to represent the velocity encoding of a 4D flow MRI acquisition. With such, the relationship between the standard deviation of the velocity field, σ v , and the SNR can be expressed as 50 : With that a corresponding σ v could be derived for the low-and high-noise configurations, respectively. Knowing this, noise was applied distributed over all voxels of the in silico flow images, following a truncated Gaussian distribution 51 , with data truncation applied above two standard deviations to avoid spurious noise outliers.
In-vivo catheter and 4D flow MRI data acquisition. For the validation of νWERP, data was obtained during a combined cardiac MRI and cardiac catheterisation procedure (XMR). Both cardiac MRI imaging and invasive pressure assessments were performed under the same general anaesthesia within a timeframe of 30 minutes to one hour. In all cases heart rate and blood pressures were constantly monitored and remained stable throughout the procedure. Invasive pressures were obtained using pullback procedures from the ventricle into the ascending aorta, proximal descending aorta, diaphragmatic aorta and the femoral artery in steps using a fluid filled catheter (4 F GL braided catheter, Terumo Medical Corp., Somerset, NJ). 10 second pressure readings were obtained at each site, with multiple beats acquired throughout the respiratory cycle. Using the simultaneously acquired ECG signals, pressure measurements at different anatomical positions could be co-registered (performed by forcing signal overlap at corresponding peak R-wave positions). At each catheter position, multiple curves were acquired throughout the respiratory cycle. Instantaneous relative pressure differences were then computed by subtracting pressures in the proximal measurement point with pressures in the more distal measurement point(s). This process was repeated to generate a family of difference curves from which mean and standard deviations were extracted to represent the entire set of catheter measurements. The position of the catheter was visually confirmed by fluoroscopy imaging, used to guide the positioning of the catheter during in-vivo measurements, and manually identified in the corresponding segmented flow field from the 4D flow MRI.
The 4D flow MRI data was acquired using a Philips Achieva system (1.5 T, average spatial resolution: 2.1 mm 3 , temporal resolution: ~31 ms), with complete flow fields reconstructed using GTFlow (GyroTools Ltd, Zurich, Switzerland), in which eddy current, field inhomogeneities, and potential aliasing were corrected for.
For the temporal alignment of 4D flow MRI data and cardiac catheterisation, this was achieved by forced overlap at peak systole, with potential cardiac phase extensions unwrapped to generate data within one given cardiac cycle only. Apart from this, the temporal sampling rate of the two measurements were kept separate, given by the used measurement equipment, respectively.
Note that all included patients underwent a clinically indicated XMR for comprehensive assessment of cardiac and vascular function for unexplained exercise intolerance in the face of a complex congenital heart disease. The subject characteristics of the 6 patients included in this study are summarised in Supplementary Table 1. As stated in the Results section, all patients participated under informed consent, with study approval granted by the Regional Ethics Committee, South East London, United Kingdom (REC, 10/H0802/65). Additionally, all methods and experiments were performed in accordance with relevant guidelines and regulations. Data analysis. The mean error, ε Δp , for the temporal νWERP estimations comparing either against available in-silico or in-vivo data is defined as: where Δp W | n is the νWERP estimate at discrete time step t n and Δp| n is the corresponding true data closest to the νWERP estimate within t n−1 < t < t n+1 , with closeness evaluated by minimising the Euclidean distance between the temporal curve of the true and νWERP pressure drop, respectively. Note that for temporal estimations over an entire cardiac cycle, T, the above is evaluated at discrete temporal acquisitions t 1 to t N where N is the number of temporal sample points. The root mean square error, ε rms , is calculated by: Using the maximum relative pressure peak, Δp max , and corresponding relative pressure peak estimated by νWERP, Δp W,max , as input. Again the νWERP estimate was identified within t n−1 < t < t n+1 , with closeness evaluated by minimising the Euclidean distance between the temporal curve of the true and νWERP pressure drop, respectively. Alternative methods for non-invasive relative pressure estimation. In this work, νWERP performance is compared to alternative methods proposed for non-invasive relative pressure estimation. Details for these methods are given in other references 24,25,32,37,52,53 . Below is a short summary of the used methods.
Bernoulli. The Bernoulli equation 52,53 , commonly used clinically as part of recommended guidelines 1 , assumes that the relative pressure, Δp, is governed solely by the advective component of the Navier-Stokes equations. By so, Δp can be derived as: With v 2 and v 1 the velocity at a single point of the inlet and outlet plane, respectively, and ρ the fluid viscosity.
Unsteady Bernoulli. Extending from the above given Bernoulli equation, the unsteady Bernoulli 37 incorporates contributions due to inertial accelerations, with the pressure drop defined as The influence of the inertial acceleration is then added by evaluating the temporal derivatives of the flow along the trace, S, of an infinitesimal particle seeded into the flow. For simplicity, this streamline is often assumed to follow the centreline of the investigated vessel.
Work Energy Relative Pressure (WERP). The Work Energy Relative Pressure (WERP) 32 is the predecessor to the proposed νWERP. The theory of WERP is identical to that provided in Theoretical derivation of νWERP. However, a key differences lies in that no virtual field is introduced, but rather the work-energy of the acquired flow field v is evaluated. By so, Δp is defined as: Compared to νWERP, all energy terms above now relate to v. Identically, Q is the flow of v over the defined inlet/outlet planes, with divergent behaviour seen when Q → 0.

Data Availability
All data generated or analysed during this study are included in this published article (and its Supplementary Information files).