Computational fluid dynamics simulations of blood flow regularized by 3D phase contrast MRI

Phase contrast magnetic resonance imaging (PC-MRI) is used clinically for quantitative assessment of cardiovascular flow and function, as it is capable of providing directly-measured 3D velocity maps. Alternatively, vascular flow can be estimated from model-based computation fluid dynamics (CFD) calculations. CFD provides arbitrarily high resolution, but its accuracy hinges on model assumptions, while velocity fields measured with PC-MRI generally do not satisfy the equations of fluid dynamics, provide limited resolution, and suffer from partial volume effects. The purpose of this study is to develop a proof-of-concept numerical procedure for constructing a simulated flow field that is influenced by both direct PC-MRI measurements and a fluid physics model, thereby taking advantage of both the accuracy of PC-MRI and the high spatial resolution of CFD. The use of the proposed approach in regularizing 3D flow fields is evaluated. The proposed algorithm incorporates both a Newtonian fluid physics model and a linear PC-MRI signal model. The model equations are solved numerically using a modified CFD algorithm. The numerical solution corresponds to the optimal solution of a generalized Tikhonov regularization, which provides a flow field that satisfies the flow physics equations, while being close enough to the measured PC-MRI velocity profile. The feasibility of the proposed approach is demonstrated on data from the carotid bifurcation of one healthy volunteer, and also from a pulsatile carotid flow phantom. The proposed solver produces flow fields that are in better agreement with direct PC-MRI measurements than CFD alone, and converges faster, while closely satisfying the fluid dynamics equations. For the implementation that provided the best results, the signal-to-error ratio (with respect to the PC-MRI measurements) in the phantom experiment was 6.56 dB higher than that of conventional CFD; in the in vivo experiment, it was 2.15 dB higher. The proposed approach allows partial or complete measurements to be incorporated into a modified CFD solver, for improving the accuracy of the resulting flow fields estimates. This can be used for reducing scan time, increasing the spatial resolution, and/or denoising the PC-MRI measurements.

CFD is an alternative that has been used to predict flow patterns in various vascular geometries, including intracranial aneurysms [10], the thoracic aorta [11], and the carotid bifurcation, both in models [12][13][14][15] and in vivo [16]. The equations describing Newtonian fluid flow are solved numerically for specified boundary and initial condition data. Such approach provides arbitrarily high spatial and temporal resolution, and is in principle capable of estimating flow fields for arbitrarily complex vessel geometries. Absolute hemodynamic parameter estimates can be obtained directly from the high-resolution flow fields produced by CFD, obviating the need for data smoothing or interpolation schemes.
The accuracy of conventional CFD routines hinges on many modeling assumptions that are not strictly true for in vivo vascular flow, including rigid vessel walls and uniform blood viscosity. Indeed, the proper choice of the underlying physics model is itself an open research question [10]. CFD predictions have so far shown variable agreement with PC-MRI measurements [10,11,13,15], and the applicability of CFD to robust flow estimation is still being debated.
The use of fluid mechanics techniques for improving PC-MRI data is an active research field. Several algorithms from the literature use regularizations based on curl and divergence of the velocity field, which are associated with the irrotationality and incompressibility characteristics of the fluid flow, respectively. These include algorithms capable of improving streamlines [31,32], and also algorithms for denoising the PC-MRI data [33][34][35][36].
However, the solutions found by those methods do not necessarily satisfy the Navier-Stokes momentum equation. Other authors integrated point-based measurements within a CFD solver, by adding a "force" term to the momentum equations that is proportional to the difference between predicted and measured velocities for a given grid point [37][38][39]. They used such an approach to integrate, respectively, Doppler-ultrasound velocity measurements and cerebral aneurysm blood flow MRI data into a CFD solver.
More recently, a method to accelerate 4D cardiac flow MRI using CFD simulations was proposed. The image model was generated by integrating numerical blood flow simulations (calculated using openFOAM [40]) into the MRI image-reconstruction algorithm [41]. However, this approach can not guarantee that the fluid physics model (momentum and continuity equations) is satisfied by the reconstructed velocity map.
Conventional CFD uses medical imaging data only to specify the vessel geometry and the flow at the inlet and outlet boundaries, or other previously known initial and boundary condition data, and uses the assumed fluid physics model to find the solution in the interior of the calculation domain [28]. The goal of this work is to develop a more general, flexible and easy to implement numerical framework for harnessing additional PC-MRI velocity measurements to construct a more robust and potentially more accurate CFD-based solution, considering PC-MRI data as ground truth. The proposed method is able to make use of full (or incomplete) PC-MRI measurements of one or more velocity components within the entire 3D volume. This is achieved through generalized Tikhonov regularization [42], obtaining a numerical solution that is close enough to the measured flow data; satisfies the fluid physics equations; reduces noise; and, in the clinical environment, can be used to reduce scan time.
Finally, this work is presented as a proof of concept of the CFD-MRI combined solver. All simulations herein were made using the finite volume method and SIMPLER algorithm in Cartesian grids, with unrealistic assumptions about the blood flow model, such as rigid walls and Newtonian viscosity. Nevertheless, the optimal numerical solution proposed in this work is general enough to be implemented: for any type of discretization method, such as finite differences, finite volume or finite element; for steady or unsteady flow; and, for any realistic physics model, such as non-Newtonian viscosity, elastic vessel walls, or slightly compressible flow. Even in the most realistic model, the discretization (finite differences, finite volume or finite elements) of the nonlinear set of differential equations produces a large and sparse system of linear equations, that forms the basis of the proposed numerical solution. The feasibility of the proposed approach is demonstrated on data from the carotid bifurcation of one healthy volunteer, and from a pulsatile carotid flow phantom. Two implementations of the regularized computational solution were evaluated and compared: one using only the PC-MRI data corresponding to the main velocity component (z axis); and another, using PC-MRI data corresponding to all three velocity components.

