A computer simulation model for Doppler ultrasound signals from pulsatile blood flow in stenosed vessels

Computers in Biology and Medicine 42 (2012) 906–914

Contents lists available at SciVerse ScienceDirect

Computers in Biology and Medicine journal homepage: www.elsevier.com/locate/cbm

A computer simulation model for Doppler ultrasound signals from pulsatile blood flow in stenosed vessels Lian Gao a, Yufeng Zhang a,n, Kexin Zhang b, Guanghui Cai a, Junhua Zhang a, Xinling Shi a a b

Department of Electronic Engineering, Information School, Yunnan University, Kunming, Yunnan 650091, China Cardiovascular Department, The Second Affiliated Hospital of Kunming Medical College, Kunming, Yunnan 650031, China

a r t i c l e i n f o

a b s t r a c t

Article history: Received 11 November 2011 Accepted 3 July 2012

A computer simulation model based on an analytic flow velocity distribution is proposed to generate Doppler ultrasound signals from pulsatile blood flow in the vessels with various stenosis degrees. The model takes into account the velocity field from pulsatile blood flow in the stenosed vessels, sample volume shape and acoustic factors that affect the Doppler signals. By analytically solving the Navier– Stokes equations, the velocity distributions of pulsatile blood flow in the vessels with various stenosis degrees are firstly calculated according to the velocity at the axis of the circular tube. Secondly, power spectral density (PSD) of the Doppler signals is estimated by summing the contribution of all scatterers passing through the sample volume grouped into elemental volumes. Finally, Doppler signals are generated using cosine-superposed components that are modulated by the PSD functions that vary over the cardiac cycle. The results show that the model generates Doppler blood flow signals with characteristics similar to those found in practice. It could be concluded that the proposed approach offers the advantages of computational simplicity and practicality for simulating Doppler ultrasound signals from pulsatile blood flow in stenosed vessels. & 2012 Elsevier Ltd. All rights reserved.

Keywords: Doppler ultrasound signal Computer simulation model Pulsatile flow Stenosed vessel

1. Introduction The stenosis, the plaque developed to the point where it significantly narrows the vessels, is one of the most serious forms of arterial disease. The clinical studies suggest that the stenosis is asymptomatic for decades and usually found in the highly curved arterial segment, vascular bifurcation, and the confluence of medium and large arteries. Once a stenosis is formed, it may disrupt the normal blood flow to the extent depending upon the degree of the deposition of substances in the arterial lumen, even leading to total vessel blockage in some instances [1] further influencing the development of the disease and arterial deformity [2]. Early and accurate detection of the artery stenosis is one of the most important diagnoses of cardiovascular diseases and has attracted the interest of many investigators. The Doppler ultrasound technique is one of the most important non-invasive diagnostic methods for cardiovascular disease by detecting and quantifying the status of pulsatile blood flow in vessels [3]. Under ideal uniform sampling conditions, the power in a particular frequency band of the Doppler spectrum is proportional to the volume of blood moving with velocities that produce frequencies in that band, and therefore the Doppler

n

Corresponding author. Tel.: þ86 0871 5031301; fax: þ86 0871 5031301. E-mail address: [email protected] (Y. Zhang).

0010-4825/$ - see front matter & 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compbiomed.2012.07.002

power spectrum should have the same shape as the velocity distribution plot for the flow in the vessel. Diagnostic information (mean velocity, maximum velocity, spectral broadening indices, skewness coefficients and turbulence indices) for clinical judgments are then derived from the spectrograms (the Doppler power spectrum varying with time) of the Doppler blood flow signals. Because of the time-varying nonstationary nature of the Doppler signals from the abnormal blood flow finding in stenosed vessels, such as vortices and turbulence, it is difficult to accurately estimate the spectrogram of Doppler blood flow signals. As a commonly used experimental method for evaluating and comparing the performance of novel Doppler blood flow signal processing approaches, simulation experiments were performed based on the data obtained from computer simulation models [4–8]. Because the theoretical velocity distribution or spectrogram of a simulated Doppler blood flow signal is known, it is more objective and accurate for the evaluation of experimental results based on simulation models by a comparison with the theoretical one. Moreover, simulation studies are less expensive and faster than those in vivo, as well as provide insights into physiological processes of pulsatile blood flow in the vessels with different stenosis degrees. Several models [9–16] for computing the velocity distribution, and then some of them further simulating the ultrasound signals from the steady [9] or pulsatile [10–16] blood flow in the stenosed vessels, have been reported. By assuming the blood to

L. Gao et al. / Computers in Biology and Medicine 42 (2012) 906–914

