Abstract

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.

1. Introduction

Vortex flow imaging has recently attracted much attention owing to a strong correlation between intraventricular flow pattern and heart function [13]. 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-of-plane 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.

2. Methods

2.1. Problem Statement

This section describes a mathematical framework of inverse problem of recovering the velocity distribution in LV by using color Doppler data. Let denote the velocity of blood flow at position 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 is governed by the incompressible Navier-Stokes equations inside the time-varying LV region : where , , and are the density, viscosity, and pressure of the blood, respectively. To determine 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 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 is necessary for computing blood flow reliably.

In this paper, we use the color Doppler data , as the additional information of one component of velocity, of the formwhere 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.

2.2. Reconstruction Model

We assume that the color Doppler data provides the first component of velocity (to be precise, ). Using the inner product of and the momentum equation in (1) with (2), we have The above system of equations can be expressed as the following matrix form: where the relation between and is given by

Remark 1. According to (4), if it is possible to measure subsidiary Doppler data along a direction , the velocity 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 and that satisfy (4), we consider an iterative procedure as follows: for the th iteration, we find and such that where and satisfy Note that proper boundary condition is needed to guarantee the uniqueness of in (8). To describe boundary conditions for , , and , 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 , respectively. Thus, the proposed reconstruction model can be rewritten as, with the boundary conditions where and satisfy Here, the velocity components and 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: , , , and , where is the Kronecker product and

Using the above discrete operators, the following discrete system for , , and can be derived: Here, is the solution to where , , , and . Note that diagonal entries of and in (13) are the derivatives of measured data . Since the linear system in (13) may not be diagonally dominant and lead to severe numerical instability, we consider a minimization problem for with a regularization term as follows: 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 and . The minimizer of (15) is given by where . During the iterative procedure, we used a stopping criterion as

2.3. Initial Guess

For better convergence in solving the minimization problem in (15), we determine a proper initial guess . For vorticity , we assume that projection of the vorticity in a direction parallel to the imaging plane passing through two valves is negligible: For computational simplicity, the third component of is set to zero: Here, and are determined by the Doppler data and divergence-free condition. From (19) and (20), we have Substituting in (21) yieldsWe consider two combinations of and directional derivatives of (22) and such that and . 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 and . Therefore, we obtain the initial guess by solving (25).

2.4. 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 in [13]: for each th tracking point , we estimate its velocity by minimizing where are neighborhood pixels of and is the regularization parameter and The first term in (26) represents local motions corresponding to the movements of speckle pattern in B-mode image based on Lucas-Kanade method, while the velocity is regularized by the second term. From , for , 2D LV contour at time is interpolated so that for .

2.5. 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. FreeCAD® is used to integrate the contours into 3D domain at each frame and COMSOL® Multiphysics is used to import the reconstructed 3D domains. The volume of the LV domain and the area of LV wall are then computed at each time step. We denote the variations of volume and area by and , respectively. Here, we consider as an average speed at all surface points of the reconstructed LV. The fluid-structure interaction (FSI) model in COMSOL Multiphysics is used to solve (1), assuming that the nodes of the LV boundary move along the outward normal direction with the average speed at each time step. The corresponding velocity 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 with constant pressure are applied at the inlet valve , while zero-velocity conditions are applied at the outlet valve : During the systole phase, the boundary conditions in (31) are imposed in the opposite way: In this study, we set  kg/m3 and  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].

3. Results and Discussion

3.1. 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 vorticity fields on planes A and B are shown in Figures 5(a)5(d) and Figures 5(e)5(h), respectively. The velocity fields are illustrated with velocity vectors, while the colored contours in Figures 5(a)5(h) represent the out-of-plane vorticity fields. Here, warm (such as red) and cool (such as blue) colors correspond to clockwise and counterclockwise rotating vortex patterns, respectively. In addition, Figures 5(i)5(l) indicate the time-varying LV volumes marked by a red dot at selected time frames in the cardiac cycle. Since we only consider ten samples from echocardiographic images over one cardiac cycle, the volume curve does not capture early and atrial waves. At the early stage of the diastole phase, Figures 5(a) and 5(e) indicate that two counterclockwise rotating vortices formed near the inlet. Later, the two vortices move toward the apex while maintaining their shape, and the left vortex starts to interact with the lateral LV wall, as shown in Figures 5(b) and 5(f). Due to the interaction, the left clockwise rotating vortex disappears in Figures 5(c) and 5(g), while the counterclockwise rotating vortex occupies the area near the apex and leads to making a larger swirling pattern. In the systole phase, the counterclockwise rotating vortex moves to the right LV wall and the corresponding blood flow heads toward the outlet valve, as shown in Figures 5(d) and 5(h).

3.2. 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 from to in order to optimize the regularization parameter in (15). The reconstruction error is determined based on -norm error between the reconstructed velocity fields and the forward data : where is the number of time steps. Figure 6 indicates the minimum error at , which is used to reconstruct the blood flow in LV.

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 is obtained by solving (25) at each time step. The reconstructed solution is obtained by repeatedly updating the minimizer in (17) until the stopping criterion (18) is satisfied. Figure 7 illustrates the performance of the proposed reconstruction model. During the diastole phase in Figures 7(i)7(k), the present reconstruction model clearly captures the counterclockwise rotating flow patterns as well as incoming velocity vector fields heading to the apex, as shown in Figures 7(a)7(c) and 7(e)7(g). Moreover, Figures 7(d) and 7(h) clearly show the movement of the counterclockwise rotating vortex to the right LV wall for the systole phase in 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 -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.

4. 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.

Competing Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

The authors thank Siemens Korea R&D center for ultrasound data support related to our work. This work was supported by National Research Foundation (NRF) of Korea grants funded by the Korean government (MSIP) (NRF-20151009350).