Inverse Problem for Color Doppler Ultrasound-Assisted Intracardiac Blood Flow Imaging

For the assessment of the left ventricle (LV), echocardiography has been widely used to visualize and quantify geometrical variations of LV. However, echocardiographic image itself is not sufficient to describe a swirling pattern which is a characteristic blood flow pattern inside LV without any treatment on the image. We propose a mathematical framework based on an inverse problem for three-dimensional (3D) LV blood flow reconstruction. The reconstruction model combines the incompressible Navier-Stokes equations with one-direction velocity component of the synthetic flow data (or color Doppler data) from the forward simulation (or measurement). Moreover, time-varying LV boundaries are extracted from the intensity data to determine boundary conditions of the reconstruction model. Forward simulations of intracardiac blood flow are performed using a fluid-structure interaction model in order to obtain synthetic flow data. The proposed model significantly reduces the local and global errors of the reconstructed flow fields. We demonstrate the feasibility and potential usefulness of the proposed reconstruction model in predicting dynamic swirling patterns inside the LV over a cardiac cycle.


Introduction
Vortex flow imaging has recently attracted much attention owing to a strong correlation between intraventricular flow pattern and heart function [1][2][3]. Possible clinical indices of cardiac functions can be obtained by characterizing and quantifying the vorticity of intraventricular blood flow. There are several studies to compute and quantify the blood flow pattern inside the left ventricle (LV) by using ultrasound imaging scanner. Echo particle image velocimetry (E-PIV) is representative of the commonly used method [4,5], which tracks the speckle pattern of blood flow. However, the E-PIV is hardly a complete noninvasive method because it requires the intravenous injection of contrast agent.
As an alternative approach to reconstruct the velocity of blood flow, color flow ultrasound-based techniques have been proposed. Color flow ultrasound is also called C-mode images, color Doppler images, color Doppler data, or color Doppler ultrasound. Color flow images reflect the velocity components projected on the direction of ultrasound beam propagation [6]. They are represented as the colors of red and blue in a region of interest (ROI) box overlapped on echo images. In general, red and blue colors indicate the velocity components coming toward and receding from scanning probe, respectively. Among the color flow ultrasound-based techniques, one method is to reconstruct blood flow by using the assumption of 2D divergence-free blood flow [7,8]. It decomposes each 2D velocity vectors into radial component obtained from color flow data and unknown angular component and computes the unknown component of blood flow velocity under the 2D divergence-free condition. However, it has limitation of oversimplification that ignores out-ofplane flows. Another method is to use the beams of multiple directions of color Doppler ultrasound. It reconstructs the 2D or 3D velocities of blood flow by using color Doppler data acquired from the beams in two or three different directions [9,10]. However, the registration of multiple imaging planes by the multiple directional beams is a very challenging issue in practical environment. Recently, we proposed a 2D incompressible Navier-Stokes model to reconstruct intraventricular flows using color Doppler data and LV boundaries extracted from echocardiography data [11]. The model was designed to cope with out-of-plane blood flows for the imaging plane by introducing a mass source term of a source-sink distribution. Although the 2D model seemed to preserve the global kinetic energy and vortex strength during cardiac cycle, the predicted velocity and vorticity fields indicated that the model did not precisely capture the location and shape of dynamic vortex patterns.
In this paper, we propose a mathematical framework in Figure 1 for reconstructing the blood flow in LV using 3D echocardiographic images and the color Doppler intensity data. The framework includes a formulation of inverse problem for undetermined velocity and pressure fields, LV boundary tracking, and a forward simulation procedure to generate synthetic flow data. Under the assumption that the high frame-rate acquisition of 3D echocardiographic data is available in the whole LV region, the color Doppler data is directly embedded in the incompressible Navier-Stokes equations, along with 3D LV boundary reconstruction. An interpolation technique using multiple planar LV contours extracted from the echocardiographic images is applied to obtain the boundary conditions on the 3D LV boundaries. We perform the forward simulation with the LV boundary data from real intensity images. Based on the simulation results, one-directional velocity component of the flow data is obtained for synthetic color Doppler data. Using the synthetic data, intracardiac blood flows inside LV are reconstructed by solving the inverse model. Finally, we demonstrate the robustness of the proposed model in visualizing time-dependent vortex patterns over one cardiac cycle and quantifying local and global errors for the reconstructed flow fields.
x z y Figure 2: Description of 3D LV domain (Ω ) and the coordinate system. The arrows indicate the LV wall movements and flow directions over the inlet and outlet valves at the diastole and systole phases.

