Polydisperse Versus Monodisperse Microbubbles: A Simulation Study for Contrast-Enhanced Ultrasound Imaging

Abstract

Objective

Contrast-enhanced ultrasound (CEUS) presents distinct advantages in diagnostic echography. Utilizing microbubbles (MBs) as conventional contrast agents enhances vascular visualization and organ perfusion, facilitating real-time, non-invasive procedures. There is a current tendency to replace traditional polydisperse MBs with novel monodisperse formulations in an attempt to optimize contrast enhancement and guarantee consistent behavior and reliable imaging outcomes. This study investigates the contrast enhancement achieved using various-sized monodisperse MBs and their influence on non-linear imaging artifacts observed in traditional CEUS.

Methods

To explore the differences between monodisperse and polydisperse populations without excessive experimentation, numerical simulations are employed for delivering precise, objective and expeditious results. The iterative non-linear contrast source (INCS) method has previously demonstrated efficacy when simulating ultrasound propagation in large populations in which each bubble has individual properties and several orders of multiple scattering are significant. Therefore, this method is employed to realistically simulate both monodisperse and polydisperse MBs.

Results

Our findings in CEUS imaging indicate that scattering from resonant monodisperse MBs is 11.8 dB stronger than scattering from polydisperse MBs. Furthermore, the amplitude of non-linear imaging artifacts downstream of the monodisperse population is 19.4 dB stronger compared with polydisperse suspension.

Conclusion

Investigating the impact of multiple scattering on polydisperse populations compared with various monodisperse suspensions has revealed that monodisperse MBs are more effective contrast agents, especially when at resonance. Despite the strong signal-to-noise ratio of monodisperse populations, imaging artifacts caused by non-linear wave propagation are also enhanced, resulting in further mis-classification of MBs as tissue.

Introduction

Achieving superior deep tissue imaging of blood vessels with ultrasound remains a challenge in medical diagnostics but contrast-enhanced imaging, particularly using microbubbles (MBs), has emerged as a promising solution [ , ]. These gas-filled microspheres, stabilized with a lipid or protein shell, enhance blood contrast for improved organ and lesion visualization. MBs, characterized by their small size, biocompatibility and vascular navigability, resonate in the ultrasound frequency range (1–10 MHz). Their efficient sound scattering in both fundamental and harmonic modes, driven by substantial acoustic impedance differences with their surroundings and highly non-linear oscillatory behavior [ , ], enhances image quality. As ultrasound waves propagate through a resonant MB suspension, they undergo non-linear distortion due to non-linear MB scattering influenced by size, shell characteristics, ultrasound pressure and frequency [ ]. Due to these properties, MBs have also been demonstrated to be efficient contrast agents in various applications in addition to contrast-enhanced ultrasound (CEUS), such as ultrasound localization microscopy [ ].

As a drawback, wave distortion can extend beyond an MB suspension, leading to the mis-identification of tissues as MBs and diminishing the specificity of CEUS imaging [ ]. Narrowing the size distribution of the MB population could be a way to provide improved acoustic scattering, reduce imaging artifacts and enhance scattering homogeneity. Historically, polydisperse MBs of varying size distributions (typical radii 0.5–15 µm) have been standard in ultrasound contrast imaging [ , ]. Recent technological breakthroughs have introduced the possibility of using monodisperse, i.e. , uniformly sized, MBs [ ]. Studies have highlighted the superiority of monodisperse MBs [ ], offering enhanced predictability, improved acoustic performance and clearer imaging signals [ , ]. Nevertheless, we believe that it is important to shed more light on the effect of monodisperse MBs as contrast agents for deep vessel imaging, particularly for the generation of clearer echoes and reducing imaging artifacts.

The use of computational tools is an efficient way to perform comprehensive investigations without performing extensive measurements. Initially, studies focused on the collective behavior of bubbly media for marine applications [ , ]. Effective medium theory facilitated 1-D computational studies for both monodisperse [ , ] and polydisperse [ , ] MB suspensions in medical ultrasound, including high-intensity focused ultrasound [ ]. Previous models successfully captured non-linear ultrasound propagation through uniform MB distributions in 2-D using iterative schemes [ , ]. Challenges have arisen when coupling the non-linear dynamics of multiple MBs in 3-D realistic simulations due to the complexity of the coupled Rayleigh-Plesset equation [ ]. Another difficulty is when the number of polydisperse MBs is small and the use of averaged quantities becomes questionable. Various computational methods have been explored to understand the dynamics between polydisperse and monodisperse MB populations. Among these, the iterative non-linear contrast source (INCS) method has demonstrated efficacy in simulating bubble cloud behavior in a 3-D domain when excited either by a plane wave or a focused beam. This method enables the generation and comparison of echoes produced by dense monodisperse MB populations, taking into account multiple scattering [ ]. This is crucial for optimizing CEUS applications and reducing the need for excessive experimentation.

