Blood flow imaging by optimal matching of computational fluid dynamics to 4D‐flow data

Three‐dimensional, time‐resolved blood flow measurement (4D‐flow) is a powerful research and clinical tool, but improved resolution and scan times are needed. Therefore, this study aims to (1) present a postprocessing framework for optimization‐driven simulation‐based flow imaging, called 4D‐flow High‐resolution Imaging with a priori Knowledge Incorporating the Navier‐Stokes equations and the discontinuous Galerkin method (4D‐flow HIKING), (2) investigate the framework in synthetic tests, (3) perform phantom validation using laser particle imaging velocimetry, and (4) demonstrate the use of the framework in vivo.


| INTRODUCTION
Three-dimensional, time-resolved, three-directional MR flow imaging (4D-flow) is a powerful method for investigation of cardiovascular physiology 1,2 and pathophysiology. [3][4][5][6][7][8] Current 4D-flow acquisition schemes provide comprehensive information on hemodynamics, but the tradeoff between image quality and scan time still limits widespread use in patient studies and the clinic.
Recent advances in compressed sensing 9 has enabled 4D-flow with reduced scan times and improved image quality [10][11][12][13] by introducing a priori information (ie, knowledge that the image is sparse in some transformed domain). Recent studies have explored using the Navier-Stokes equations as a source of a priori information for blood flow imaging by merging 4D-flow and computational fluid dynamics (CFD). [14][15][16][17][18][19][20][21][22][23] The Navier-Stokes equations give a complete physical description of fluid flow and provides powerful a priori information, which may be used to increase resolution or reduce scan times.
Current methods for matching CFD to 4D-flow typically consist of three parts: an accurate CFD simulation with free parameters to optimize, a metric to measure the difference between CFD and 4D-flow, and an efficient strategy to update the parameters to minimize the difference metric. Most studies to date use CFD simulations with low order of accuracy, such as particle methods, finite differences, or low-order finiteelement methods. [14][15][21][22][23][24] Furthermore, the difference metric is usually not designed to model the 4D-flow measurement process. [14][15][16][20][21][22][23] Finally, efficient optimization can be limited by a lack of efficient computation of the gradient of the difference metric. 20,22 Funke et al 16 used an adjoint CFD solver to compute gradients for efficient optimization. However, they used a low-order CFD method and a difference metric that does not take the 4D-flow measurement process into account.
We have developed a new optimization-driven simulationbased framework for efficient matching of high-order accurate CFD to 4D-flow data, called the 4D-flow High-resolution Imaging with a priori Knowledge Incorporating the Navier-Stokes equations and the discontinuous Galerkin method (4D-flow HIKING). This framework uses a high-order CFD solver, a difference metric modeling the 4D-flow measurement process, and an efficient optimization strategy based on gradient computations through the adjoint CFD equations. The HIKING framework therefore uses the velocity field from acquired 4D-flow images in a segmented anatomic structure as input to an iterative CFD optimization, with the aim to improve temporal and spatial resolution with higher accuracy and precision. This can be applied to 4D-flow data from any sequence and scanner, as it takes place after MRI reconstruction.
Therefore, this study aims to (1) present the 4D-flow HIKING framework, (2) investigate the performance of the framework in synthetic test cases, (3) validate the framework using laser particle imaging velocimetry (PIV) in a phantom setup, and (4) demonstrate the use of the framework in vivo.