be represented by Herschel–Bulkley fluid, Sankar and Hemalatha [10] derived an analytic solution based on a perturbation method to analyze the flow velocity and wall shear stress distribution under the effects of pulsatility, stenosis and non-Newtonian behavior of blood. The other models [9,11–16], instead of using analytically described flow behavior, represented complex blood movement from velocity fields obtained by using computational fluid dynamics (CFD), which is a series of numerical solutions of the nonlinear Navier–Stokes equations by dividing the complex 3D geometry of blood vessels into a large number of small, regular elements and utilizing the finite-element method (FEM) [9,11–15], or the finite-volume method (FVM) [16] to acquire the pressure and velocity at each node. By further acoustic modeling of blood as a collection of point scatterers, result pulse wave Doppler signals [9,11–13] or RF-signals [14–16] could be efficiently retrieved using an existing ultrasound simulation model [9,11–14] or a commercial software (FIELD II software) [15–16]. Therefore, it is shown that the computer simulation models based on the CFD velocity distribution present flexibility and realism to a certain extent, and have been extensively studied. However, compare with the models based on an analytic flow velocity distribution, the CFD models require more complicated boundary and initial conditions for solving nonlinear equations, as well as high-performance computational hardware and long computational time, which limits its wide application. In this paper, a computer simulation model based on an analytic flow velocity distribution is proposed to generate Doppler ultrasound signals from pulsatile blood flow in the vessels with various stenosis degrees. The model takes into account the velocity field from pulsatile blood flow in the stenosed vessels, sample volume shape and acoustic factors that affect the Doppler signals. By solving the Navier–Stokes equations, the velocity distributions of pulsatile blood flow in the vessels with various stenosis are firstly calculated according to the velocity at the axis of the circular tube. Secondly, the power spectral density (PSD) of the Doppler ultrasound signals is estimated by the overalldistribution nonparametric estimation method. Finally, Doppler ultrasound signals are generated using cosine-superposed components that are modulated by a PSD function varying over the cardiac cycle. The paper is organized as follows: Section 2 presents the method of the proposed simulation model, including the governing equations for the pulsatile blood flow velocity distribution, the solutions for the blood velocity, the overall-distribution nonparametric estimation method for the spectrogram calculation and the simulation procedure for Doppler signals. Section 3 deals with the detailed description of experiments on simulation of Doppler pulsatile blood flow signals in normal vessels and vessels with various stenosis degrees of 20%, 40% and 70% (by diameter), followed by results, discussions and conclusions.

2. Method

Fig. 1. Geometric model of the stenosis.

artery in the non-stenotic region, 2x0 is the axial length of the stenosis and d is the measure of the degree of the stenosis. 2.2. The blood flow velocity distribution It was assumed that the blood flow in the stenotic cosineshaped arterial segment is pulsatile, axisymmetric, where the flowing blood is treated to be an incompressible Newtonian fluid. This investigation related to pulsatile blood flow in stenosed vessels with various stenosis degrees is based on the use of the linearized Navier–Stokes equations. The governing equations for the x and r components of momentum, together with the equation of continuity, in the cylindrical coordinate system can be written as [17] @u @x

þ 1r @ðrvÞ @r ¼ 0

@u 1 @p Z @2 u 1 @u ¼ þ þ @t r @x r @r 2 r @r

A cosine-shaped vessel segment with an axially symmetric stenosis is modeled as a rigid tube with a circular cross section. Let x-axis be taken along the axis of the artery while r is the radial coordinates. The geometry of the stenosis is shown in Fig. 1 and described as    d px RðxÞ ¼ R0 1 ð1Þ , x A ½x0 ,x0  1 þ cos x0 2R0 where R(x) denotes the radius of cosine-shaped arterial segment in the constricted region, R0 is the constant radius of the straight

!

@p @p ¼ ¼0 @r @y

ð2Þ

where u and v are the axial and the radial velocity components, respectively, p designates the pressure, r is the density, Z represents the dynamic viscosity of blood, and t is time. The boundary conditions are prescribed, where the velocity boundary conditions on the vessel wall are taken as vjr ¼ R ¼ 0;

ujr ¼ R ¼ 0

ð3Þ

and at the same time the condition on the symmetry axis is  @u ¼0 ð4Þ @r r ¼ 0 Because the transport of blood in the circulatory system depends on the periodic pumping action of the heart producing a pressure gradient throughout the artery, the blood flow can reasonably be considered to be periodic, and any periodic function can be expanded in a Fourier series. Suppose that f(t) is a periodic function of time, in this case, we can attempt to decompose f(t) as a sum of complex exponential functions: ! 1 X f ðtÞ ¼ F 0 þ Re F i ejoi t F0 ¼ Fi ¼

2.1. The stenosis model

907

i¼1 R 1 T f ðtÞdt ði ¼ T 0 R 2 T jwi t dt T 0 f ðtÞe

1,2,3, . . . :Þ

ð5Þ

hence, some flow parameters, such as pressure, pressure gradient and velocity, can be expressed as the Fourier series in the following complex exponential form [17] ! N X U i ejoi t u ¼ U 0 þRe i¼1

v ¼ V 0 þ Re

N X

! joi t

V ie

i¼1 @p @x

¼ C 0 þRe

N X i¼1

ð6Þ !

C i ejoi t

908

L. Gao et al. / Computers in Biology and Medicine 42 (2012) 906–914

where oi ¼ oi ¼ i 2Tp denotes angular frequency, T denotes the cardiac cycle, Re() is a function pffiffiffiffiffiffiffi which can get the real part of a complex number, and j ¼ 1 . The number of Fourier series expansion equations is finite as u,v and @p=@x are physiological flow parameters, and let it be N. Assuming that the pressure gradient is determined, it is clear that in (6), C0 and Ci are the determined quantities, and U0, Ui, V0 and Vi are the undetermined quantities to be computed. Substituting (6) into (2), and using the boundary conditions (3) and (4), the following equation can be obtained [18,19]. 2

U0 ¼

R ðy2 y21 ÞC 0 4Z

ð7Þ

U n0 ¼

C n0 R20 4Z

2 3 C ni R20 4 1  15 Ui ¼ jZa2i J j3=2 ai 0

4ZU n0

39 ,8 2 < = 1 2 15 R0 4  C i ¼ jZa2i U i 3=2 : ; J0 j ai

#

3=2