Blood flow model
The general model for fluid motion in 3D Euclidian space is given by the Navier-Stokes momentum equation [43]: where ρ is the fluid density, � ν = (u, v, w) is the flow velocity vector (u, v, and w are the velocity components associated with spatial axes x, y, and z, respectively), t is time, ∇ is the gradient operator, p is the hydrodynamic pressure, τ is the stress tensor, and b represents the body forces acting on the fluid during the flow. The stress tensor represents the momentum transferred in virtue of the molecular motions and interactions within the fluid. It is a function of the scalar invariants of the strain rate tensor ê = (1/2) ∇� ν + (∇� ν) T , where T denotes the transpose operation. For an incompressible fluid, it can be written as τ = −µ(ê)ê, where the scalar µ(ê) is the generalized Newtonian viscosity for a given ê.
In this work, blood is modeled as a Newtonian, incompressible and isothermal fluid, with constant viscosity µ and constant density ρ. We are also assuming that there are no body forces acting on the blood flow. Then, the simplification of the general momentum equation, Eq. (1), provides our blood model equation [43]: where is the Laplacian operator.
Since there are no sources of blood inside an artery, the flow field must also satisfy mass conservation [43], which is expressed by the continuity equation:

SIMPLER algorithm
Equations (2) and (3) must be solved for the unknown scalar field variables u, v, w, and p. Those equations are non-linear and coupled, and attempting to solve them directly in one step is a formidable, if not impossible, task.
The semi-implicit method for pressure-linked equations revised (SIMPLER) algorithm [44] is a well-known and established numerical routine for solving the momentum and continuity equations, subject to given boundary and initial conditions. It belongs to a class of algorithms capable of solving the non-linear coupled fluid dynamic equations, which also includes the SIMPLE, SIMPLEC, and PISO algorithms [45]. For our purposes, SIMPLER's major advantage is that it does not require an initial guess for the pressure field; instead an initial estimate of the velocity field is used.
The discretization of the momentum equation, Eq. (2), forms the basis of the iterative CFD routine, yielding three linear systems. Let N be the total number of grid points in the discrete 3D calculation domain, i.e., in a rectangular grid, N = N x · N y · N z , where N x , N y , and N z represent the number of grid points along the x, y, and z axes, respectively. Then, for the n-th iteration, we have: where S u,n−1 , S v,n−1 , and S w,n−1 are N × N square hepta-diagonal sparse matrices, each containing previous iteration information about all three velocity components, as well as the values of the density and viscosity constants; the three N × 1 column vectors u n , v n , and w n store the current iteration of the u, v, and w velocity component values, respectively, associated with all grid points in the 3D calculation domain (Fig. 1); and each of the constant N × 1 column vectors f u,n−1 , f v,n−1 , and f w,n−1 contains previous iteration information about all three velocity components, current iteration pressure difference values, and the physical parameters ρ and µ.
(4) S u,n−1 u n = f u,n−1 The main difficulty when attempting to predict the flow of an incompressible fluid is that there is no equation for pressure. Hence, a discretized pressure equation is deduced applying Eqs. (4), (5), and (6) in the discretized continuity equation, giving rise to a discretization of the following Poisson equation: where δt is the time step, i.e. the time increment between iterations (this is further discussed in the next section). For a given time instant t = t i , convergence is achieved when ∇ · � ν = 0. The pressure field is updated at each iteration based on the current velocity field estimate, and hence does not appear explicitly in Eqs. (4), (5), and (6). It is important to note that u, v, and w values must be defined on regular grids, staggered by half a grid spacing (along the three directions) with respect to the grid on which pressure values are defined. This is to avoid non-physical wiggle solutions for the pressure and velocity fields [44].
The steps of the SIMPLER algorithm, for a given time instant t = t i , are summarized in Algorithm 1. A detailed explanation of the discretization of the Navier-Stokes equations is provided in Refs. [44,45].  14:110 Algorithm 1 Steps of the SIMPLER algorithm, for a given time instant t = t i . [44] 1 Begin with a initial guess for the velocity components, u 0 , v 0 , and w 0 . 2 Using Eqs. (4), (5), and (6), the final velocity component estimates from the previous iteration (u n−1 , v n−1 , and w n−1 ) -or the initial guess, if this is the first iteration -and pressure field pn = 0 (pn is a N × 1 stacked column vector, as illustrated in Fig. 1; and 0 is the N × 1 null vector), find a first estimate of un, vn, and wn. 3 Using Eq. (7) and the velocity field found in step 2, find an estimate of the pressure field, pn. 4 Using Eqs. (4), (5), and (6), the final velocity component estimates from the previous iteration (u n−1 , v n−1 , and w n−1 ) -or the initial guess, if this is the first iteration -and the pressure field found in step 3, find an updated estimate of un, vn, and wn. 5 Using Eq. (7) and the velocity field found in step 4, find an updated estimate of pn. 6 Using the pressure field found in step 5, find the final estimate of the velocity field, adding α∇pn -where α is a relaxation factor -to each of the velocity component estimates (un, vn, and wn). 7 Repeat from step 2 until convergence, i.e., ∇ · = 0, where is the vector field represented by the velocity component estimates (un, vn, and wn) found in step 6. Note that correcting the pressure field after step 6 is pointless, as it will be reset to a null vector when step 2 is repeated.