| Overview
The 4D-flow HIKING framework aims to create a detailed blood-flow computational fluid dynamics (CFD) simulation that matches the patient's blood flow. The simulation contains free parameters that prescribe the time-dependent flow at the boundaries. In earlier CFD approaches, information for this is typically taken from 2D through-plane flow scans, 25,26 which means that real-world information enters the simulation model only at the inflow boundaries. Outflow conditions are typically specified by a flow split between vessels 27,28 or a Windkessel model. 29 This potentially leads to a loss of simulation fidelity away from the inflow boundaries.
To improve the simulation, 4D-flow HIKING includes data from a full volumetric 4D-flow scan to optimally and automatically determine boundary conditions. This is done by minimizing the difference between the CFD simulation and the 4D-flow data, as shown in Figure 1. Because the simulation data are defined on a high-resolution mesh and the 4D-flow data are acquired on a low-resolution grid, the comparison cannot be made directly. Therefore, the CFD data are subjected to a measurement model that is designed to imitate the spatial and temporal averaging in the 4D-flow measurement.
In detail, the 4D-flow HIKING framework uses the 4D-flow velocity data v MR and a forward model M (described in detail in the following subsection) that takes the spatial and temporal smoothing of the MRI measurement into account, to recover the high-resolution CFD data v CFD . The CFD solution is uniquely determined by the geometry and a parameter vector µ, which defines flow at inlets and outlets. Recovery of the high-resolution CFD data v CFD is then formulated as the following optimization problem ( Figure 1G): The obtained CFD solution is optimal in the sense that it is the closest CFD approximation to the 4D-flow data, while also considering the 4D-flow measurement process in the forward model M. This formulation does not include tunable parameters to weigh constraints and data, which simplifies the optimization.
A flowchart of the algorithm to minimize Equation 1 is given in Figure 2. The process is governed by the Interior Point OPTimizer software. 30 The optimizer (Figure 2A) is initialized with an estimate of the parameter vector µ, and then starts the forward CFD solver ( Figure 2B) to compute a CFD solution v CFD . The CFD solution is combined with the 4D-flow data ( Figure 2C) to compute the value of the objective function I( ). The adjoint CFD solver is then used to compute the gradient dI( )∕d ( Figure 2B). The gradient information is used to update µ in the direction of the steepest descent of the objective function, thereby improving the fit of the CFD to the 4D-flow data. This process is repeated until convergence is attained.

| Construction of the forward model M
The forward model M projects the CFD solution to the geometry of the MRI data, while also considering spatial and temporal smoothing in 4D-flow. As previously shown, 18 spatial smoothing can be modeled by the Cartesian MRI point spread function. The temporal smoothing was modeled as a smoothed averaging in time, with the window width equal to the temporal resolution of the 4D-flow scan.
In detail, let the Cartesian physical coordinates (x, y, z) be defined by the top-left pixel in the first slice of the 4D-flow data set, and its readout, phase-encoding, and slice-encoding directions. The MRI data are indexed using voxel numbers I = 1 … N x , j = 1 … N y , and k = 1 … N z , and the timeframe number n = 1 … N t . Each voxel has a center coordinate (x ijk , y ijk , z ijk ) with the voxel size (Δx, Δy, Δz), and the temporal resolution of the MRI data is Δt.

| Temporal smoothing
Temporal averaging for each time step was modeled using a smoothed "box" function χ centered at position t 0 with the width w and smoothing parameter σ (Supporting Information Figure S1A), as follows: The factor α is used to normalize the area under the curve to 1. The temporal smoothing function w n t for each 4D-flow time step n is then defined as The parameter Δt was equal to the temporal resolution of the 4D-flow data. The smoothed shape of the box function was used to avoid discontinuities and improve numerical accuracy of integrals. The parameter σ t was set to 4 ms.

| Spatial smoothing
For Cartesian 4D-flow, 18 the spatial smoothing W ijk xyz for each voxel (i, j, k) is given in terms of the function sinc (x) = sin ( x) ∕ ( x), truncated by the smooth-box function χ as follows (Supporting Information Figure S1B): The truncation smoothing parameter σ xyz was set to σ xyz = 0.1 × min(Δx, Δy, Δz) (Supporting Information Figure S1B,C).
With these definitions, the model M for the voxel (i, j, k) at the time step n was computed by multiplying the CFD velocity field v CFD by the smoothing functions w n t and W ijk xyz and integrating as follows: The outer integral (T) was taken over the whole time interval of the 4D-flow data, and the inner integral (Ω) over the whole FOV. The integral was computed for each MRI voxel and for every time step.