R J0 ðj ai yÞ 1 C i jZa2i J 0 ðj3=2 ai y1 Þ

   1 3 3=2 3=2 3=2 R3 4@ 1 J 0 j ai y yA dC i J1 j ai y1 J 1 j ai y 0 5    y1 C i Vi ¼  þ 3=2 2 dx jZa2i j ai J 0 j3=2 ai y1 J0 2 j3=2 ai y1

ð10Þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffi where y ¼ ðr=R0 Þ, y1 ¼ ðRðxÞ=R0 Þ, ai ¼ R0 roi =Z, i ¼ 1, 2, 3, . . ., y10 ¼ ðdy1 ðxÞ=dxÞ, and J0 and J1 are Bessel functions of the first kind of order zero and first, respectively. Substituting (7) and (9) into (6), the velocity distributions of pulsatile blood flow in the vessels with various stenosis degrees are obtained and written as 2  3=2 3 0 1 N 2 J0 j ai y X  R2 2 R 4  15C i ejoi t A uðy,t Þ ¼ y y1 2 C 0 þRe@ 4Z jZa2i J j3=2 ai y i¼1 0

1

ð11Þ Obviously, the pressure gradient (i.e. by the given coefficients C0 and Ci) is used to totally determine the axial velocity components. So it is necessary to calculate the pressure gradient in the following. For the stenotic arterial segment, the pressure gradient is changing continuously along the axis of the varying-area artery. Using the first boundary conditions of (3), (8) and (10) can be transformed to the following equations: ð12Þ 



J 21 j3=2 ai y1 dC i 2y10   Ci ¼ 0 þ dx y1 J j3=2 a y J j3=2 a y i 1 2 i 1 0

n

ð18Þ

ð9Þ

20

dC 0 4y10 þ C0 ¼ 0 dx y1

ð17Þ

R20

n

"

ð16Þ

Eqs. (15) and (16) can be rearranged in order to obtain the coefficients (C0 and Ci) of the pressure gradient at the site x¼x* of the artery:

ð8Þ

and 2

ð15Þ

n

C n0 ¼ 

  R y dC 0 4y1 y10 C 0 þ ð2y21 y2 Þ V0 ¼ 14Z dx 2

Ui ¼

segment of the stenotic artery by taking y ¼0 and y1 ¼1 into (7) and (9), and then the following equations can be obtained:

By solving the differential Eq. (12), C0 can be expressed as  Z  4y0 C 0 ¼ C exp  dx y1 1 ð19Þ ¼C 4 y1 According to the boundary condition given by (17), (19) can be rearranged as ð20Þ

The expression of Ci can be achieved as follows by solving the differential Eq. (13): ! Z 2J 2 1 ðj3=2 ai y1 Þ dy1 C i ¼ C exp  ð21Þ y1 J0 ðj3=2 ai y1 ÞJ 2 ðj3=2 ai y1 Þ Letting y2 ¼ j3=2 ai y1 , Eq. (21) can be rewritten as ! Z 2J2 1 ðy2 Þ dy2 C i ¼ C exp  y2 J0 ðy2 ÞJ 2 ðy2 Þ

Jn1 ðxÞ þ J n þ 1 ðxÞ ¼

2n J ðxÞ x n

ð23Þ

and for n ¼1, (23) can be reduced to

ð13Þ

x ðJ ðxÞ þ J2 ðxÞÞ 2 0

ð24Þ

by inserting (24) into (22), Ci is obtained by

 Z  2J 1 ðy2 Þ y22 ðJ0 ðy2 Þ þ J2 ðy2 ÞÞ C i ¼ C exp  dy2 y2 J 0 ðy2 ÞJ2 ðy2 Þ

 Z  y2 J 1 ðy2 ÞJ 0 ðy2 Þþ y2 J 1 ðy2 ÞJ 2 ðy2 ÞJ 0 ðy2 ÞJ 2 ðy2 Þ þ J 0 ðy2 ÞJ 2 ðy2 Þ ¼ C exp  dy2 y2 J 0 ðy2 ÞJ 2 ðy2 Þ  Z    y2 J 1 ðy2 ÞJ0 ðy2 ÞJ 0 ðy2 ÞJ 2 ðy2 Þ J ðy ÞJ ðy Þ y2 J 1 ðy2 ÞJ 2 ðy2 Þ ¼ C exp  þ 0 2 2 2  dy2 y2 J 0 ðy2 ÞJ 2 ðy2 Þ y2 J 0 ðy2 ÞJ 2 ðy2 Þ y2 J0 ðy2 ÞJ 2 ðy2 Þ   Z   y2 J 1 ðy2 ÞJ 2 ðy2 Þ 1 J ðy Þ ¼ C exp  þ  1 2 dy2 y2 J 2 ðy2 Þ y2 J 0 ðy2 Þ

  ¼ C exp  log½y2 J 0 ðy2 Þ2J1 ðy2 Þ þ logðy2 Þlog½J 0 ðy2 Þ ¼C

J 0 ðy2 Þ y2 2 J 0 ðy2 Þ2y2 J 1 ðy2 Þ

ð25Þ

Using the boundary condition (18), (25) takes the following form:

i¼1

where U n0 and U ni can be determined through the centerline velocity u(x*,0,t) at the certain position of the upstream faraway

ð22Þ

To solve Eq. (22), one of the recurrence relations of the Bessel function of the first kind is introduced as

J1 ðxÞ ¼