Algorithm implementation
On constructing the numerical solution of the unsteady Navier-Stokes equation, we assume that a velocity field and the boundary conditions at a given time instant t = t 0 are known. For this initial set of data, the numerical solution for the next time step t 1 = t 0 + δt is constructed, and it converges toward the solution when the continuity equation is satisfied (step 7 in Algorithm 1). Starting with the solution for t 1 , the same iterative procedure is repeated to obtain the solution for t 2 = t 1 + δt, and so forth. In this manner, a time-dependent flow field is computed.
In unsteady flow simulations of the Navier-Stokes equations by implicit numerical routines, the nature of the transient SIMPLER iterative procedures is equivalent to steady-state SIMPLER calculations applied, until convergence, for each time instant [45]. In other words, solving transient problems using SIMPLER is equivalent to solving successive steady-state problems. Moreover, steady-state calculations may be interpreted as pseudo-transient solutions with spatially-varying time steps [45]. In other words, the steady-state solutions are, in practice, unsteady solutions, considering a virtual time step δt with fixed boundary conditions and initial data. This approach has been used by many recent authors [13,16,[37][38][39]. While this was the numerical strategy adopted in this work, the approach proposed here can also be used for unsteady flow predictions.
The steady-state solution � ν ∞ corresponding to a given cardiac phase (a temporal frame within the cardiac cycle) is calculated using the MRI-measured inlet and outlet velocities for that cardiac phase. CFD calculations begin with an initial guess for ν, and simulations are carried forward in time until convergence, i.e.: given a suficiently small ε (�·� denotes vector magnitude) and suficiently small time step δt. Note that, here, time t is a simulation-only parameter, and is unrelated to time instants within the cardiac cycle. It is also not related with the iteration steps (n) of the SIMPLER algorithm (Algorithm 1), since the entire algorithm-with multiple iterations until convergence criterion ∇ · � ν = 0 is satisfied-is performed at each time instant t, until the convergence criterion shown in Eq. 8 is satisfied. At this point, � ν ∞ is obtained. If multiple cardiac phases were to be reconstructed, then � ν ∞ would be independently calculated for each cardiac phase.
Our implementation of the SIMPLER algorithm was validated with the bidimensional lid-driven cavity flow problem, known in the literature as a benchmark for testing CFD algorithms [46][47][48]. All algorithms were implemented in Matlab (The MathWorks, Inc., Natick, MA, USA). Linear systems were solved using the biconjugate gradients stabilized method.

Proposed numerical solution
In this paper, we solve for a simulated velocity field, , that is close enough to the MRI-measured vector field � ν mri = (u mri , v mri , w mri ), and satisfies the fluid dynamics equations, Eqs. (2) and (3).
Let M be the total number of voxels in the reconstructed ν mri 3D velocity field, i.e., where M x , M y , and M z represent the number of voxels along the x, y, and z axes, respectively. Consider u mri , v mri , and w mri as the stacked M × 1 column vectors with the PC-MRI measurements. Since the numerical solution of the Navier-Stokes continuity-equation system is based on the solution of linear systems, we propose that our numerical optimal solution is obtained by minimizing, for each velocity component, at iteration n, the following equations: The first term on the right hand side of Eqs. (9)-(11) is related to the numerical solution of the Navier-Stokes continuity equations, and the second term is related to the comparison between the numerical solution and the PC-MRI velocity field. The S matrices and f vectors are defined in Eqs. (4)- (6), but note that we dropped the "n − 1" subscripts for simplicity and clarity; these are updated by velocity and pressure values calculated in the previous iteration. Coefficients u , v , and w are regularization factors, which weight consistance with PC-MRI data against conformance with the momentum equations.
Matrices Ŵ u , Ŵ v , and Ŵ w are of size M × N, and model the blurring effects due to finite k-space coverage in PC-MRI (this is further discussed below), while adjusting the number of points on the CFD grid in order to allow a comparison between ν n and ν mri . In this approach, the number of grid points in the CFD and MRI grids are not necessarily the same; we can use a finer grid in CFD than in MRI, for example. The optimal solutions for Eqs. (9)-(11) are straightforward [42], and given by To understand the construction of the Ŵ matrices, consider that in the absence of noise, artifacts, and distortions, the MRI-measured vector field � ν mri = (u mri , v mri , w mri ) is a blurred version of the true vector field � ν = (u, v, w). For the u component, for example, we can write: where * denotes convolution, and blurring kernel ψ u (x, y, z) is the point-spread function associated with the k-space coverage that was used when measuring u mri . Similarly, we can write v mri = v * ψ v and w mri = w * ψ w . If all three PC-MRI velocity components are measured using the same k-space coverage, then ψ u = ψ v = ψ w . For a 3DFT acquisition, these spatial blurring kernels are equal to where δx, δy and δz are the spatial resolutions of ν mri along the x, y, and z axis, respectively.
We want the CFD-estimated vector field � ν ∞ = (u ∞ , v ∞ , w ∞ ) to be an accurate representation of the true vector field ν. If this is so, then we should expect u mri ≈ u ∞ * ψ u , v mri ≈ v ∞ * ψ v , and w mri ≈ w ∞ * ψ w . The discretization of these equations yields three linear systems. Then, using the same notation introduced earlier, for the nth iteration of the CFD algorithm, we can write: The coefficients of Ŵ u , Ŵ v , and Ŵ w are calculated from ψ u , ψ v , and ψ w , respectively. If all three PC-MRI velocity components are measured using the same k-space coverage, and reconstructed onto identical grids, then The MRI-guided CFD estimate corresponding to one cardiac phase was calculated as a steady-state solution � ν ∞ . All three components of the PC-MRI velocity field ν mri measured at the z positions at the boundaries of the calculation domain were used as inlet and outlet boundary conditions for that cardiac phase. Note that this steady state solution � ν ∞ is the closest fit in the least-squares sense to the direct PC-MRI measurements that satisfy both momentum equation (Eq. 2) and continuity equation (Eq. 3). This is guaranteed by the fact that the optimal solutions Eqs. (12)- (14) are solved in each iteration of the SIMPLER algorithm (steps 2 and 4, in Algorithm 1), and by the convergence criterion (step 7).
In each of our experiments, all three PC-MRI velocity components were measured using the same k-space coverage, and reconstructed onto identical grids. In the phantom experiments, we used the same grid size for both ν mri and � ν ∞ , because the phantom data were measured with high spatial resolution. In these experiments, ν mri was reconstructed without zero-padding, i.e., onto δx × δy × δz voxels, and the CFD grid points were defined at the center of each of ν mri 's voxels. Hence, Ŵ u = Ŵ v = Ŵ w was defined as an identity matrix. In the in vivo experiments, ν mri was reconstructed using 2-fold zeropadding along each of the spatial axes, since the data was acquired with low spatial resolution. Then, Ŵ was an N × N symmetric matrix, with coefficients calculated from the point spread function ψ(x, y, z), defined in Eq. (16). This infinite support point spread function was truncated by multiplication with the box function where rect (w) = 1, if −1 ≤ w ≤ 1, and rect (w) = 0, otherwise.

Experimental setup: phantom demonstration
PC-MRI data of a pulsatile carotid flow phantom (Phantoms by Design, Inc., Bothell, WA) (Fig. 2) were obtained with high spatial resolution and high signal-to-noise ratio, from four time-resolved 3DFT FGRE image volumes (three acquired each with a velocity encoding bipolar gradient on one of the three axes, and one without a bipolar gradient). The scan parameters were: 0.5 × 0.5 × 1.0 mm 3 spatial resolution; The through-slab (z) axis was oriented along the S/I direction. The phantom's pulse cycle was set to 60 bpm. Only the temporal frame corresponding to peak flow was reconstructed. PC-MRI velocity component maps u mri , v mri and w mri were calculated using data from all channels of the receive coil array. The lumen was segmented by manually outlining the vessel borders from a stack of 2D axial images, obtained from the reconstructed 3D volume. A few voxels presented phase-wrap artifacts; these voxels were manually identified and their velocities were corrected by adding 2π to their values.
The combined solver calculations assumed fluid viscosity of µ = 0.005 Pa s and density of ρ = 1100 kg/m 3 (these values were provided by the phantom manufacturer). Calculations were performed with time step δt = 0.1 ms on a Cartesian grid of 0.5 × 0.5 × 1.0 mm 3 voxel size.
The CFD simulation domain was rectangular, of size 32.5 × 9.0 × 41.0 mm 3 . Each iteration required about 10 seconds of computation time on an Intel Core i7 processor running at 2.8 GHz.
Three simulated steady-state velocity fields � ν ∞ were obtained: 1. using the conventional SIMPLER algorithm, i.e., not using the PC-MRI data to constrain the CFD solution ( ν mri was used only as inlet and outlet velocities for the geometry); 2. using the velocity component associated with the main flow axis (z) measured with PC-MRI (w mri ) to constrain the CFD solution (u and v components were determined solely from the fluid physics model); and 3. using all three velocity components measured with PC-MRI (u mri , v mri , and w mri ) to constrain the CFD solution. The first approach is equivalent to making w = u = v = 0; in the second approach, we used w = 1 and u = v = 0; in the third approach, we used u = v = w = 1 . All three approaches used all three components of the PC-MRI velocity field ν mri measured at the z positions at the boundaries of the calculation domain as inlet and outlet boundary conditions. The number of iterations until convergence for the above simulations was 89, 40 and 5 iterations, respectively.

Experimental setup: in vivo demonstration
PC-MRI data of the carotid bifurcation of one healthy volunteer were obtained from four time-resolved 3DFT FGRE image volumes (three acquired each with a velocity encoding bipolar gradient on one of the three axes, and one without a bipolar gradient). The scan parameters were: 1.0 × 1.0 × 2.5 mm 3 spatial resolution; FOV 7.5 × 12.0 × 36.0 cm 3 ; TR 7.0 ms; flip angle 15 • ; temporal resolution 56 ms; VENC 160 cm/s; 7 min per scan; 1 NEX. The data were acquired on a GE Signa 3T EXCITE HD system (40 mT/m and 150 T/m/s max gradient amplitude and slew rate), with a 4-channel neck receive coil array. The through-slab (z) axis was oriented along the S/I direction. The institutional review board of the University of Southern California approved the imaging protocols. The subject was screened for MRI risk factors and provided informed consent in accordance with institutional policy.
Only the cardiac phase corresponding to peak flow was reconstructed. PC-MRI velocity component maps u mri , v mri and w mri were calculated using data from only one channel of the receive coil array. Residual linear velocity offsets in each velocity component map (e.g., due to eddy-currents) were removed by performing a linear fit to manually defined 3D regions containing only stationary tissue. The lumen was segmented by manually outlining the vessel borders from a stack of 2D axial images, obtained from the reconstructed 3D volume.
The combined solver calculations assumed blood viscosity µ = 0.0032 Pa s and density of ρ = 1060 kg/m 3 [49]. Calculations were performed with time step δt = 0.25 ms on a Cartesian grid of 0.50 × 0.50 × 1.25 mm 3 voxel size. The CFD simulation domain was rectangular, of size 30 × 37 × 125 mm 3 (the PC-MRI data was cropped to match this grid size). Each iteration required about 180 s of computation time on an Intel Core i7 processor running at 2.8 GHz.
Three simulated steady-state velocity fields � ν ∞ were obtained, using the same three approaches used in the phantom experiment. The number of iterations until convergence for the simulations was 1058, 190 and 6, respectively.

Quantitative evaluation
The CFD-simulated velocity fields were quantitatively compared with the PC-MRI measurements by means of the signal-to-error ratio (SER). The SER measures the ratio between the energy of the signal and the energy of the estimation error. We considered the PC-MRI velocity field, � ν mri = (u mri , v mri , w mri ), as our ground-truth "signal"; consequently, the estimation error is the vector difference between the CFD-estimated velocity field, � ν ∞ = (u ∞ , v ∞ , w ∞ ), and the ground-truth field, ν mri . Thus, the SER is calculated (in decibels) as: where integers i, j, and k represent grid-point indexes along the x, y, and z axes, respectively. Similarly, the SER was also calculated individually for each of the velocity components, as: Using these SER values, the three CFD approaches-pure CFD, CFD driven by one PC-MRI velocity component, and CFD driven by all three PC-MRI velocity componentswere quantitatively evaluated and compared.

Evaluation of denoising properties
Under our hypothesis, CFD simulations provide a smooth, noise-free flow field. Therefore, we expect that the proposed approach can be used as a denoising mechanism for PC-MRI flow assessment. In order to verify the denoising effects of the combined solver, we added zero-mean Gaussian noise with standard deviation 8 cm/s to the phantom's measured velocity field, ν mri . This noisy flow field was used to constrain the CFD calculations, using the approach in which all three velocity components measured with PC-MRI are used. In this experiment, we used u = v = w = ; and four different values of were evaluated: 5 × 10 −9 , 5 × 10 −8 , 5 × 10 −7 , and 5 × 10 −6 . The SER between the proposed approach and the original PC-MRI measurements was calculated, and compared with the SER of the noisy flow field. The pure CFD approach, in which the noisy PC-MRI data are used only as inlet and outlet velocities for the geometry, was also evaluated (this is equivalent to making = 0).
The phantom data was used in this denoising experiment, because it was acquired using 9 NEX-which results in high signal-to-noise ratio (SNR), while the in vivo data was acquired using only 1 NEX. The noise levels on the phantom's measured velocity components-estimated as the standard deviation in regions of uniform mean velocity-are lower than 3 cm/s; while the SNR of the magnitude images exceeds 26 dB. The velocity-to-noise ratio (VNR) [50,51] for the u and w components reach 28 and 31 dB, respectively (the VNR for the v component was not calculated, because v is approximately null over the entire geometry).
Finally, in order to justify our denoising experiment, we analyze the noise distribution in PC-MRI images. We note that, from a maximum likelihood perspective, Eqs. (9)- (11) assume that the PC-MRI data is degraded by Gaussian noise. Under certain conditions, one can prove that velocity field noise in PC-MRI satisfies a zero-mean Gaussian distribution [52]. Therefore, the additive noise acting on the velocity fields can be assumed to be Gaussian distributed [27,52]. Hence, the proposed minimization is well-suited in terms of the MR noise distribution. Figure 3 shows a qualitative velocity-map comparison between the PC-MRI phantom measurements and the three simulated results. The PC-MRI velocity field does not satisfy the continuity equation, since its divergence is nonzero within the lumen (Fig. 3a). The pure CFD solution produced a velocity field that satisfies the physical model, but is considerably smoother than the PC-MRI measurements (Fig. 3b). Using one MRI-measured velocity component (w mri ) to guide the CFD simulation resulted in a solution that is qualitatively more similar to the MRI-measured field, while still satisfying the continuity and momentum equations (Fig. 3c). Even better agreement was achieved when all three MRI-measured velocity components (u mri , v mri , and w mri ) were used to guide the CFD simulation (Fig. 3d). These improvements can be also appreciated on a vector field visualization of the flow field over the entire tridimensional volume (Fig. 4). Table 1 shows the SER between the phantom experiment results from each of the three CFD approaches, relative to the MRI-measured velocity field. The approach using only one PC-MRI velocity component (w mri ) to drive the CFD calculations (labeled "CFD + 1D") provided a 11.09 dB improvement in SER relative to pure CFD (labeled simply "CFD") with respect to the w component; however, the improvement was only 1.81 dB when considering all three components. The approach using all three PC-MRI velocity components (u mri , v mri , and w mri ) to drive the CFD calculations (labeled "CFD + 3D") provided a 6.56 dB improvement in SER relative to pure CFD when considering all three components, while still providing a 8.02 dB improvement when considering only the w component.

Phantom demonstration
Note that, for each of the CFD approaches, the SER was lower for the u and v components than it was for w (Table 1). This can be explained by the fact that the same VENC was used for measuring all three PC-MRI velocity components (which implies similar noise levels), but the velocities along the z axis are considerably higher than those along the x and y axes (the energy of w mri was 11.42 dB higher than that of u mri , and 15.72 dB higher than that of v mri ). As a consequence, the SNR of w mri is substantially higher than that of u mri and v mri . Thus, even if the energy of the absolute error between the CFDestimated velocities and the MRI-measured velocities was the same for the three components, SER w would be higher than SER u and SER v . Also, note that CFD approaches provide a smooth, (ideally) noise-free velocity field, but the SER was calculated with respect to a noisy PC-MRI field. This means that the denoising properties of the CFD approaches could actually hurt the SER, especially if noise levels are relatively high when compared to the velocity values (this is the case for the u and v components). However, denoising is a desirable feature, so this does not necessarily indicate an unwanted result. Figure 5 and Table 2 show the results of the experiment in which the denoising properties of the proposed approach were evaluated. CFD calculations were constrained by all three PC-MRI velocity components, with added Gaussian noise. The combined  Table 1 Signal-to-error ratio (in dB) between the phantom experiment results from each of the three CFD approaches-pure CFD (labeled "CFD"); CFD driven by one PC-MRI velocity component (labeled "CFD + 1D"); and CFD driven by all three PC-MRI velocity components (labeled "CFD + 3D"), relative to the MRI-measured velocity field

Table 2 Signal-to-error ratio (in dB) between noisy and original PC-MRI measurements; and between the MRI-guided CFD estimates and the original PC-MRI measurements
Additive zero-mean Gaussian noise with standard deviation of σ = 8 cm/s was used. CFD estimates were obtained using the combined solver with different values of the weight parameter , constrained by the noisy PC-MRI measurements. All three PC-MRI velocity components were used to guide the CFD calculations. SER values were calculated according to Eqs.    v, and w), for an axial slice at the bifurcation of the carotid flow phantom: a PC-MRI; b PC-MRI with added Gaussian noise (σ = 8 cm/s); c pure CFD solution using noisy PC-MRI data as inlet and oulet velocities; d CFD guided by the noisy PC-MRI data, with = 5 × 10 −9 ; e CFD guided by the noisy PC-MRI data, with = 5 × 10 −8 ; f CFD guided by the noisy PC-MRI data, with = 5 × 10 −7 ; g CFD guided by the noisy PC-MRI data, with = 5 × 10 −6 solver improved the SER of each individual velocity component, for all different weight parameters we evaluated ( Table 2). The CFD solution constrained by PC-MRI using = 5 × 10 −8 (Fig. 5e) provides a velocity field that is less noisy and visually more similar to the measured PC-MRI velocity field (Fig. 5a) than the pure CFD solution obtained using the noisy PC-MRI velocity field as boundary data (Fig. 5c). Using smaller values of (Fig. 5d) results in solutions that are closer to the pure CFD solution (Fig. 5c).
These results illustrate the potential of the proposed numerical framework also as a denoising technique for PC-MRI. Figure 6 provides a qualitative velocity-map comparison between the PC-MRI in vivo measurements and the three simulated results. As in the phantom experiment, the PC-MRI velocity field does not satisfy the continuity equation, since its divergence is nonzero within the lumen (Fig. 6a). The pure CFD solution produced a velocity field that satisfies the physical model, but differs considerably from the PC-MRI measurements (Fig. 6b). Using one MRI-measured velocity component (w mri ) to guide the CFD simulation resulted in a solution that is qualitatively more similar to the MRI-measured field, while still satisfying the continuity and momentum equations (Fig. 6c). Even better agreement was achieved when all three MRI-measured velocity components (u mri , v mri , and w mri ) were used to guide the CFD simulation (Fig. 6d). These improvements can be also appreciated on a vector field visualization of the flow field over the entire tridimensional volume (Fig. 7). Table 3 shows the SER between the in vivo experiment results from each of the three CFD approaches, relative to the MRI-measured velocity field. The approach using only one PC-MRI velocity component (w mri ) to drive the CFD calculations (labeled "CFD + 1D") provided a 9.79 dB improvement in SER relative to pure CFD (labeled simply "CFD") with respect to the w component; however, the SER was only 0.76 dB higher when considering all three components, since the agreement with PC-MRI for the v component was made worst. The approach using all three PC-MRI velocity components (u mri , v mri , and w mri ) to drive the CFD calculations (labeled "CFD + 3D") improved the SER for all three components, by 0.24, 0.59, and 6.00 dB, respectively; as a result, there was a 2.15 dB overall improvement in SER relative to pure CFD, and a 1.39 dB improvement relative to CFD guided only by w mri . As in the phantom experiments, the SER was lower for the u and v components than it was for w, for all three CFD approaches. This can again be explained by the fact that the same VENC was used for all three velocity components, while the velocities along the z axis are considerably higher than those along the x and y axes (the energy of w mri was 29.76 dB higher than that of u mri , and 27.09 dB higher than that of v mri ), and is a direct (and desirable) consequence of the denoising properties of the proposed approach, as previously discussed.