| Computational fluid dynamics implementation
Blood flow was modeled using the compressible Navier-Stokes equations with the isentropic assumption (constant entropy; flow is internally adiabatic and reversible), which provably approaches the incompressible Navier-Stokes equations in the low-velocity limit. 31 In our simulations, flows were nearly incompressible (less than 1% density variation); therefore, the isentropic assumption accurately models incompressible flow without explicitly dealing with the incompressibility constraint, which would add significant analytical and computational complexity in the method. Details on the mathematical formulation and discretization are presented in the Supporting Information.
The method was implemented in the in-house 3DG software, 32 a high-order multiphysics solver based on an efficient parallel implementation of the Discontinuous Galerkin (DG) method that allows for unstructured meshes of tetrahedral and hexahedral elements, implements a range of high-order implicit time-integrators, and uses automated symbolic differentiation to evaluate Jacobians and sensitivities.
To evaluate the objective function in Equation 1, the integral in Equation 5 is computed in a solver-consistent manner. 33 Specifically, the basis functions and quadrature rule used in the DG method were used to compute the spatial integral. The temporal integral was converted into an ordinary differential equation and integrated with the same diagonally implicit Runge-Kutta (DIRK) scheme used in the forward solver, which ensures that the truncation error of the objective function and governing equations match exactly. The gradient of the objective function was computed using a previously described adjoint method. 33 The adjoint equations of the fully discrete CFD method are solved (linearized equations solved backward in time), and the adjoint variables are used to compute the gradient dI( )∕d of the objective function I( ). Therefore, the cost to evaluate the objective gradient is approximately that of one standard (primal, nonlinear) CFD simulation and one adjoint linearized solution. In practice, the computational cost of the linearized solution may be somewhat larger than the primal simulation due to practical considerations such as file input and output (I/O) and evaluating partial derivatives. The adjoint method and solver-consistent approximation of the objective function were implemented in C++ within the in-house 3DG framework. 33 The optimization problem (Equation 1) was solved using Interior Point OPTimizer. 30 The solver uses a limitedmemory quasi-Newton approximation of the objective function Hessian to compute a search direction to update the parameter vector μ. A line search is used to determine step sizes, to sufficiently reduce the objective function. Therefore, each adjoint solution (gradient evaluation) was followed by multiple primal solutions to determine an appropriate step size. Interior Point OPTimizer uses the Karush-Kuhn-Tucker conditions to end the optimization; however, it turns out to be more useful for the present application to stop when the relative change ε in μ is smaller than a threshold (taken here as < 0.01) for three consecutive iterations. For the in vivo experiment, a relative reduction of the objective function I( ) of less than 0.01 for two consecutive iterations was used.
Synthetic test cases were run using 16 threads on a computation server with 256 GB of memory and 20 cores (dual Intel Xeon E5-2680). Flow phantom simulations were run on the Edison system at the National Energy Research Scientific Computation Center using 10 nodes, with two 12-core Intel "Ivy Bridge" processors (2.4 GHz) and 64 GB memory per node, for a total of 240 cores. The in vivo computations were performed using 24 cores on a 36-core server (dual Intel Xeon Gold 6154) with 384 GB memory.

| Synthetic test cases
To investigate the performance of the method in a controlled setting with a known reference solution (synthetic test cases), a 2D synthetic numerical flow was constructed. A rectangular domain ( Figure 3A) with dimensions 180 × 175 mm 2 was defined, with a rectangular inlet of dimensions 37.5 × 25 mm 2 . The mesh used rectangular elements of order P = 3 with resolution h = 6.25 mm. A time window of 2.0 seconds was used with 80 time steps.
The inflow boundary condition was defined by polynomial interpolation on three nodes over the inlet (Supporting Information Figure S2A), with parameters µ 1 , µ 2 , and µ 3 . The temporal parametrization was a Gaussian with center µ 4 and SD µ 5 (Supporting Information Figure S2B).
A reference solution was computed using the forward CFD solver, with a ground-truth parameter vector µ* = (26.25 cm/s, 35 cm/s, 26.25 cm/s, 0.5 s, 0.2 s). Kinematic viscosity was set to ν = 3.1·10 −6 m 2 /s. The Reynolds number based on the inflow dimensions (25 mm), and peak inflow velocity (35 cm/s) was Re ≈ 2900. Three numeric MRI data sets were constructed with reduced resolution and Gaussian noise added to the velocity (Table 1): (1) 4D-flow with spatial and temporal resolution following a recent consensus statement 34 and previously published noise levels, 35 (2) reduced spatial and temporal resolution (5 mm and 100 ms, respectively), and (3) higher noise levels. The 4D-flow HIKING framework was then used to optimize the flow field by optimizing Equation 1 with respect to the parameter vector µ = (µ 1 , µ 2 , µ 3, µ 4, µ 5 ). The initial guess for the parameter vector was µ 0 = (17.5 cm/s, 17.5 cm/s, 17.5 cm/s, 0 s, 1 s). The final flow estimation error was computed as the RMS (averaged L 2 norm) of the difference between the optimized velocity field and the initial reference, with the mean taken over space and time.