where J2 is the Bessel function of the first kind of order two. In order to obtain the solutions of Eqs. (12) and (13), appropriate boundary conditions are needed. The pressure gradient at the upstream faraway segment of the stenotic artery is determined by using the local axis velocity. Currently, using ultrasound Doppler technology, one can measure precisely the centerline (y¼0) velocity waveform u(x*,0,t) in the stenotic arterial segment upstream at certain position without the stenosis (for example, x¼ x* ¼ 5R0, where R(x*)¼R0 and y1 ¼1). The arterial segment at this position could be considered as a vessel without stenosis. Its Fourier series can be written as ! N X

 u xn ,0,t ¼ U n0 þ Re U ni ejoi t ð14Þ

1 y41

C 0 ¼ C n0

Ci ¼ Cin

J0 ðj3=2 ai Þ ðj

3=2

2

ai Þ J0 ðj

3=2

ai Þ2j

3=2

3=2

ai J 0 ðj

J0 ðy2 Þ 2 J ðy Þ2y J ðy Þ y ai Þ 2 0 2 2 1 2

ð26Þ

L. Gao et al. / Computers in Biology and Medicine 42 (2012) 906–914

909

Finally, by substituting C0 and Ci into (11), axial velocity u(y,t) from pulsatile blood flow in stenosed vessels is obtained.

contribution of all scatterers with the same flow velocity passing through the sample volume.

2.3. The power spectral density for the velocity distribution

2.4. The Doppler ultrasound signal simulation

When computing the PSD of the blood Doppler signal with a velocity distribution, the cross section artery is divided into a set of Cartesian grids, and then, the sample volume (SV) is also divided into many elemental volumes (EV) along the arterial axis. The division of the vessel for the PSD computation is shown in Fig. 2. The velocities of scatterers in each cirque area are presumed to be equal if the separation along the radius is small enough. The neighbor scatterers with similar velocities are grouped into an EV. The PSD is calculated by the overalldistribution nonparametric estimation method [20]. The probability mass associated with an EV region R is given by Z p ¼ pðxÞdx ð27Þ

After the PSD from a certain sample volume is obtained, we simulate Doppler ultrasound signals z(t) from pulsatile blood flow in vessels with various stenosis degrees. The method proposed by Mo and Cobbold [21] is used to simulate the Doppler signals z(t) from this sample volume:

R

where p(x) is overall probability density function. We collect N scatterers from p(x), and each has a probability P of falling within R. Then the number K of scatterers that lie inside R will be distributed according to the binomial distribution: P K ¼ C KN PK ð1PÞNK

ð28Þ

The mean of scatterers falling inside R is NP. For large N, this distribution will be sharply peaked around the mean: K  NP

R

where V is the volume of R. Eqs. (29) and (30) combined gives the overall probability density: K NV

ð31Þ

As N is lager enough, the overall probability p(x) approaches uniform distribution (1/V), and then the number of scatterers in an EV, which contributes the total PSD from this EV, is proportional to its volume. Thus, the volume of a single EV at a radial and axial position can be used to simply represent the PSD from this EV, where the frequency shift is calculated from the velocity of this sub-volume position by f ðtÞ ¼

M X

am exp½jð2pf m t þ jm Þ

ð33Þ

m¼1

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

 where f m ¼ mð1=2Þ Df , am ¼ 2Sz ðf m ÞDf ym , Df ¼ f max =M, M ¼ 2f max T. The ym and jm are independent random variables that account for the random nature of the signal. The random phase jm is uniformly distributed over [0,2p], ym is a w-squared random variable with two degree of freedom (w22 ), S(t,fm) is the spectrogram of Doppler signals at this position with the maximum frequency of fmax. The sinusoidal amplitudes am are Rayleigh distributed, and the number M of sinusoid functions is determined by the signal duration T and its maximum frequency fmax. An excellent approximation for the number of sinusoids needed estimate by M ¼ 2f max T

ð34Þ

ð29Þ

If we also assume that R is sufficiently small, then we have Z ð30Þ P ¼ pðxÞdx  pðxÞV

pðxÞ 

zðtÞ 

2uðy,tÞf 0 cosðyÞ c

ð32Þ

where y is the angle between the ultrasonic sound beam and the direction of blood flow, u(y,t) is the magnitude of the velocity of target, f0 is the frequency of transmitted ultrasound and c is the average velocity of ultrasound in soft tissue. The power of the Doppler signals can be conceived as a summation of the

3. Experiments With a view to examining the validity of the simulation model under consideration, a specific numerical illustration has been taken up to carry out numerical computation. The basic equations and the numerical methods mentioned above are applied to generate pulsatile blood flow in the normal vessel and vessel with different stenosis degrees of 20%, 40% and 70% at the positions x ¼1R0, 2R0 and 3R0 downstream of the point of the maximum stenosis, respectively. All the above simulation and analysis are conducted in Matlab 7.8 (The Mathworks, Inc., Natick, MA, USA). Parameters used in the simulation model are chosen to closely resemble those in the realistic Doppler ultrasound measurement system. Computer computations have been carried out using the following parameter values: the constant radius of the straight artery in non-stenotic region R0 ¼5 mm, stenosis length of the vessels 2x0 ¼8R0, blood density r ¼1050 kg/m3, blood viscosity Z ¼3.5  10  3 Pa s, cardiac cycle T¼ 1 s, the number of Fourier series expansion equations N ¼14, the Doppler transmitting

0.5

Element Volume u*(m/s)

0.4

0.3

0.2

0.1 0

0.2

0.4

0.6

0.8

1

t(s) Fig. 2. Division of the vessel for the PSD computation.

Fig. 3. Axial centerline velocity waveform u*(t) (at x ¼  5R0).

910

L. Gao et al. / Computers in Biology and Medicine 42 (2012) 906–914