Problem Statement.
This section describes a mathematical framework of inverse problem of recovering the velocity distribution in LV by using color Doppler data. Let k(x, ) denote the velocity of blood flow at position x and time , and let Ω be a 3D LV region at time , as shown in Figure 2. Assuming that blood flow is an incompressible Newtonian fluid, the velocity k is governed by the incompressible Navier-Stokes equations inside the time-varying LV region Ω : Computational and Mathematical Methods in Medicine   3 where , , and are the density, viscosity, and pressure of the blood, respectively. To determine k in (1) uniquely, we need to impose proper boundary conditions of LV wall motion involving the inlet and outlet conditions over the inlet and outlet valves that are denoted by Γ and Γ in Figure 2, respectively. Since the velocity k obtained by solving (1) is sensitive to boundary conditions, it is very important to extract accurate LV wall motion. However, it is almost impossible to capture 3D LV wall motion accurately in practical environment. Therefore, supplementary information of k is necessary for computing blood flow reliably.
In this paper, we use the color Doppler data D, as the additional information of one component of velocity, of the form where a is the ultrasound beam direction. Based on (1) and (2), our reconstruction model assimilates velocity components with the additional information. Details of the proposed reconstruction model are explained in the following section.

Reconstruction Model.
We assume that the color Doppler data provides the first component of velocity (to be precise, a = [1 0 0] ). Using the inner product of a and the momentum equation in (1) with (2), we have (3) The above system of equations can be expressed as the following matrix form: where the relation between and k is given by Remark 1. According to (4), if it is possible to measure subsidiary Doppler data D 2 along a direction a 2 ̸ = a, the velocity k can be explicitly expressed by ) .
Hence, if two pieces of linearly independent Doppler data are available, then the inverse problem can be simply solved without the information of LV boundary Ω .
To solve the velocity components V 2 and V 3 that satisfy (4), we consider an iterative procedure as follows: for the th iteration, we find V +1 in Ω , where k = (D, V 2 , V 3 ) and +1 satisfy Note that proper boundary condition is needed to guarantee the uniqueness of in (8). To describe boundary conditions for V +1 2 , V +1 3 , and +1 , we first simplify LV boundary near the valves by considering linearly connected parts to the valves in the LV domain. Since it is still challenging to accurately capture the valve geometry in B-mode images, we focus on the LV wall motion, while we ignore the motion of valves in this paper. Then, we impose proper boundary conditions on the simplified boundary. Considering the simplified LV domain Ω , the corresponding boundaries near valves and physical LV wall are defined as Γ and Γ = Ω \ Γ , 4 Computational and Mathematical Methods in Medicine respectively. Thus, the proposed reconstruction model can be rewritten as, in Ω , with the boundary conditions where k and +1 satisfy Here, the velocity components V wall 2 and V wall 3 are estimated from LV boundary segmentation. Moreover, the Neumann boundary condition is imposed to allow the reconstructed flow to pass through Γ .
The 3D domain is discretized by × × node points, where , , and are the numbers of nodes along , , and directions, respectively, with uniform grid size ℎ. we define discrete differential operators based on the finite difference method as follows: where ⊗ is the Kronecker product and Using the above discrete operators, the following discrete system for ( , , ) ] can be derived: Here, p +1 is the solution to where . Note that diagonal entries of C and C in (13) are the derivatives of measured data D. Since the linear system in (13) may not be diagonally dominant and lead to severe numerical instability, we consider a minimization problem for k +1 = [k +1 2 , k +1 3 ] with a regularization term as follows: Computational and Mathematical Methods in Medicine 5 where is the regularization parameter and The first term in (15) is a penalty term which makes a solution satisfy the proposed reconstruction model, while the second term mitigates the ill-posedness of the diagonal matrices C and C . The minimizer of (15) is given by where = × × . During the iterative procedure, we used a stopping criterion as