| Physical flow phantom validation
The flow phantom is shown in Figure 4 (previously described in detail 36 ). In the main tank, a nozzle with diameter = 25 mm is used to shape the pulsatile inflow from a custombuilt pump into a vortex ring, a 3D and pulsatile yet stable, well-defined, and reproducible flow field. 37,38 Vortex rings appear naturally in the left ventricle of the human heart and are physiologically relevant in diastolic dysfunction and heart failure, 2,39 and are therefore appropriate as a test case for the phantom validation. Figure 4E,F shows the computational mesh used, containing 50 465 tetrahedra ("50k mesh"). The maximum element size in the part of the mesh with significant flow was 6 mm. High-order curved elements (P = 3, 4th-order accurate) were used to accurately represent the circular shape of the inflow nozzle ( Figure 4F). A refined mesh (246 996 tetrahedra, "250k mesh") was used to generate a high-resolution flow field with the final optimized parameter set for visualization purposes (typical element size = 2.8 mm). Geometry and meshing was performed using the software Gmsh 3.0.3 40 with the constructive solid geometry module in the OpenCASCADE module (Open Cascade, Guyancourt, France).
Phantom 4D-flow data were acquired on a 1.5T Philips Achieva scanner (Philips Healthcare, Best, The Netherlands). Spatial resolution was 3 × 3 × 3 mm 3 and temporal resolution 50.4 ms. Sequence parameters are given in Table 1. To demonstrate the potential of 4D-flow HIKING to compute high-resolution flow from low-resolution data, the acquired 4D-flow data were downsampled by a factor of 2 in space and time. The resulting input data to the framework had 6 × 6 × 6 mm 3 voxels and a temporal resolution of 100 ms.
To provide a reference flow field, PIV was performed using a LaVision (Göttingen, Germany) system as previously described. 36 The system consisted of a 10-Hz dual-pulse  The 4D-flow HIKING framework was applied with the downsampled 4D-flow velocity images as the input data. A cropped slice (20 × 13 voxels) in the symmetry plane of the nozzle was extracted from the full 4D-flow volume to accelerate the computation. All three flow directions and all timeframes of the downsampled 4D-flow data were used. The inflow in the circular nozzle was described as a plug flow (constant over the nozzle inlet) with a generalized Gaussian profile over time, resulting in four free parameters that were optimized: the peak velocity of the flow (µ 1 ), the position of the Gaussian (µ 2 ), the width of the Gaussian (µ 3 ), and the shape parameter β (µ 4 , Supporting Information Figure S2B).
The parameters were initialized to match the through-plane flow in the phantom nozzle: µ 0 = 32.7 cm/s, 0.22 s, 0.125 s, 3.47. Kinematic viscosity was set to ν = 1.0·10 −6 m 2 /s, consistent with water at approximately 20°C. Mesh independence was investigated by running the simulation on both the 50k mesh and 250k mesh.

