Computational and Mathematical Methods in Medicine
Volume 2016 (2016), Article ID 6371078, 10 pages
http://dx.doi.org/10.1155/2016/6371078

Research Article

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

Jaeseong Jang,¹ Chi Young Ahn,² Jung-Il Choi,¹ and Jin Keun Seo¹

¹Department of Computational Science and Engineering, Yonsei University, 50 Yonsei-ro, Seodaemun-gu, Seoul 03722, Republic of Korea
²Division of Integrated Mathematics, National Institute for Mathematical Sciences, 70 Yuseong-daero 1689 beon-gil, Yuseong-gu, Daejeon 34047, Republic of Korea

Received 27 February 2016; Accepted 28 April 2016

Academic Editor: Yuhai Zhao

Copyright © 2016 Jaeseong Jang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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 [1–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-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.

Figure 1: Overall framework for LV blood flow reconstructions based on echocardiographic images.

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.

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.

In this paper, we use the color Doppler data , as the additional information of one component of velocity, of the form

where 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) yields

We 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 3: Multiple views of 3D cardiac volume images. The left column shows coronal, sagittal, and axial views in order, and right column represents an image for 3D volume rendering.

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.

Figure 4: LV model for numerical experiments. The model was reconstructed from ultrasound images.

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/m³ 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).

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

Figure 5: Synthetic intraventricular velocity fields and out-of-plane vorticity component on imaging planes A and B in the first and second rows, respectively, along with volume curves in the third row. On each volume curve, a red dot indicates the corresponding phase for the velocity fields and vorticity contours.

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.

Figure 6: Error change depending on the regularization parameter .

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.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

Figure 7: Reconstructed intraventricular velocity and out-of-plane vorticity fields on imaging planes A and B in the first and second rows, respectively, along with volume curves in the third row. On each volume curve, a red dot indicates the corresponding phase for the velocity and vorticity fields.

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.

Table 1: Comparison of -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 , respectively. The third and fourth columns represent normalized -norm in global energy estimates of velocity and vorticity, respectively. Here, , , ,  , and . Note that the superscripts and denote the forward and reconstructed data, respectively.

0.05
16.7%
62.4%
0.4%
4.1%

0.15
11.8%
20.4%
0.2%
3.0%

0.25
18.7%
29.3%
2.9%
6.2%

0.35
15.9%
26.0%
3.1%
4.8%

0.45
16.8%
27.6%
0.3%
4.5%

0.55
20.9%
30.9%
1.4%
4.9%

0.65
26.8%
38.6%
2.0%
6.7%

0.75
29.9%
45.1%
5.0%
6.9%

0.85
32.2%
50.4%
0.4%
4.1%

0.95
32.2%
48.6%
0.4%
4.1%

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

References

G.-R. Hong, G. Pedrizzetti, G. Tonti et al., “Characterization and quantification of vortex flow in the human left ventricle by contrast echocardiography using vector particle image velocimetry,” JACC: Cardiovascular Imaging, vol. 1, no. 6, pp. 705–717, 2008.

S. Cimino, G. Pedrizzetti, G. Tonti et al., “In vivo analysis of intraventricular fluid dynamics in healthy hearts,” European Journal of Mechanics, B/Fluids, vol. 35, pp. 40–46, 2012.

K. Lampropoulos, W. Budts, A. Van de Bruaene, E. Troost, and J. P. Van Melle, “Visualization of the intracavitary blood flow in systemic ventricles of Fontan patients by contrast echocardiography using particle image velocimetry,” Cardiovascular Ultrasound, vol. 10, article 5, 2012.

H. Gao, P. Claus, M.-S. Amzulescu, I. Stankovic, J. D'Hooge, and J.-U. Voigt, “How to optimize intracardiac blood flow tracking by echocardiographic particle image velocimetry? Exploring the influence of data acquisition using computer-generated data sets,” European Heart Journal Cardiovascular Imaging, vol. 13, no. 6, pp. 490–499, 2012.

H. Gao, B. Heyde, and J. D'hooge, “3D intra-cardiac flow estimation using speckle tracking: a feasibility study in synthetic ultrasound data,” in Proceedings of the IEEE International Ultrasonics Symposium (IUS '13), pp. 68–71, Prague, Czech Republic, July 2013.

D. H. Evans and W. N. McDicken, Doppler Ultrasound: Physics, Instrumentation, and Signal Processing, John Wiley & Sons, Chichester, UK, 2000.

D. Garcia, J. C. del Álamo, D. Tanné et al., “Two-dimensional intraventricular flow mapping by digital processing conventional color-doppler echocardiography images,” IEEE Transactions on Medical Imaging, vol. 29, no. 10, pp. 1701–1713, 2010.

S. Ohtsuki and M. Tanaka, “The flow velocity distribution from the doppler information on a plane in three-Dimensional flow,” Journal of Visualization, vol. 9, no. 1, pp. 69–82, 2006.

M. Arigovindan, M. Sühling, C. Jansen, P. Hunziker, and M. Unser, “Full motion and flow field recovery from echo Doppler data,” IEEE Transactions on Medical Imaging, vol. 26, no. 1, pp. 31–45, 2007.

A. Gomez, K. Pushparajah, J. M. Simpson, D. Giese, T. Schaeffter, and G. Penney, “A sensitivity analysis on 3D velocity reconstruction from multiple registered echo Doppler views,” Medical Image Analysis, vol. 17, no. 6, pp. 616–631, 2013.

J. Jang, C. Y. Ahn, K. Jeon et al., “A reconstruction method of blood flow velocity in left ventricle using color flow ultrasound,” Computational and Mathematical Methods in Medicine, vol. 2015, Article ID 108274, 14 pages, 2015.

K. E. Leung and J. G. Bosch, “Automated border detection in three-dimensional echocardiography: principles and promises,” European Heart Journal-Cardiovascular Imaging, vol. 11, no. 2, pp. 97–108, 2010.

C. Y. Ahn, “Robust myocardial motion tracking for echocardiography: variational framework integrating local-to-global deformation,” Computational and Mathematical Methods in Medicine, vol. 2013, Article ID 974027, 14 pages, 2013.

Y. Cheng, H. Oertel, and T. Schenkel, “Fluid-structure coupled CFD simulation of the left ventricular flow during filling phase,” Annals of Biomedical Engineering, vol. 33, no. 5, pp. 567–576, 2005.

X. Zheng, J. H. Seo, V. Vedula, T. Abraham, and R. Mittal, “Computational modeling and analysis of intracardiac flows in simple models of the left ventricle,” European Journal of Mechanics B: Fluids, vol. 35, pp. 31–39, 2012.

J. H. Seo, V. Vedula, T. Abraham et al., “Effect of the mitral valve on diastolic flow patterns,” Physics of Fluids, vol. 26, no. 12, Article ID 12190, 2014.