Initial Guess.
For better convergence in solving the minimization problem in (15), we determine a proper initial guess k 0 . For vorticity fl ∇ × k, we assume that projection of the vorticity in a direction d parallel to the imaging plane passing through two valves is negligible: For computational simplicity, the third component of d is set to zero: Here, 1 and 2 are determined by the Doppler data D and divergence-free condition. From (19) and (20), we have Substituting V 1 = D in (21) yields We consider two combinations of and directional derivatives of (22) and the second equation of (3) ⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟⏟ such that − / (22) + 1 ( / )( * ) and / (22) + 1 ( / )( * ). Then, we have the following two equations: with boundary conditions Using the discrete operators in Section 2.2, we obtain the following linear system: where L = (L ⊗I ⊗I )+(I ⊗L ⊗I ) and L = I ⊗ I ⊗ L . Therefore, we obtain the initial guess k 0 = (V 0 2 , V 0 3 ) by solving (25).

Estimation of Boundary Condition: LV Contour
Segmentation. The proposed reconstruction model in (9) requires time-dependent boundary conditions defined in (10). In this study, for the 3D velocity boundary conditions on the LV wall, time-varying LV boundaries are tracked and local displacements of the LV wall are estimated. Similar to the LV boundary tracking in 3D echocardiographic data [12], we generate 3D LV surface by interpolating planar LV contours tracked on multiple 2D echocardiographic imaging planes, where each 2D LV contour is tracked by LV boundary tracking method on each imaging plane. Among various 2D LV tracking methods, we use the LV tracking method combined with the Lucas-Kanade method and a constraint formulated by the global deformation of nonrigid heart motion proposed 6 Computational and Mathematical Methods in Medicine in [13]: for each th tracking point x ( ) = ( ( ), ( )), we estimate its velocity k = x / by minimizing . . .

Setting for Forward Simulation.
For the forward simulation, we use a series of real 3D LV volume data, acquired by a Siemens ACUSON SC2000 imaging system with a 4Z1c probe. The acquisition conditions are given by the transmitting center frequency 2.8 MHz, the mechanical index 0.92, and the thermal index 0.46. The dataset consists of 10 frames within the whole cardiac cycle. The multiple slices are set to sagittal and coronal planes, as shown in Figure 3, and we apply a 2D LV contour tracking method independently of the sagittal and coronal planes. Although ultrasound images including echo and color flow Doppler are obtained from real-time measurements, the scanning time for ultrasound images is bounded below due to dependence on the sound speed. Lately, the development of ultrafast imaging system is ongoing in some research groups. Such a system is expected to have a frame-rate higher than 1,000 fps (for echo images)  using plane wave beam-forming techniques. Unfortunately, we do not have such an imaging system yet but have real 3D echocardiography data. Figure 4 illustrates 3D LV surface reconstructed by simplifying valve geometry into the pipe-type structure and integrating the LV contours tracked on the two orthogonal imaging planes A and B, which correspond to the coronal and sagittal planes in Figure 3, respectively. In the 3D LV model, the lengths of two pipes corresponding to mitral and aorta valves are set to be similar to the axial dimension of the LV so that the pressure is developed well from the end of the pipes to the interior of LV region. The cross sections of mitral and aorta pipes are modeled by elliptical and circular shapes, respectively.
We use commercial software programs for convenience in performing numerical simulations of the forward problem. Multiphysics is used to solve (1), assuming that the nodes of the LV boundary move along the outward normal direction n with the average speed at each time step. The corresponding velocity k at the nodes of the LV boundary are defined as The Neumann boundary condition for pressure is applied at the nodes: Since the aorta valve is closed but the mitral valve remains open during the diastole phase, we consider that no viscous stress conditions for k with constant pressure are applied at the inlet valve Γ , while zero-velocity conditions are applied at the outlet valve Γ : n ⋅ (∇k + ∇k ) = 0, During the systole phase, the boundary conditions in (31) are imposed in the opposite way: In this study, we set = 1050 kg/m 3 and = 0.00316 Pa⋅s for the density and viscosity of the blood flow, respectively [14]. Stroke volume is about 70 mL, heart rate is 60 per minute [15], and ratio of the diastole to one cardiac cycle is set to 0.6 [16].