| In vivo proof of concept
One healthy subject was scanned using a 7T MRI system (Philips Achieva 7T, Philips Healthcare, Best, The Netherlands) at the National 7T Facility in Lund, Sweden. The study was approved by the local ethical review board in Lund, Sweden. A 4D-flow scan was performed in a transversaloblique orientation to cover the distal left internal carotid artery and proximal middle cerebral artery. A time-resolved gradientecho angiographic image was acquired as a high-resolution anatomical reference, to be used for segmentation and meshing. Furthermore, a reference 2D flow measurement was acquired in the first segment (M1) of the left middle cerebral artery ( Figure 1B). Sequence parameters are given in Table 2. The left internal carotid artery, middle cerebral artery, anterior cerebral artery, and three major branches were segmented from the angiographic image stack using ITK-SNAP. 41 The segmentation was then imported into ICEM CFD v19.2 (ANSYS, Canonsburg, PA) to generate tetrahedral meshes ( Figure 5A-C). Mesh independence was investigated by running simulations on two grids: one with 11 645 elements and 150 time steps (11k mesh), and one with 18 856 elements and 300 time steps (18k mesh). Voxels located within the CFD mesh were used for the optimization. Voxels with a 3D velocity higher than 1 m/s were excluded to limit the influence of outliers on the optimization.  The mesh contained six boundaries with flow: one inlet (internal carotid artery) and five outlets ( Figure 5B,C). For each flow boundary, the flow was spatially constant (plug flow), and the temporal parametrization consisted of six elements with a parabolic curve shape in each element, giving 12 parameters per boundary ( Figure 5D) (ie, 72 parameters to be optimized in total). Velocity at vessel walls was set to zero (no-slip condition). Parameters were initialized from the 4D-flow data at the midpoint of each boundary. Kinematic viscosity was set to ν = 3.5·10 −6 m 2 /s. Convergence was assumed when the objective function I ( ) decreased less than 1% in three consecutive iterations.

| Statistical methods
Agreement was quantified using linear regression, R 2 values, and modified Bland-Altman analysis with the laser PIV data on the x-axis. Bias between laser PIV data and 4D-flow HIKING optimized flow is presented as mean ± SD for all analyzed voxels. For the phantom validation, differences in mean and SD of the error with respect to the PIV reference were tested using an unpaired t-test (for mean error) and an F-test (for differences in error SD). Mesh independence was analyzed by comparing point-wise velocities using Bland-Altman plots. For statistical tests, P-values less than 0.05 were considered statistically significant. Figure 3D-F shows 4D-flow HIKING optimizations for all three 2D synthetic test cases. Qualitatively, the optimized flow fields are similar to the reference CFD solution. Quantitative accuracy was high, with errors of 0.22, 0.14, and 0.24 cm/s, respectively (corresponding to 0.6%, 0.4%, and 0.7% of the peak inflow velocity). The optimization of test case 1 required 81 minutes of wall-clock computation time (22 core hours). Case 2 required 64 minutes (17 core hours), and case 3 used 76 minutes total (20 core hours). Figures 6 and 7 show phantom results, with better agreement with laser PIV for 4D-flow HIKING compared with the downsampled 4D-flow data used as input to the framework, both in the horizontal (mean difference −0.05 vs − 1.11 cm/s, P < .001; SD 1.86 vs 4.26 cm/s, P < .001) and vertical directions (mean 0.05 vs −0.04 cm/s, P = .29; SD 1.36 vs 3.95 cm/s, P < .001). Additionally, 4D-flow HIKING results showed better agreement with laser PIV than the original non-downsampled 4D-flow data in the horizontal (mean −0.05 vs −0.36 cm/s, P < .001; SD 1.86 vs 2.25 cm/s, P < .001) and vertical directions (mean 0.05 vs 0.02 cm/s, P = .56; SD 1.36 vs 2.06 cm/s, P < .001). Figure 7B,C shows convergence of the phantom validation parameters. The optimization converged after 7 hours and 19 minutes of wall-clock time (1756 core hours). The final parameter vector was µ* = 37.3 cm/s, 239 ms, 110 ms, 3.46. Results on the 50k and 250k meshes showed strong agreement. The error was 0.02 ± 0.93 cm/s and −0.04 ± 0.57 cm/s in the horizontal and vertical directions, respectively (Supporting Information Figures S3 and S4).  Figure 8B). After 4D-flow HIKING optimization, the flow curves show better agreement (flow rate error 3.5%, Figure 8C). Computation time was 5 hours and 45 minutes of wall-clock time (138 core hours).