frequency f0 ¼9000 Hz, the ultrasound velocity c¼1540 m/s, and the angle between the ultrasound beam and the blood flow direction y ¼301. The cross section of the vessel is divided by a 100  100 Cartesian grid, and the SV with the length 10 mm is divided into the 60 equal elements along the arterial axis. The blood flow is segmented into 50 cirques with the equal thickness along the radius. According to this model, the velocity distribution of the blood flow in the vessels with various stenosis degrees are computed using Eqs. (2)–(26) with inputs of the centerline (y¼0) blood velocity waveform u(x*,0,t) in the stenotic arterial segment upstream at the position (x¼x* ¼  5R0), which is given in precedence and shown in Fig. 3. Then the Doppler signals are simulated using Eq. (33) according to the PSD, which is estimated by Eqs. (31,32) from the blood flow velocity distribution. Finally the Doppler signals are simulated and the time-varying Doppler spectrograms are computed using a short time Fourier transform (STFT) with a temporal window of 10 ms to compare with the PSDs for performance evaluation.

4. Results and discussions

Table 1 Maximum blood pressure gradient pgmax and minimum blood pressure gradient pgmin in the normal vessel and the vessel with the stenosis degrees of 20%, 40% and 70%. Pressure gradient

pgmax (kPa/m) pgmin (kPa/m)

Stenosis degree 0% (Normal vessel)

20%

40%

70%

5.9439  2.5585

9.7326  3.8531

18.6851  6.2417

104.7027  9.6646

6 4

4

frequency(kHz)

pressure gradient(kPa/m)

Fig. 3 shows the centerline (y¼0) velocity waveform u(x*,0,t) in the stenotic arterial segment upstream at the position (x¼x*¼ 5R0), which is measured as the ensemble averaged mean velocity waveform from 50 realizations of a normal common carotid artery (the carotid artery without stenosis) by using a portable Doppler ultrasound system (KJ2V2U, Nanjing KeJin Industral Limted, Nanjing, China) and a 4-MHz transducer. According to the Fourier expansion coefficients of the centerline velocity waveform, the Fourier expansion coefficients of pressure gradient at the different positions of the artery are determined, and then the pressure gradients for the normal and stenosed vessels are obtained. As shown in Fig. 4(a) and Fig. 5, the pressure gradients from pulsatile blood flow in the normal vessel (Fig. 4(a)) and the vessel with

different stenosis degrees of 20% (Fig. 5(a)), 40% (Fig. 5(b)) and 70% (Fig. 5(c)) at different positions (x ¼ 4R0 ,3R0 , . . . ,3R0 and 4R0) along the stenosed segment are demonstrated for comparison. It does not lose generality to choose the pressure gradient in the normal vessel for discussing its change with time as an example. Just as shown in Fig. 4(a), when blood is ejected from the heart during systole, pressure in the artery rises, and during the early phase of diastole, it falls and even presents a negative value, and after that it propagates as a pulse wave. It also follows (see Fig. 5) that for a certain time, the axisymmetrical pressure gradient reaches its maximum along the axis where maximum stenosis is reached (x¼0), meanwhile, the pressure decreases along the radius, and then with the increasing distance downstream from the tightly-stenotic site along the axial direction. Also, as shown in Fig. 5(a)–(c), we find that the maximum and minimum values of pressure gradient become more larger with the increase of the stenosis degree, and those corresponding numerical value (the maximum value pgmax and minimum value pgmin of pressure gradient in the normal vessel and the vessel with different stenosis degrees) obtained from these figures are all listed in Table 1, which indicates that the all values of pgmax and pgmin are increasing with reduction of the vessel radius. Table 1 also shows that the values of pgmax and pgmin in the normal artery are smaller than those in stenosed arteries as larger radius of normal arteries results in the smaller impedance in the vessel generally. Especially, for the highest stenosis degree (70%) listed in

2 0 -2

3 2 1 0 -1

-4 0

0.2

0.4 t(s)

0.6

0.8

1

0

0.2

0.4

0.6

0.8

time(s)

Fig. 4. Pressure gradient (a), axial velocity distribution (b) and its PSD (c) of the blood flow from a normal common carotid artery.

Fig. 5. Pressure gradients from pulsatile blood flow in the vessels with different stenosis degrees of 20% (a), 40% (b) and 70% (c).

1

L. Gao et al. / Computers in Biology and Medicine 42 (2012) 906–914

the Table 1, the value of pgmax is 1661.5%, and the value of pgmin is 277.7% more than those in the normal vessel, respectively. Meanwhile, the pgmax and pgmin with the highest stenosis degree (70%) are 975.8% and 150.8% more than those with the lowest stenosis degree (20%), respectively. These experimental results about the pressure changing in the normal vessel and the vessels with various stenosis degrees shown in Fig. 4(a), Fig. (5) and Table 1 match the clinical observation [22]. The variation of axial velocity in the normal vessel and vessel with various stenosis degrees of 20%, 40% and 70% at different axial positions x ¼1R0, 2R0 and 3R0 over a cardiac cycle are shown in Figs. 4(b) and 6, respectively. Fig. 4(b) shows the axial velocity profile in the normal vessel. Note that, for this case, the pulsatile velocity profiles have no longer been monotonous profile at any certain time. The pulsatile velocity at the axis reaches its maximum, but during the quite big radial region, axial velocity along vessel radius changes flatly, and it falls rapidly and reaches to zero only near the vessel wall (where y ¼1). In addition, the blood flow velocity distribution for a certain radial position is also in a pulse pattern in response to the pulsatile pressure gradient (Fig. 4(a)). From Fig. 6, it can be clearly found that the changes of the pulsatile velocity during a cardiac cycle for a certain radial position in the stenosed vessels are the same as those in the normal vessel; the changes of the pulsatile velocity with radial