The aim of this numerical study is to investigate the efficacy of monodisperse and polydisperse populations when used as contrast agents for deep tissue imaging. More precisely, this article discusses an extension of the INCS method to simulate the behavior of a population of polydisperse scatterers.

Furthermore, the effectiveness of the extended INCS method is illustrated by simulating the multiple scattering that occurs inside a population of polydisperse MBs, each with individual properties represented by its own Marmottant model [ ]. INCS is based on an iterative scheme for computing the scattered acoustic signals [ , ]. Numerically, the accuracy of the final result improves after each iteration. In a physical sense, each iteration adds an extra order of multiple scattering corresponding to an additional path of wave propagation.

First, in the ‘Inclusion of a polydisperse MB population’ section below, the fundamental theory behind the INCS method is explained, followed by its extension with the introduction of polydisperse point scatterers. In the section ‘Configurations used in the simulations,’ configurations for the numerical experiments are discussed. Next, in the ‘Numerical results’ section, the results from the numerical simulations for each different test case are presented. Concluding remarks are given in the ‘Conclusion’ section.

Inclusion of a polydisperse MB population

Linear field

The linear pressure field generated by an external source in a linear, lossless, homogeneous acoustic background medium is described by the wave equation (eqn [ 1 ]):

𝑐−20𝜕2𝑝(𝑥,𝑡)𝜕𝑡2−𝛻2𝑝(𝑥,𝑡)=𝑆pr(𝑥,𝑡)
c0−2∂2p(x,t)∂t2−∇2p(x,t)=Spr(x,t)

Here, x [m] is the Cartesian position vector and t [s] is the time. Furthermore, p ( x,t ) [Pa] is the acoustic pressure and 𝑐0=1/𝜌0𝜅0‾‾‾‾‾√
c 0 = 1 / ρ 0 κ 0
[m/s] is the small signal sound speed in the background medium, where ρ ₀ [kg⋅m ⁻³ ] is the mass density and κ ₀ [Pa ⁻¹ ] is the compressibility. The symbol ∇ ²indicates the Laplacian operator. The acoustic field is generated by the primary source term, S _pr (eqn [ 2 ]):

𝑆pr(𝑥,𝑡)=𝜌0𝜕𝑞(𝑥,𝑡)𝜕𝑡−𝛻·𝑓(𝑥,𝑡)
Spr(x,t)=ρ0∂q(x,t)∂t−∇·f(x,t)

where q ( x ,t ) [s ^-1 ] is the volume density of the volume injection rate and f ( x ,t ) [N/m ³ ] is the volume force density of the external source. A source with a plane aperture, e.g. , a phased array transducer, can be represented by a pressure jump in the particle velocity or pressure. This can formally be described by a primary source term with a Dirac delta function [ ], but in our numerical implementation of INCS, the source was modeled as a boundary condition at the location of the source aperture.

Non-linear field due to contrast agents

In medical ultrasound, non-linearities arising from contrast media can have a significant impact on the propagation of the acoustic signals. To incorporate any phenomena that affect the pressure field, we extended eqn (1) with a contrast source term, S _cs (eqn [ 3 ]):

𝑐−20𝜕2𝑝𝜕𝑡2−𝛻2𝑝=𝑆pr+𝑆cs(𝑝).
c0−2∂2p∂t2−∇2p=Spr+Scs(p).

With this approach, multiple contrast sources can be accommodated that represent global non-linear effects [ , ], attenuation [ , ], non-homogeneous medium properties [ ] or local non-linear effects [ ]. In contrast-enhanced imaging, the non-linear oscillatory behavior of the MBs influences the pressure field. To include the contribution of a population of a total number N of MBs, each was described as a point scatterer and the source term was written as eqn (4) [ ]:

𝑆cs(𝑥,𝑡)=∑𝑁𝑖=1𝑆𝑀𝐵𝑖=𝜌0∑𝑁𝑖=1𝑑2𝑉(𝑖)(𝑥(𝑖)sc,𝑡)𝛿(𝑥−𝑥(𝑖)sc)
Scs(x,t)=∑i=1NSMBi=ρ0∑i=1Nd2V(i)(xsc(i),t)δ(x−xsc(i))

