Med Biol Eng Comput DOI 10.1007/s11517-014-1193-3
ORIGINAL ARTICLE
Development and feasibility study of a two‑dimensional ultrasonic‑measurement‑integrated blood flow analysis system for hemodynamics in carotid arteries Takaumi Kato · Kenichi Funamoto · Toshiyuki Hayase · Shusaku Sone · Hiroko Kadowaki · Tadashi Shimazaki · Takao Jibiki · Koji Miyama · Lei Liu
Received: 25 April 2014 / Accepted: 25 August 2014 © International Federation for Medical and Biological Engineering 2014
Abstract Prevention and early detection of atherosclerosis are critical for protection against subsequent circulatory disease. In this study, an automated two-dimensional ultrasonic-measurement-integrated (2D-UMI) blood flow analysis system for clinical diagnosis was developed, and the feasibility of the system for hemodynamic analysis in a carotid artery was revealed. The system automatically generated a 2D computational domain based on ultrasound color Doppler imaging and performed a UMI simulation of blood flow field to visualize hemodynamics in the domain. In the UMI simulation, compensation of errors was applied by adding feedback signals proportional to the differences between Doppler velocities by measurement and computation while automatically estimating the cross-sectional average inflow velocity. The necessity of adjustment of the feedback gain was examined by analyzing blood flow in five carotid arteries: three healthy, one sclerosed, and one stenosed. The same feedback gain was generally applicable for the 2D-UMI simulation in all carotid arteries, depending on target variables. Thus, the present system was shown to be versatile in the sense that the parameter
T. Kato · S. Sone Graduate School of Biomedical Engineering, Tohoku University, Sendai 980‑8579, Japan K. Funamoto (*) · T. Hayase Institute of Fluid Science, Tohoku University, 2‑1‑1 Katahira, Aoba‑ku, Sendai 980‑8577, Japan e-mail: [email protected] H. Kadowaki Graduate School of Engineering, Tohoku University, Sendai 980‑8579, Japan T. Shimazaki · T. Jibiki · K. Miyama · L. Liu GE Healthcare Japan, Hino 191‑8503, Japan
is patient independent. Moreover, the possibility of a new diagnostic method based on the hemodynamic information obtained by the 2D-UMI simulation, such as a waveform of the cross-sectional average inflow velocity and wall shear stress distributions, was suggested. Keywords Hemodynamics · Ultrasound color Doppler imaging · Computational fluid dynamics · Measurementintegrated simulation · Carotid artery
1 Introduction Circulatory disease is a major cause of death, its mortality rate being comparable with that of cancer. Since atherosclerosis is a key factor leading to the disease, much attention has been paid to its prevention and early detection [21]. It is well known that hemodynamic stresses, such as wall pressure and wall shear stress (WSS) [20], are closely related to the development and progression of atherosclerosis and subsequent circulatory disease. In order to utilize hemodynamic information for the development of a therapeutic strategy, various metrics concerning blood flow and hemodynamic stresses have been proposed [13] and their thresholds for pathological events have been discussed. Thus, acquisition of accurate, detailed information on hemodynamics is critical for advanced diagnosis of circulatory disease, but none of the present methodologies can completely provide such information. For instance, ultrasound Doppler measurement provides only a one-dimensional velocity component (Doppler velocity) of blood flow, which is a projection of the three-dimensional (3D) velocity vector in the ultrasound beam direction, with a compact system noninvasively [18, 25], but hemodynamic stresses cannot be measured. Phase-contrast magnetic resonance imaging (PC
13
MRI) may be used to obtain information on 3D velocity vectors of blood flow, as well as blood vessel configuration, enabling calculation of WSS distribution [9]. However, the temporal and spatial resolutions of PC MRI are limited, and its accuracy is insufficient to quantitatively evaluate hemodynamics. On the other hand, a numerical simulation of blood flow in a patient-specific vessel geometry reconstructed from medical images reveals the blood flow field and hemodynamic stresses in detail [1, 23]. However, it is difficult to correctly specify upstream boundary conditions and physiological parameters of the blood and the vessel wall [11, 15, 16, 26], and therefore, the computed result may differ from the real in vivo hemodynamics. In addition, validation of the computational result is difficult. As a solution for individual limitations in measurement and computation, the authors have previously proposed ultrasonic-measurement-integrated (UMI) simulation, in which artificial body forces proportional to the differences between measured and computed Doppler velocities are fed back to the numerical simulation, and the computational result converges to the real blood flow field [4]. Fundamental knowledge of the effects of feedback was obtained using two-dimensional UMI (2D-UMI) simulation of blood flow in a thoracic aneurysm [5]. In addition, efficiency and other characteristics of three-dimensional UMI (3D-UMI) simulation of blood flow were also investigated by numerical experiments [6, 7]. In spite of recent improvements in computational hardware, 3D simulation of blood flow still requires a significant amount of computational time, and therefore, a 3D-UMI blood flow analysis system is not realistic for clinical use. In contrast, a 2D-UMI blood flow analysis system could be clinically useful since it is convenient to use and provides results in a shorter time. However, such a 2D-UMI simulation system for clinical diagnosis is not yet available. Moreover, in former studies of UMI simulations, blood vessel morphology was extracted manually and the blood flow volume was assumed based on references, resulting in consideration under a semi-patient-specific condition. In this study, an automated 2D-UMI blood flow analysis system for clinical diagnosis was developed, and the feasibility of the system for hemodynamic analysis in a carotid artery was examined. The system automatically generated a 2D computational domain based on ultrasound color Doppler imaging and performed a UMI simulation of blood flow field to visualize hemodynamics in the domain. In the UMI simulation, compensation of errors was applied by adding feedback signals proportional to the differences between Doppler velocities by measurement and computation while automatically estimating the cross-sectional average inflow velocity. If it is necessary to adjust the feedback gain, which is a variable parameter significantly influencing the computational accuracy of the UMI simulation, the 2D-UMI blood flow analysis system
13
Med Biol Eng Comput
lacks the versatility for hemodynamic analysis. Hence, the versatility of the system for hemodynamic analysis in the carotid artery was examined by investigating effects of the feedback gain on the analysis results of five carotid arteries, including three healthy arteries, one sclerosed artery, and one stenosed artery. Then, the possibility of a new diagnostic method based on the hemodynamic information obtained by the 2D-UMI simulation was also investigated.
2 Methods 2.1 System configuration A block diagram and a snapshot of the 2D-UMI blood flow analysis system and the definition of the coordinate systems are shown in Fig. 1. The system is equipped with the following functions required for clinical diagnosis: (1) an automatic generation function of the 2D computational domain based on ultrasound color Doppler imaging, (2) a UMI simulation function for analysis of the blood flow field with compensation of errors by application of feedback signals proportional to the differences between measured and computed Doppler velocities, (3) an automatic estimation function of the cross-sectional average inflow velocity, and (4) a function for visualization of analysis results. In the system, ultrasound color Doppler imaging of hemodynamics is first carried out with conventional ultrasound diagnostic imaging equipment. The ultrasound color Doppler data are imported into a computer in DICOM format, and a longitudinal 2D blood vessel shape is automatically extracted from the data. A 2D-UMI simulation of blood flow is then performed by applying feedback signals, while the cross-sectional average inflow velocity is adjusted to reduce the errors between measured and computed Doppler velocities in the feedback domain. After a convergent result is obtained, the adverse effect of the feedback on the pressure field in the UMI simulation is compensated and the blood flow field (velocity vector and pressure distribution) and hemodynamic parameters such as WSS are visualized. In this study, ultrasound diagnostic imaging equipment (LOGIQ 7, GE Healthcare Japan, Tokyo, Japan) equipped with a linear ultrasound probe (12L, Fc = 8.2 MHz, GE Healthcare Japan, Tokyo, Japan) was employed for ultrasound color Doppler imaging. To perform the UMI simulation, a custom-built export function of raw data of ultrasound color Doppler was installed in the ultrasonic equipment. The computations were performed using one central processing unit (CPU) at 1.60 GHz in a server (Altix 3700 Bx2, SGI Japan, Tokyo, Japan). It should be noted that the computations can be performed with a general-purpose personal computer or a built-in computer in ultrasound diagnostic imaging equipment.
Med Biol Eng Comput
(a)
(b)
(c)
Fig. 1 The 2D-UMI blood flow analysis system: a schematic diagram, b overview of the system, and c coordinate systems in the system. Blood flow in a carotid artery of subject ID #1 is visualized in (b)
2.2 Analysis methods In spite of the three-dimensionality of blood flow, the 2D-UMI blood flow analysis system dealt with a cross section of a blood flow field along the axis of a blood vessel as a convenient clinical system. Deformation of the blood vessel was ignored for simplicity. 2.2.1 Extraction of blood vessel shape The shape of a blood vessel was extracted based on color data of Doppler velocity and on gray-scale data of tissue morphology, both of which were simultaneously provided in ultrasound color Doppler imaging. In a color Doppler image, the color indication of blood flow often extended through the deep-side blood vessel wall (the lower side wall in the ultrasound image) due to dispersion of ultrasound. Locations with nonzero time-averaged Doppler velocities over a measurement period tended to show a blood vessel lumen larger than the actual one. Hence, it was necessary to modify the blood vessel shape with reference to the gray-scale tissue information. The time-averaged gray-scale data in the measurement period were binarized after enhancement of the blood vessel wall with a differentiation filter in the ultrasound beam direction (vertical to the blood vessel). The blood vessel lumen was then determined by means of a logical product of both
binarized data of time-averaged color and gray-scale information. The coordinate system was defined as the Y and X directions toward the normal ultrasound beam direction and the head, respectively, as shown in Fig. 1c. To facilitate hemodynamic analysis, the xy coordinate system was used in the computation, rotating the XY coordinate system at an angle of β (counterclockwise is positive) to align the centerline of the extracted blood vessel shape with the x-direction. The ultrasound beam angle, θ, had a positive value clockwise against the Y-direction, and the net angle of the ultrasound beam, α, was expressed as the summation of θ and β. 2.2.2 2D‑UMI simulation The outline of the previously reported UMI simulation is described here [3]. The governing equations of the 2D-UMI blood flow simulation are the 2D Navier–Stokes equations for incompressible, viscous fluid flow, and the pressure equation as follows: ∂u + (u · ∇)u = µ�u − ∇p + f, ρ (1) ∂t
�p = −∇ · ρ(u · ∇)u + ∇ · f,
(2)
where u = (u, v) is the velocity vector, p is the pressure, t is the time, ρ is the density, μ is the viscosity, and f represents
13
Med Biol Eng Comput
a feedback signal to make the computational result approach the measurement data (real flow). The feedback signals, f, are applied as artificial body forces defined by the differences between Doppler velocities Vm and Vc on a cross section, which are measured by ultrasound color Doppler imaging and evaluated based on computed velocity vectors:
f=
Vc −Kv∗
− Vm U
ρU 2 L
· b,
(3)
where K*v is the feedback gain (nondimensional), U is the characteristic velocity, L is the characteristic length, and b is the unit vector along the ultrasound beam. Note that the special case with K*v = 0 is an ordinary simulation without feedback. The application of the feedback signals makes the velocity field converge to the real blood flow, satisfying the equation of continuity. However, it may adversely affect the pressure field in the case that the feedback signal does not satisfy the divergence-free condition in Eq. (2) [3]. Hence, the pressure field is modified by adding a compensation term calculated with the irrotational part of the feedback signal after a convergent result of the pressure field is obtained. The governing equations explained above were discretized and solved by a method similar to the SIMPLER method [17]. The convective terms were discretized with the reformulated QUICK scheme [10], and the time derivative terms were discretized with the first Euler implicit scheme. Linear algebraic equations were solved using a modified strongly implicit (MSI) scheme [22]. Feedback signals in the UMI simulation were added at the computational grid points in the feedback domain defined within the computational domain. Treatments were applied to detect the major measurement errors, such as aliasing, wall filter, and lack of data, and to reset the feedback signal to zero at the feedback point [8]. For evaluation of the computational accuracy, the space-averaged error of the Doppler velocity, e(t), in a monitoring domain Ωe was defined as follows:
e(t) =
1 |Vc − Vm | , NΩe uref
(4)
Ωe
where uref is the representative velocity for normalization and NΩe is the total number of the grid points in the monitoring domain Ωe. In addition, the time–space-averaged error, eave, was calculated by averaging the space-averaged error, e(t), in four cardiac cycles after a convergent solution was obtained. 2.2.3 Automatic estimation of the cross‑sectional average inflow velocity The exact boundary conditions are usually unknown for UMI blood flow simulation. A correct setting of the
13
cross-sectional average inflow velocity, uin(t), minimizes the error between computation and measurement [11, 26]. Consequently, the cross-sectional average inflow velocity uin(t) for a numerical simulation was estimated by means of the golden section search [19], which minimized the summation of the errors of the computed Doppler velocities Vc from the measured ones Vm in a monitoring domain Ωq defined in a computational domain: |Vc − Vm |. min uin (5) Ωq
A parabolic velocity profile with the estimated cross-sectional average inflow velocity was applied at the upstream boundary. In addition, free-flow condition (du/dx = 0 and dv/dy = 0) was applied at the downstream boundary, and the no-slip condition (u = 0) was assumed on the vessel wall. 2.2.4 Visualization of analysis results Application software for visualization of analysis results on Windows OS was programmed with Visual C++ by using Microsoft Visual Studio (Microsoft Japan, Tokyo, Japan). Velocity vectors at arbitrary locations obtained by the 2D-UMI simulation were visualized on a color Doppler image, in which Doppler velocities were contoured with warm or cold colors according to their positive or negative values, respectively. In addition, for clinical use, the above-mentioned visualized results were superimposed on the corresponding ultrasound B-mode image, and WSS was indicated with a color display on bands aligned with the blood vessel walls. 2.3 Analysis conditions In order to confirm the versatility of the system for hemodynamic analysis in the carotid artery, objectives were blood flow in five carotid arteries, including three healthy arteries, one sclerosed artery, and one stenosed artery, in five patients as summarized in Table 1. The study was approved by the local ethics committee, and informed consent was obtained from all patients. Clinical ultrasound color Doppler imaging was conducted following the ethics guidelines of GE Healthcare Japan Corporation. In the measurement, a longitudinal cross section was measured along the blood vessel axis as much as possible so as not to cause a large difference of flow volume among the axial positions. Ninety-one ultrasonic scan lines were acquired with an ultrasound beam angle θ of 20° or −20° from the ultrasound probe, and measurement data at 104 sampling points on each ultrasonic scan line were measured. The lateral and longitudinal pixel spacing was 303 and 182 μm, respectively. The angle of blood vessel β, net beam angle α,
Med Biol Eng Comput Table 1 Information on carotid arteries, ultrasonic measurement, and computation
ID
1
2
3
4
5
Age Sex Position Condition
53 Male Right Healthy 5.72
56 Male Left Healthy 5.72
77 Male Right Healthy 5.99
76 Female Right Sclerosed 7.16
71 Male Right Stenosed 5.68
-20 -5.2 -25.2
20 13.8 33.8
20 9.8 29.8
-20 -3.3 -23.3
20 11.8 31.8
Δt (ms) Nx × Ny
5.0 5.4 83.6 69 × 104
5.0 5.4 83.6 68 × 104
5.0 4.4 59.2 68 × 104
5.0 4.4 59.2 69 × 104
5.0 4.4 59.2 68 × 104
Δx, Δy (μm)
277, 164
227, 130
251, 146
287, 171
240, 138
Inlet diameter, D (mm) Beam angle θ (°) Angle of blood vessel β (°) Net beam angle α (= θ + β) (°) Working frequency (MHz) PRF (kHz)
working frequency, pulse repetition frequency (PRF), and temporal resolution Δt in the ultrasound imaging are summarized in Table 1. As for analysis conditions in the 2D-UMI simulation, the density ρ and viscosity μ of the blood were assumed to be 1.00 × 103 kg/m3 and 4.0 × 10−3 Pa s, respectively, within a normal range. The spatial resolution (Δx, Δy) of the computational grid was set to be that of the original ultrasound color Doppler imaging. The computational grid points are summarized in Table 1. The computational time increment was set to be identical to the temporal resolution Δt of the ultrasound imaging. The domain between 1/8 and 7/8 of the computational domain from the upstream boundary was defined as the feedback domain [5], and feedback signals, f, were added at all the grid points in the domain. In the estimation of the cross-sectional average inflow velocity, the monitoring domain Ωq was set so as to be identical to the feedback domain, and the value of uin which satisfied Eq. (5) was searched for in an adequately wide range between −0.1 m/s and 0.5 m/s. All variables were non-dimensionalized with characteristic values. The inlet diameter, D, was used as the characteristic length, L, an expedient value of 0.1 m/s was used as the characteristic velocity, U, and the density of blood, ρ, was used as the other characteristic value. The representative time-averaged velocity in a common carotid artery of 0.39 m/s was used as uref in the evaluation of the computational accuracy of Eq. (4) [24]. From here on, the same symbols are used for both dimensional and non-dimensional values for the sake of simplicity. The convergence criteria were set as a compromise between computational time and accuracy: 1.00 × 10−2 (nondimensional) for the residuals in the blood flow analysis, and 1.00 × 10−3 (1.00 × 10−4 m/s) for the estimation of the cross-sectional average inflow velocity. In order to confirm the versatility of the system for hemodynamic analysis in the carotid artery, the necessity
of adjustment of the feedback gain was examined. For the five carotid arteries, including healthy and diseased ones, the 2D-UMI simulation was performed by changing the feedback gain from 0 (ordinary simulation) to 100 in increments of 10, and was also performed with the feedback gains of 200, 300, 400, 500, and 1,000. Then, effects of the feedback gain on analysis results were investigated.
3 Results An example of the display visualizing an analysis result obtained by the present 2D-UMI blood flow analysis system is shown in Fig. 1b. A contour of Doppler velocities and velocity vectors are displayed, superimposed on the ultrasound B-mode image, and WSS distributions are visualized on bands along the blood vessels. A movie can be displayed to observe the variation of blood flow and hemodynamic stresses according to the heartbeat. Moreover, it is possible to compare the computational results with the measurement data obtained by the ultrasound diagnostic imaging equipment on the left-hand side. In order to confirm the versatility of the system, hemodynamic analysis in five carotid arteries, including three healthy arteries, one sclerosed artery, and one stenosed artery, was performed with various feedback gains. The measured color Doppler images at the peak flow of each subject (ID #1–ID #5) are shown in the left column of Fig. 2. Blood flows from right to left in all images, but Doppler velocities are colored with warm or cold color based on differences of the ultrasound beam angle θ (see Table 1). Subjects of ID #1 through ID #3 (Fig. 2a–c) are healthy carotid arteries; the bulge at the downstream of subject ID #2 (Fig. 2b) is the carotid bifurcation, and subject ID #3 (Fig. 2c) is a region extending from a carotid bifurcation into the internal carotid artery. On the other
13
Fig. 2 Measured color Doppler image at the peak flow (the left column) and computational results of the corresponding color Doppler image with velocity vectors and the Doppler velocity error distribution by the ordinary simulation (K*v = 0, the middle column) and the UMI simulation (K*v = 500, the right column) of each subject: a ID #1 (t = 0.92 s), b ID #2 (t = 2.09 s), c ID #3 (t = 3.55 s), d ID #4 (t = 3.02 s), and e ID #5 (t = 3.44 s). Blood flows from the right side to the left side. The computational domain of each subject is indicated with a two-way arrow in the measurement image
Med Biol Eng Comput
(a)
(b)
(c)
(d)
(e)
hand, subject ID #4 (Fig. 2d) is a common carotid artery with moderate atherosclerosis, the color Doppler image of which contains aliasing as indicated by an arrow. Subject ID #5 is a region extending from a carotid bifurcation into
13
the internal carotid artery with a stenosis. The color Doppler image at the stenosis shows a mosaic pattern (see a dotted circle region) since the hemodynamics are complicated by fast blood flow velocity.
Med Biol Eng Comput Fig. 3 Variations of the space-averaged error of Doppler velocity in the feedback domain, e, in the ordinary simulation (K*v = 0) and the UMI simulations (K*v = 200 and 500) of each subject: a ID #1, b ID #2, c ID #3, d ID #4, and e ID #5
(a)
(b)
(c)
(d)
(e)
The velocity vectors and color Doppler images obtained by the ordinary simulation (K*v = 0) and the UMI simulation with K*v = 500 at the same time step are represented in the upper middle and upper right, respectively, of each group in Fig. 2. The computational domains are indicated in the corresponding measured color Doppler images on the left. In the ordinary simulations in the middle column, the blood flow maintains the parabolic velocity profile applied at the upstream boundary since the carotid arteries are almost straight except for the stenosed artery in Fig. 2e. In contrast, the UMI simulation in the right column shows complicated velocity profiles in the feedback domain (the domain between 1/8 and 7/8 of the computational domain) even in the healthy carotid arteries. The errors of Doppler velocities between the measurement and each computation are shown in the lower monochrome images in Fig. 2, in which black indicates small errors. Error against the measurement data exists in the whole computational domain in the ordinary simulation, whereas it is suppressed in the feedback domain in the UMI simulations. Somewhat large errors are observed outside the feedback domain and at the locations of aliased data (see the arrow and the dotted circle in Fig. 2d, e).
Periodic variations of the space-averaged error of Doppler velocity in the feedback domain, e(t), the estimated cross-sectional average inflow velocity, uin(t), and the spaceaveraged WSS, τsave(t), in four cardiac cycles are compared between the ordinary simulation (K*v = 0) and the UMI simulations with K*v = 200 and 500 in Figs. 3, 4 and 5. Note that the cross-sectional average inflow velocity was successfully obtained after 21-time iterative computations of the golden section method at each time step, and the value of τsave was calculated by averaging WSS acting on both the shallow-side and deep-side (the upper and lower sides in an ultrasound image) walls of the blood vessel. The spaceaveraged error of Doppler velocity becomes small as the feedback gain increases in all subjects, implying improvement of the computational accuracy (Fig. 3). On the other hand, the estimated cross-sectional average inflow velocity is less sensitive to the feedback gain except for the subject with a stenosis of ID #5 (Fig. 4). The space-averaged WSS increases as the feedback gain increases (Fig. 5). Variations of the time–space-averaged error of Doppler velocity in the feedback domain, eave, the maximum value of the cross-sectional average inflow velocity, umax, the
13
Fig. 4 Variations of the cross-sectional average inflow velocity, uin, by the ordinary simulation (K*v = 0) and the UMI simulations (K*v = 200 and 500) of each subject: a ID #1, b ID #2, c ID #3, d ID #4, and e ID #5
Med Biol Eng Comput
(a)
(b)
(c)
(d)
(e)
time–space-averaged WSS, τtsave, and the average computational time for one time step, T, with the feedback gain are shown in Fig. 6, being normalized by the corresponding values of the ordinary simulation (K*v = 0). The values used for the normalization are summarized in Table 2. For each subject, the normalized time–space-averaged error of the Doppler velocity, eave/eave0, sharply decreases to half of that of the ordinary simulation in K*v ≤ 200 and converges with increasing feedback gain (Fig. 6a). The maximum value of the cross-sectional average inflow velocity, umax/umax0, shows little variation in K*v ≥ 500 (Fig. 6b). On the other hand, the time–space-averaged WSS, τtsave/τtsave0, increases with feedback gain, and the change of τtsave/τtsave0 becomes small in K*v ≥ 500 (Fig. 6c). The average computational time for one time step, T/T0, almost linearly increases with feedback gain for all subjects (Fig. 6d).
4 Discussion An automated 2D-UMI blood flow analysis system for clinical diagnosis was successfully developed. The system was
13
equipped with the following four functions: (1) an automatic generation function of the 2D computational domain based on ultrasound color Doppler imaging, (2) a UMI simulation function for analysis of the blood flow field with compensation of errors by application of feedback signals proportional to the differences between Doppler velocities by measurement and computation, (3) an automatic estimation function of the cross-sectional average inflow velocity, and (4) a function for visualization of analysis results. The versatility of the system for hemodynamic analysis in the carotid artery was proved, showing more accurate results than the ordinary simulation for blood flow in all carotid arteries treated in this study. Moreover, the possibility of a new diagnostic method based on the hemodynamic information obtained by the 2D-UMI simulation, such as a waveform of the cross-sectional average inflow velocity and WSS distributions, was also suggested. Effects of the feedback gain on parameters obtained by the present system are different. Although an increase of the feedback gain decreases the error of Doppler velocity against the measurement (Fig. 6a) as well as increasing the WSS (Fig. 6c), it hardly influences the estimation of
Med Biol Eng Comput Fig. 5 Variations of the spaceaveraged WSS, τsave, by the ordinary simulation (K*v = 0) and the UMI simulations (K*v = 200 and 500) of each subject: a ID #1, b ID #2, c ID #3, d ID #4, and e ID #5
(a)
(b)
(c)
(d)
(e)
Fig. 6 Variations of parameters with the feedback gain in the UMI simulation: a the time– space-averaged error of Doppler velocity in the feedback domain, eave, b the maximum value of the cross-sectional average inflow velocity, umax, c the time–space-averaged WSS, τtsave, in the four cardiac cycle, and d the average computational time, T, required for one time step. Each result was normalized with the corresponding result by the ordinary simulation (K*v = 0) in Table 2
(a)
(b)
(c)
(d)
13
Med Biol Eng Comput
Table 2 Parameters in the ordinary simulation (K*v = 0): the time– space-averaged error of Doppler velocity in the feedback domain, eave0, the maximum value of the cross-sectional average inflow velocID
1
2
3
ity, umax0, the time–space-averaged WSS, τtsave0, in the four cardiac cycle, and the average computational time, T0, required for one time step 4
5
Ave.
SD
eave0 (–)
0.075
0.101
0.065
0.064
0.101
0.081
0.019
umax0 (m/s) τtsave0 (Pa)
0.405
0.338
0.183
0.358
0.229
0.303
0.093
0.339
0.279
0.199
0.495
1.063
0.475
0.346
T0 (s)
13.4
12.7
10.8
12.2
20.2
13.9
3.7
the cross-sectional average inflow velocity (Fig. 6b). WSS is sensitive to slight changes in velocity since it is calculated from the spatial gradient of the velocity. Thereby, the variations of the error of Doppler velocity and the WSS become small with a feedback gain of more than about 200 and 500, respectively. The variations of the parameters are independent of the net beam angle of the ultrasound (see Table 1). For clinical use of the 2D-UMI blood flow analysis system, it is preferable to obtain an analysis result in a shorter time. The feedback gain should not be increased without restriction since the computational time almost linearly increases (Fig. 6d). Comprehensively considering the effects of the feedback gain on the analysis parameters and the computational time, the UMI simulation with a feedback gain between 200 and 500 can provide an appropriate result with better computational accuracy than the ordinary simulation. Comparison of the results of hemodynamic analysis with the same feedback gain between healthy and diseased carotid arteries shows a difference of feedback effects. The application of feedback in healthy carotid arteries reduces the error of Doppler velocity more than such application in diseased ones (Fig. 6a). The reduction rate of the time– space-averaged error of Doppler velocity by the feedback is smaller in the diseased carotid arteries than that in the healthy ones, though it becomes less than 50 % in all cases. This phenomenon could be derived from the complicated blood flow field in the diseased carotid artery and might be useful for diagnosis of carotid arteries: Diseased carotid arteries could be detected based on the reduction rate of the error. As observed in Fig. 2, the diseased carotid arteries have more complicated velocity distributions than the healthy ones, resulting in difficulty of reproduction of the hemodynamics. Moreover, the values of WSS in healthy carotid arteries by the ordinary simulation (K*v = 0) fluctuate close to zero, but those by the UMI simulations increase and show remarkable variations (Fig. 5a–c). In contrast, the values of WSS in diseased carotid arteries show fluctuations with maximum values of several Pascal, and the differences of the results between the ordinary simulation and the UMI simulations are not large (Fig. 5d, e). This reflects changes in the velocity field of blood flow near the
13
vessel wall owing to the feedback. For instance, the crosssectional average inflow velocity in a healthy carotid artery steeply changes in the acceleration and deceleration phases (Fig. 4a), and the velocity profile has a rectangular shape (Fig. 2a). However, the cross-sectional average inflow velocity in a carotid artery with atherosclerosis gradually decreases in the deceleration phase (Fig. 4d), and the velocity profile is approximately parabolic (Fig. 2d). Thus, the differences of variation and profile of blood flow velocity depending on the conditions of blood vessels cause differences in WSS distribution. This implies that the present system has the potential to detect atherosclerosis. It is expected that analysis by the 2D-UMI blood flow analysis system for more patients would clarify the relationships between circulatory diseases, such as atherosclerosis and stenosis, and hemodynamics, and help establish criteria to diagnose them. However, several tasks remain to be considered. Variability of the results obtained by different medical personnel, and usefulness of the system for follow-up diagnosis should be investigated. More importantly, a real blood flow has three-dimensionality [2, 12, 14, 25], and therefore, it is necessary to evaluate the effect of the three-dimensionality of the blood flow on the result of 2D-UMI simulation. Moreover, the equidistant orthogonal grid used in the system introduces inaccuracies in the reproducibility of the blood vessel configuration. Usage of a boundary-fitted or an unstructured grid system merits future consideration. Furthermore, neglect of wall motion influences the accuracy of the analysis result, especially near the wall, and thus, a moving boundary problem in the 2D-UMI blood flow analysis system should also be considered in future work.
5 Conclusion In this paper, a 2D-UMI blood flow analysis system was developed to easily obtain hemodynamic information. In order to confirm the versatility of the system for hemodynamic analysis in the carotid artery, the necessity of adjustment of the feedback gain of a variable parameter was examined by investigating effects of the feedback gain on
Med Biol Eng Comput
analysis results for five carotid arteries: three healthy, one sclerosed, and one stenosed. It was shown that the same feedback gain was generally applicable for the 2D-UMI simulation in all carotid arteries, yielding more accurate results with increasing feedback gain than those by the ordinary simulation without feedback. Thus, the present system is versatile in the sense that the parameter is patient independent. Moreover, the possibility of a new diagnostic method based on the hemodynamic information obtained by the 2D-UMI simulation, such as a waveform of the cross-sectional average inflow velocity and WSS distributions, was suggested. Acknowledgments Part of this study was supported by Grantin-Aid for Scientific Research (B) (24360064). All computations were performed using the supercomputer system at the Advanced Fluid Information (AFI) Research Center, Institute of Fluid Science, Tohoku University.
References 1. Filipovic N, Teng ZZ, Radovic M, Saveljic I, Fotiadis D, Parodi O (2013) Computer simulation of three-dimensional plaque formation and progression in the carotid artery. Med Biol Eng Comput 51:607–616 2. Ford MD, Xie YJ, Wasserman BA, Steinman DA (2008) Is flow in the common carotid artery fully developed? Physiol Meas 29:1335–1349 3. Funamoto K, Hayase T (2013) Reproduction of pressure field in ultrasonic-measurement-integrated simulation of blood flow. Int J Numer Methods Biomed Eng 29:726–740 4. Funamoto K, Hayase T, Shirai A, Saijo Y, Yambe T (2005) Fundamental study of ultrasonic-measurement-integrated simulation of real blood flow in the aorta. Ann Biomed Eng 33:415–428 5. Funamoto K, Hayase T, Saijo Y, Yambe T (2006) Numerical study on variation of feedback methods in ultrasonic-measurement-integrated simulation of blood flow in the aneurysmal aorta. JSME Int J, Ser C 49:144–155 6. Funamoto K, Hayase T, Saijo Y, Yambe T (2008) Numerical experiment for ultrasonic-measurement-integrated simulation of three-dimensional unsteady blood flow. Ann Biomed Eng 36:1383–1397 7. Funamoto K, Hayase T, Saijo Y, Yambe T (2009) Numerical experiment of transient and steady characteristics of ultrasonicmeasurement-integrated simulation in three-dimensional blood flow analysis. Ann Biomed Eng 37:34–49 8. Funamoto K, Hayase T, Saijo Y, Yambe T (2011) Numerical analysis of effects of measurement errors on ultrasonicmeasurement-integrated simulation. IEEE Trans Biomed Eng 58:653–663 9. Gelfand BD, Epstein FH, Blackman BR (2006) Spatial and spectral heterogeneity of time-varying shear stress profiles in the carotid bifurcation by phase-contrast MRI. J Magn Reson Imaging 24:1386–1392
10. Hayase T, Humphrey JAC, Greif R (1992) A consistently formulated QUICK scheme for fast and stable convergence using finitevolume iterative calculation procedures. J Comput Phys 98:108–118 11. Hoi Y, Wasserman BA, Lakatta EG, Steinman DA (2010) Carotid bifurcation hemodynamics in older adults: effect of measured versus assumed flow waveform. J Biomech Eng 132:071006 12. Kamenskiy AV, Dzenis YA, Mactaggart JN, Desyatova AS, Pipinos II (2011) In vivo three-dimensional blood velocity profile shapes in the human common, internal, and external carotid arteries. J Vasc Surg 54:1011–1020 13. Lee SW, Antiga L, Steinman DA (2009) Correlations among indicators of disturbed flow at the normal carotid bifurcation. J Biomech Eng 131:061013 14. Manbachi A, Hoi Y, Wasserman BA, Lakatta EG, Steinman DA (2011) On the shape of the common carotid artery with implications for blood velocity profiles. Physiol Meas 32:1885–1897 15. Morbiducci U, Gallo D, Ponzini R, Massai D, Antiga L, Montevecchi FM, Redaelli A (2010) Quantitative analysis of bulk flow in image-based hemodynamic models of the carotid bifurcation: the influence of outflow conditions as test case. Ann Biomed Eng 38:3688–3705 16. Moyle KR, Antiga L, Steinman DA (2006) Inlet conditions for image-based CFD models of the carotid bifurcation: is it reasonable to assume fully developed flow? J Biomech Eng 128:371–379 17. Patankar SV (1980) Numerical heat transfer and fluid flow. Hemisphere Pub. Corp, Washington DC 18. Poepping TL, Rankin RN, Holdsworth DW (2010) Flow patterns in carotid bifurcation models using pulsed doppler ultrasound: effect of concentric vs. eccentric stenosis on turbulence and recirculation. Ultrasound Med Biol 36:1125–1134 19. Press WH, Flannery BP, Teukolsky SA, Vetterling WT (1986) Numerical recipes: the art of scientific computing. Cambridge University Press, Cambridge, pp 274–282 20. Reneman RS, Hoeks APG (2008) Wall shear stress as measured in vivo: consequences for the design of the arterial system. Med Biol Eng Comput 46:499–507 21. Santhiyakumari N, Rajendran P, Madheswaran M, Suresh S (2011) Detection of the intima and media layer thickness of ultrasound common carotid artery image using efficient active contour segmentation technique. Med Biol Eng Comput 49:1299–1310 22. Schneider GE, Zedan M (1981) A modified strongly implicit procedure for the numerical solution of field problems. Numer Heat Transfer 4:1–19 23. Steinman DA (2002) Image-based computational fluid dynam ics modeling in realistic arterial geometries. Ann Biomed Eng 30:483–497 24. The Joint Committee of “The Japan Academy of Neurosonology” and “The Japan Society of Embolus Detection and Treatment” on Guideline for Neurosonology (2006) Carotid ultrasound examination. Neurosonology 19:49–69 25. Tortoli P, Michelassi V, Bambi G, Guidi F, Righi D (2003) Interaction between secondary velocities, flow pulsation and vessel morphology in the common carotid artery. Ultrasound Med Biol 29:407–415 26. Wake AK, Oshinski JN, Tannenbaum AR, Giddens DP (2009) Choice of in vivo versus idealized velocity boundary conditions influences physiologically relevant flow patterns in a subject-specific simulation of flow in the human carotid bifurcation. J Biomech Eng 131:021013
13
ORIGINAL ARTICLE
Development and feasibility study of a two‑dimensional ultrasonic‑measurement‑integrated blood flow analysis system for hemodynamics in carotid arteries Takaumi Kato · Kenichi Funamoto · Toshiyuki Hayase · Shusaku Sone · Hiroko Kadowaki · Tadashi Shimazaki · Takao Jibiki · Koji Miyama · Lei Liu
Received: 25 April 2014 / Accepted: 25 August 2014 © International Federation for Medical and Biological Engineering 2014
Abstract Prevention and early detection of atherosclerosis are critical for protection against subsequent circulatory disease. In this study, an automated two-dimensional ultrasonic-measurement-integrated (2D-UMI) blood flow analysis system for clinical diagnosis was developed, and the feasibility of the system for hemodynamic analysis in a carotid artery was revealed. The system automatically generated a 2D computational domain based on ultrasound color Doppler imaging and performed a UMI simulation of blood flow field to visualize hemodynamics in the domain. In the UMI simulation, compensation of errors was applied by adding feedback signals proportional to the differences between Doppler velocities by measurement and computation while automatically estimating the cross-sectional average inflow velocity. The necessity of adjustment of the feedback gain was examined by analyzing blood flow in five carotid arteries: three healthy, one sclerosed, and one stenosed. The same feedback gain was generally applicable for the 2D-UMI simulation in all carotid arteries, depending on target variables. Thus, the present system was shown to be versatile in the sense that the parameter
T. Kato · S. Sone Graduate School of Biomedical Engineering, Tohoku University, Sendai 980‑8579, Japan K. Funamoto (*) · T. Hayase Institute of Fluid Science, Tohoku University, 2‑1‑1 Katahira, Aoba‑ku, Sendai 980‑8577, Japan e-mail: [email protected] H. Kadowaki Graduate School of Engineering, Tohoku University, Sendai 980‑8579, Japan T. Shimazaki · T. Jibiki · K. Miyama · L. Liu GE Healthcare Japan, Hino 191‑8503, Japan
is patient independent. Moreover, the possibility of a new diagnostic method based on the hemodynamic information obtained by the 2D-UMI simulation, such as a waveform of the cross-sectional average inflow velocity and wall shear stress distributions, was suggested. Keywords Hemodynamics · Ultrasound color Doppler imaging · Computational fluid dynamics · Measurementintegrated simulation · Carotid artery
1 Introduction Circulatory disease is a major cause of death, its mortality rate being comparable with that of cancer. Since atherosclerosis is a key factor leading to the disease, much attention has been paid to its prevention and early detection [21]. It is well known that hemodynamic stresses, such as wall pressure and wall shear stress (WSS) [20], are closely related to the development and progression of atherosclerosis and subsequent circulatory disease. In order to utilize hemodynamic information for the development of a therapeutic strategy, various metrics concerning blood flow and hemodynamic stresses have been proposed [13] and their thresholds for pathological events have been discussed. Thus, acquisition of accurate, detailed information on hemodynamics is critical for advanced diagnosis of circulatory disease, but none of the present methodologies can completely provide such information. For instance, ultrasound Doppler measurement provides only a one-dimensional velocity component (Doppler velocity) of blood flow, which is a projection of the three-dimensional (3D) velocity vector in the ultrasound beam direction, with a compact system noninvasively [18, 25], but hemodynamic stresses cannot be measured. Phase-contrast magnetic resonance imaging (PC
13
MRI) may be used to obtain information on 3D velocity vectors of blood flow, as well as blood vessel configuration, enabling calculation of WSS distribution [9]. However, the temporal and spatial resolutions of PC MRI are limited, and its accuracy is insufficient to quantitatively evaluate hemodynamics. On the other hand, a numerical simulation of blood flow in a patient-specific vessel geometry reconstructed from medical images reveals the blood flow field and hemodynamic stresses in detail [1, 23]. However, it is difficult to correctly specify upstream boundary conditions and physiological parameters of the blood and the vessel wall [11, 15, 16, 26], and therefore, the computed result may differ from the real in vivo hemodynamics. In addition, validation of the computational result is difficult. As a solution for individual limitations in measurement and computation, the authors have previously proposed ultrasonic-measurement-integrated (UMI) simulation, in which artificial body forces proportional to the differences between measured and computed Doppler velocities are fed back to the numerical simulation, and the computational result converges to the real blood flow field [4]. Fundamental knowledge of the effects of feedback was obtained using two-dimensional UMI (2D-UMI) simulation of blood flow in a thoracic aneurysm [5]. In addition, efficiency and other characteristics of three-dimensional UMI (3D-UMI) simulation of blood flow were also investigated by numerical experiments [6, 7]. In spite of recent improvements in computational hardware, 3D simulation of blood flow still requires a significant amount of computational time, and therefore, a 3D-UMI blood flow analysis system is not realistic for clinical use. In contrast, a 2D-UMI blood flow analysis system could be clinically useful since it is convenient to use and provides results in a shorter time. However, such a 2D-UMI simulation system for clinical diagnosis is not yet available. Moreover, in former studies of UMI simulations, blood vessel morphology was extracted manually and the blood flow volume was assumed based on references, resulting in consideration under a semi-patient-specific condition. In this study, an automated 2D-UMI blood flow analysis system for clinical diagnosis was developed, and the feasibility of the system for hemodynamic analysis in a carotid artery was examined. The system automatically generated a 2D computational domain based on ultrasound color Doppler imaging and performed a UMI simulation of blood flow field to visualize hemodynamics in the domain. In the UMI simulation, compensation of errors was applied by adding feedback signals proportional to the differences between Doppler velocities by measurement and computation while automatically estimating the cross-sectional average inflow velocity. If it is necessary to adjust the feedback gain, which is a variable parameter significantly influencing the computational accuracy of the UMI simulation, the 2D-UMI blood flow analysis system
13
Med Biol Eng Comput
lacks the versatility for hemodynamic analysis. Hence, the versatility of the system for hemodynamic analysis in the carotid artery was examined by investigating effects of the feedback gain on the analysis results of five carotid arteries, including three healthy arteries, one sclerosed artery, and one stenosed artery. Then, the possibility of a new diagnostic method based on the hemodynamic information obtained by the 2D-UMI simulation was also investigated.
2 Methods 2.1 System configuration A block diagram and a snapshot of the 2D-UMI blood flow analysis system and the definition of the coordinate systems are shown in Fig. 1. The system is equipped with the following functions required for clinical diagnosis: (1) an automatic generation function of the 2D computational domain based on ultrasound color Doppler imaging, (2) a UMI simulation function for analysis of the blood flow field with compensation of errors by application of feedback signals proportional to the differences between measured and computed Doppler velocities, (3) an automatic estimation function of the cross-sectional average inflow velocity, and (4) a function for visualization of analysis results. In the system, ultrasound color Doppler imaging of hemodynamics is first carried out with conventional ultrasound diagnostic imaging equipment. The ultrasound color Doppler data are imported into a computer in DICOM format, and a longitudinal 2D blood vessel shape is automatically extracted from the data. A 2D-UMI simulation of blood flow is then performed by applying feedback signals, while the cross-sectional average inflow velocity is adjusted to reduce the errors between measured and computed Doppler velocities in the feedback domain. After a convergent result is obtained, the adverse effect of the feedback on the pressure field in the UMI simulation is compensated and the blood flow field (velocity vector and pressure distribution) and hemodynamic parameters such as WSS are visualized. In this study, ultrasound diagnostic imaging equipment (LOGIQ 7, GE Healthcare Japan, Tokyo, Japan) equipped with a linear ultrasound probe (12L, Fc = 8.2 MHz, GE Healthcare Japan, Tokyo, Japan) was employed for ultrasound color Doppler imaging. To perform the UMI simulation, a custom-built export function of raw data of ultrasound color Doppler was installed in the ultrasonic equipment. The computations were performed using one central processing unit (CPU) at 1.60 GHz in a server (Altix 3700 Bx2, SGI Japan, Tokyo, Japan). It should be noted that the computations can be performed with a general-purpose personal computer or a built-in computer in ultrasound diagnostic imaging equipment.
Med Biol Eng Comput
(a)
(b)
(c)
Fig. 1 The 2D-UMI blood flow analysis system: a schematic diagram, b overview of the system, and c coordinate systems in the system. Blood flow in a carotid artery of subject ID #1 is visualized in (b)
2.2 Analysis methods In spite of the three-dimensionality of blood flow, the 2D-UMI blood flow analysis system dealt with a cross section of a blood flow field along the axis of a blood vessel as a convenient clinical system. Deformation of the blood vessel was ignored for simplicity. 2.2.1 Extraction of blood vessel shape The shape of a blood vessel was extracted based on color data of Doppler velocity and on gray-scale data of tissue morphology, both of which were simultaneously provided in ultrasound color Doppler imaging. In a color Doppler image, the color indication of blood flow often extended through the deep-side blood vessel wall (the lower side wall in the ultrasound image) due to dispersion of ultrasound. Locations with nonzero time-averaged Doppler velocities over a measurement period tended to show a blood vessel lumen larger than the actual one. Hence, it was necessary to modify the blood vessel shape with reference to the gray-scale tissue information. The time-averaged gray-scale data in the measurement period were binarized after enhancement of the blood vessel wall with a differentiation filter in the ultrasound beam direction (vertical to the blood vessel). The blood vessel lumen was then determined by means of a logical product of both
binarized data of time-averaged color and gray-scale information. The coordinate system was defined as the Y and X directions toward the normal ultrasound beam direction and the head, respectively, as shown in Fig. 1c. To facilitate hemodynamic analysis, the xy coordinate system was used in the computation, rotating the XY coordinate system at an angle of β (counterclockwise is positive) to align the centerline of the extracted blood vessel shape with the x-direction. The ultrasound beam angle, θ, had a positive value clockwise against the Y-direction, and the net angle of the ultrasound beam, α, was expressed as the summation of θ and β. 2.2.2 2D‑UMI simulation The outline of the previously reported UMI simulation is described here [3]. The governing equations of the 2D-UMI blood flow simulation are the 2D Navier–Stokes equations for incompressible, viscous fluid flow, and the pressure equation as follows: ∂u + (u · ∇)u = µ�u − ∇p + f, ρ (1) ∂t
�p = −∇ · ρ(u · ∇)u + ∇ · f,
(2)
where u = (u, v) is the velocity vector, p is the pressure, t is the time, ρ is the density, μ is the viscosity, and f represents
13
Med Biol Eng Comput
a feedback signal to make the computational result approach the measurement data (real flow). The feedback signals, f, are applied as artificial body forces defined by the differences between Doppler velocities Vm and Vc on a cross section, which are measured by ultrasound color Doppler imaging and evaluated based on computed velocity vectors:
f=
Vc −Kv∗
− Vm U
ρU 2 L
· b,
(3)
where K*v is the feedback gain (nondimensional), U is the characteristic velocity, L is the characteristic length, and b is the unit vector along the ultrasound beam. Note that the special case with K*v = 0 is an ordinary simulation without feedback. The application of the feedback signals makes the velocity field converge to the real blood flow, satisfying the equation of continuity. However, it may adversely affect the pressure field in the case that the feedback signal does not satisfy the divergence-free condition in Eq. (2) [3]. Hence, the pressure field is modified by adding a compensation term calculated with the irrotational part of the feedback signal after a convergent result of the pressure field is obtained. The governing equations explained above were discretized and solved by a method similar to the SIMPLER method [17]. The convective terms were discretized with the reformulated QUICK scheme [10], and the time derivative terms were discretized with the first Euler implicit scheme. Linear algebraic equations were solved using a modified strongly implicit (MSI) scheme [22]. Feedback signals in the UMI simulation were added at the computational grid points in the feedback domain defined within the computational domain. Treatments were applied to detect the major measurement errors, such as aliasing, wall filter, and lack of data, and to reset the feedback signal to zero at the feedback point [8]. For evaluation of the computational accuracy, the space-averaged error of the Doppler velocity, e(t), in a monitoring domain Ωe was defined as follows:
e(t) =
1 |Vc − Vm | , NΩe uref
(4)
Ωe
where uref is the representative velocity for normalization and NΩe is the total number of the grid points in the monitoring domain Ωe. In addition, the time–space-averaged error, eave, was calculated by averaging the space-averaged error, e(t), in four cardiac cycles after a convergent solution was obtained. 2.2.3 Automatic estimation of the cross‑sectional average inflow velocity The exact boundary conditions are usually unknown for UMI blood flow simulation. A correct setting of the
13
cross-sectional average inflow velocity, uin(t), minimizes the error between computation and measurement [11, 26]. Consequently, the cross-sectional average inflow velocity uin(t) for a numerical simulation was estimated by means of the golden section search [19], which minimized the summation of the errors of the computed Doppler velocities Vc from the measured ones Vm in a monitoring domain Ωq defined in a computational domain: |Vc − Vm |. min uin (5) Ωq
A parabolic velocity profile with the estimated cross-sectional average inflow velocity was applied at the upstream boundary. In addition, free-flow condition (du/dx = 0 and dv/dy = 0) was applied at the downstream boundary, and the no-slip condition (u = 0) was assumed on the vessel wall. 2.2.4 Visualization of analysis results Application software for visualization of analysis results on Windows OS was programmed with Visual C++ by using Microsoft Visual Studio (Microsoft Japan, Tokyo, Japan). Velocity vectors at arbitrary locations obtained by the 2D-UMI simulation were visualized on a color Doppler image, in which Doppler velocities were contoured with warm or cold colors according to their positive or negative values, respectively. In addition, for clinical use, the above-mentioned visualized results were superimposed on the corresponding ultrasound B-mode image, and WSS was indicated with a color display on bands aligned with the blood vessel walls. 2.3 Analysis conditions In order to confirm the versatility of the system for hemodynamic analysis in the carotid artery, objectives were blood flow in five carotid arteries, including three healthy arteries, one sclerosed artery, and one stenosed artery, in five patients as summarized in Table 1. The study was approved by the local ethics committee, and informed consent was obtained from all patients. Clinical ultrasound color Doppler imaging was conducted following the ethics guidelines of GE Healthcare Japan Corporation. In the measurement, a longitudinal cross section was measured along the blood vessel axis as much as possible so as not to cause a large difference of flow volume among the axial positions. Ninety-one ultrasonic scan lines were acquired with an ultrasound beam angle θ of 20° or −20° from the ultrasound probe, and measurement data at 104 sampling points on each ultrasonic scan line were measured. The lateral and longitudinal pixel spacing was 303 and 182 μm, respectively. The angle of blood vessel β, net beam angle α,
Med Biol Eng Comput Table 1 Information on carotid arteries, ultrasonic measurement, and computation
ID
1
2
3
4
5
Age Sex Position Condition
53 Male Right Healthy 5.72
56 Male Left Healthy 5.72
77 Male Right Healthy 5.99
76 Female Right Sclerosed 7.16
71 Male Right Stenosed 5.68
-20 -5.2 -25.2
20 13.8 33.8
20 9.8 29.8
-20 -3.3 -23.3
20 11.8 31.8
Δt (ms) Nx × Ny
5.0 5.4 83.6 69 × 104
5.0 5.4 83.6 68 × 104
5.0 4.4 59.2 68 × 104
5.0 4.4 59.2 69 × 104
5.0 4.4 59.2 68 × 104
Δx, Δy (μm)
277, 164
227, 130
251, 146
287, 171
240, 138
Inlet diameter, D (mm) Beam angle θ (°) Angle of blood vessel β (°) Net beam angle α (= θ + β) (°) Working frequency (MHz) PRF (kHz)
working frequency, pulse repetition frequency (PRF), and temporal resolution Δt in the ultrasound imaging are summarized in Table 1. As for analysis conditions in the 2D-UMI simulation, the density ρ and viscosity μ of the blood were assumed to be 1.00 × 103 kg/m3 and 4.0 × 10−3 Pa s, respectively, within a normal range. The spatial resolution (Δx, Δy) of the computational grid was set to be that of the original ultrasound color Doppler imaging. The computational grid points are summarized in Table 1. The computational time increment was set to be identical to the temporal resolution Δt of the ultrasound imaging. The domain between 1/8 and 7/8 of the computational domain from the upstream boundary was defined as the feedback domain [5], and feedback signals, f, were added at all the grid points in the domain. In the estimation of the cross-sectional average inflow velocity, the monitoring domain Ωq was set so as to be identical to the feedback domain, and the value of uin which satisfied Eq. (5) was searched for in an adequately wide range between −0.1 m/s and 0.5 m/s. All variables were non-dimensionalized with characteristic values. The inlet diameter, D, was used as the characteristic length, L, an expedient value of 0.1 m/s was used as the characteristic velocity, U, and the density of blood, ρ, was used as the other characteristic value. The representative time-averaged velocity in a common carotid artery of 0.39 m/s was used as uref in the evaluation of the computational accuracy of Eq. (4) [24]. From here on, the same symbols are used for both dimensional and non-dimensional values for the sake of simplicity. The convergence criteria were set as a compromise between computational time and accuracy: 1.00 × 10−2 (nondimensional) for the residuals in the blood flow analysis, and 1.00 × 10−3 (1.00 × 10−4 m/s) for the estimation of the cross-sectional average inflow velocity. In order to confirm the versatility of the system for hemodynamic analysis in the carotid artery, the necessity
of adjustment of the feedback gain was examined. For the five carotid arteries, including healthy and diseased ones, the 2D-UMI simulation was performed by changing the feedback gain from 0 (ordinary simulation) to 100 in increments of 10, and was also performed with the feedback gains of 200, 300, 400, 500, and 1,000. Then, effects of the feedback gain on analysis results were investigated.
3 Results An example of the display visualizing an analysis result obtained by the present 2D-UMI blood flow analysis system is shown in Fig. 1b. A contour of Doppler velocities and velocity vectors are displayed, superimposed on the ultrasound B-mode image, and WSS distributions are visualized on bands along the blood vessels. A movie can be displayed to observe the variation of blood flow and hemodynamic stresses according to the heartbeat. Moreover, it is possible to compare the computational results with the measurement data obtained by the ultrasound diagnostic imaging equipment on the left-hand side. In order to confirm the versatility of the system, hemodynamic analysis in five carotid arteries, including three healthy arteries, one sclerosed artery, and one stenosed artery, was performed with various feedback gains. The measured color Doppler images at the peak flow of each subject (ID #1–ID #5) are shown in the left column of Fig. 2. Blood flows from right to left in all images, but Doppler velocities are colored with warm or cold color based on differences of the ultrasound beam angle θ (see Table 1). Subjects of ID #1 through ID #3 (Fig. 2a–c) are healthy carotid arteries; the bulge at the downstream of subject ID #2 (Fig. 2b) is the carotid bifurcation, and subject ID #3 (Fig. 2c) is a region extending from a carotid bifurcation into the internal carotid artery. On the other
13
Fig. 2 Measured color Doppler image at the peak flow (the left column) and computational results of the corresponding color Doppler image with velocity vectors and the Doppler velocity error distribution by the ordinary simulation (K*v = 0, the middle column) and the UMI simulation (K*v = 500, the right column) of each subject: a ID #1 (t = 0.92 s), b ID #2 (t = 2.09 s), c ID #3 (t = 3.55 s), d ID #4 (t = 3.02 s), and e ID #5 (t = 3.44 s). Blood flows from the right side to the left side. The computational domain of each subject is indicated with a two-way arrow in the measurement image
Med Biol Eng Comput
(a)
(b)
(c)
(d)
(e)
hand, subject ID #4 (Fig. 2d) is a common carotid artery with moderate atherosclerosis, the color Doppler image of which contains aliasing as indicated by an arrow. Subject ID #5 is a region extending from a carotid bifurcation into
13
the internal carotid artery with a stenosis. The color Doppler image at the stenosis shows a mosaic pattern (see a dotted circle region) since the hemodynamics are complicated by fast blood flow velocity.
Med Biol Eng Comput Fig. 3 Variations of the space-averaged error of Doppler velocity in the feedback domain, e, in the ordinary simulation (K*v = 0) and the UMI simulations (K*v = 200 and 500) of each subject: a ID #1, b ID #2, c ID #3, d ID #4, and e ID #5
(a)
(b)
(c)
(d)
(e)
The velocity vectors and color Doppler images obtained by the ordinary simulation (K*v = 0) and the UMI simulation with K*v = 500 at the same time step are represented in the upper middle and upper right, respectively, of each group in Fig. 2. The computational domains are indicated in the corresponding measured color Doppler images on the left. In the ordinary simulations in the middle column, the blood flow maintains the parabolic velocity profile applied at the upstream boundary since the carotid arteries are almost straight except for the stenosed artery in Fig. 2e. In contrast, the UMI simulation in the right column shows complicated velocity profiles in the feedback domain (the domain between 1/8 and 7/8 of the computational domain) even in the healthy carotid arteries. The errors of Doppler velocities between the measurement and each computation are shown in the lower monochrome images in Fig. 2, in which black indicates small errors. Error against the measurement data exists in the whole computational domain in the ordinary simulation, whereas it is suppressed in the feedback domain in the UMI simulations. Somewhat large errors are observed outside the feedback domain and at the locations of aliased data (see the arrow and the dotted circle in Fig. 2d, e).
Periodic variations of the space-averaged error of Doppler velocity in the feedback domain, e(t), the estimated cross-sectional average inflow velocity, uin(t), and the spaceaveraged WSS, τsave(t), in four cardiac cycles are compared between the ordinary simulation (K*v = 0) and the UMI simulations with K*v = 200 and 500 in Figs. 3, 4 and 5. Note that the cross-sectional average inflow velocity was successfully obtained after 21-time iterative computations of the golden section method at each time step, and the value of τsave was calculated by averaging WSS acting on both the shallow-side and deep-side (the upper and lower sides in an ultrasound image) walls of the blood vessel. The spaceaveraged error of Doppler velocity becomes small as the feedback gain increases in all subjects, implying improvement of the computational accuracy (Fig. 3). On the other hand, the estimated cross-sectional average inflow velocity is less sensitive to the feedback gain except for the subject with a stenosis of ID #5 (Fig. 4). The space-averaged WSS increases as the feedback gain increases (Fig. 5). Variations of the time–space-averaged error of Doppler velocity in the feedback domain, eave, the maximum value of the cross-sectional average inflow velocity, umax, the
13
Fig. 4 Variations of the cross-sectional average inflow velocity, uin, by the ordinary simulation (K*v = 0) and the UMI simulations (K*v = 200 and 500) of each subject: a ID #1, b ID #2, c ID #3, d ID #4, and e ID #5
Med Biol Eng Comput
(a)
(b)
(c)
(d)
(e)
time–space-averaged WSS, τtsave, and the average computational time for one time step, T, with the feedback gain are shown in Fig. 6, being normalized by the corresponding values of the ordinary simulation (K*v = 0). The values used for the normalization are summarized in Table 2. For each subject, the normalized time–space-averaged error of the Doppler velocity, eave/eave0, sharply decreases to half of that of the ordinary simulation in K*v ≤ 200 and converges with increasing feedback gain (Fig. 6a). The maximum value of the cross-sectional average inflow velocity, umax/umax0, shows little variation in K*v ≥ 500 (Fig. 6b). On the other hand, the time–space-averaged WSS, τtsave/τtsave0, increases with feedback gain, and the change of τtsave/τtsave0 becomes small in K*v ≥ 500 (Fig. 6c). The average computational time for one time step, T/T0, almost linearly increases with feedback gain for all subjects (Fig. 6d).
4 Discussion An automated 2D-UMI blood flow analysis system for clinical diagnosis was successfully developed. The system was
13
equipped with the following four functions: (1) an automatic generation function of the 2D computational domain based on ultrasound color Doppler imaging, (2) a UMI simulation function for analysis of the blood flow field with compensation of errors by application of feedback signals proportional to the differences between Doppler velocities by measurement and computation, (3) an automatic estimation function of the cross-sectional average inflow velocity, and (4) a function for visualization of analysis results. The versatility of the system for hemodynamic analysis in the carotid artery was proved, showing more accurate results than the ordinary simulation for blood flow in all carotid arteries treated in this study. Moreover, the possibility of a new diagnostic method based on the hemodynamic information obtained by the 2D-UMI simulation, such as a waveform of the cross-sectional average inflow velocity and WSS distributions, was also suggested. Effects of the feedback gain on parameters obtained by the present system are different. Although an increase of the feedback gain decreases the error of Doppler velocity against the measurement (Fig. 6a) as well as increasing the WSS (Fig. 6c), it hardly influences the estimation of
Med Biol Eng Comput Fig. 5 Variations of the spaceaveraged WSS, τsave, by the ordinary simulation (K*v = 0) and the UMI simulations (K*v = 200 and 500) of each subject: a ID #1, b ID #2, c ID #3, d ID #4, and e ID #5
(a)
(b)
(c)
(d)
(e)
Fig. 6 Variations of parameters with the feedback gain in the UMI simulation: a the time– space-averaged error of Doppler velocity in the feedback domain, eave, b the maximum value of the cross-sectional average inflow velocity, umax, c the time–space-averaged WSS, τtsave, in the four cardiac cycle, and d the average computational time, T, required for one time step. Each result was normalized with the corresponding result by the ordinary simulation (K*v = 0) in Table 2
(a)
(b)
(c)
(d)
13
Med Biol Eng Comput
Table 2 Parameters in the ordinary simulation (K*v = 0): the time– space-averaged error of Doppler velocity in the feedback domain, eave0, the maximum value of the cross-sectional average inflow velocID
1
2
3
ity, umax0, the time–space-averaged WSS, τtsave0, in the four cardiac cycle, and the average computational time, T0, required for one time step 4
5
Ave.
SD
eave0 (–)
0.075
0.101
0.065
0.064
0.101
0.081
0.019
umax0 (m/s) τtsave0 (Pa)
0.405
0.338
0.183
0.358
0.229
0.303
0.093
0.339
0.279
0.199
0.495
1.063
0.475
0.346
T0 (s)
13.4
12.7
10.8
12.2
20.2
13.9
3.7
the cross-sectional average inflow velocity (Fig. 6b). WSS is sensitive to slight changes in velocity since it is calculated from the spatial gradient of the velocity. Thereby, the variations of the error of Doppler velocity and the WSS become small with a feedback gain of more than about 200 and 500, respectively. The variations of the parameters are independent of the net beam angle of the ultrasound (see Table 1). For clinical use of the 2D-UMI blood flow analysis system, it is preferable to obtain an analysis result in a shorter time. The feedback gain should not be increased without restriction since the computational time almost linearly increases (Fig. 6d). Comprehensively considering the effects of the feedback gain on the analysis parameters and the computational time, the UMI simulation with a feedback gain between 200 and 500 can provide an appropriate result with better computational accuracy than the ordinary simulation. Comparison of the results of hemodynamic analysis with the same feedback gain between healthy and diseased carotid arteries shows a difference of feedback effects. The application of feedback in healthy carotid arteries reduces the error of Doppler velocity more than such application in diseased ones (Fig. 6a). The reduction rate of the time– space-averaged error of Doppler velocity by the feedback is smaller in the diseased carotid arteries than that in the healthy ones, though it becomes less than 50 % in all cases. This phenomenon could be derived from the complicated blood flow field in the diseased carotid artery and might be useful for diagnosis of carotid arteries: Diseased carotid arteries could be detected based on the reduction rate of the error. As observed in Fig. 2, the diseased carotid arteries have more complicated velocity distributions than the healthy ones, resulting in difficulty of reproduction of the hemodynamics. Moreover, the values of WSS in healthy carotid arteries by the ordinary simulation (K*v = 0) fluctuate close to zero, but those by the UMI simulations increase and show remarkable variations (Fig. 5a–c). In contrast, the values of WSS in diseased carotid arteries show fluctuations with maximum values of several Pascal, and the differences of the results between the ordinary simulation and the UMI simulations are not large (Fig. 5d, e). This reflects changes in the velocity field of blood flow near the
13
vessel wall owing to the feedback. For instance, the crosssectional average inflow velocity in a healthy carotid artery steeply changes in the acceleration and deceleration phases (Fig. 4a), and the velocity profile has a rectangular shape (Fig. 2a). However, the cross-sectional average inflow velocity in a carotid artery with atherosclerosis gradually decreases in the deceleration phase (Fig. 4d), and the velocity profile is approximately parabolic (Fig. 2d). Thus, the differences of variation and profile of blood flow velocity depending on the conditions of blood vessels cause differences in WSS distribution. This implies that the present system has the potential to detect atherosclerosis. It is expected that analysis by the 2D-UMI blood flow analysis system for more patients would clarify the relationships between circulatory diseases, such as atherosclerosis and stenosis, and hemodynamics, and help establish criteria to diagnose them. However, several tasks remain to be considered. Variability of the results obtained by different medical personnel, and usefulness of the system for follow-up diagnosis should be investigated. More importantly, a real blood flow has three-dimensionality [2, 12, 14, 25], and therefore, it is necessary to evaluate the effect of the three-dimensionality of the blood flow on the result of 2D-UMI simulation. Moreover, the equidistant orthogonal grid used in the system introduces inaccuracies in the reproducibility of the blood vessel configuration. Usage of a boundary-fitted or an unstructured grid system merits future consideration. Furthermore, neglect of wall motion influences the accuracy of the analysis result, especially near the wall, and thus, a moving boundary problem in the 2D-UMI blood flow analysis system should also be considered in future work.
5 Conclusion In this paper, a 2D-UMI blood flow analysis system was developed to easily obtain hemodynamic information. In order to confirm the versatility of the system for hemodynamic analysis in the carotid artery, the necessity of adjustment of the feedback gain of a variable parameter was examined by investigating effects of the feedback gain on
Med Biol Eng Comput
analysis results for five carotid arteries: three healthy, one sclerosed, and one stenosed. It was shown that the same feedback gain was generally applicable for the 2D-UMI simulation in all carotid arteries, yielding more accurate results with increasing feedback gain than those by the ordinary simulation without feedback. Thus, the present system is versatile in the sense that the parameter is patient independent. Moreover, the possibility of a new diagnostic method based on the hemodynamic information obtained by the 2D-UMI simulation, such as a waveform of the cross-sectional average inflow velocity and WSS distributions, was suggested. Acknowledgments Part of this study was supported by Grantin-Aid for Scientific Research (B) (24360064). All computations were performed using the supercomputer system at the Advanced Fluid Information (AFI) Research Center, Institute of Fluid Science, Tohoku University.
References 1. Filipovic N, Teng ZZ, Radovic M, Saveljic I, Fotiadis D, Parodi O (2013) Computer simulation of three-dimensional plaque formation and progression in the carotid artery. Med Biol Eng Comput 51:607–616 2. Ford MD, Xie YJ, Wasserman BA, Steinman DA (2008) Is flow in the common carotid artery fully developed? Physiol Meas 29:1335–1349 3. Funamoto K, Hayase T (2013) Reproduction of pressure field in ultrasonic-measurement-integrated simulation of blood flow. Int J Numer Methods Biomed Eng 29:726–740 4. Funamoto K, Hayase T, Shirai A, Saijo Y, Yambe T (2005) Fundamental study of ultrasonic-measurement-integrated simulation of real blood flow in the aorta. Ann Biomed Eng 33:415–428 5. Funamoto K, Hayase T, Saijo Y, Yambe T (2006) Numerical study on variation of feedback methods in ultrasonic-measurement-integrated simulation of blood flow in the aneurysmal aorta. JSME Int J, Ser C 49:144–155 6. Funamoto K, Hayase T, Saijo Y, Yambe T (2008) Numerical experiment for ultrasonic-measurement-integrated simulation of three-dimensional unsteady blood flow. Ann Biomed Eng 36:1383–1397 7. Funamoto K, Hayase T, Saijo Y, Yambe T (2009) Numerical experiment of transient and steady characteristics of ultrasonicmeasurement-integrated simulation in three-dimensional blood flow analysis. Ann Biomed Eng 37:34–49 8. Funamoto K, Hayase T, Saijo Y, Yambe T (2011) Numerical analysis of effects of measurement errors on ultrasonicmeasurement-integrated simulation. IEEE Trans Biomed Eng 58:653–663 9. Gelfand BD, Epstein FH, Blackman BR (2006) Spatial and spectral heterogeneity of time-varying shear stress profiles in the carotid bifurcation by phase-contrast MRI. J Magn Reson Imaging 24:1386–1392
10. Hayase T, Humphrey JAC, Greif R (1992) A consistently formulated QUICK scheme for fast and stable convergence using finitevolume iterative calculation procedures. J Comput Phys 98:108–118 11. Hoi Y, Wasserman BA, Lakatta EG, Steinman DA (2010) Carotid bifurcation hemodynamics in older adults: effect of measured versus assumed flow waveform. J Biomech Eng 132:071006 12. Kamenskiy AV, Dzenis YA, Mactaggart JN, Desyatova AS, Pipinos II (2011) In vivo three-dimensional blood velocity profile shapes in the human common, internal, and external carotid arteries. J Vasc Surg 54:1011–1020 13. Lee SW, Antiga L, Steinman DA (2009) Correlations among indicators of disturbed flow at the normal carotid bifurcation. J Biomech Eng 131:061013 14. Manbachi A, Hoi Y, Wasserman BA, Lakatta EG, Steinman DA (2011) On the shape of the common carotid artery with implications for blood velocity profiles. Physiol Meas 32:1885–1897 15. Morbiducci U, Gallo D, Ponzini R, Massai D, Antiga L, Montevecchi FM, Redaelli A (2010) Quantitative analysis of bulk flow in image-based hemodynamic models of the carotid bifurcation: the influence of outflow conditions as test case. Ann Biomed Eng 38:3688–3705 16. Moyle KR, Antiga L, Steinman DA (2006) Inlet conditions for image-based CFD models of the carotid bifurcation: is it reasonable to assume fully developed flow? J Biomech Eng 128:371–379 17. Patankar SV (1980) Numerical heat transfer and fluid flow. Hemisphere Pub. Corp, Washington DC 18. Poepping TL, Rankin RN, Holdsworth DW (2010) Flow patterns in carotid bifurcation models using pulsed doppler ultrasound: effect of concentric vs. eccentric stenosis on turbulence and recirculation. Ultrasound Med Biol 36:1125–1134 19. Press WH, Flannery BP, Teukolsky SA, Vetterling WT (1986) Numerical recipes: the art of scientific computing. Cambridge University Press, Cambridge, pp 274–282 20. Reneman RS, Hoeks APG (2008) Wall shear stress as measured in vivo: consequences for the design of the arterial system. Med Biol Eng Comput 46:499–507 21. Santhiyakumari N, Rajendran P, Madheswaran M, Suresh S (2011) Detection of the intima and media layer thickness of ultrasound common carotid artery image using efficient active contour segmentation technique. Med Biol Eng Comput 49:1299–1310 22. Schneider GE, Zedan M (1981) A modified strongly implicit procedure for the numerical solution of field problems. Numer Heat Transfer 4:1–19 23. Steinman DA (2002) Image-based computational fluid dynam ics modeling in realistic arterial geometries. Ann Biomed Eng 30:483–497 24. The Joint Committee of “The Japan Academy of Neurosonology” and “The Japan Society of Embolus Detection and Treatment” on Guideline for Neurosonology (2006) Carotid ultrasound examination. Neurosonology 19:49–69 25. Tortoli P, Michelassi V, Bambi G, Guidi F, Righi D (2003) Interaction between secondary velocities, flow pulsation and vessel morphology in the common carotid artery. Ultrasound Med Biol 29:407–415 26. Wake AK, Oshinski JN, Tannenbaum AR, Giddens DP (2009) Choice of in vivo versus idealized velocity boundary conditions influences physiologically relevant flow patterns in a subject-specific simulation of flow in the human carotid bifurcation. J Biomech Eng 131:021013
13