911

location y for a certain time are also similar to those of the normal vessel. However, owing to the much sharper decrease of the pressure gradient, the pulsatile velocity falls rapidly, and even results in negative velocity only near vessel wall (where y¼y1). We can see from Fig. 6 that for a certain axial position x, the increase of the stenosed degree results in the raise of maximum velocities for both the forward and reverse blood flows, and this influence is much more obvious when the stenosed degree reaches 70%, especially at axial position x¼1R0 shown in Fig. 6(c). However, for a comparatively small stenosed degree (20%), the influence of stenotic arteries on velocity distribution is weaker. Besides, for a certain stenosed degree, with the increase of the distance downstream from the tightly-stenosed position, the influence of stenosed artery on velocity distribution decreases. For convenience of the quantitative discussion about the influence of stenosis on the blood flow velocities, the maximum forward flow velocity uf_max and the maximum reverse flow velocity ur_max in the normal vessel and the vessel with the stenosis degrees of 20%, 40% and 70% at the certain position x¼1R0 are listed in Table 2. Compared with those of the normal vessel, the maximum velocity of forward flow increases by 649.0% while stenosis degree is 70%, 142.8 % while stenosis degree is 40%, and 47.6% while stenosis degree is 20%; the maximum velocity of the reverse flow increases by 2333.0% while stenosis degree is

Fig. 6. Axial velocity distribution from pulsatile blood flow in the vessel with different stenosis degrees of 20% (a, d, g), 40% (b, e, h) and 70% (c, f, i) at different axial positions of x ¼1R0(a, b, c), x¼ 2R0(d, e, f) and x ¼3R0(g, h, i).

912

L. Gao et al. / Computers in Biology and Medicine 42 (2012) 906–914

70%, 662.2% while stenosis degree is 40% and -2.7% while stenosis degree is 20%. This quantitative analytical conclusion is consistent with the theoretical one that the maximum velocities of the forward and reverse blood flows increase rapidly in the serious stenotic positions, but the change of the maximum forward flow velocity is smaller and the reverse flow velocity is hardly changed in the mild stenotic position [22]. Figs. 4(c) and 7 show the PSDs of the SV with the length 5 mm centered on x¼1R0, 2R0 and 3R0 in the normal vessel and vessels with the stenosis degrees of 20%, 40% and 70%, respectively. From these figures, it can be found that the PSDs of the blood flow are varied in a pulsatile pattern over a carotid cycle in response to the variation of the blood flow velocity. The frequency shift increases rapidly at the early phase of cardiac cycle and reaches its Table 2 Maximum forward flow velocity uf_max and the maximum reverse flow velocity ur_max in the normal vessel and the vessel with the stenosis degrees of 20%, 40% and 70% at the certain position x¼ 1R0. Velocity

Stenosis degree

uf_max (m/s) ur_max (m/s)

0% (Normal vessel)

20%

40%

70%

0.4885  0.0185

0.7209  0.0180

1.1859  0.1410

3.6590  0.4501

maximum value near time t¼ 0.2 s. When the stenosis degree varies and the axial centered position keeps unchanged, characteristic parameters of the Doppler blood flow spectrogram, the values of the bandwidth, the maximum frequency and the components of the negative frequency, become larger and larger with the increase of the stenosis degree. This indicates that the presence of the stenoses makes the bandwidth, the maximum frequency and the components of the negative frequency of Doppler blood flow signals bigger than the instance of the normal vessel. Moreover, with the increase of stenosis degree, the bandwidth of Doppler signals is broadened, and then the width of the spectral window becomes narrowing. Particularly, the change of the characteristic parameters for the vessel with the stenosis degree more than 40% appears more obviously than those less than 40%. Under the certain stenosed degree, the effect of stenosis on the characteristic parameters of spectrograms is stronger if the distance from the tightly-stenosed position is smaller, and the effect is weaker if the distance is larger. These experimental results are in agreement with the characteristics of the clinical spectrograms of the Doppler ultrasound blood flow signals from the normal common carotid artery and the artery with various stenosis degrees [22]. Finally, the Doppler ultrasound signals from pulsatile blood flow in vessels with various stenosis degrees are simulated by using the method of cosine-superposed components that are

Frequency(kHz)

Frequency(kHz)

5 3 1

Frequency(kHz)

14 7

10 6 2 -2

-1 0

0.2

0.4 0.6 Time(s)

0.8

1

0

0.2

0.4 0.6 Time(s)

0.8

56 48 40 32 24 16 8 0 -8

1

0

0.2

0.4 0.6 Time(s)

0.8

1

0

0.2

0.4 0.6 Time(s)

0.8

1

0

0.2

0.4 0.6 Time(s)

0.8

1

5 3 1

8

Frequency(kHz)

Frequency(kHz)

Frequency(kHz)

7 6 4 2 0 -2

-1 0

0.2

0.4 0.6 Time(s)

0.8

1

0

0.2

0.4 0.6 Time(s)

0.8

1

0

0.2

0.4 0.6 Time(s)

0.8

1

5 1

9

5 3 1 -1

-1

9

1

Frequency(kHz)

Frequency(kHz)

Frequency(kHz)

3

13

-3

7 5

17

7 5 3 1 -1

0

0.2

0.4 0.6 Time(s)