where V ⁽ ⁱ⁾ [m ³ ] is the volume of the i th MB,

𝑥(𝑖)sc
x sc ( i )
is the position vector of its center and

𝛿(𝑥−𝑥(𝑖)sc)=𝛿(𝑥−𝑥(𝑖)sc)𝛿(𝑦−𝑦(𝑖)sc)𝛿(𝑧−𝑧(𝑖)sc)
δ ( x − x sc ( i ) ) = δ ( x − x sc ( i ) ) δ ( y − y sc ( i ) ) δ ( z − z sc ( i ) )
[m ^-3 ] is the 3-D Dirac delta distribution. Each scatterer’s volume depends on the bubble radius, R , as a function of time, which in our case was calculated by solving the Marmottant equation [ , ].

In the case of a population of monodisperse MBs, the equilibrium radius, R ₀, is the same for all the scatterers, whereas for a polydisperse distribution, each scatterer has its own R ₀ , 𝑅(𝑖)0
R 0 ( i )
.

Configurations used in the simulations

Simulation of pressure fields

Incident field and contrast domain

To study the influence of different populations of MBs on a propagating plane wave, we considered the computational domain and the domain for the contrast media, as depicted in Figure 1 a. This configuration is used in sections ‘Comparison of INCS and effective medium theory’ and ‘Plane wave: monodisperse vs polydisperse populations’ for INCS validation and the comparison between different populations, respectively. Computational domain dimensions X × Y × Z = 20 mm × 20 mm × 30 mm were used. The scatterers were placed in a domain with dimensions X × Y × Z =15 mm × 15 mm × 4.444 mm, resulting in 1 mL volume. These configuration choices were made to simplify the comparison between polydisperse and monodisperse populations.

The incident pressure field is a plane wave being generated at z = 0 and propagating in the positive z-direction. A plane wave is used to initially excite all the scatterers by an incident wave that has the same pressure amplitude everywhere. The temporal signature of the incident pressure is given as (eqn [ 5 ]):

𝑝(𝑡)=𝑃0exp[−(𝑡−𝑇𝑑𝑇𝑤/2)2]sin[2𝜋𝑓0(𝑡−𝑇𝑑)],
p(t)=P0exp[−(t−TdTw/2)2]sin[2πf0(t−Td)],

where

𝑇𝑤=3/𝑓0
T w = 3 / f 0
is the width and

𝑇𝑑=6/𝑓0
T d = 6 / f 0
is the delay of a Gaussian envelope with a duration of

12/𝑓0
12 / f 0
, where

𝑓0=1
f 0 = 1
MHz is the center frequency. Therefore, this is a narrow-band pulse. The peak pressure is P ₀ = 200 kPa. The scatterers were embedded in water with a density of ρ ₀ = 1060 kg/m ³ and a speed of sound of c ₀ = 1482 m/s. In the considered situations, water had negligible losses and non-linear effects were hardly noticeable. Therefore, we assumed that the embedding medium was lossless and linear. A sampling frequency of 18 MHz was used as the basis for the discretization of the spatiotemporal domain [ ].

Configuration for validation

To validate INCS, we compared our results with those following effective medium theory. Analytical expressions describing the effective behavior of a population of isotropic linear scatterers (LSs) were derived from Foldy [ , ]. A similar validation was used in a previous publication for a monodisperse population of scatterers [ ], but here we considered a polydisperse population. For INCS implementation, we assumed that the contrast source term for each LS is given by eqn (6) :

𝑆sc(𝑥,𝑡)=−𝑓(𝑅0)𝑉0𝜌0𝜌1𝑐21𝜕2𝑝(𝑥sc,𝑡)𝜕𝑡2𝛿(𝑥−𝑥sc),
Ssc(x,t)=−f(R0)V0ρ0ρ1c12∂2p(xsc,t)∂t2δ(x−xsc),

where R ₀ is the equilibrium radius, V ₀ is its initial volume, ρ ₁ is the density of mass of the gas inside the LS, c ₁ is the speed of sound inside the LS and f ( R ₀ ) is the polydispersity co-efficient given by eqn (7) :

𝑓(𝑅0)=𝐴0(𝑅0/𝑅0,ref)𝛾.
f(R0)=A0(R0/R0,ref)γ.

The constant A ₀ was used to adjust the scattering strength if necessary and γ was the polydispersity scale parameter to control the scattering distribution of the population.

In the case of a plane wave excitation as in eqn (5) , the scattered pressure was given by eqn (8) [ ]: