CC BY
46Electronic Journal «Technical Acoustics» http://webcenter.ru/~eeaa/ejta/
2005, 23
Riccardo Leonardi1, Elisabetta Ronchieri1,2*
1 Department of Electrical System and Automation,University Of Pisa Via Diotisalvi, 2, 56100 Pisa, Italy
2 INFNCNAF, Viale Berti Pichat 6/2, 40127Bologna, Italy
Scattering of a distribution of encapsulated gas bubbles, a computer model
Received 17.05.2005, published 08.08.2005
The interest in using suspensions of polymeric encapsulated microbubbles as contrast agents in medical ultrasound imaging has increased during the last few years. Most of the knowledge on the non-linear mechanical and acoustical behaviour of encapsulated microbubbles is the result of empirical measurements, while the theoretical study of these phenomena, in particular the echo of a bubble distribution, is still an open research area. This paper presents a computer simulation program that calculates the global echo of a bubble distribution according to a modified Rayleigh-Plesset model for each bubble and by summing the contribution of each bubble to highlight the higher harmonics. Thanks to this approach, the study of the effects of the parameters of both the bubble and transmitted acoustic impulse on the higher harmonics contained in the echo is simplified. The results of this study are fundamental to the design of the microbubble or the echographic system.
INTRODUCTION
The echographic systems insonify human tissues by emitting acoustic impulses and receiving the corresponding echo. Extracting information about a specific tissue can be difficult, as the global echo consists of signals backscattered by any insonified tissue. Under these circumstances, contrast agents can improve the performance of the echographic systems by increasing the resolution and highlighting the vascularization, leading to a better diagnosis of pathologies like haemorrhages, blood clots, and partial localized reduction of the blood flux.
Using gas bubbles as acoustic contrast agents has some interesting advantages. Because of the contained gas, bubbles are more compressible than the other particles of similar size, such as the blood particles. This means that bubbles have a higher acoustical scattering cross section than the other particles. The echo from the bubbles can be distinguished from that of the surrounding tissues because the bubbles to tissue ratio is greater at the second harmonic than it is at the fundamental signal frequency [1]. In additions, bubbles’ non-linear behaviour can generate echoes at higher harmonics, while the tissue has virtually zero third harmonic with respect to the incoming acoustic impulse [1]. Hence, the global echo received by the
* contact person, e-mail: elisabetta.ronchieri@cnaf.infn.it
transducer can be filtered to extract only the higher harmonics, removing the effect of the tissue in order to obtain the spatial distribution of the bubbles. Using the echographic techniques the vascularization can be visualized permitting the diagnosis of vascular pathologies. The closer the bubbles are to the acoustic axis of the transducer, the higher the power of the impulse insonifying the bubbles and the more enhanced the non-linear behaviour of the bubbles are. Hence, by filtering the global echo at higher harmonics, the contrasts reduce the sidelobe artefacts and the lateral resolution of the acoustic system is improved [2].
In order to avoid the dissolution of the gas in the vascular system, the bubbles can be encapsulated, for instance by using polymeric shells. Unfortunately, the non-linear behaviour of the contrast agents is strongly influenced by the geometric and mechanical characteristics of the shell (e.g. the radius, the thickness and the shear parameters). Consequently, the study of the bubbles depends on the shell characterization. The interest in using suspensions of encapsulated microbubbles as contrast agents in medical ultrasound imaging has increased during the last few years. Some manufacturers have already developed instruments having the broad bandwidth transducers required to acquire the higher harmonics [1] and second-harmonic analysis is already used in the Medical UltraSound Contrast Imaging.
Most of the knowledge of the non-linear mechanical and acoustical behaviour of the encapsulated microbubbles has been produced by empirical measurements. The physical aspects of these phenomena have been the subject of several studies and a number of mathematical models for single bubble behaviour [3, 4, 5, 6, 7, 8, 9] are already available. Some interesting simulation software [10, 11] have been developed, but most of them, for example the BubbleSim simulator [11], only study a single bubble. In particular, the BubbleSim simulates the response of an encapsulated micrometer-sized bubble exposed to an ultrasound pulse, including the radial oscillation mode and the shell parameters. It is a valid instrument to predict the contrast bubble’s oscillation and the resulting echo for a specified driving ultrasound pulse. Other existing simulators, as the one described in [2], study bubbles distributions without considering the shell. The aim of the work presented in this paper is to reproduce and study more realistic scenarios by using the new BubbleDistributionEcho simulator. The 3D spatial distribution of encapsulated bubbles and the dispersion of the parameters values for all the bubbles in the distribution are introduced. The time behaviour of the global distribution echo is calculated by taking into account the non-linear behaviour of the bubbles and the spectrogram is produced.
This paper is organized as follows. Section 1 describes the bubble theory on which the present analysis is based, shows a brief overview of the solution of the non-linear equation for the motion of a bubble and summarizes the theory of the propagation process. Section 2 outlines the BubbleDistributionEcho simulator: the possible research applications, the parameters for the simulator, the calculations, and the simulated results. Section 3 compares results for a single bubble from two simulators - BubbleDistributionEcho and BubbleSim -by using the same parameters. Section 4 shows the results from our simulation on several different bubble distributions. Finally, the paper summarizes the results and a discussion of future research activities.
1. THEORY
This section presents the theory on which this work is based. The single bubble model and its solution are discussed. In addition, the theory of the propagating process used for calculating the global echo of the bubbles distribution is described.
1.1. The Rayleigh-Plesset bubble model
The simulator is based on the non-linear resolution of the Rayleigh-Plesset differential equation, modified for gas bubbles enclosed in a solid and incompressible viscoelastic shell [6]. The following model and assumptions [6, 7, 8] are applied.
We consider a spherical bubble. We assume the content of the bubble to be homogeneous and both the inside temperature and pressure to be always uniform. At equilibrium, the pressure in the gas inside the bubble is assumed to be equal to the hydrostatic pressure in the surrounding fluid. This means that the tension in the shell can be considered constant at equilibrium. The shell is considered to be solid, incompressible and viscoelastic. The mechanical parameters of the shell are considered to be constant within the range of frequency of this study. The shell is considered to be thin compared to the particle radius. The shell is assumed to reduce the surface tensions at the gas-shell and shell-liquid interfaces, so that surface tension can be neglected [6]. The temperature in the surrounding medium (for instance, the blood), is assumed to be a simple constant, as the temperature gradients are considered to be negligible, and the uniform heating of the liquid due to internal heat sources or radiation is not considered. We consider the medium to be described as an incompressible Newtonian fluid. The fluid density and dynamic viscosity are considered to be constant in time and uniform in space. The incident pressure is assumed to be known and controlled. It is the acoustic impulse generated by the transducer of the biomedical system.
Under these assumptions, the non-linear reduced equation of motion for the bubble is [6]
pf I rr + 3 r 2 J = Po
- Pi (t) - 4ßf — - 12^^-12G5^ ll - ^ I, (1)
r r r r I r I
where Pf is the density of the surrounding fluid, r(t) is the outer bubble radius, pO is the
atmospheric pressure, re is the outer bubble radius at the equilibrium, k is the polytrophic gas
exponent, pt (t) is the incoming acoustic pressure, is the viscosity of the surrounding
fluid, jUS is the shear viscosity of the shell, dSe is the shell thickness, and GS is the shear
modulus of the shell. It contains bubble radius r(t) as the only variable.
The polytrophic exponent of the gas is calculated by using the following formula (according to the equation (25) from reference [6]):
f d \ 1
K, - 4GS — V r, I
(2)
k =
where KB is the bulk modulus of the bubbles [Pa]
1.2. An alternative form of the non-linear bubble motion equation
The equation for the motion of the bubble boundary can be solved by applying the theory of the ordinary differential equations (ODE's) [13]. There are many different methods for solving ODE's at different levels of complexity. The single-bubble model was numerically solved with a fourth-order Runge Kutta self-starting method [13]. It represents an appropriate compromise between the competing requirements of a low truncation error per step and a low computational cost per step. This method is reasonably simple and robust and is a good general candidate for numerical solution of differential equations when combined with an intelligent adaptive step-size routine.
According to the theory of the ODE, we are able to solve the non-linear equation for the motion of the bubble (1) by representing it in the form r and specifying the initial conditions of r and r :
r (0) = re, r(0) = 0.
At this point, we can make the following substitutions:
(3)
(4)
r = y^
r = y 2 ,
which transforms the non-linear second-order equation in the following equivalent system of first-order equations:
y1 = y2 >
y =-
1
Pfy,
( f \3k \
' * ' - 1 yi
- Pi (t )-4^f — -
-- 12gM yi y{ yi yi
2 ( r e
1-------e
yi
(5)
For this system n = 2 and we require two boundary conditions in order to define the solution. These conditions must specify values of the dependent variables at certain points. We have the following initial conditions which enable us to integrate the equations numerically from the point t = 0 to some specified end-point:
(6)
y:(0) = re,
y2(0) = 0.
Table 1 reports the values used as initial conditions. The shear modulus of the shell G,,
the shear viscosity of the shell /us, the bulk modulus KB and the ratio —^— values shown in
d
Se
Table 1 are those for a contrast agent produced by Nycomed (Nycomed Amersham, Oslo, Norway) [6].
Table 1. Initial Conditions
Parameters Values Unit
yi = r re (e §- 10"6) [m]
y 2 = r 0 [m/s]
Pf 1050 [ Kg / m3]
Vf 0.004 [Paxs]
Us 0.39 [Paxs]
Gs 10.6 [MPa]
Po 101325 [Pa]
KB 2.5 [MPa]
re dSe 20 dimensionless
1.3. The propagating process
This section describes the theory and the assumptions about the transmitting process on which the present work is based [8, 14]. The spherical propagation without absorption is assumed for both the transducer transmitted impulse and bubble echoes. The frequency-dependent absorption due to the tissue has typically heavy consequences on the received echo, nevertheless in this work it has been neglected [2] (this will be included in future work). It is assumed that the bubbles do not affect the propagation of both the transmitted impulse and echoes of other bubbles. This is a valid assumption for the typical bubbles density and scattering strengths in this specific application [2], even if the multiple scattering can not be always neglected. The following calculations are executed for any bubble.
The transit time of the received impulse Tt generated by the transducer is
T = -, (7)
c
where d is the distance between the centre of the active surface of the transducer and the
bubble [m], c is the sound speed in the surrounding medium [m/s].
The transmitting directivity function of the rectangular transducer is
sini sin(a) — 1 sini sin(3) — 1 Ds(a, 3) = -±--------------3^—------------------------------------------------------------------^, (8)
sin(a) ~x sin(3) ~X where a and 3 are respectively the azimuth and the elevation orientation angles [rad], b and h are the corresponding dimensions of the active surface [m], A is the wave length of the
sound in the surrounding medium at the signal central frequency [m].
The incident pressure on the bubble is obtained by considering a spherical propagation without absorption. Hence, the power of the propagating impulse does not change if calculated on the total surface of concentric spheres starting from the sound source. By applying the following algebra, the incident pressure is calculated as a function of the acoustic signal power of the transmitted impulse, the transducer directivity factor, the transmitting directivity and the medium characteristics.
The maximum instantaneous a-directional intensity Ia is
G
Ia = ^ ’ (9)
where G is the maximum instantaneous acoustic signal power of the transmitted impulse [W] (setting the transmitter acoustic power should be more useful than setting the transmitted pressure magnitude for the target users of the BubbleDistributionEcho simulator).
The maximum instantaneous directional intensity Id on the acoustic axis inversely
depends on the transducer directivity factor Fd, see [14]:
1- = Ft , (1°)
d
A
d = 4nA - 2Ay[A + 2A2 ’ (11)
where A is the active surface of the transducer [m2 ].
The maximum instantaneous incident intensity It and the corresponding magnitude of the incident pressure Pi are
I= IdDs(a, 3)2, (12)
Pi ¡i 2pfc , (13)
where pf is the density of the surrounding medium [ Kg/m3].
The incident pressure is used for calculating the time behaviour of the bubble radius r(t) by solving the modified Rayleigh-Plesset model (see Section 1.2). The time behaviour of the backscattered pressure pb (t) calculated at 1 [m] from the bubble is (according to equation (2.7) from reference [8])
Pb (t) = Pf
F = r 2r,
where db is the distance from the bubble (in this case, 1 [m]).
(14)
The maximum instantaneous backscattered intensity Ib at 1 [m] and the corresponding maximum instantaneous backscattered power Gb [W] at 1 [m] are
(15)
(16)
where Pb is the magnitude of the backscattered pressure pb(t) [Pa].
By considering the spherical propagation without absorption and assuming that both the transmitting and receiving directivities of a rectangular transducer have the same value, the maximum instantaneous intensity I r received by the transducer and the magnitude of the corresponding received pressure Pr are
The transit time of the backscattered echo is assumed to be the same of the incident impulse on the bubble. Finally, the backscattered signals of all the bubbles are summed together according to the corresponding round trip transit time in order to obtain the total echo received by the acoustic transducer.
2. BUBBLE DISTRIBUTION ECHO SIMULATOR
The BubbleDistributionEcho simulator has been developed by using Matlab [15]. It has been designed and coded in order to be easily scalable and reusable. The BubbleDistributionEcho simulator matches the needs of the designers of medical echographic systems, producers of echographic contrast agents, and medical researchers. The following paragraphs describe the parameters, the possible calculations and the results from the use of the current release of the simulator.
2.1. Parameters
The parameters are grouped into the following areas: electro-acoustic transducer, surrounding medium, bubble properties, bubble distribution, acoustic impulse, and signal processing. Table 2 summarizes the most important simulator parameters, which are described in the following paragraphs.
2.1.1. Electro-Acoustic Transducer Parameters
Currently the only transmitter-receiver supported is the simple rectangular transducer, represented by specifying the 3D coordinates of the centre of the active face, the elevation and the azimuth orientation angles and the dimensions of the active surface. It is possible to simulate the same bubbles distribution by changing the position and the orientation of the transducer in order to study how this can affect the global echo. The effects of different directivities can be simulated by changing the dimensions of the active surface.
(17)
(18)
2.1.2. Surrounding Medium Parameters
The surrounding medium is the fluid where the bubbles are suspended (for example the blood), and is characterized by the sound speed c [m/s] (1550 for the blood), the density pf,
the viscosity , and the equilibrium pressure pO. The volumetric portion of tissue
containing the bubble distribution is represented by its dimensions (e.g. a cube having a side of 0.05 [m]).
Table 2 Simulator Parameters
Areas Parameters Meanings Unit
Surrounding Medium c sound speed [m/sec]
Pf fluid density [ Kg / m3]
Vf fluid viscosity [Paxs]
Po equilibrium pressure [Pa]
Bubble Properties Vs shear viscosity of the shell [Paxs]
Gs shear modulus of the shell [MPa]
r bubble radius [m]
re bubble radius at the equilibrium [m]
dSe shell thickness at the equilibrium [m]
fa y z) coordinates [m]
KB bulk modulus [MPa]
Acoustic Impulse a sign dimensionless
f central frequency [Hz]
g maximum instantaneous signal power [W]
N total number of hemicycles in the impulse dimensionless
2.1.3. Bubble Properties Parameters
Each bubble is characterized by its mechanical, geometrical and spatial parameters that are the shear viscosity of the shell /us, the shear modulus of the shell GS, the bubble outer radius
r, the bubble outer radius at equilibrium re, the shell thickness at equilibrium dSe, and the coordinates (x, y, z). The inner gas is represented by specifying the bulk modulus KB. The statistical distribution of the bubbles radii values is supported by setting the minimum and maximum radius values.
2.1.4. Bubble Distribution Parameters
The current simulator supports the cubic crystal, the fan and the random distributions. The cubic and the fan distributions are not realistic but allow rapid, controllable and easily understandable analysis, as shown in Figures 1 and 2. The random distribution is more realistic and can be used for comparing simulating and in-vitro experimental results. For any side of the cubic crystal distribution it is possible to specify the number of bubbles and the distance between two adjacent ones.
Fig. 1. A 3D cubic crystal distribution of bubbles having 5 bubbles per side. Bubbles are represented by the small blue spheres, while the transducer is the large black sphere
Fig. 2. A fan distribution of bubbles on the XZ plane. Bubbles are represented by the small blue spheres, while the transducer is the large black sphere
The cubic crystal can be used for representing cubic, parallelepiped and planar distributions that are useful for highlighting the solid angle of the transducer directivity, as shown in Figure 3. In the fan distribution it is possible to represent a portion of the area enclosed by two concentric circles by specifying both the minimum and maximum radius, the total angle, and the angular and radial distances between adjacent bubbles. The fan distribution is useful for highlighting the planar directivity of the transducer, as shown in Figure 4. The random distribution can be uniform (see Figure 5) or normal: if uniform, minimum and maximum values can be set for each coordinate; if normal, the mean value and the variance can be set for each coordinate. The analysis of bubble distributions during an insonification is simplified since the insonified bubbles are represented by red diamonds and the non-insonified ones by blue spheres (see Figure 3, Figure 4 and Figure 5).
2.1.5. Transmitted Acoustic Impulse Parameters
Currently the transmitted acoustic impulse supported by the simulator is a truncated sequence of sine cycles modulated with a low-frequency sine wave according to the following formula:
pt (t) = aPt sin(2nf t)window , (21)
window =
0 t < 0, t > T,
sin i^~j 0 < t < T,
(22)
T = N-
2f
where a is the initial sign, f is the central frequency [Hz], Pt is the maximum instantaneous pressure magnitude [Pa], and N is the total number of hemicycles in the impulse. The initial sign of the pressure specifies if starting with decompression (a = -1) or compression (a = 1). The window can be enabled or disabled.
Fig. 3. A 3D cubic crystal distribution with bubbles insonified by the acoustic impulse emitted by a transducer. It is a parallelepiped having 5 bubbles on Y and Z sides and 10 bubbles on X side (i.e., the acoustic axis of the transducer). The insonified bubbles are represented by red diamonds (see Section 2.2), while noninsonified bubbles are the small blue spheres, and the transducer is the large black sphere
Fig. 4. A possible fan distribution with bubbles insonified by the acoustic impulse emitted by a transducer. The insonified bubbles are represented by red diamonds (see Section 2.2), while non-insonified bubbles are the small blue spheres, and the transducer is the large black sphere
Fig. 5. 1000 bubbles uniformly
distributed in 0.01x0.005x0.005 [ m3] and insonified by the acoustic impulse emitted by a transducer. The insonified bubbles are represented by red diamonds (see Section 2.2), while non-insonified bubbles are the small blue spheres, and the transducer is the large black sphere
The sign parameter is utilized for studying how the bubble behaviour and its backscattered echo change whether the initial phase of the incident impulse is compression or decompression. The maximum instantaneous pressure magnitude Pt is calculated according
to the maximum instantaneous signal power G (see Section 1.3). The frequency and the number of hemicycles taken into account during the simulation are respectively 5 [MHz] and 10.
2.1.6. Signal Processing Parameters
The signal processing of the global echo received by the transducer is done by calculating its spectrogram, i.e. the time behaviour of the spectra of the acoustic signal. It is obtained by applying the discrete-time Fourier transform to the total echo by using partially overlapping sliding windows. The parameters that can be set to calculate the spectrogram are the length of the Fourier transform (that determines the frequencies at which the discrete-time Fourier transform is computed), the sampling frequency, the width of the window (the Hanning window is used [16]), and the number of samples by which the sections overlap.
2.2. Calculating
The sequence of calculations currently executed by the simulator is as follows. The computation starts by generating the specified 3D spatial distribution of bubbles. Both the mechanical, geometric and spatial parameters are initialized for any bubble. The transit time of the impulse generated by the transducer, the transmitting directivity and the incident pressure are calculated for each bubble (see Section 1.3). Bubbles contributing to the calculation of the global echo are chosen among all the bubbles by applying a configurable threshold to the value of the corresponding incident pressure. Only bubbles having the maximum value of the incident pressure over this threshold are taken into account. This "filtering" has been introduced in order to have a configurable reduction of the CPU and memory usage. Obviously, by using a zeroed threshold, all the bubbles will be involved in the calculation of the echo.
The radius and backscattered pressure time behaviours are calculated for any bubble by solving the Rayleigh-Plesset equation discussed previously. The equation is numerically solved by using the ode45 function from Matlab [17], implementing a fourth-order Runge Kutta method (see Section 1.2).
The transit time, the attenuation and the receiving directivity are calculated for each bubble echo received by the transducer (see Section 1.3). All the bubbles backscattered signals are summed together obtaining the total echo received by the acoustic transducer. The spectrogram of the total signal is calculated by applying the specgram function from Matlab [18].
2.3. Outputs
Currently, the BubbleDistributionEcho simulator is able to provide several types of visual outputs, including the 3D visualization of the transducer and the distribution of bubbles, the time behaviours of both the radius and echo of any bubble, the time behaviour of the total echo, and the 2D colour map specifying the spectrogram of the total echo. Figure 6(a) and Figure 6(b) show examples of the spectrogram 2D colour map.
(a) (b)
Fig. 6. These figures show the spectrograms of a simulated echo of a distribution of bubbles,
i.e., the time behaviour of the frequency components of the echo. (a) The region in the white box is zoomed in (b), where the higher harmonics are highlighted
3. VALIDATION TEST
Comparative tests have been executed by using the existing BubbleSim software by Hoff [10, 11]. These tests successfully confirmed the reliability of the BubbleDistributionEcho in relation to the simulation of the non-linear mechanical and acoustical behaviour of a single encapsulated bubble. The echo backscattered by the encapsulated bubble has been calculated at 1 [m] from the bubble geometric centre by setting the same values both in BubbleSim and BubbleDistributionEcho simulators. We simulated a gas bubble encapsulated in a thin polymeric shell from Nycomed [6].
Table 3. Simulator Parameters
Parameters Values Unit
Liquid Model Rayleigh-Plesset
Gas Model Isothermal
Liquid Blood
Pulse Envelope Rectangular
Pulse Amplitude 1.5 [MPa]
Pulse No. of Cycles 5
Pulse Center Freq. 5 [MHz]
Sample Rate 160 [MHz]
Bubble Radius 1 [pm]
Shell Thickness 50 [nm]
Shell Shear Modulus 10.6 [MPa]
Shell Viscosity 0.39 [Paxs]
Different sets of values of the parameters have been tested. For the sake of simplicity, only one of the tests executed is reported in what follows. Table 3 summarizes the values used as inputs for the two simulators during the experiments. Figure 7(a) and Figure 7(b) show the backscattered echoes calculated by using the Bubblesim and BubbleDistributionEcho simulators when applying the settings specified in Table 3. It is important to claim that only a qualitative agreement between the two simulators is obtained. A quantitative comparison of the two will be part of the future work.
(a) (b)
Fig. 7. The comparison of the backscattered echo of the same bubble (see Table 3) simulated by using the BubbleSim (solid line) and the BubbleDistributionEcho (dotted line). It also is an example of the non-linear response of the microbubbles when insonified.
(b) is a portion of (a)
4. RESULTS
Tests done at different amplitudes of the incident pulse confirmed the fact that the fraction of power associated with higher harmonics of the backscattered echo increases as the power of the incident impulse increases [2]. This behavior is explained by the fact that the non-linear behavior of the bubble is enhanced by increasing the incident power. In order to avoid bubble destruction, the power of the incident pressure cannot increase freely, so the simulations have been executed by monitoring the magnitude of the incident pressure on each bubble of the distribution, that was always less than 1 [MPa], a typical diagnostic imaging peak pressure [2]. In addition, simulations performed for different shell thickness values suggested that the portion of power associated with higher harmonics of the backscattered echo increases as the ratio between the bubble radius and the shell thickness increases (see Figure 8). The values in Table 1 were used and the ratio between the bubble radius and the shell thickness was changed in the range [15, 35]. This range was inspired from shell thickness values of real contrast agents, as the Nycomed one reported in Table 1.
The presence of the second and higher harmonics in the received echo has been demonstrated by simulating randomly distributed bubbles with statistical distribution of the radii values. Figure 9 shows the spectrogram of the global echo produced by a population of
1000 bubbles uniformly distributed in a volume of 0.01x0.005x0.005 [m3] (see Figure 5) by applying a uniform distribution of the bubbles radii in the range from 1 [pm] to 2 [pm]. This is a valid basis for the comparison with in-vitro experimental measurements.
(a) w (b)
(c)
Fig. 8. Different spectrograms obtained by using the BubbleDistributionEcho simulator for a bubble distribution. The same distribution is simulated only by changing the ratio between the bubble radius and the shell thickness (i.e., 15 in (a), 25 in (b) and 35 in (c)). All the bubbles in the distribution are identical. Figures (b) and (c) point out that the portion of the power associated to the higher harmonics of the backscattered echo increases as the ratio between the bubble radius and the shell thickness
increases
Fig. 9. The spectrogram of the global echo produced by a population of 1000 bubbles uniformly distributed in a volume of
0.01x0.005x0.005 [ m3] (see Figure 5) by applying a uniform distribution of the bubbles radii in the range from 1 [pm] to
2 [pm]. The presence of the second and higher harmonics is highlighted
The simulating experiments described in this paper have been executed by using a standard Pentium IV machine (2.8 [GHz], 512 [MB] RAM). The BubbleDistributionEcho time
consumption for a distribution of few bubbles (e.g. 1000) is negligible, but it increases almost linearly with the number of the bubbles in the distribution because the Rayleigh-Plesset equation is resolved for each bubble. When executing a sequence of simulations of a realistic distribution of contrast agents in the tissues (i.e. a possible usage of the BubbleDistributionEcho simulator), a single machine might not be performant enough. Future works will investigate this, also taking into account the Grid technology [19] that provides a very good means to add additional computing power for large scale simulation runs. A typical example for biomedical applications using this technology can be found in [20, 21].
CONCLUSIONS AND FUTURE ENHANCEMENTS
In this paper, we presented new software for simulating the acoustic backscattered echo of distributions of microbubbles and their non-linear mechanical and acoustical behaviour called BubbleDistributionEcho. The acoustic contrast agents used in the Medical Ultrasound Imaging Systems have been introduced, in particular the encapsulated bubbles. The Rayleigh-Plesset differential equation and the assumptions for the present simulation study of the bubble behaviour have been discussed, and the calculation used for obtaining the global distribution echo has been described, together with settings, features and outputs of the simulator. In addition, the BubbleDistributionEcho simulator has been compared to the BubbleSim simulator, showing that consistent results have been obtained when considering one bubble. Finally, preliminary but interesting results have been obtained, in particular the correlation between the shell thickness and the non-linear behaviour of the bubble, and the evidence of the second and higher harmonics by simulating randomly distributed bubbles with statistical distribution of the radii values.
This work sets the basis for the next activities related to the BubbleDistributionEcho simulator and for interesting future research applications. One of the major activities will be the comparison with experimental results. The future developing phases are focused on the following new aspects: the statistical dispersion of the bubble parameters (in the present work we studied the statistical dispersion of the radii values. Future works will address other bubble parameters, e.g. the shell thickness); the extraction of statistics from groups of simulations; the extraction of bubbles information (e.g. the 3D distribution) by filtering the global echo; more realistic 3D distribution of the bubbles in order to simulate the 3D vascularization of the human tissues; the study of other commercial acoustic contrast agents; different types of echographic transducers (currently only the rectangular transducer is taken into account); different types of realistic echographic impulses; the mechanical and acoustical interaction with the surrounding tissues; other differential equations that model the time behaviour of the bubble; performance analysis by simulating realistic bubbles distribution also in a Grid environment.
ACKNOWLEDGEMENTS
This research was supported by the Liguori Foundation of Livorno, Italy. The authors would like to thank Alessandro Mangoni and Henk Van Hasselt for their suggestions during the developing of the simulator. The authors are also grateful to Sergio Andreozzi, Andrea
Caiti, Tiziana Ferrari, Charles Loomis, Andrea Sciaba for their comments during the
preparation of this paper. The authors would like to thank referees for their very careful
comments and helpful suggestions which improved this note significantly.
REFERENCES
1. M. J. Monaghan. Contrast Echocardiography - Current clinical applications. The Eighth European Symposium on Ultrasound Contrast Imaging, 2003.
2. C. T. Chin, P. N. Burns. Predicting the acoustic response of a microbubble population for contrast imaging in medical ultrasound. Ultrasound in Med. & Biol., 26(8), 1293-1300, 2000.
3. A. Prosperetti. Bubble phenomena in sound fields: Part one. Ultrasonics 22, 69-77, 1984.
4. A. Prosperetti. Bubble phenomena in sound fields: Part two. Ultrasonics 22, 115-124, 1984.
5. T. G. Leighton. The Acoustic Bubble. Academic Press, London, 1994.
6. L. Hoff, Per C. Sontum, J. M. Hovem. Oscillations of polymeric microbubbles: Effect of the encapsulating shell. J. Acoust. Soc. Am., 107(4), 2272-2280, 2000.
7. L. Hoff. Acoustic characterization of Contrast Agents for Medical Ultrasound Imaging. Kluwer Academic Publishers, Dordrecth, The Netherlands, 2001.
8. C. E. Brennen. CAVITATION AND BUBBLE DYNAMICS. Oxford University Press, 1995.
9. C. S. Clay, H. Medwin. Acoustical Oceanography, principles and applications. Wiley-Interscience Publication, New York, 1997.
10. L. Hoff. Modelling and Characterization of Ultrasound Contrast Agent. The Leading Edge in Diagnostic Ultrasound in Atlantic City, NJ, 8 May 2001. http://home.online.no/~fam.hoff/Bubblesim/AbstractLeadingEdge.PDF.
11. L. Hoff. Ultrasound Contrast Bubble Simulation. Bubblesim. 2001.
http://home.online.no/~fam.hoff/Bubblesim/Readme.pdf.
12. P. Mazzoldi, M. Nigro, C. Voci. Fisica, Volume I, Meccanica - Termodinamica. EdiSES, Napoli.
13. E. Hairer, S. P. Narsett, G. Wanner. Solving Ordinary Differential Equations I: Nonstiff Problems. 2nd ed., Springer-Verlag, Berlin, 1993.
14. G. Pazienza. Fondamenti della Localizzazione Sottomarina. Poligrafico Accademia Navale, 1990.
15. L. Tsang, J. A. Kong, Kung-Hau Ding. Scattering of Electromagnetic Waves: Theories and Applications. John Wisley & Sons, Inc., 2000.
16. A. V. Oppenheim, R. W. Schafer. Discrete-Time Signal Processing. Englewood Cliffs, NJ: Prentice-Hall, 1989.
17. L. F. Shampine, M. W. Reichelt. The MATLAB ODE Suite, to appear in SIAM Journal on Scientific Computing, vol. 18-1, 1-22, 1997.
18. P. A. Lynn, W. Fuerst, B. Thomas. Introductory Digital Signal Processing with Computer Applications. Paperback - 494 pages, 2nd Ed edition (29 April, 1998) John Wiley and Sons; ISBN: 0471976318.
19. G. Avellino, et al. The DataGrid Workload Management System: Challenges and Results. Journal of Grid Computing, vol. 2, N°4, 353-367, 2004.
20. J. Montagnat, et al. Medical Images Simulation, Storage, and Processing on the European DataGrid Testbed. Journal of Grid Computing, vol. 2, N°4, 387-400, 2004.
21. N. Jacq, et al. Grid as a bioinformatic tool. Parallel Computing, vol. 30, 1093-1107, 2004.