0.8

1

Fig. 7. PSD from pulsatile blood flow in the vessels with different stenosis degrees of 20% (a, d, g), 40% (b, e, h) and 70% (c, f, i) at different axial positions of x¼ 1R0(a, b, c), x¼ 2R0(d, e, f) and x ¼3R0(g, h, i).

L. Gao et al. / Computers in Biology and Medicine 42 (2012) 906–914

modulated by the computed PSDs. In order to visually evaluate the performance of the proposed simulation model, the spectrograms of the simulated Doppler blood flow signals are computed by using the STFT with a temporal window of 10 ms for comparison with the corresponding theoretical ones. Fig. 8 shows the estimated spectrograms in the vessel with stenosis degrees of 20%, 40% and 70% at the axial positions of x¼ 1R0, 2R0 and 3R0. It could be found from Fig. 8 that the characteristics of the estimated spectrograms from the simulated signals are in accordance with those of the corresponding PSDs, namely, increase of the stenosis in the vessel causes the proportional increase of forward and reverse directional blood flow velocity, broadening of the bandwidth, and decrease of the spectral window. Additionally, the profiles of the estimated spectrograms from the simulated signals (Fig. 8) are nearly as same as the corresponding theoretical ones (Fig. 7), except the granular pattern, known as Doppler speckles, appearing in the estimated spectrograms, which indicates that the synthesized signals have very similar speckle patterns as the clinical Doppler blood flow signals. However, due to the several limitations, which are axisymmetric flow, rigid stenosed vessel, and idealized stenotic geometry, imposed for analytically solving the velocity field, the proposed model is less flexible and realistic than those based on

the CFD velocity distribution. Moreover, we simply regard the blood flow in the arteries as a homogeneous and incompressible Newtonian fluid in the simulation model, and omit the effects of the radial velocity, the wall shear stress, the wall pressure and viscoelasticity, as well as the effect of viscoelastic tissue on pulsatile flow in arteries. In fact, the periodically pulsatile blood flow in the artery may cause the circumferential and axial motion of the elastic blood vessel and in turn the oscillation of the vessel affects that of the blood flow [23]. With considering more effect factors, the simulation model could be improved in future to accord with reality better. In the experiments, the spectrograms of the simulated Doppler blood flow signals are estimated by using the STFT for the model performance evaluation. It has been known that there are limitations for the STFT to perform time frequency analysis, namely, there is an obvious tradeoff between the poor spectral resolution introduced by short data windows and the spectral broadening that arises from nonstationary characteristics of the signal when using longer data windows [24]. Obvious errors between the estimated (Fig. 8) and theoretical (Fig. 7) spectrograms could be found, especially for the segments of the Doppler ultrasound signal from blood flow having a relative wide bandwidth and changing rapidly with time.

14

3 1

Frequency(kHz)

Frequency(kHz)

Frequency(kHz)

7 5

10 6 2 -2

-1 0

0.2

0.4 0.6 Time(s)

0.8

1

0

0.2

913

0.4 0.6 Time(s)

0.8

56 48 40 32 24 16 8 0 -8

1

0

0.2

0.4 0.6 Time(s)

0.8

1

0

0.2

0.4 0.6 Time(s)

0.8

1

0

0.2

0.4 0.6 Time(s)

0.8

1

5 3 1

6 4 2 0 -2

-1 0.2

0.4 0.6 Time(s)

0.8

1

Frequency(kHz)

5 3 1

0.2

0.4 0.6 Time(s)

0.8

0

0.2

0.4 0.6 Time(s)

0.8

1

13 9 5 1

1

7

9

5

7

3 1 -1

-1

17

-3 0

Frequency(kHz)

0

Frequency(kHz)

8

Frequency(kHz)

Frequency(kHz)

Frequency(kHz)

7

5 3 1 -1

0

0.2

0.4 0.6 Time(s)

0.8

1

Fig. 8. STFT-based spectrograms of a simulated Doppler signal from pulsatile blood flow in the vessels with different stenosis degrees of 20% (a, d, g), 40% (b, e, h) and 70% (c, f, i) at different axial positions of x¼ 1R0(a, b, c), x ¼2R0(d, e, f) and x¼ 3R0(g, h, i).

914

L. Gao et al. / Computers in Biology and Medicine 42 (2012) 906–914

5. Conclusions A computer simulation model based on the analytic velocity field, sample volume shape and acoustic factors that affect the Doppler signals is proposed to generate Doppler ultrasound signals from pulsatile blood flow in the vessels with various stenosis degrees in this study. By analytically solving the Navier–Stokes equations, the velocity distributions of pulsatile blood flow in the vessels with various stenosis are firstly calculated according to the velocity at the axis of the circular tube. Secondly, power spectral density (PSD) of the Doppler signals is estimated by summing the contribution of all scatterers passing through the sample volume grouped into elemental volumes. Finally, Doppler signals are generated using cosinesuperposed components that are modulated by a PSD function that varies over the cardiac cycle. The experimental results show that the proposed model provides the good performance for generating Doppler blood flow signals with characteristics in accordance with those found from the clinical normal common carotid artery and arteries with various stenosis degrees. It could be concluded that the proposed approach offers the advantages of computational simplicity and practicality for simulating Doppler ultrasound signals from pulsatile blood flow in stenosed vessels, and could be easily used as a preliminary experimental tool for the analysis of Doppler blood flow signals.