| In vivo proof of concept
Strong agreement was found between simulations on the two meshes (Supporting Information Figures S5 and S6), with low errors in all three flow directions (right-left: 0.02 ± 2.19 cm/s, anterior-posterior: 0.07 ± 1.05 cm/s, feet-head: 0.13 ± 1.26 cm/s).

| DISCUSSION
This paper describes a new framework for optimizationdriven simulation-based matching of high-resolution CFD flow data to 4D-flow by incorporating the Navier-Stokes equations as a priori information, called 4D-flow HIKING. Synthetic test cases show the feasibility of the proposed formulation to recover high-resolution flow from noisy data with low resolution. Phantom validation using PIV shows excellent accuracy. Furthermore, data from a healthy volunteer show that the method is feasible in vivo.

| Relation to previous studies
Previous studies 42,43 have used flow MRI data to inform CFD simulations. However, the flow MRI data typically enters only as boundary conditions at the inlets, which means that the full potential of 4D-flow is not used. Additionally, because only the inflow boundary conditions are constrained, accuracy of the computed flow further downstream in the flow is not enforced. In contrast, our framework enables the use of a full volumetric 4D-flow velocity scan to optimize both inflow and outflow conditions to fit the CFD flow field in the whole region of interest. In this way, CFD modeling assumptions are automatically and optimally guided by the full 4D-flow data set, thereby leading to more realistic CFD data.
Rispoli et al 18 previously demonstrated high-resolution flow optimization in the carotid bifurcation using an MRI signal model and the SIMPLER 44 CFD solver on a Cartesian grid. In contrast, our work uses a fully unstructured grid, potentially enabling local refinement in regions of interest (eg, near walls). Furthermore, they use a weighting approach between the 4D-flow data and CFD solver, which leads to solutions that approximately fulfill the discretized Navier-Stokes equations. In this way, the relative influence of the MRI data or the CFD solver can be adjusted. However, adding such a weight factor to produce approximate CFD solutions is not feasible in our framework due to analytical and computational reasons. Instead, our framework searches for the exact CFD solution with the boundary conditions that give the best agreement with the existing 4D-flow data.
Funke et al 16 applied a method similar to 4D-flow HIKING to the human aorta, including adjoint-based gradient computation. However, they did not use a forward model of the MRI measurement process, opting instead to linearly interpolate the 4D-flow data onto the computational mesh nodes. Furthermore, they used a low-order CFD solver, which limits numerical efficiency. Guerra et al 17 approached the CFD data-matching problem using optimal control theory, which after discretization led to a series of quadratic optimization problems, in which approximate gradients are computed in each step. In contrast, 4D-flow HIKING computes exact gradients of the CFD solution. Goenezen et al 20 used velocity data from optical coherence tomography to perform CFD optimization of blood flow in chick embryo hearts. They used a heuristic algorithm to compute the inlet pressure level that resulted in a best-fit velocity field to the Doppler data. A similar strategy was used by Mohd Adib et al 22 by adjusting the pressure levels of outlets until a good fit to the 4D-flow data was found. The main limitation with this approach was that the number of parameter combinations to test grows quickly with the number of parameters. In contrast, the adjoint-based formulation used here can efficiently find the optimal solution using gradient descent, even for many parameters as shown in our in vivo proof of concept.
Furthermore, several groups used a CFD solver to compute a library of flow fields for a given geometry, which are then combined in a filtering process to create the optimized flow. 19,21,23 This constitutes a different approach compared with purely optimization-based methods such as 4D-flow HIKING. Bakhshinejad et al 21 showed good performance in denoising simulated 4D-flow data in synthetic tests, and Gaidzik et al 23 showed good accuracy in a physical phantom setup with PIV as the reference standard. Further studies are needed to elucidate the relative strengths and weaknesses of these filtering methods compared with 4D-flow HIKING.

| Performance and accuracy
The synthetic test cases were designed to investigate the ability of the 4D-flow HIKING framework to recover a known flow field from MRI-like data. Good agreement with the ground truth was found in all three cases. This suggests that 4D-flow HIKING can not only be used to increase spatial and temporal resolution in a given 4D-flow data set, but also to accelerate scans by acquiring images at lower spatial and temporal resolutions, or by using acceleration techniques that result in a higher noise level.
The laser PIV phantom validation showed good agreement for the optimized CFD data, indicating that the forward model with spatial and temporal smoothing used in the optimization is appropriate for Cartesian 4D flow. Future studies may further investigate an improved forward model for additional performance and accuracy. Interestingly, the optimized 4D-flow HIKING data, which used only the downsampled 4D-flow data as input, had a better agreement with the laser PIV than the original high-resolution 4D-flow data. This highlights the potential of simulation-based imaging to both accelerate existing 4D-flow acquisitions and to improve accuracy.
For in vivo data, 4D-flow HIKING optimization improved agreement between the CFD model and the 2D-flow reference data. This shows the potential of 4D-flow HIKING to improve in vivo imaging, even in a setting with multiple inlets and outlets in arbitrary geometries, and many parameters to optimize. Because our method merges 4D flow and CFD, which are methods with separate error sources, limitations and assumptions, more comprehensive validation studies are needed to investigate the performance of the framework in vivo.
Previous studies have shown that in vivo blood-flow simulations can be sensitive to in-plane (secondary) flow, helicity, F I G U R E 7 Accuracy improvement using 4D-flow HIKING in phantom data (A) and plots of parameter convergence (B,C). A, Improvement in accuracy of the 4D-flow HIKING optimization compared with the 4D-flow data. B,C, The change in optimization parameters over the 12 iterations. ***P < .001 for mean bias; † † † P < .001 for SD of bias and spatial variations in the inflow boundary conditions. [45][46][47][48] As our present framework uses 1-dimensional "plug" flow at the boundaries, adding these features may improve accuracy in future in vivo applications. Using our adjoint-based optimization framework, these features of the inflow profile can be automatically optimized to fit the available 4D-flow data. The fitting process can further benefit from accurate initialization using optimized mappings from 4D-flow data to the CFD mesh, which conserves the flow rate and velocity field characteristics. 49 Furthermore, Bozzi et al used Monte Carlo methods to evaluate the sensitivity of aortic CFD to boundary conditions. They found that boundary condition uncertainties can be amplified in the CFD solution and influence both velocity and wall shear stress data. 50 These effects can be further investigated using accelerated uncertainty quantification methods. 51

| Future use of the 4D-flow HIKING framework
Optimization-driven, simulation-based flow imaging, using 4D-flow HIKING, generally has several potential applications: (1) to increase spatial and temporal resolution, in which the Navier-Stokes equations are used to compute blood flow at a higher resolution than possible to acquire with 4D flow; (2) to accelerate scans, in which a low-resolution 4D-flow data set is acquired in a fast scan and simulation-based imaging can then be used to recover an intermediate-resolution flow; and (3) to compute quantitative parameters not apparent in the 4D-flow data set (eg, wall shear stress in the aorta, carotid arteries, intracranial arteries) to predict aneurysm formation and rupture 52 and stenotic pressure loss, such as in coarctation of the aorta. 53

| Limitations
This study does not include adaptations of the mesh near walls to capture sharp velocity gradients in the boundary layers. Such mesh modifications may improve accuracy. In vivo blood flow was approximated as Newtonian, which has previously been shown to be appropriate. 54 Furthermore, vessel walls were assumed to be rigid, and making them compliant may improve simulation fidelity.
Although the adjoint-based gradient computations enable efficient optimization, 4D-flow HIKING is currently limited by long computation times (several hours for our in vivo proof of concept). This may be alleviated by an optimized implementation of the adjoint CFD, or using reduced-order models that provably converge to the full optimization problem. 51

| CONCLUSIONS
The presented framework for optimization-driven, simulation-based flow imaging, called 4D-flow HIKING, shows good performance in both synthetic test cases and phantom validation using laser PIV. In vivo results show promise for future applications with higher resolution, shorter scan times, and accurate quantification of physiological parameters.