Forward Simulation of Blood Flow in LV.
We perform the forward numerical simulation using the FSI model in COMSOL Multiphysics in order to obtain 3D synthetic flow data inside the moving LV domain. Figure 5 shows synthetic intraventricular velocity and out-of-plane vorticity fields on two orthogonal imaging planes. The velocity and

Reconstruction of Blood Flow in LV.
Before proceeding further, we investigate the errors of the reconstructed velocity fields that are assessed with respect to the change of 2 from 10 −7 to 10 −1 in order to optimize the regularization parameter in (15). The reconstruction error is determined based on 2norm error between the reconstructed velocity fields k and the forward data k : where is the number of time steps. Figure 6 indicates the minimum error at 2 = 10 −6 , which is used to reconstruct the blood flow in LV. Table 1: Comparison of 2 -norm errors for the velocity field recovered by the proposed reconstruction model and synthetic velocity field. The first and second columns represent normalized pointwise error of velocity and vorticity in 2 , respectively. The third and fourth columns represent normalized 2 -norm in global energy estimates of velocity and vorticity, respectively. Here, Note that the superscripts and denote the forward and reconstructed data, respectively. We consider one-directional velocity component of the synthetic flow data from the forward simulation as color Doppler data. Based on the velocity component, the velocity fields in LV are reconstructed by solving the minimization problem in (15). For better convergence, the initial guess k 0 is obtained by solving (25) at each time step. The reconstructed solution k is obtained by repeatedly updating the minimizer in (17) until the stopping criterion (18) is satisfied.  Figure 7(l). It is worth noting that the proposed model provides accurate reconstructions of the blood flow in a strong vortex region, while the model seems to smooth out small-scale vortex patterns due to the regularization. The maximum magnitude of local velocity error is less than 0.16 m/s and 0.29 m/s during the diastole and systole phases, respectively. Note that the maximum speeds in the diastole and the systole cycle are 0.62 m/s and 1.2 m/s, respectively. This implies that the maximum error scaled by the maximum speed in LV is around 25% during the whole cardiac cycle.
For further quantitative comparison, we investigate reconstruction errors of the velocity and vorticity fields based on the local 2 -norm errors and the differences in the global energy estimation, as listed in Table 1. The local velocity errors are less than 21% during the diastole phase, while the errors are slightly increased up to 32% during the systole phase. In the diastole phase except for the early stage, the local vorticity errors are less than 30%. However, the vorticity errors are increased up to 50% in the rest of the cardiac cycle, because the blood flow in LV exhibits weak vorticity patterns, relative to the diastole phase. Although the local reconstruction errors are not ignorable, the differences of the global energy estimates between the reconstructed flow fields and reference data are less than 7% both for the velocity and for vorticity fields. Thus, the pointwise errors at the local regions may not affect the major features of vortex flows. Furthermore, compared to the previous 2D reconstruction model [11], Table 1 confirms that the proposed 3D model provides better reconstruction performance in terms of local and global errors of the velocity and vorticity fields.

Conclusions
We proposed a mathematical framework involving a reconstruction model for three-dimensional blood flow inside LV. This framework consists of the extraction of time-varying LV boundaries from multiple echocardiographic images, the forward simulation of the blood flow, and the reconstruction model that combines 3D incompressible Navier-Stokes equations with one-direction velocity component from the synthetic flow data (or color Doppler data) from the forward simulation (or measurement). Assuming that an ultrasound imaging device provides color Doppler data in the whole 3D LV region, we formulated an inverse problem to reconstruct the blood flows, directly embedding the Doppler data in the Navier-Stokes equations and incorporating with the extracted LV boundaries. To generate synthetic Doppler data, we performed the forward simulation of the blood flow using the FSI model and extracted one-directional velocity component of the synthetic flow data. The blood flow inside LV was reconstructed by solving the inverse problem with a proper initial guess of unknown velocity components. Similar to the forward simulation, time-evolving vortex patterns over a cardiac cycle were clearly found in the reconstructed velocity and vorticity fields obtained from the proposed model. We also quantified the reconstruction errors based on the local and global differences between the reconstructed and synthetic flow data. Compared to the previous reconstruction model [11], the proposed model significantly improves the performance of the reconstruction of the blood flow. Through the numerical simulation, we demonstrated the feasibility and potential usefulness of the proposed model in reconstructing the intracardiac blood flows inside the LV.