Conflict of interest statement The work ‘‘A Computer Simulation Model for Doppler Ultrasound Signals from Pulsatile Blood Flow in Stenosed Vessels’’ is supported by the Grant (60861001) from the National Natural Science Foundation of China and the Grant (2009CD016) from the Yunnan Natural Science Foundation. However, such supports do not influence the research work presented here. We certify here that we do not have any financial and personal relationships with other people or organizations that could inappropriately influence our work presented in this manuscript.

Acknowledgments This work was supported by the Grant (60861001) from the National Natural Science Foundation of China and the Grant (2009CD016) from the Yunnan Natural Science Foundation. References [1] T. Dahi, J. Bang, A. Ushakova, S. Lydersen, H.O. Myhre, Parameters describing motion in carotid artery plaques from ultrasound examination: a reproducibility study, Ultrasound Med. Biol. 30 (2004) 1133–1143.

[2] J.M. Romero, R.H. Ackerman, N.A. Dault, M.H. Lev, Noninvasive evaluation of carotid artery stenosis: indications, strategies, and accuracy, Neuroimaging Clin. N. Am. 15 (2005) 351–365. [3] D.H. Evans, W.N. McDicken, Doppler Ultrasound: Physics, Instrumentation and Signal Processing, Wiley, Chichester, 2000. [4] Q. Tao, Y. Wang, P.J. Fish, W. Wang, J.C. Cardoso, The wall signal removal in Doppler ultrasound systems based on recursive PCA, Ultrasound Med. Biol. 30 (2004) 369–379. [5] Y. Zhang, Y. Wang, W. Wang, B. Liu, Doppler ultrasound signal denoising based on wavelet frames, IEEE Trans. Ultrason. Ferroelectr. Freq. Control 48 (2001) 709–716. [6] Y. Zhang, N. Su, Z. Li, Z. Gou, Q. Chen, Y. Zhang, Assessment of arterial distension based on continuous wave Doppler ultrasound with an improved Hilbert-huang processing, IEEE Trans. Ultrason. Ferroelectr. Freq. Control 57 (2010) 203–213. [7] Y. Zhang, Y. Gao, L. Wang, J. Chen, X. Shi, The removal of wall components in Doppler ultrasound signals by using the empirical mode decomposition algorithm, IEEE Trans. Biomed. Eng. 54 (2007) 1631–1642. [8] D.S. Sankar, A two-fluid model for pulsatile flow in catheterized blood vessels, Int. J. Non-Linear Mech. 44 (2009) 337–351. [9] X. Fang, Y. Wang, W.Q. Wang, Doppler ultrasound signals simulation from vessels with various stenosis degrees, Ultrasonics 44 (Suppl. 1) (2006) e173–e177. [10] D.S. Sankar, K. Hemalatha, Pulsatile flow of Herschel–Bulkley fluid through stenosed arteries—a mathematical model, Int. J. Non-Linear Mech. 41 (2006) 979–990. [11] M. Khoshniat, M.L. Thorne, T.L. Poepping, S. Hirji, D.W. Holdsworth, D.A. Steinman, Real-time numerical simulation of Doppler ultrasound in the presence of nonaxial flow, Ultrasound Med. Biol. 31 (2005) 519–528. [12] M. Khoshniat, Real-Time Virtual Doppler Ultrasound, Thesis, University of Western Ontario, Canada, 2004. [13] H. Oung, F. Forsberg, Doppler ultrasound simulation model for pulsatile flow with nonaxial components, Ultrason. Imaging 18 (1996) 157–172. [14] K.W. Ferrara, V.R. Algazi, A theoretical and experimental analysis of the received signal from disturbed blood flow, IEEE Trans. Ultrason. Ferroelectr. Freq. Control 41 (1994) 172–184. [15] S. Balocco, O. Basset, J. Azencot, P. Tortoli, C. Cachard, 3D dynamic model of healthy and pathologic arteries for ultrasound technique evaluation, Med. Phys. 35 (2008) 5440–5450. [16] A. Swillens, L. Lovstakken, J. Kips, H. Torp, P. Segers, Ultrasound simulation of complex flow velocity fields based on computational fluid dynamics, IEEE Trans. Ultrason. Ferroelectr. Freq. Control 56 (2009) 546–556. [17] Z.R. Liu, Cardiovascular Hydrodynamics, Fudan University Presses, Shanghai, 1986 (in Chinese). [18] B. Liu, B. Guo, Z. Liu, Velocity profile of periodic oscillatory blood flow in circular tube, Chin. Q. Mech. 23 (2002) 15–22. (in Chinese). [19] H. Sun, Z. Liu, Velocity and shear stress of periodic oscillatory blood flow in vessel with stenosis, Chin. Q. Mech. 23 (2002) 148–156. (in Chinese). [20] Z.Q. Bian, X.G. Zhang, Pattern Recognition, Tsinghua University Press, 2000, pp. 1–60. [21] L.Y.L. Mo, R.S.C. Cobbold, Speckle in continuous wave Doppler ultrasound spectra: a simulation study, IEEE Trans. Ultrason. Ferroelectr. Freq. Control 33 (1986) 747–753. [22] W.C. Machey, D. Teso, Comprehensive Vascular and Endovascular Surgery, second ed., Mosby, USA, 2009. [23] C. Wang, Z. Liu, Analysis of oscillatory flow in consideration of a plasma layer in arterial stenoses, Appl. Math. Mech. (English Edition) 17 (1996) 1127–1135. [24] Y. Zhang, L. Xu, J. Chen, H. Ma, X. Shi, Correction for broadening in Doppler blood flow spectrum estimated using wavelet transform, Med. Eng. Phys. 28 (2006) 596–603.