Discussion
The proposed methodology uses the three-dimensional SIMPLER algorithm with Cartesian uniform meshes to perform blood flow simulations under the influence of threedimensional MRI-measured velocity profiles. The combined MRI-CFD methodology attempts to correct the MRI-measured flow field, forcing it to satisfy the fluid mechanics equations. The choice of the SIMPER algorithm and Cartesian discretization was performed in order to facilitate implementation of the algorithm. We showed that the proposed technique provides better agreement with the PC-MRI measurements than pure CFD simulations. We also showed that this MRI-guided CFD approach can be used as a means of reducing noise in the PC-MRI measurements. It can also be used to reduce computation time: when all three MRI-measured velocity components are used to guide the CFD simulations, only a few iterations are required for convergence, i.e., for finding the flow field that is the most similar to the PC-MRI measurements (in the least squares sense) while satisfying the fluid mechanics equations.  Table 3 Signal-to-error ratio (in dB) between the in vivo experiment results from each of the three CFD approaches-pure CFD (labeled "CFD"); CFD driven by one PC-MRI velocity component (labeled "CFD + 1D"); and CFD driven by all three PC-MRI velocity components (labeled "CFD + 3D"), relative to the MRI-measured velocity field SER values were calculated according to Eqs.   14:110 We believe that the proposed method can also be used as a technique for reducing scan time. The proposed methodology allows the use of only part of the velocity measurements obtained with MRI to guide computational solutions by appropriately choosing the Ŵ matrices in Eqs. (12)- (14). In this paper, we used identical grids for MRI and CFD. Therefore, the Ŵ matrices were square with size N × N. In each of the following suggested approaches, at least one of the Ŵ u , Ŵ v or Ŵ w matrices may be of size M × N , with M < N. Possible approaches for reducing scan time include: (1) acquiring the MRI data with reduced spatial resolution, and using the MRI-guided CFD simulation to improve the spatial resolution; (2) measuring portion(s) of the volume (e.g., carotid inlet and outlets) with full spatial resolution, while measuring the rest of the volume with reduced spatial resolution; (3) acquiring only a few slices along the bifurcation, and using the MRI-guided CFD simulation to fill in the gaps; (4) measuring one or two velocity components with full spatial resolution, while measuring the other component(s) with reduced spatial resolution; (5) measuring one or two velocity components across the entire volume, while measuring the other component(s) for only a few slices; and (6) any combination of these approaches. We plan to explore these ideas in future studies.
It is well known that blood is a non-Newtonian fluid, therefore its viscosity is not uniform. While there exists many constitutive models and studies in the literature regarding the non-Newtonian rheological properties of blood [43,53,54], no gold-standard constitutive model exists, and the assumption of constant whole blood viscosity is common practice [13,16,33,37,55,56]. Therefore, we used the Newtonian blood flow model (constant whole blood viscosity), which greatly simplified the implementation.
It is often desirable to estimate the wall shear rate at the carotid bifurcation. Determining wall shear rates from CFD simulation results would require using highly refined meshes around the neighborhood of the vessel wall. Our implementation of the SIM-PLER algorithm does not allow using spatially-varying grid spacing, and the mesh is uniform all over the integration domain. Using a very fine grid over the entire integration domain is possible, in principle. However, this would drastically increase the computational complexity, and could make the proposed methodology impractical in a clinical environment, for example.
Finally, the proposed approach does not take the effects of vessel wall elasticity into consideration. While the pulsatile carotid flow phantom we used has rigid tube walls, the human carotid vessel wall is generally elastic. While there exists blood flow CFD simulation methods that incorporate elastic wall effects [57,58], the assumption of rigid vessel walls is another common practice [13,16,33,37,55]. The SIMPLER algorithm implemented in this work uses a Cartesian uniform mesh, and does not allow the use of the elastic wall models.
Both wall shear stress calculations and elastic vessel walls could be properly addressed with an implementation using finite elements-which would allow an easier adaptation of the mesh near the wall and also allow simulating the effects of fluid-structure interaction. Despite the fact that the implementation of a finite-element solver would be substantially more complex than our implementation, all the proposed methodology described in this proof of concept paper can still be applied in the same fashion, since the problem of solving the set of differential equations in a finite element discretization is replaced by a sparse system of linear equations similar to the ones obtained in this work.
In future works, this MRI-guided CFD methodology for unsteady flow will be implemented on FreeFem++, 1 a partial differential equation solver capable of solving the Navier-Stokes equations using finite elements [59]. FreeFem++ allows access to the linear systems that can be modified in order to make use of the principles introduced here. With this approach, we expect more general, higher-quality solutions, allowing the calculation of biomarkers, such as wall shear stress, since FreeFem++ can handle different types of triangular finite elements, large varieties of linear system solvers, and automatic mesh adaptation. Note that there are other free software that could be used to reproduce this methodology, such as openFOAM, 2 which uses finite volume discretization. Both FreeFem++ and openFOAM allow modifications of the linear systems, and have CFD solvers already implemented, which could facilitate the implementation of the methodology proposed in this study [40,59].
All the assumptions and simplifications disscussed above (Newtonian blood flow, imprecise geometry, non-compliant walls) contribute non-linearly to the differences observed between CFD solutions and MRI-measured velocity fields. Using the MRImeasured velocity field to constrain the CFD solution indirectly addresses these simplifications, and provides a more realistic CFD solution. However, we are still unable to clearly identify the main factors responsible for the disagreement between MRI-measured and CFD-computed velocity fields. Nevertheless, an implementation using a more robust solver, as proposed above, could improve on these limitations and further reduce the gap between MRI-measured and CFD-computed results.

Conclusion
We have proposed a framework for obtaining flow field estimates that are influenced by both PC-MRI measurements and a fluid physics model. The results showed that the proposed technique provides better agreement with the PC-MRI measurements than pure CFD simulations, and has reduced computation time (faster convergence). MRI-guided CFD can be used to correct the MRI-measured flow field, forcing it to satisfy the fluid mechanics equations. It can also be used as a means of reducing noise in the PC-MRI measurements, and has potential as a method for reducing scan time.
The proposed framework offers a general approach to in vivo blood flow assessment, that is complementary to improvements in PC-MRI acquisition and reconstruction techniques, and can be applied to the study and diagnosis of a broad range of cardiovascular flow mapping applications.