Microbubble physics, contrast imaging and the ULM signal source
A contrast microbubble is a gas spring with a liquid mass, tuned to megahertz by its size and stiffened by its shell. Everything you do in CEUS and ULM follows from that oscillator: why bubbles outshine tissue, how pulse sequences separate them, how time–intensity curves measure perfusion, and why ULM must trade bubble density against acquisition time.
You will be able to
- Derive the Minnaert resonance, correct it for surface tension, heat exchange and a lipid shell, and quote \(f_0R_0\) with its assumptions
- Integrate the Rayleigh–Plesset equation with a Marmottant shell (RK4) and read harmonics, subharmonics and compression-only behaviour from the emitted spectrum
- Predict what pulse inversion, amplitude modulation and their combination keep and reject, including the residues from tissue harmonics, motion and amplitude errors
- Fit bolus and flash-replenishment time–intensity curves in the right (linear) domain
- Size the ULM concentration from the PSF, and estimate how long capillaries take to fill
Assumes
- CEUS practice (low-MI imaging, agents, perfusion curves)
- Module 06 · Eigendecomposition & SVD: SVD clutter filtering; Module 12 · Estimation theory & the Cramér–Rao bound: ULM localization bounds
- Module 15 · Wave physics, nonlinear propagation and harmonic imaging: harmonics and nonlinear propagation; Module 22 · Acoustic output, bioeffects and regulatory measurement: MI and derating
28.1The agent: gas core, shell, size, and why compressibility wins
Clinical agents are suspensions of gas microbubbles, 1–10 µm in diameter, stabilized by a monolayer of phospholipid or a thin shell of denatured albumin. The gas is a heavy, poorly soluble fluorinated gas because an uncoated micron bubble does not last: the Laplace pressure \(2\sigma/R\) (144 kPa at \(R\) = 1 µm in water) pushes gas into solution, and the Epstein–Plesset estimate for a 1 µm air bubble in air-saturated water is of order milliseconds. Low solubility and diffusivity of the core gas and a shell that cancels most of the surface tension stretch that to minutes.
| Agent (US label) | Gas | Shell | Mean diameter | Bubbles per mL |
|---|---|---|---|---|
| Lumason / SonoVue | SF6 | phospholipid (DSPC, DPPG) with palmitic acid | 1.5–2.5 µm | 1.5–5.6 × \(10^{8}\) |
| Definity | C3F8 (perflutren) | phospholipid | 1.1–3.3 µm | up to 1.2 × \(10^{10}\) |
| Optison | C3F8 | human albumin | 3.0–4.5 µm | 5–8 × \(10^{8}\) |
| Sonazoid (not US-approved) | C4F10 | phospholipid | see the regional label | |
Numbers are from the US prescribing information (mean diameters are number-weighted; the acoustic response is dominated by the larger bubbles, see the cards below). A 2.4 mL Lumason dose at \(3\times 10^{8}\) per mL diluted in about 5 L of blood gives of order \(10^{5}\) bubbles per mL (140 per µL) at first pass, before lung filtering and destruction. That number reappears in 28.8.
Why a bubble is so bright. A sphere of radius \(a \ll \lambda\) in a fluid scatters through two channels: its compressibility contrast drives a monopole (volume pulsation) and its density contrast drives a dipole (rigid-body oscillation). For an incident plane wave the Rayleigh cross-section is (28.1)
- \(\sigma_s\)
- scattering cross-section, m²
- \(k\)
- wavenumber \(2\pi/\lambda\), 1/m
- \(a\)
- sphere radius, m
- \(\kappa,\ \kappa_s\)
- compressibility of the fluid and of the sphere, 1/Pa
- \(\rho,\ \rho_s\)
- density of the fluid and of the sphere, kg/m³
The ratio bubble to rigid sphere of the same size is \((4/7)(15\,400)^{2} \approx 1.4\times 10^{8}\), and against a tissue inhomogeneity with 5 % compressibility and density contrast it is about \(10^{11}\). That is far below resonance. Near resonance the bubble gains another factor \(1/\delta^{2}\) (28.2, (28.5)). One bubble is therefore visible on its own in a beamformed image, which is the physical basis of ULM.
Tissue scatters because it is a little stiffer or denser than its surroundings. A bubble scatters because it is about 15 000 times softer: the same incident pressure makes it change volume 15 000 times more than the water it replaces, and a pulsating volume is a monopole source.
28.2The bubble as a mass on a gas spring: Minnaert and beyond
Mass. Incompressible radial flow around a sphere has velocity \(u(r) = R^{2}\dot R/r^{2}\). Its kinetic energy is \(\tfrac12\int \rho u^{2}\, 4\pi r^{2}\, dr = \tfrac12 (4\pi R^{3}\rho)\, \dot R^{2}\), so the effective mass is \(m = 4\pi R_0^{3}\rho\), three times the mass of the displaced liquid. Spring. A polytropic gas gives \(p_g = P_0(R_0/R)^{3\kappa}\), so a radial displacement \(x\) changes the pressure by \(-3\kappa P_0 x/R_0\) over the area \(4\pi R_0^{2}\): the stiffness is \(s = 12\pi\kappa P_0 R_0\); hence (28.2).
- \(\omega_0,\ f_0\)
- angular and ordinary resonance frequency, rad/s and MHz
- \(R_0\)
- equilibrium bubble radius, µm
- \(\kappa\)
- polytropic exponent of the gas (dimensionless)
- \(P_0\)
- ambient pressure, Pa
- \(\rho\)
- liquid density, kg/m³
- \(s,\ m\)
- stiffness, N/m, and effective mass, kg
Three corrections matter at micron sizes, and all three are larger than the difference between 3.26 and 3.29.
- Surface tension. The gas sits at \(P_0 + 2\sigma/R_0\), and the tension term itself responds to \(R\). Linearizing gives \(\rho\omega_0^{2}R_0^{2} = 3\kappa P_0 + 2(3\kappa - 1)\sigma/R_0\). A free 1 µm air bubble resonates at 4.74 MHz, not 3.29 (+44 %); at 2 µm, 2.04 versus 1.64 MHz.
- Heat exchange. The gas is adiabatic only if the thermal diffusion length \(\sqrt{D_{\mathrm{th}}/2\omega}\) is much smaller than \(R_0\). For air at 3 MHz it is 0.73 µm, comparable to the radius, so \(\kappa\) lies between 1 and \(\gamma\) and heat conduction adds a thermal damping (Prosperetti's analysis). For contrast agents the polytropic exponent is usually taken close to the heat-capacity ratio of the heavy gas (≈ 1.07–1.1), which is already near 1.
- Shell. In the Marmottant model the elastic shell makes the surface tension radius dependent, \(\sigma(R) = \chi(R^{2}/R_b^{2} - 1)\). Its slope adds stiffness (28.3):
- \(\sigma_0\)
- \(\sigma(R_0)\), shell surface tension at rest, N/m
- \(\chi\)
- shell elasticity, N/m
- \(\rho,\ \omega_0,\ R_0,\ \kappa,\ P_0\)
- as in (28.2)
Damping and scattering. Three losses act on the linear oscillator. With \(\delta = b/\omega\) the dimensionless damping (\(Q = 1/\delta\) at resonance) (28.4)–(28.5):
- \(\delta\)
- dimensionless damping \(b/\omega\) (\(Q = 1/\delta\) at resonance)
- \(c\)
- speed of sound in the liquid, m/s
- \(\mu\)
- liquid viscosity, Pa·s
- \(\kappa_s\)
- shell dilatational viscosity, kg/s
- \(\omega\)
- angular drive frequency, rad/s
- \(R_0,\ \rho\)
- as in (28.2)
- \(\delta_{\mathrm{rad}},\ \delta_{\mathrm{vis}},\ \delta_{\mathrm{shell}}\)
- radiation, liquid-viscosity and shell-viscosity damping
- \(\sigma_s(f)\)
- scattering cross-section at frequency \(f\), m²
- \(f_0\)
- resonance frequency, MHz
- \(\delta\)
- total damping
- \(R_0\)
- equilibrium bubble radius, µm (as in (28.2))
- \(f\)
- driving frequency, MHz
| \(R_0\) | Model | \(f_0\) (MHz) | \(\delta_{\mathrm{rad}}\) | \(\delta_{\mathrm{vis}}\) | \(\delta_{\mathrm{shell}}\) | Q | \(\sigma_s/\pi R_0^{2}\) at \(f_0\) |
|---|---|---|---|---|---|---|---|
| 1 µm | free air | 4.74 | 0.020 | 0.134 | 0 | 6.5 | 168 |
| 1 µm | lipid shell | 7.95 | 0.034 | 0.080 | 0.401 | 1.9 | 15 |
| 2 µm | free air | 2.04 | 0.017 | 0.078 | 0 | 10.5 | 440 |
| 2 µm | lipid shell | 2.99 | 0.025 | 0.053 | 0.133 | 4.7 | 89 |
| 3 µm | free air | 1.28 | 0.016 | 0.055 | 0 | 14.0 | 779 |
| 3 µm | lipid shell | 1.72 | 0.022 | 0.041 | 0.069 | 7.6 | 231 |
Water: \(\rho\) = 1000 kg/m³, \(c\) = 1480 m/s, \(\mu\) = 1 mPa·s, \(\sigma_w\) = 0.072 N/m. Shell: \(\kappa\) = 1.07, \(\chi\) = 0.5 N/m, \(\sigma_0\) = 0.02 N/m, \(\kappa_s = 5\times 10^{-9}\) kg/s, inside the range fitted to phospholipid agents (Marmottant et al. used \(\chi\) = 1 N/m and \(\kappa_s = 1.5\times 10^{-8}\) kg/s for SonoVue). Fitted shell parameters differ between studies and between models; do not mix \(\chi\) from one model with \(\kappa_s\) from another.
Resonance explorer: size, gas, tension and shell
Approx. modelSmall-amplitude linear oscillator: Minnaert (exact for its assumptions), Laplace-corrected free bubble and linearized Marmottant shell; thermal damping omitted; rigid sphere from the Rayleigh formula. Shell values are illustrative.Resonant radii at clinical frequencies (1.5–3 µm at 2–5 MHz with a shell) sit inside the agent size distributions, which is no accident: agents are sized to pass the pulmonary capillaries (diameters a few µm) and to resonate in the abdominal and cardiac band. High-frequency preclinical imaging (15–40 MHz) is far above resonance for most of the population, so the echo is weaker and more linear per bubble, and nonlinear contrast modes lose much of their advantage unless the agent is size-selected.
"\(f_0 \approx 3/R_0\)" is for a free air bubble with no surface tension. For a 1 µm free bubble the tension correction is +44 %, and for a shelled bubble the shell term can triple the stiffness. Quote which model and which constants you used.
28.3The Rayleigh–Plesset equation and shell models
Large oscillations need the nonlinear equation. Integrate the radial momentum equation of the incompressible flow \(u = R^{2}\dot R/r^{2}\) from the bubble wall to infinity (where the pressure is \(P_0 + P_{\mathrm{ac}}(t)\)), and use the normal-stress balance at the wall, \(p_{\mathrm{liquid}}(R) = p_g - 2\sigma/R - 4\mu\dot R/R\), to obtain (28.6)–(28.8):
- \(R,\ \dot R,\ \ddot R\)
- bubble radius (m), wall velocity (m/s) and acceleration (m/s²)
- \(p_g(R)\)
- gas pressure, Pa
- \(\sigma(R)\)
- radius-dependent surface tension, N/m
- \(\mu\)
- liquid viscosity, Pa·s
- \(\kappa_s\)
- shell dilatational viscosity, kg/s
- \(P_{\mathrm{ac}}(t)\)
- acoustic drive pressure, Pa
- \(\kappa\)
- polytropic exponent of the gas, as in (28.2)
- \(\rho,\ P_0,\ c\)
- liquid density, ambient pressure, sound speed in the liquid
- \(t\)
- time, s
- \(p_g(R)\)
- gas pressure inside the bubble of radius \(R\), Pa
- \(P_0,\ \sigma_0,\ R_0,\ \kappa\)
- as in (28.2) and (28.3)
- \(R\)
- instantaneous bubble radius, µm
- \(R < R_b\)
- buckled
- \(R_b \le R \le R_r\)
- elastic
- \(R > R_r\)
- ruptured
- \(R_b,\ R_r\)
- buckling and rupture radii, µm
- \(\sigma_w\)
- surface tension of water, N/m (0.072)
- \(R\)
- instantaneous bubble radius, µm
- \(\chi\)
- shell elasticity, N/m, as in (28.3)
- \(\sigma(R)\)
- shell surface tension, N/m
Earlier shell models treat the coating as a thin viscoelastic layer of fixed properties: de Jong and colleagues added a shell stiffness and friction term to the linear oscillator (1992), and Church (1995) derived a Rayleigh–Plesset form for a shell of finite thickness with shear modulus and viscosity. They describe albumin and hard shells well at small amplitude. Marmottant's model adds what lipid monolayers do at large amplitude: buckle in compression, rupture in expansion.
What the probe receives. The bubble is a monopole. Far from it (\(r \gg R\)) the velocity potential \(\varphi = -R^{2}\dot R/r\) gives (28.9)
- \(p_s\)
- emitted (scattered) pressure, Pa
- \(r\)
- distance from the bubble centre, m
- \(V\)
- bubble volume, m³
- \(\rho,\ R,\ \dot R,\ \ddot R\)
- as in (28.6)
- \(t\)
- time, s
- \(d^2/dt^2\)
- second time derivative
Numerical integration. State \((R, \dot R)\), classical fourth-order Runge–Kutta, 200 steps per acoustic period, with each step split further whenever \(R\) would move more than 0.5 % (the violent collapses). For a free 2 µm bubble at 2 MHz and 200 kPa, \(R_{\mathrm{min}}/R_0\) = 0.4163 at 200 steps per period and 0.415 at 400 and 800: the default is converged. The demo runs the same solver as m28_rayleigh_plesset.m.
Rayleigh–Plesset solver (RK4): radius, emitted pressure and spectrum
SimulationRK4 Rayleigh–Plesset with Marmottant shell, JS in page (same as m28_rayleigh_plesset.m); first-order radiation damping, polytropic gas, no thermal damping or lipid shedding; illustrative shell \(\chi\) = 0.5 N/m, \(\sigma_0\) = 0.02 or 0 N/m, \(\kappa_s = 5\times 10^{-9}\) kg/s, \(\kappa\) = 1.07.28.4Nonlinear signatures: harmonics, subharmonics, compression-only behaviour, destruction
Harmonics. The gas spring is asymmetric: \((R_0/R)^{3\kappa}\) stiffens without bound in compression and softens in expansion, and the liquid inertia makes the collapse faster than the growth. Near resonance the radial excursion is amplified by \(Q\), so a bubble becomes nonlinear at incident pressures where tissue propagation is still nearly linear. In the demo a free 2 µm bubble driven at its resonance has a second harmonic of −13 dB already at 10 kPa (MI 0.007; 20-cycle burst); the tissue second harmonic after 40 mm at 100 kPa and 2.5 MHz is −28 dB (Module 15).
Subharmonics and ultraharmonics. Above a threshold the oscillation doubles its period and emits at \(f/2\) (and \(3f/2\), \(5f/2\)). For a bubble the threshold is lowest when the drive is near twice the resonance, because the \(f/2\) response then falls on the resonance. For the free 2 µm bubble at \(2f_0\) the model gives \(f/2\) at −43, −30, −13 and −4 dB for 50, 100, 150 and 200 kPa: a threshold between about 100 and 175 kPa. A 1 µm bubble, with almost twice the viscous damping, stays below −25 dB up to 200 kPa. Tissue propagation creates no subharmonic, so subharmonic imaging has a very high contrast-to-tissue ratio, and because the subharmonic amplitude of some agents depends on ambient pressure it has been studied for noninvasive pressure estimation. Phospholipid bubbles show subharmonics at lower pressures than the free-bubble model predicts; single-bubble optical measurements found subharmonic thresholds as low as 5 kPa near twice the resonance, far below viscoelastic-shell predictions, and attribute them to the buckling regime of the shell (Sijl et al., JASA 2010).
Compression-only behaviour. A lipid bubble near the buckling radius has no tension in compression (\(\sigma = 0\), soft) but a rising tension in expansion (stiff). It therefore compresses much more than it expands, an effect observed optically by de Jong et al. (UMB 2007). In the model at 3 MHz the buckled bubble's compression/expansion ratio is 3.0 at 10 kPa and 2.5 at 50 kPa, against 1.1–1.4 for the elastic one. The asymmetry is an even-order nonlinearity at very low pressure: exactly what pulse inversion and amplitude modulation detect.
Destruction. As pressure rises the bubble sheds lipid, fragments, loses gas by diffusion, or collapses inertially. Rayleigh–Plesset cannot model fragmentation, but \(R_{\mathrm{max}}/R_0 \ge 2\) is the classic rough criterion for inertial collapse (Flynn); the analysis that led to MI (Apfel and Holland, UMB 1991) used a related collapse-temperature threshold of 5000 K. Clinical practice follows the same logic: real-time CEUS runs at low MI, below about 0.3 and often 0.05–0.1 on modern systems (EFSUMB guidance), to keep bubbles intact and nonlinear; a few high-MI "flash" frames destroy them on purpose for replenishment studies. MI is defined from the derated peak rarefactional pressure of the beam (Module 22), not from the pressure at a given bubble, so the bubble's local pressure depends on depth, attenuation and focusing.
28.5Contrast-specific pulse sequences
Write the echo of a scatterer as a power series in the transmit amplitude \(a\): \(e(a) = a_1 a\,p + a_2 a^{2}p^{2} + a_3 a^{3}p^{3} + \dots\), where the \(a_n\) are operators (filters) and \(p^{n}\) carries energy near \(n f_0\) (and, for odd \(n\), back at \(f_0\)). Linear tissue has only \(a_1\). A multi-pulse sequence transmits several scaled or inverted copies and adds the echoes with weights that cancel \(a_1\):
| Sequence | Transmit amplitudes | Receive sum | What survives | Band of the residue |
|---|---|---|---|---|
| Pulse inversion (PI) | +1, −1 | \(e(1) + e(-1)\) | \(2a_2p^{2}\) (all even orders) | \(2f_0\) (and DC) |
| Amplitude modulation (AM, power modulation) | 1, ½ | \(e(1) - 2e(\tfrac12)\) | \(\tfrac12 a_2p^{2} + \tfrac34 a_3p^{3}\) | \(2f_0\) and \(f_0\) (from \(p^{3}\)) |
| AM + PI (CPS-type) | 1, −½ | \(e(1) + 2e(-\tfrac12)\) | \(\tfrac32 a_2p^{2} + \tfrac34 a_3p^{3}\) | \(f_0\) and \(2f_0\) |
PI (Simpson, Chin and Burns, IEEE TUFFC 1999) needs no filter to separate harmonic from fundamental, so it keeps the full transducer bandwidth, but its signal lives at \(2f_0\), where attenuation is higher. AM keeps the odd-order "nonlinear fundamental", which propagates and is received at \(f_0\): better penetration, at the price of an exact half-amplitude transmit (systems commonly fire half of the elements to get it). Vendors combine amplitude and phase steps over three or more pulses (the "contrast pulse sequence" modes of commercial scanners); the algebra is always the same: choose weights that null \(a_1\) and keep the orders you want.
Three things leave tissue in a contrast image, and the demo models each:
- Nonlinear propagation. The incident pulse at depth already contains a tissue-generated second harmonic. To second order \(p_{\mathrm{inc}} = a\,p_1 + a^{2}\,\frac{\beta x}{2\rho c^{3}}\,\frac{d(p_1^{2})}{dt}\), which survives PI (it is even in \(a\)) and is attenuated four times in the half pulse of AM. At 100 kPa, 2.5 MHz, \(x\) = 40 mm, \(\beta\) = 4.5 the CW ratio \(B_2/B_1 = \sigma/2\) is −28 dB.
- Motion. A scatterer moving at \(v\) between pulses shifts by \(\tau = 2v/(c\cdot\mathrm{PRF})\). For a linear scatterer the PI residue is \(\lvert 1 - e^{-j2\pi f\tau}\rvert = 2\lvert\sin(\pi f\tau)\rvert\) of the echo: 10 mm/s at PRF 5 kHz and 2.5 MHz gives −27.8 dB, at 5 MHz −21.8 dB.
- Amplitude and phase errors. AM needs the half pulse to be exactly half: a gain error \(e\) (dB) leaves \(\lvert 1 - 10^{e/20}\rvert\) of the linear echo, −38.7 dB for 0.1 dB and −24.5 dB for 0.5 dB.
Pulse inversion, amplitude modulation and AM + PI: what survives for a bubble and for tissue
SimulationRK4 bubble echo (same shell as Demo 28.2) and a point tissue scatterer insonified by the second-order plane-wave incident pulse; FFT band energies; JS (same as m28_contrast_pulse_sequences.m)."PI removes tissue" is true only for linear tissue at rest. Near the focus at moderate MI, tissue harmonics from propagation survive PI exactly like bubble harmonics do; the remedies are lower MI, AM (its tissue-harmonic residue is a quarter of PI's, 12 dB lower) or a sequence that uses the bubble's odd orders. And every multi-pulse method is a high-pass filter in slow time: motion between pulses sets the floor.
28.6Perfusion quantification: low-MI imaging and time–intensity curves
Real-time low-MI imaging follows the agent through the microcirculation. Two acquisition schemes give two kinds of time–intensity curve (TIC) from a region of interest:
- Bolus. The ROI intensity rises and falls as the bolus passes. Models: lognormal (below), gamma variate \(I(t) = K(t - t_0)^{\alpha}e^{-(t - t_0)/b}\) (peak at \(t_0 + \alpha b\), mean transit \((\alpha + 1)b\) after \(t_0\)), local density random walk. Descriptors: peak enhancement (PE), time to peak (TTP), rise time, wash-in rate (WiR, maximum slope), area under the curve (AUC), mean transit time (MTT).
- Flash-replenishment. During a constant infusion, a few high-MI frames destroy the bubbles in the slice; the low-MI frames that follow record the refill. Wei et al. (Circulation 1998) modelled it as \(I(t) = A\left(1 - e^{-\beta t}\right)\): the plateau \(A\) is proportional to the microvascular blood volume, the rate \(\beta\) to the mean velocity divided by the elevational thickness of the destroyed slice (the slice must be refilled from its sides), and \(A\beta\) to flow.
- \(I(t)\)
- linear (power) intensity in the ROI
- \(\mathrm{AUC}\)
- area under the curve
- \(t_0\)
- arrival time, s
- \(\mu,\ s\)
- lognormal location (of \(\ln\) s) and shape
- \(I_{\mathrm{bg}}\)
- background intensity
- \(t\)
- time after injection, s
- \(\ln,\ \exp\)
- natural logarithm, exponential
- \(\mathrm{MTT},\ \mathrm{TTP}\)
- mean transit time and time to peak, s
- \(\mathrm{PE}\)
- peak enhancement
- \(\mu,\ s,\ t_0,\ \mathrm{AUC}\)
- as in (28.10)
- \(e\)
- Euler’s number
Fit linear intensity. All of these models ((28.10)–(28.11)) describe a quantity proportional to bubble concentration, which is the linear (power) intensity, not the log-compressed display. A monotone map preserves the peak time but changes the shape: fitting the same lognormal model to 40 dB log-compressed gray levels doubles the apparent MTT (×2.2 in m28_tic_perfusion.m) and moves TTP by 5 %. Quantification software either works on linear data or inverts the display map. Within a small ROI the intensity also carries speckle: averaging \(N_s\) independent speckle cells leaves a multiplicative fluctuation of \(1/\sqrt{N_s}\) (14 % for 50 cells).
Other sources of variability are larger than the fitting error: attenuation and shadowing by upstream bubbles, nonlinear propagation, bubble destruction by the imaging frames themselves, depth-dependent pressure, injection technique and cardiac output (reviewed by Tang et al., Interface Focus 2011). Ratios to a reference region and fixed system settings matter more than the choice of model.
TIC fitting: bolus and flash-replenishment, linear versus log-compressed
SimulationMonte Carlo speckle (gamma, shape \(N_s\)) on lognormal or Wei model curves, fitted by Levenberg–Marquardt, JS (same as m28_tic_perfusion.m). The true curves are models, not patient data.28.7Separating bubbles from tissue in slow time
Ultrafast ULM usually separates bubbles from tissue with the spatiotemporal SVD filter (Module 06; Demené et al., IEEE TMI 2015): stack \(N_t\) frames as a Casorati matrix, discard the largest singular components (tissue: spatially coherent, slowly varying, 30–60 dB brighter), keep the intermediate ones (moving bubbles), discard the smallest (noise). It exploits motion, not bubble physics. Two consequences matter for ULM design:
- Slow bubbles look like tissue. A bubble moving axially at \(v\) has Doppler \(2vf/c\): 6.5 Hz at 1 mm/s and 5 MHz, the same Doppler as tissue moving at 1 mm/s, which is ordinary physiological motion. The SVD cut that removes tissue removes capillary bubbles too. Nonlinear sequences (AM, PI) remove tissue by physics and do not care about velocity, at the price of two or three transmits per frame and motion sensitivity (28.5).
- Bubbles decorrelate. Between frames a bubble moves \(v/\mathrm{FR}\) (10 µm at 10 mm/s and 1 kHz), its amplitude changes as it crosses the beam and its echo flickers with its radial dynamics. Tracking needs displacements per frame well below the PSF size; very fast bubbles at low frame rates break the linking, and very slow ones sit inside the clutter band.
28.8ULM: the bubble as a point source, and the concentration trade-off
A point source with an unknown, flickering amplitude. A 2 µm bubble is \(\lambda/150\) of a 5 MHz wavelength, so its image is the system PSF times a complex amplitude. Localization precision is then set by the PSF width and the energy-to-noise ratio, \(\sigma_x \ge w_x/\sqrt{\mathrm{ENR}}\) (Module 12). The bubble physics enters through the amplitude: \(\sigma_s\) varies by orders of magnitude across the size distribution at a given frequency (28.1, 28.2), so the ENR, and therefore the precision, is a distribution, dominated by near-resonant and large bubbles. Choosing \(f\) for ULM trades PSF width (\(\propto 1/f\)) against bubble response (resonant population at lower \(f\)) and attenuation.
The density trade-off. Localization needs isolated bubbles. If bubbles in the image plane are a Poisson process with areal density \(n\), and two bubbles closer than \(k\) FWHM (elliptical metric) cannot be localized separately, then (28.12)
- \(P_{\mathrm{iso}}\)
- probability that a bubble has no neighbour within the exclusion area
- \(n\)
- areal bubble density in the image, per mm²
- \(A\)
- exclusion area, mm²
- \(k\)
- separation needed, in FWHM (elliptical metric)
- \(F_x,\ F_z\)
- lateral and axial FWHM, mm
- \(e\)
- Euler’s number
- \(n^{*}\)
- density that maximizes the isolated-bubble yield, per mm²
The density in the image is \(n = C\cdot(\text{slice thickness})\cdot(\text{vascular fraction})\). With \(C\) = 140 per µL (the first-pass estimate in 28.1), a 0.5 mm slice and 5 % blood volume: \(n\) = 3.5 per mm², \(nA\) = 0.78, \(P_{\mathrm{iso}}\) = 0.46. A clinical bolus at 5 MHz is already near the overlap optimum on average, and well above it in large vessels, where the local density is twenty times higher.
The time trade-off. A capillary is mapped only when bubbles pass through it. With diameter \(d\) and speed \(v\) it receives \(C\cdot(\pi d^{2}/4)\cdot v\) bubbles per second: for \(d\) = 8 µm, \(v\) = 1 mm/s and \(C\) = 140 per µL, one bubble every 142 s. Lowering \(C\) to reduce overlap lengthens this wait in proportion. This is the argument of Hingot et al. (Sci. Rep. 2019): microvascular flow, not frame rate, sets the acquisition time for full capillary maps, typically minutes.
ULM: bubble density, PSF overlap and capillary waiting time
SimulationPoisson Monte Carlo of bubble positions with a complex Gaussian PSF (JS); \(P_{\mathrm{iso}} = e^{-nA}\) and the waiting time are exact for that model; uniform density, no vessels.Safety and the MI budget in ULM
ULM needs the same bubbles to survive hundreds to thousands of frames, so it runs at low MI, and plane-wave sequences with no focus keep the peak pressure low. The time-averaged exposure is another matter: 0.5–1 kHz compounded frame rates raise \(I_{\mathrm{SPTA}}\) and TI compared with conventional imaging (Module 22). Contrast agents lower the pressure threshold for cavitation-related bioeffects; animal studies have shown microvascular effects when bubbles are destroyed at high MI (for example Skyba et al., Circulation 1998). The design rule is to keep the MI at the bubble below the destruction threshold of the agent and to report MI, frame rate, duration and agent dose together.
Where ULM innovation lives
| Bottleneck | Physics from this module | Levers people pull |
|---|---|---|
| Acquisition time (minutes) | Capillary flux \(C\cdot Q\); \(n^{*} = 1/A\) | Localize at higher density (overlap-tolerant or sparse estimators, learned localizers), 3D or smaller PSF to shrink \(A\) |
| Slow-flow detection | SVD separates by Doppler; 1 mm/s ↔ 6.5 Hz at 5 MHz | Nonlinear sequences (AM, PI) before or instead of SVD; adaptive cut-offs |
| Precision spread | \(\sigma_s(R, f)\) spans orders of magnitude; ENR follows it | Frequency matched to the size distribution, size-selected agents, amplitude-weighted maps |
| Elevation overlap | \(A\) in 2D is a projection of a 3D slice | Matrix and row–column arrays (Module 19); volumetric ULM |
| Tracking failures | Displacement \(v/\mathrm{FR}\) versus PSF; amplitude flicker | Higher frame rate, motion models, Kalman or graph-based linking |
| Exposure | MI at the bubble, \(I_{\mathrm{SPTA}}\) at kHz frame rates | Low-MI diverging or plane waves, duty-cycle limits, coded transmits (Module 10) |
28.9Visual summary
Thirteen pictures that carry the module. Each figure is computed from a formula or from the RK4 solver in the text (water constants of 28.2; illustrative shell \(\chi\) = 0.5 N/m, \(\sigma_0\) = 0.02 N/m, \(\kappa_s = 5\times 10^{-9}\) kg/s unless stated); the badge on each card says which. "Approx. model" here means the small-amplitude linear oscillator or, for the PI and AM cards, a memoryless polynomial echo.
28.10How experts think about it
An expert pictures the bubble as a damped oscillator whose spring is gas plus shell and whose mass is three times the displaced liquid: resonance from size and stiffness, Q from shell viscosity, brightness from compressibility. Nonlinearity is the asymmetry of that spring (and of the lipid shell's buckling), and every contrast sequence is a choice of weights that cancels the linear term. For ULM the same expert thinks in densities and fluxes: the PSF sets how many bubbles fit in a frame, the vasculature sets how fast new ones arrive.
The radial flow carries three times the displaced liquid's mass; the gas gives a stiffness \(12\pi\kappa P_0R_0\). Then \(\omega_0^{2} = 3\kappa P_0/(\rho R_0^{2})\) (28.2): 3.29 MHz·µm for air (\(\kappa\) = 1.4, 101.3 kPa, \(\rho\) = 998), 3.26 with 1 bar and \(\rho\) = 1000, 2.78 isothermal. A 2 µm bubble: 1.64 MHz; 1 µm: 3.29 MHz.
Linearizing the Laplace term gives \(\rho\omega_0^{2}R_0^{2} = 3\kappa P_0 + 2(3\kappa - 1)\sigma/R_0\). For air with \(\sigma\) = 0.072 N/m the correction is +44 % at 1 µm, +24 % at 2 µm and +10 % at 5 µm. Heat exchange works the other way: for an isothermal gas the Minnaert value drops 16 %. At MHz the thermal length (0.73 µm for air at 3 MHz) is comparable to \(R_0\), so neither limit is exact.
The shell term grows as \(1/R_0\) while the gas term is constant, so the shift is largest for small bubbles: at 2 µm, \(4\chi/R_0\) = 1.0 MPa against \(3\kappa P_0\) = 0.33 MPa (28.3), and \(f_0\) rises from 1.44 to 2.99 MHz (with \(\sigma_0\) = 0.02 N/m). At 1 µm the shelled bubble resonates at 7.95 MHz. This is why agents of nominal 1–3 µm still answer at 2–8 MHz.
Radiation damping (\(\omega R_0/c \approx 0.02\)) is small; liquid viscosity scales as \(1/R_0^{2}\) and shell viscosity as \(1/R_0^{3}\) (28.4). With \(\kappa_s = 5\times 10^{-9}\) kg/s the shell term is 0.13 at 2 µm and 0.40 at 1 µm. A low \(Q\) broadens the response, so a polydisperse agent responds over a wide band, and it shortens ring-down (good for axial resolution).
Both scale as \(k^{4}a^{6}\) (28.1), but the bubble's monopole term is \(((\kappa_s - \kappa)/\kappa)^{2} \approx 2.4\times 10^{8}\) against 7/4 for a rigid fixed sphere. The resonance formula (28.5) gives the same low-frequency limit, \(4\pi R_0^{2}(f/f_0)^{4}\). Because of the \(R^{6}\) law the largest bubbles of a population dominate the backscatter far below resonance.
The resonant response is limited only by damping. A free 2 µm bubble (\(\delta\) = 0.095) scatters like a disc of radius \(21\,R_0\) = 42 µm; the coated one (\(\delta\) = 0.21) like \(9.4\,R_0\). Above resonance \(\sigma_s \to 4\pi R_0^{2}\), so the advantage over tissue shrinks at high frequency for a given size.
From \(\varphi = -R^{2}\dot R/r\) and \(p = -\rho\,\partial\varphi/\partial t\) in the far field (28.9). Linear and at resonance the pressure ratio is \((R_0/r)/\delta = 2.1\times 10^{-3}\) at 1 cm for a free 2 µm bubble. In the nonlinear regime the collapse makes \(\ddot R\) spike, which is why the emitted pressure is a train of short pulses rich in harmonics while \(R(t)\) still looks smooth.
Tissue harmonics grow with pressure and distance, \(B_2/B_1 \approx \beta\omega P x/(2\rho c^{3})\). A resonant bubble's harmonic content is already high at small pressure because \(Q\) amplifies the excursion, and the lipid shell's buckling adds an even-order nonlinearity at 10 kPa. In the figure (RK4, coated 2 µm bubble at 3 MHz) H2 stays between −15 and −21 dB from 10 to 100 kPa, while tissue at 40 mm and the same frequency rises from −47 to −27 dB.
Period doubling is a parametric instability: the \(f/2\) response is strongest when it falls on the bubble's own resonance. The figure sweeps \(f/f_0\) at 150 kPa for the free 2 µm bubble: the \(f/2\) line peaks at \(f/f_0 \approx 1.9\) (the large-amplitude resonance sits slightly below the linear one) and is 30–40 dB lower at 1.6 and 2.4. Because tissue propagation generates no \(f/2\), a subharmonic image is bubble-only.
With \(\sigma(R) = 0\) in compression and \(\chi(R^{2}/R_b^{2} - 1)\) in expansion (28.8) the spring is soft one way and stiff the other. In the RK4 model at 3 MHz a buckled 2 µm bubble compresses 3.0 times further than it expands at 10 kPa and 2.5 times at 50 kPa; the elastic bubble only 1.1–1.4 times. The asymmetry is an even-order (\(a_2\)) response that PI and AM detect at very low MI.
Substitute \(a = \pm 1,\ \tfrac12,\ -\tfrac12\) into \(e(a) = \sum a_n a^{n} p^{n}\) and apply the receive weights. The \(a_1\) term cancels in all three. Only AM-type sequences keep the cubic term, whose energy at \(f_0\) propagates with fundamental-band attenuation: that is why AM penetrates deeper than PI for the same bubbles.
Two echoes of a linear scatterer that has moved by \(v/\mathrm{PRF} = c\tau/2\) between pulses differ by a delay \(\tau\); inverted and summed they leave \(\lvert 1 - e^{-j2\pi f\tau}\rvert\). The floor rises 6 dB per doubling of velocity or frequency and falls 6 dB per doubling of PRF. Cardiac and abdominal motion of several mm/s limits CTR more than the sequence algebra does.
After the flash, bubbles enter the destroyed slice of elevational thickness \(E\) from its sides at the mean velocity, so the refill time is of order \(E/v\); the exponential is the fit used by Wei et al. In m28_tic_perfusion.m, \(A\) = 50, \(\beta = 0.8\ \mathrm{s^{-1}}\) and 50 speckle cells give \(A\) to 3.6 % and \(\beta\) to 12 % (one standard deviation): \(\beta\) is always the noisier parameter.
Perfusion models describe concentration, which is proportional to linear intensity. The log map compresses the peak and lifts the tail, so the same model fits a much wider curve. TTP survives (a monotone map keeps the peak in place) but MTT, AUC and wash-in rate do not. In the figure the normalized log-compressed curve (dashed) is still at 43 % of its range when the linear curve has fallen to 10 % of its peak.
\(\frac{d}{dn}\left(n\,e^{-nA}\right) = 0\) gives \(n^{*} = 1/A\) and a maximum of \(1/(eA)\) isolated bubbles per unit area (28.12): 1.65 per mm² at 5 MHz (\(F_x = \lambda\), \(F_z = 0.75\lambda\), \(k\) = 1). Above \(n^{*}\) adding bubbles reduces the useful count. Because \(A \propto \lambda^{2}\), the optimum density grows as \(f^{2}\): 9 times higher at 15 MHz than at 5 MHz.
The flux through a vessel is concentration times volumetric flow. An 8 µm capillary at 1 mm/s carries 0.05 nL/s, so at 140 bubbles per µL it sees 0.007 bubbles per second. The frame rate does not change this; only the concentration and the vasculature do, and the concentration is capped by \(n^{*}\) (rule above). Hence ULM acquisitions of minutes for capillary-level maps.
The singular components that capture tissue also capture anything with the same slow-time spectrum. Bubbles slower than a few mm/s fall into the tissue subspace and are cut with it. Nonlinear sequences reject tissue by its linearity instead, so they are the natural complement to SVD for slow flow, at the cost of extra transmits per frame and motion residue.
The image of a sub-wavelength scatterer is the PSF; only its complex amplitude depends on the bubble. That amplitude follows \(\sigma_s(R, f)\), which in the shell model of 28.2 spans 35 dB across radii of 1–4 µm at 3 MHz and 18 dB at 5 MHz, so localization precision is a distribution. Reporting a single "\(\lambda/10\)" hides that spread.
28.11Cheat sheet
| Concept | Formula | Ultrasound use |
|---|---|---|
| Minnaert resonance | \(f_0R_0 = \sqrt{3\kappa P_0/\rho}\,/\,2\pi \approx\) 3.29 MHz·µm (air, 1 atm) | First estimate of resonant size |
| With surface tension | \(\rho\omega_0^{2}R_0^{2} = 3\kappa P_0 + 2(3\kappa - 1)\sigma/R_0\) | Free micron bubbles |
| Marmottant shell (elastic) | \(\rho\omega_0^{2}R_0^{2} = 3\kappa P_0 + 2(3\kappa + 1)\sigma_0/R_0 + 4\chi/R_0\) | Lipid agents, spectroscopy fits |
| Damping | \(\delta = \omega R_0/c + 4\mu/(\rho\omega R_0^{2}) + 4\kappa_s/(\rho\omega R_0^{3})\) | Q, bandwidth, ring-down |
| Scattering cross-section | \(\sigma_s = 4\pi R_0^{2}/[((f_0/f)^{2} - 1)^{2} + \delta^{2}]\) | Backscatter, attenuation of agents |
| Rayleigh small sphere | \((4\pi/9)k^{4}a^{6}[(\Delta\kappa/\kappa)^{2} + (1/3)(3\Delta\rho/(2\rho_s + \rho))^{2}]\) | Bubble versus tissue versus rigid |
| Rayleigh–Plesset (Marmottant) | \(\rho(R\ddot R + 1.5\dot R^{2}) = p_g(1 - 3\kappa\dot R/c) - P_0 - 2\sigma(R)/R - 4\mu\dot R/R - 4\kappa_s\dot R/R^{2} - P_{\mathrm{ac}}\) | Nonlinear bubble response |
| Emitted pressure | \(p_s = \rho(R^{2}\ddot R + 2R\dot R^{2})/r\) | Echo of one bubble |
| Laplace pressure | \(\Delta p = 2\sigma/R\) (144 kPa at 1 µm, water) | Dissolution, why shells and heavy gases |
| Pulse inversion | \(e(1) + e(-1) = 2a_2p^{2} + \dots\) | Harmonic contrast, full bandwidth |
| Amplitude modulation | \(e(1) - 2e(\tfrac12) = \tfrac12 a_2p^{2} + \tfrac34 a_3p^{3}\) | Nonlinear fundamental, penetration |
| Motion residue | \(2\lvert\sin(\pi f\tau)\rvert\), \(\tau = 2v/(c\cdot\mathrm{PRF})\) | CTR floor of multi-pulse modes |
| Lognormal TIC | \(\mathrm{MTT} = e^{\mu + s^{2}/2}\), \(\mathrm{TTP} = t_0 + e^{\mu - s^{2}}\) | Bolus perfusion |
| Flash-replenishment | \(I = A(1 - e^{-\beta t})\), flow \(\propto A\beta\) | Infusion perfusion |
| ULM isolation | \(P_{\mathrm{iso}} = e^{-nA}\), \(A = \pi k^{2}F_xF_z\), \(n^{*} = 1/A\) | Choosing concentration |
| Capillary flux | bubbles/s \(= C\cdot\pi d^{2}v/4\) | ULM acquisition time |
| Doppler of a bubble | \(f_D = 2vf/c\) | SVD cut-off versus slow flow |
28.12Interview and qualifying-exam questions
I1Derive the resonance frequency of a gas bubble in water. What is it for \(R_0\) = 2 µm? · screen
Radial flow \(u = R^{2}\dot R/r^{2}\) has kinetic energy \(\tfrac12(4\pi R^{3}\rho)\dot R^{2}\): effective mass \(4\pi R_0^{3}\rho\). Polytropic gas: stiffness \(12\pi\kappa P_0R_0\). \(\omega_0^{2} = 3\kappa P_0/(\rho R_0^{2})\), \(f_0R_0 \approx\) 3.3 MHz·µm for air, so 1.64 MHz at 2 µm. Then I would add surface tension (2.04 MHz for a free 2 µm air bubble) and say which \(\kappa\) I assumed. Follow-up: "And for a lipid-coated agent?" Add \(4\chi/R_0\): with \(\chi\) = 0.5 N/m and a perfluorocarbon-like \(\kappa\) about 3 MHz.
I2Why do microbubbles scatter so much more than tissue or red cells of similar size? · screen
Rayleigh scattering from a small sphere has a compressibility (monopole) term \(((\kappa_s - \kappa)/\kappa)^{2}\) and a density (dipole) term. Gas is about 15 000 times more compressible than water, so the monopole term is about \(2.4\times 10^{8}\), against order \(10^{-3}\) for a tissue inhomogeneity of a few percent: some \(10^{11}\) per particle of equal size. Near resonance there is another factor of order \(1/\delta^{2}\) (about 20–200 for real bubbles). Follow-up: "Then why not use bigger bubbles?" They do not pass the lung capillaries, and \(\sigma_s\) above resonance saturates at \(4\pi R_0^{2}\).
I3Write the Rayleigh–Plesset equation with a Marmottant shell and explain each term. How would you integrate it? · onsite
\(\rho(R\ddot R + 1.5\dot R^{2}) = p_g(R)(1 - 3\kappa\dot R/c) - P_0\) \(- 2\sigma(R)/R - 4\mu\dot R/R - 4\kappa_s\dot R/R^{2} - P_{\mathrm{ac}}(t)\) (28.6): liquid inertia; polytropic gas with first-order radiation damping; ambient pressure; radius-dependent tension (0 buckled, \(\chi(R^{2}/R_b^{2} - 1)\) elastic, \(\sigma_w\) ruptured); liquid viscosity; shell viscosity; drive. Integrate \((R, \dot R)\) with RK4 or an adaptive RK45, 100–200 steps per period with refinement during collapse, and check convergence by halving the step (\(R_{\mathrm{min}}\) should not move). Follow-up: "Why not Euler?" Collapse is stiff; explicit Euler needs orders of magnitude more steps and loses energy accuracy.
I4Compare pulse inversion and amplitude modulation. What does each cancel and keep, and where does the signal live? · onsite
With \(e(a) = \sum a_n a^{n}p^{n}\): PI sums \(e(1) + e(-1) = 2a_2p^{2}\), the even orders, at \(2f_0\). AM forms \(e(1) - 2e(\tfrac12) = \tfrac12 a_2p^{2} + \tfrac34 a_3p^{3}\), which includes an odd-order nonlinear signal in the fundamental band: better penetration. AM needs an exact half amplitude (0.1 dB error leaves −38.7 dB of tissue); PI needs exact inversion. Both are hurt by motion: \(2\lvert\sin(\pi f\tau)\rvert\) with \(\tau = 2v/(c\cdot\mathrm{PRF})\). Follow-up: "Which removes tissue harmonics better?" AM, because the half pulse carries a quarter of the tissue harmonic; PI keeps it entirely.
I5In a contrast mode you see bright tissue around the focus before the agent arrives. Explain and fix. · onsite
Nonlinear propagation: near the focus the tissue-generated second harmonic is largest (\(B_2/B_1 \approx \beta\omega P x/2\rho c^{3}\), −28 dB at 100 kPa, 2.5 MHz, 40 mm), and it is even in amplitude, so PI keeps it. Fixes: lower MI (the residue falls 6 dB per halving of \(P\), the bubble signal much less), move the focus deeper or use unfocused transmits, switch to AM or AM + PI, and check for motion residue and transmit asymmetry. Follow-up: "How would you measure it?" Image a tissue-mimicking phantom without agent in contrast mode and plot the residue against MI.
I6What is compression-only behaviour, and why does it matter for low-MI imaging? · research
A phospholipid bubble near its buckling radius has no surface tension in compression and a rising one in expansion (Marmottant). It compresses several times further than it expands (de Jong et al. 2007 observed it optically; the model gives a ratio of 3 at 10 kPa). That asymmetry is a strong even-order response at very low pressure, which is why PI and AM work at MI below 0.1, and it varies from bubble to bubble with the initial surface tension. Follow-up: "Consequence for quantification?" The echo is not proportional to concentration times incident pressure, so calibrations at one MI do not transfer to another.
I7How would you measure perfusion with flash-replenishment, and what limits its accuracy? · onsite
Constant infusion, low-MI imaging, a few high-MI frames to clear the slice, then fit \(I(t) = A(1 - e^{-\beta t})\) on linearized intensity: \(A \propto\) blood volume, \(\beta \propto\) velocity/slice thickness, \(A\beta \propto\) flow. Limits: incomplete destruction, bubble destruction by the imaging frames (biases \(\beta\) upward), attenuation and shadowing, receiver saturation of the plateau, beam elevation profile, speckle noise (\(\beta\) is the noisiest parameter, 12 % in the MATLAB example). Follow-up: "Why not fit the log-compressed curve?" The model is linear in concentration; the log map changes the shape (MTT doubles in the bolus example).
I8Why is the subharmonic attractive for contrast imaging, and where is its threshold lowest? · research
Tissue propagation does not generate \(f/2\), so a subharmonic image has a very high CTR. Period doubling has a threshold that is lowest for a drive near \(2f_0\), where \(f/2\) hits the resonance; for a free 2 µm bubble the model gives \(f/2\) at −30 dB at 100 kPa and −4 dB at 200 kPa. Costs: long bursts (poor axial resolution), the threshold behaviour, and the need for a transducer with bandwidth at \(f\) and \(f/2\). Follow-up: "What else is it used for?" Ambient-pressure estimation, because the subharmonic amplitude of some agents depends on hydrostatic pressure.
I9How do you choose the microbubble concentration for ULM? · research
Model bubbles in the image as a Poisson process with areal density \(n\); a bubble is usable if no other lies within \(k\) FWHM: \(P_{\mathrm{iso}} = e^{-nA}\), \(A = \pi k^{2}F_xF_z\) (28.12). The isolated count \(n\,e^{-nA}\) peaks at \(n^{*} = 1/A\) (\(P_{\mathrm{iso}} = 1/e\)). At 5 MHz with \(F_x = \lambda\), \(F_z = 0.75\lambda\): 4.5 per mm². Convert \(n = C\) × slice × vascular fraction, then remember that large vessels are locally much denser. Follow-up: "And the price of a lower concentration?" Capillaries are sampled more slowly in proportion (next question).
I10Why does ULM take minutes to map capillaries, and how would you shorten it? · research
A capillary of diameter \(d\) and speed \(v\) receives \(C\cdot\pi d^{2}v/4\) bubbles per second: about one per 140 s for 8 µm, 1 mm/s and 140 per µL. Frame rate does not help; the concentration is capped by overlap at \(n^{*} = 1/A\). Ways out: localize at higher densities (overlap-tolerant or sparsity-based estimators, learned localizers), shrink \(A\) (higher frequency, 3D removes elevation overlap), accept a coarser map for capillaries (density rather than tracks), or separate subpopulations (for example by acoustic signature). Follow-up: "What did Hingot et al. show?" That microvascular flow sets the compromise between resolution and acquisition time.
I11Your SVD-filtered ULM data miss the slow capillary flow seen in histology. Why, and what would you change? · research
The SVD cut separates by slow-time spectrum and spatial coherence. A bubble at 1 mm/s has a Doppler of 6.5 Hz at 5 MHz, inside the tissue-motion band, so it is projected out with tissue. Options: an adaptive cut chosen from the spatial-similarity of singular vectors, lower cut with better motion correction, longer ensembles for finer frequency resolution, or a nonlinear sequence (AM) that removes tissue regardless of velocity, accepting the motion residue and the extra transmits. Follow-up: "How would you validate?" A flow phantom with known slow flow, or a microchannel, reporting detection rate versus velocity.
28.13MATLAB lab
Base MATLAB only (no toolboxes); each script runs as-is in MATLAB R2019b+ and GNU Octave 8. The RK4 integrator is written out in each script (no ode45). Paste into the editor and press Run.
m28_rayleigh_plesset.mRK4 Rayleigh–Plesset with a Marmottant shell: radius, emitted pressure and spectrum for six cases · change R0, f, P, shell state · runtime ≈ 15 s (Octave)
%% m28_rayleigh_plesset.m -- Rayleigh-Plesset (Marmottant form) with a hand-written RK4
% Purpose : integrate the radial dynamics of a free gas bubble and of a lipid-coated
% bubble (Marmottant buckling/rupture shell), compute the emitted pressure
% p_s = rho (R^2 R'' + 2 R R'^2) / r and its spectrum, and read off the
% harmonic, subharmonic and compression-only signatures.
% Model : rho (R R'' + 3/2 R'^2) = Pg0 (R0/R)^(3g) (1 - 3g R'/c) - P0 - 2 sig(R)/R
% - 4 mu R'/R - 4 ks R'/R^2 - Pac(t)
% sig(R) = 0 (R < Rb), chi (R^2/Rb^2 - 1) (Rb <= R <= Rr), sigma_w (R > Rr).
% Change : R0, f, P (cases table below), chi, ks, s0 (surface tension at rest).
% Expect : free 2 um bubble (f0 = 2.04 MHz) at 2 MHz, 200 kPa: Rmax/R0 = 2.30,
% 2nd harmonic -10.6 dB; driven at 2 f0 (4.08 MHz) the subharmonic is
% -30 dB at 100 kPa and -4 dB at 200 kPa (above threshold); a buckled
% shell (s0 = 0) compresses 2.5 times further than it expands.
% Provenance: simulation (this script). Shell values chi = 0.5 N/m, s0 = 0.02 N/m, ks = 5e-9 kg/s,
% g = 1.07 are illustrative (inside the range fitted to lipid agents, not one agent).
% Runtime : about 15 s in Octave (RK4 in an interpreted loop), a few s in MATLAB.
clear; close all; rng(1);
%% Liquid and gas constants (water at 20 C)
rho = 1000; % kg/m^3
c = 1480; % m/s
mu = 1e-3; % Pa s
P0 = 101325; % Pa
sw = 0.072; % N/m, water surface tension
rObs = 0.01; % m, observation distance for the emitted pressure
%% Cases: {label, R0 (m), shell (0/1), gamma, s0 (N/m), f (Hz or 0 = 2 f0), P (Pa)}
chi = 0.5; ks = 5e-9; % shell elasticity (N/m) and dilatational viscosity (kg/s)
cases = { 'free, f = 2 MHz, 50 kPa', 2e-6, 0, 1.40, 0, 2e6, 50e3;
'free, f = 2 MHz, 200 kPa', 2e-6, 0, 1.40, 0, 2e6, 200e3;
'free, f = 2 f0, 100 kPa', 2e-6, 0, 1.40, 0, 0, 100e3;
'free, f = 2 f0, 200 kPa', 2e-6, 0, 1.40, 0, 0, 200e3;
'shell elastic s0=0.02, 3 MHz, 100 kPa', 2e-6, 1, 1.07, 0.02, 3e6, 100e3;
'shell buckled s0=0, 3 MHz, 50 kPa', 2e-6, 1, 1.07, 0, 3e6, 50e3 };
ncyc = 30; % cycles in the tone burst (sin^2 ramps of one cycle)
nper = 200; % RK4 steps per acoustic period (split further when R moves fast)
nc = size(cases, 1);
res = zeros(nc, 7);
keep = cell(nc, 1);
fprintf('%-40s %7s %6s %6s %7s %7s %7s %8s\n', 'case', 'f0 MHz', 'Rmax', 'Rmin', 'H2 dB', 'H3 dB', 'f/2 dB', 'ps,max Pa');
for ic = 1:nc
R0 = cases{ic, 2}; shell = cases{ic, 3}; g = cases{ic, 4}; s0 = cases{ic, 5};
%% Shell law and equilibrium gas pressure
if shell
Rb = R0 / sqrt(1 + s0 / chi); Rr = Rb * sqrt(1 + sw / chi);
sig = @(R) (R >= Rb & R <= Rr) .* chi .* (R.^2 / Rb^2 - 1) + (R > Rr) .* sw;
sig0 = s0; kss = ks;
K = 3*g*P0 + 2*(3*g + 1)*sig0/R0 + 4*chi/R0; % linear stiffness (elastic state)
else
sig = @(R) sw + 0*R; sig0 = sw; kss = 0;
K = 3*g*P0 + 2*(3*g - 1)*sw/R0;
end
f0 = sqrt(K / rho) / (2*pi*R0);
f = cases{ic, 6}; if f == 0, f = 2 * f0; end
P = cases{ic, 7};
Pg0 = P0 + 2 * sig0 / R0;
T = ncyc / f; tr = 1 / f;
env = @(t) (t > 0 & t < T) .* sin(pi/2 * min(1, min(t / tr, (T - t) / tr))).^2;
pac = @(t) P * env(t) .* sin(2*pi*f*t);
acc = @(t, R, U) ((Pg0 * (R0 ./ R).^(3*g) .* (1 - 3*g*U/c) - P0 - 2*sig(R)./R ...
- 4*mu*U./R - 4*kss*U./R.^2 - pac(t)) / rho - 1.5*U.^2) ./ R;
%% RK4 integration
dt = 1 / (f * nper); n = round((T + 4/f) / dt);
t = (0:n) * dt; R = zeros(1, n + 1); U = zeros(1, n + 1); A = zeros(1, n + 1);
R(1) = R0; A(1) = acc(0, R0, 0); Rk = R0; Uk = 0; tk = 0;
for i = 1:n
m = max(1, min(64, ceil(abs(Uk) * dt / Rk / 0.005))); % sub-steps when R moves > 0.5 % per step
h = dt / m;
for s = 1:m
k1R = Uk; k1U = acc(tk, Rk, Uk);
k2R = Uk + h/2*k1U; k2U = acc(tk + h/2, Rk + h/2*k1R, Uk + h/2*k1U);
k3R = Uk + h/2*k2U; k3U = acc(tk + h/2, Rk + h/2*k2R, Uk + h/2*k2U);
k4R = Uk + h*k3U; k4U = acc(tk + h, Rk + h*k3R, Uk + h*k3U);
Rk = Rk + h/6*(k1R + 2*k2R + 2*k3R + k4R);
Uk = Uk + h/6*(k1U + 2*k2U + 2*k3U + k4U);
tk = tk + h;
end
R(i+1) = Rk; U(i+1) = Uk; A(i+1) = acc(tk, Rk, Uk);
end
ps = rho * (R.^2 .* A + 2 * R .* U.^2) / rObs; % emitted pressure at rObs (Pa)
%% Spectrum of the steady part (30 % to 95 % of the burst), Hann window, 8x zero padding
idx = find(t > 0.3*T & t < 0.95*T); x = ps(idx) - mean(ps(idx)); L = numel(x);
w = 0.5 - 0.5*cos(2*pi*(0:L-1)/(L-1));
nf = 2^nextpow2(8*L); X = abs(fft(x .* w, nf)); fa = (0:nf-1) / (nf*dt);
lev = @(ff) max(X(max(1, round(ff/fa(2)) - 2):round(ff/fa(2)) + 4));
L1 = lev(f);
res(ic, :) = [f0/1e6, max(R)/R0, min(R)/R0, 20*log10(lev(2*f)/L1), 20*log10(lev(3*f)/L1), 20*log10(lev(f/2)/L1), max(ps)];
fprintf('%-40s %7.2f %6.2f %6.2f %7.1f %7.1f %7.1f %8.1f\n', cases{ic, 1}, res(ic, :));
keep{ic} = struct('t', t, 'R', R / R0, 'ps', ps, 'f', f, 'fa', fa(1:nf/2), 'X', 20*log10(X(1:nf/2)/L1 + 1e-12));
end
%% Compression-only index for the buckled shell: (R0 - Rmin)/(Rmax - R0)
fprintf('\nCompression/expansion ratio (R0 - Rmin)/(Rmax - R0): elastic %.2f, buckled %.2f\n', ...
(1 - res(5, 3)) / (res(5, 2) - 1), (1 - res(6, 3)) / (res(6, 2) - 1));
%% Step-size check on the most violent case (free, 2 MHz, 200 kPa): already converged at 200/period
fprintf('Step check: Rmin/R0 = %.4f at %d steps/period (compare 0.415 at 400 and 800 steps/period)\n', res(2, 3), nper);
%% Figures
figure('Position', [50 50 900 640]);
sel = [1 2 4 6];
for k = 1:4
s = keep{sel(k)};
subplot(4, 2, 2*k - 1); plot(s.t * 1e6, s.R, 'b'); grid on;
ylabel('R / R_0'); title(cases{sel(k), 1});
if k == 4, xlabel('time (\mus)'); end
subplot(4, 2, 2*k); plot(s.fa / s.f, s.X, 'r'); grid on; xlim([0 4]); ylim([-80 5]);
ylabel('dB re f'); title('emitted spectrum');
if k == 4, xlabel('frequency / drive frequency'); end
end
What to look for: the free 2 µm bubble (\(f_0\) = 2.04 MHz) at 2 MHz grows to \(R_{\mathrm{max}}/R_0\) = 1.59 at 50 kPa and 2.30 at 200 kPa, with the 2nd harmonic at −7.9 and −10.6 dB; driven at \(2f_0\) the subharmonic is −30 dB at 100 kPa and −4.3 dB at 200 kPa; the elastic shell (\(f_0\) = 2.99 MHz) at 100 kPa compresses 1.3 times further than it expands, the buckled one 2.5 times at 50 kPa; the step check prints \(R_{\mathrm{min}}/R_0\) = 0.4163, converged.
m28_resonance_shell.mMinnaert constants, tension and shell corrections, RK4 ring-down check of the linear formula, damping and cross-sections · change \(\chi\), \(\kappa_s\), \(\sigma_0\), \(\kappa\) · runtime ≈ 2 s
%% m28_resonance_shell.m -- Minnaert, surface tension, shell stiffness, damping, cross-section
% Purpose : (1) evaluate the Minnaert constant f0*R0 for the usual assumptions;
% (2) add surface tension and a Marmottant shell to the linear resonance;
% (3) verify the linear formula by a small-amplitude RK4 ring-down;
% (4) compute damping (radiation, liquid viscosity, shell viscosity) and the
% scattering cross-section, and compare with a rigid sphere.
% Linear : rho w0^2 R0^2 = 3 g P0 + 2(3g - 1) sigma/R0 (free bubble)
% rho w0^2 R0^2 = 3 g P0 + 2(3g + 1) s0/R0 + 4 chi/R0 (shell, elastic state)
% sigma_s = 4 pi R0^2 / [((f0/f)^2 - 1)^2 + delta^2],
% delta = w R0/c + 4 mu/(rho w R0^2) + 4 ks/(rho w R0^3)
% Change : chi, ks, s0, gamma, the radius list.
% Expect : f0*R0 = 3.29 MHz um (air, 1 atm, adiabatic), 3.26 at 1 bar, 2.78 isothermal;
% a free 1 um air bubble resonates at 4.74 MHz (not 3.29) because of surface
% tension; a 2 um shelled bubble (g = 1.07) at 2.99 MHz versus 1.44 MHz
% (Minnaert, same gas) and 1.75 MHz (free, with tension); ring-down
% frequencies agree with sqrt(w0^2 - (b/2)^2) to better than 0.1 %.
% Provenance: closed-form linear formulas checked by an RK4 ring-down. Shell values chi = 0.5 N/m, s0 = 0.02 N/m, ks = 5e-9 kg/s,
% g = 1.07 are illustrative (inside the range fitted to lipid agents, not one agent).
% Runtime : about 2 s.
clear; close all; rng(1);
rho = 1000; c = 1480; mu = 1e-3; P0 = 101325; sw = 0.072;
chi = 0.5; ks = 5e-9; s0 = 0.02; % shell: N/m, kg/s, N/m
%% 1. Minnaert constant (f0 * R0 in m/s = MHz * um)
mc = @(g, p, rr) sqrt(3 * g * p / rr) / (2*pi);
fprintf('Minnaert f0*R0 (MHz um): air 1 atm adiabatic %.3f | 1 bar, rho 1000: %.3f | isothermal %.3f | g = 1.07 (C3F8/C4F10-like) %.3f | SF6 g = 1.095: %.3f\n', ...
mc(1.4, P0, 998), mc(1.4, 1e5, 1000), mc(1.0, P0, rho), mc(1.07, P0, rho), mc(1.095, P0, rho));
%% 2. Resonance versus radius for four models
R0s = [0.5 1 1.5 2 3 4 5] * 1e-6;
fM = @(R, g) sqrt(3*g*P0/rho) ./ (2*pi*R);
fT = @(R, g) sqrt((3*g*P0 + 2*(3*g - 1)*sw./R) / rho) ./ (2*pi*R);
fS = @(R, g) sqrt((3*g*P0 + 2*(3*g + 1)*s0./R + 4*chi./R) / rho) ./ (2*pi*R);
fprintf('\n R0 (um) | Minnaert air | + tension air | free g=1.07 + tension | shell g=1.07 (MHz)\n');
for R = R0s
fprintf(' %4.1f | %6.2f | %6.2f | %6.2f | %6.2f\n', R*1e6, fM(R, 1.4)/1e6, fT(R, 1.4)/1e6, fT(R, 1.07)/1e6, fS(R, 1.07)/1e6);
end
%% 3. Ring-down check with RK4 (release from R = 1.002 R0, no drive)
fprintf('\nRing-down check (damped frequency, MHz)\n');
cfg = {'free air', 1.4, 0; 'shell', 1.07, 1};
for ic = 1:2
g = cfg{ic, 2}; shell = cfg{ic, 3};
for R0 = [1 2 3] * 1e-6
if shell
Rb = R0 / sqrt(1 + s0/chi); Rr = Rb * sqrt(1 + sw/chi);
sig = @(R) (R >= Rb & R <= Rr) .* chi .* (R.^2/Rb^2 - 1) + (R > Rr) .* sw; sig0 = s0; kk = ks;
w0 = 2*pi*fS(R0, g);
else
sig = @(R) sw + 0*R; sig0 = sw; kk = 0; w0 = 2*pi*fT(R0, g);
end
Pg0 = P0 + 2*sig0/R0;
acc = @(R, U) ((Pg0*(R0./R).^(3*g).*(1 - 3*g*U/c) - P0 - 2*sig(R)./R - 4*mu*U./R - 4*kk*U./R.^2)/rho - 1.5*U.^2)./R;
b = (4*mu/R0 + 4*kk/R0^2 + 3*g*Pg0/c) / (rho*R0); % linear damping coefficient (1/s)
wd = sqrt(w0^2 - (b/2)^2);
dt = 2*pi/w0/400; n = 2400; x = zeros(1, n); Rk = 1.002*R0; Uk = 0;
for i = 1:n
k1R = Uk; k1U = acc(Rk, Uk);
k2R = Uk + dt/2*k1U; k2U = acc(Rk + dt/2*k1R, Uk + dt/2*k1U);
k3R = Uk + dt/2*k2U; k3U = acc(Rk + dt/2*k2R, Uk + dt/2*k2U);
k4R = Uk + dt*k3U; k4U = acc(Rk + dt*k3R, Uk + dt*k3U);
Rk = Rk + dt/6*(k1R + 2*k2R + 2*k3R + k4R); Uk = Uk + dt/6*(k1U + 2*k2U + 2*k3U + k4U);
x(i) = Rk - R0;
end
zc = find(x(1:end-1) < 0 & x(2:end) >= 0); % upward zero crossings
zf = zc - x(zc) ./ (x(zc+1) - x(zc)); % sub-sample interpolation
fmeas = 1 / (mean(diff(zf(1:min(4, end)))) * dt);
fprintf(' %-8s R0 = %.0f um: formula f0 = %.3f, damped sqrt(w0^2-(b/2)^2) = %.3f, RK4 = %.3f, Q = %.1f\n', ...
cfg{ic, 1}, R0*1e6, w0/2/pi/1e6, wd/2/pi/1e6, fmeas/1e6, w0/b);
end
end
%% 4. Damping and cross-section at resonance (physical radiation damping w R0 / c)
fprintf('\nDamping at resonance and cross-section (free bubble g = 1.4 with tension; shell g = 1.07)\n');
fprintf(' R0 (um) | model | f0 MHz | d_rad d_vis d_shell | Q | sigma_s/(pi R0^2) at f0\n');
for R0 = [1 2 3] * 1e-6
for shell = 0:1
if shell, f0 = fS(R0, 1.07); kk = ks; lab = 'shell'; else, f0 = fT(R0, 1.4); kk = 0; lab = 'free '; end
w = 2*pi*f0; dr = w*R0/c; dv = 4*mu/(rho*w*R0^2); dsh = 4*kk/(rho*w*R0^3); d = dr + dv + dsh;
fprintf(' %3.0f | %s | %5.2f | %.3f %.3f %.3f | %4.1f | %6.0f\n', R0*1e6, lab, f0/1e6, dr, dv, dsh, 1/d, 4/d^2);
end
end
kappaRatio = rho*c^2/(1.4*P0);
fprintf('\nGas/water compressibility ratio (air, adiabatic): %.0f; bubble / rigid sphere cross-section far below resonance: (4/7) ratio^2 = %.2e\n', ...
kappaRatio, 4/7*kappaRatio^2);
%% 5. Figures: resonance versus radius, and cross-section versus frequency for R0 = 2 um
Rg = linspace(0.5, 5, 200) * 1e-6;
figure('Position', [50 50 900 380]);
subplot(1, 2, 1);
loglog(Rg*1e6, fM(Rg, 1.4)/1e6, 'k--', Rg*1e6, fT(Rg, 1.4)/1e6, 'b', Rg*1e6, fT(Rg, 1.07)/1e6, 'c', Rg*1e6, fS(Rg, 1.07)/1e6, 'r'); grid on;
xlabel('R_0 (\mum)'); ylabel('f_0 (MHz)'); title('Linear resonance'); xlim([0.5 5]);
legend('Minnaert, air', 'air + tension', '\gamma = 1.07 + tension', 'lipid shell', 'Location', 'southwest');
subplot(1, 2, 2);
fg = linspace(0.2, 10, 600) * 1e6; R0 = 2e-6;
for shell = 0:1
if shell, f0 = fS(R0, 1.07); kk = ks; st = 'r'; else, f0 = fT(R0, 1.4); kk = 0; st = 'b'; end
w = 2*pi*fg; d = w*R0/c + 4*mu./(rho*w*R0^2) + 4*kk./(rho*w*R0^3);
ss = 4*pi*R0^2 ./ (((f0./fg).^2 - 1).^2 + d.^2);
semilogy(fg/1e6, ss/(pi*R0^2), st); hold on;
end
k = 2*pi*fg/c; semilogy(fg/1e6, 7/9*(k*R0).^4, 'k--');
grid on; xlabel('f (MHz)'); ylabel('\sigma_s / \pi R_0^2'); title('Scattering cross-section, R_0 = 2 \mum');
legend('free (air)', 'lipid shell', 'rigid sphere', 'Location', 'southeast'); ylim([1e-10 1e4]);
What to look for: \(f_0R_0\) = 3.287 (air, 1 atm), 3.262 (1 bar), 2.775 (isothermal); the resonance table matches 28.2; the ring-down frequencies match \(\sqrt{\omega_0^{2} - (b/2)^{2}}\) to the third decimal (4.724 vs 4.723 MHz at 1 µm free); the damping table reproduces \(Q\) = 10.5 and 4.7 at 2 µm; the bubble/rigid-sphere ratio prints \(1.36\times 10^{8}\).
m28_contrast_pulse_sequences.mPI, AM and AM + PI on an RK4 bubble and a nonlinearly insonified tissue scatterer, with motion and gain error · change P, f, R0, x, v, PRF, error · runtime ≈ 6 s
%% m28_contrast_pulse_sequences.m -- PI, AM and AM+PI residuals for a bubble and for tissue
% Purpose : transmit the pulse families used in contrast imaging, compute the echo of a
% lipid-coated bubble (Rayleigh-Plesset/Marmottant, RK4) and of a linear tissue
% scatterer that sees the same nonlinearly propagated incident pulse, combine the
% echoes, and measure what survives.
% PI : e(+p) + e(-p) keeps even orders (2 a2 p^2)
% AM : e(p) - 2 e(p/2) keeps a2 p^2/2 + 3 a3 p^3/4 (odd orders land in the fundamental band)
% AMPI : e(p) + 2 e(-p/2) keeps 3 a2 p^2/2 + 3 a3 p^3/4
% Tissue : incident pulse p_inc = a p1 + a^2 (beta x / (2 rho c^3)) d(p1^2)/dt
% (second-order plane-wave solution), echo = p_inc delayed by motion 2 v k / (c PRF).
% Change : P (kPa at the bubble), f, R0, depth x, motion v, half-amplitude error errDb.
% Expect : at 100 kPa, 2.5 MHz, R0 = 2 um, 40 mm: tissue suppressed by 23.4 dB (PI, limited
% by the tissue harmonic), 35.4 dB (AM), 25.9 dB (AMPI); bubble retained within 1.5 dB
% (AM) but 7.5 dB down (PI); 10 mm/s motion at PRF 5 kHz caps AM at about 27 dB;
% without tissue nonlinearity and motion every sequence cancels tissue exactly.
% Provenance: simulation (this script). Shell values chi = 0.5 N/m, s0 = 0.02 N/m, ks = 5e-9 kg/s,
% g = 1.07 are illustrative (inside the range fitted to lipid agents, not one agent).
% Runtime : about 6 s.
clear; close all; rng(1);
rho = 1000; cL = 1480; mu = 1e-3; P0 = 101325; sw = 0.072; % liquid around the bubble
cT = 1540; rhoT = 1000; beta = 4.5; % tissue path: c, rho, 1 + B/2A
R0 = 2e-6; g = 1.07; chi = 0.5; ks = 5e-9; s0 = 0.02; % lipid-coated bubble
f = 2.5e6; N = 3; P = 100e3; x = 0.04; prf = 5000; % 3-cycle Hann pulse, 40 mm path
Rb = R0 / sqrt(1 + s0/chi); Rr = Rb * sqrt(1 + sw/chi);
sig = @(R) (R >= Rb & R <= Rr) .* chi .* (R.^2/Rb^2 - 1) + (R > Rr) .* sw;
Pg0 = P0 + 2*s0/R0;
Tp = N / f; fs = 100 * f; dt = 1/fs; n = round((Tp + 6/f) / dt) + 1; t = (0:n-1) * dt;
p1 = @(t) (t > 0 & t < Tp) .* 0.5 .* (1 - cos(2*pi*t/Tp)) .* sin(2*pi*f*t);
dp1 = @(t) (p1(t + dt/4).^2 - p1(t - dt/4).^2) / (dt/2); % d(p1^2)/dt
knl = beta * x / (2 * rhoT * cT^3); % s/Pa
fprintf('Tissue 2nd harmonic for a CW wave at this P, f, x: B2/B1 = sigma/2 = %.3f (%.1f dB)\n', ...
knl * 2*pi*f * P, 20*log10(knl * 2*pi*f * P));
scen = {'default', 0, 0, beta; 'motion 10 mm/s', 0.01, 0, beta; 'half-pulse error +0.1 dB', 0, 0.1, beta; 'linear tissue, no motion', 0, 0, 0};
nfft = 2^nextpow2(4*n); df = fs/nfft;
bandE = @(y, lo, hi) sum(abs(subsref(fft(y, nfft), struct('type', '()', 'subs', {{(ceil(lo/df):floor(hi/df)) + 1}}))).^2);
fprintf('\n%-26s %-5s | tissue: all f-band 2f-band | bubble: all f-band 2f-band | CTR gain (dB)\n', 'scenario', 'seq');
for is = 1:size(scen, 1)
v = scen{is, 2}; errDb = scen{is, 3}; knlS = knl * (scen{is, 4} > 0);
half = 0.5 * 10^(errDb/20);
amps = [1, -1, half, -half]; idxp = [0, 1, 1, 1]; % pulse amplitudes and pulse index (motion)
E_t = zeros(4, n); E_b = zeros(4, n);
for k = 1:4
a = amps(k); tau = idxp(k) * 2 * v / (cT * prf) + 1e-7;
pinc = @(tt) a * P * p1(tt - tau) + a^2 * P^2 * knlS * dp1(tt - tau);
E_t(k, :) = pinc(t);
acc = @(tt, R, U) ((Pg0*(R0./R).^(3*g).*(1 - 3*g*U/cL) - P0 - 2*sig(R)./R - 4*mu*U./R - 4*ks*U./R.^2 - pinc(tt))/rho - 1.5*U.^2)./R;
Rk = R0; Uk = 0; ps = zeros(1, n); ps(1) = rho*(R0^2*acc(0, R0, 0))/0.01;
for i = 2:n
tk = t(i-1); h = dt;
k1R = Uk; k1U = acc(tk, Rk, Uk);
k2R = Uk + h/2*k1U; k2U = acc(tk + h/2, Rk + h/2*k1R, Uk + h/2*k1U);
k3R = Uk + h/2*k2U; k3U = acc(tk + h/2, Rk + h/2*k2R, Uk + h/2*k2U);
k4R = Uk + h*k3U; k4U = acc(tk + h, Rk + h*k3R, Uk + h*k3U);
Rk = Rk + h/6*(k1R + 2*k2R + 2*k3R + k4R); Uk = Uk + h/6*(k1U + 2*k2U + 2*k3U + k4U);
ps(i) = rho*(Rk^2*acc(t(i), Rk, Uk) + 2*Rk*Uk^2)/0.01;
end
E_b(k, :) = ps;
end
W = [1 1 0 0; 1 0 -2 0; 1 0 0 2]; names = {'PI', 'AM', 'AMPI'};
refT = [bandE(E_t(1,:), 0, fs/2), bandE(E_t(1,:), 0.5*f, 1.5*f)];
refB = [bandE(E_b(1,:), 0, fs/2), bandE(E_b(1,:), 0.5*f, 1.5*f)];
for s = 1:3
rt = W(s, :) * E_t; rb = W(s, :) * E_b;
T3 = 10*log10([bandE(rt, 0, fs/2)/refT(1), bandE(rt, 0.5*f, 1.5*f)/refT(2), bandE(rt, 1.5*f, 2.5*f)/refT(1)] + 1e-30);
B3 = 10*log10([bandE(rb, 0, fs/2)/refB(1), bandE(rb, 0.5*f, 1.5*f)/refB(2), bandE(rb, 1.5*f, 2.5*f)/refB(1)]);
fprintf('%-26s %-5s | %6.1f %6.1f %7.1f | %6.1f %6.1f %7.1f | %6.1f\n', scen{is, 1}, names{s}, T3, B3, B3(1) - T3(1));
if is == 1, keepT(s, :) = rt; keepB(s, :) = rb; end
end
if is == 1, E1t = E_t(1, :); E1b = E_b(1, :); end
end
fprintf('(all: residual energy re the single-pulse echo; f-band: 0.5f-1.5f re the single echo in that band; 2f-band re the whole single echo; tissue at -300 dB means exact cancellation)\n');
%% Motion law check: PI residual of a linear scatterer = |1 - exp(-j 2 pi f tau)| = 2 |sin(pi f tau)|
v = 0.01; tau = 2*v/(cT*prf);
fprintf('\nMotion: v = 10 mm/s, PRF 5 kHz -> shift %.2f ns, 2|sin(pi f tau)| = %.1f dB (the f-band floor above)\n', tau*1e9, 20*log10(2*abs(sin(pi*f*tau))));
%% Figure: single echoes and residuals
figure('Position', [50 50 900 520]);
subplot(2, 2, 1); plot(t*1e6, E1t/max(abs(E1t)), 'k', t*1e6, keepT(1,:)/max(abs(E1t)), 'b', t*1e6, keepT(2,:)/max(abs(E1t)), 'r'); grid on;
title('tissue: single echo and residuals'); xlabel('time (\mus)'); ylabel('re single echo'); legend('single', 'PI', 'AM');
subplot(2, 2, 2); plot(t*1e6, E1b/max(abs(E1b)), 'k', t*1e6, keepB(1,:)/max(abs(E1b)), 'b', t*1e6, keepB(2,:)/max(abs(E1b)), 'r'); grid on;
title('bubble: single echo and residuals'); xlabel('time (\mus)'); ylabel('re single echo');
fa = (0:nfft/2-1)*df/1e6; sp = @(y) 20*log10(abs(fft(y, nfft)) + 1e-30);
SpT = sp(E1t); SpB = sp(E1b); mT = max(SpT); mB = max(SpB);
subplot(2, 2, 3); s1 = sp(keepT(1,:)); s2 = sp(keepT(2,:));
plot(fa, SpT(1:nfft/2) - mT, 'k', fa, s1(1:nfft/2) - mT, 'b', fa, s2(1:nfft/2) - mT, 'r'); grid on; xlim([0 10]); ylim([-80 5]);
xlabel('f (MHz)'); ylabel('dB'); title('tissue spectra');
subplot(2, 2, 4); s1 = sp(keepB(1,:)); s2 = sp(keepB(2,:));
plot(fa, SpB(1:nfft/2) - mB, 'k', fa, s1(1:nfft/2) - mB, 'b', fa, s2(1:nfft/2) - mB, 'r'); grid on; xlim([0 10]); ylim([-80 5]);
xlabel('f (MHz)'); ylabel('dB'); title('bubble spectra');
What to look for: tissue suppression of 23.4 (PI), 35.4 (AM) and 25.9 dB (AM + PI), bubble retention −7.5, −1.5 and −1.1 dB; with 10 mm/s motion every sequence's fundamental-band tissue residue sits at −27.4 to −27.6 dB, matching the printed \(2\lvert\sin(\pi f\tau)\rvert\) = −27.8 dB; a +0.1 dB half-pulse error leaves −38.6 dB in the AM fundamental band; linear tissue at rest cancels exactly (printed as −300 dB).
m28_tic_perfusion.mLognormal bolus and flash-replenishment fits with a hand-written Levenberg–Marquardt loop; log-compression bias · change Ns, truth, display range · runtime ≈ 2 s
%% m28_tic_perfusion.m -- Time-intensity curves: bolus (lognormal) and flash-replenishment fits
% Purpose : simulate ROI-averaged linear intensity with speckle (gamma-distributed, shape Ns)
% and fit (1) a lognormal bolus model and (2) the Wei flash-replenishment model
% with a hand-written Levenberg-Marquardt loop (numeric Jacobian). Then show why
% fitting log-compressed display values biases the perfusion parameters.
% Models : bolus I(t) = AUC / ((t - t0) s sqrt(2 pi)) exp(-(ln(t - t0) - m)^2 / (2 s^2)) + bg
% MTT = exp(m + s^2/2) (after t0), TTP = t0 + exp(m - s^2), PE = peak above bg
% flash I(t) = A (1 - exp(-beta t)) + bg, flow index A * beta
% Change : Ns (speckle cells in the ROI), frame rate, true parameters, DR of the display map.
% Expect : linear-domain fits recover MTT = 11.97 s and TTP = 9.98 s within a few percent
% (Monte Carlo std printed); fitting the 40 dB log-compressed values lengthens MTT
% by a factor of about 2 and the flash fit recovers A, beta to a few percent.
% Provenance: Monte Carlo speckle on model curves (lognormal, Wei), not patient data.
% Runtime : about 2 s.
clear; close all; rng(1);
%% Bolus truth and acquisition
fr = 4; t = (0:1/fr:60)'; % s, 4 frames per second, 60 s
th = [1000, log(10), 0.6, 3, 2]; % [AUC (a.u. s), m, s, t0 (s), bg]
lgn = @(p, t) (t > p(4)) .* p(1) ./ (max(t - p(4), 1e-9) * p(3) * sqrt(2*pi)) ...
.* exp(-(log(max(t - p(4), 1e-9)) - p(2)).^2 / (2 * p(3)^2)) + p(5);
Ns = 50; % independent speckle cells averaged in the ROI
noisy = @(I) I .* (-sum(log(rand(numel(I), Ns)), 2) / Ns); % mean I, std I/sqrt(Ns)
mtt = @(p) exp(p(2) + p(3)^2/2); ttp = @(p) p(4) + exp(p(2) - p(3)^2);
fprintf('Truth: MTT = %.2f s, TTP = %.2f s, peak above bg = %.1f, AUC = %.0f\n', mtt(th), ttp(th), max(lgn(th, t)) - th(5), th(1));
%% Levenberg-Marquardt, written inline (scripts cannot define local functions portably)
keepFit = cell(1, 2);
nTrial = 100; est = zeros(nTrial, 4); estLog = zeros(nTrial, 2);
DR = 40; % display dynamic range (dB)
for trial = 1:nTrial
y = noisy(lgn(th, t));
for mode = 1:2
if mode == 1
yy = y;
else % naive: fit the log-compressed gray levels
Iref = max(y); yy = 255 * min(max((10*log10(y / Iref) + DR) / DR, 0), 1);
end
[pk, ip] = max(yy); base = median(yy(1:10));
p = [ (pk - base) * 8, log(max(t(ip) - 3, 1)), 0.6, 3, base ]; % initial guess
lam = 1e-2; r = yy - lgn(p, t); cost = r' * r;
for it = 1:200
J = zeros(numel(t), 5);
for k = 1:5
dp = zeros(1, 5); dp(k) = 1e-6 * max(1, abs(p(k)));
J(:, k) = (lgn(p + dp, t) - lgn(p - dp, t)) / (2 * dp(k));
end
Hm = J' * J; g = J' * r;
stepOK = false;
while ~stepOK && lam < 1e12
d = (Hm + lam * diag(diag(Hm) + 1e-12)) \ g;
pn = p + d'; pn(3) = max(pn(3), 0.05); pn(4) = min(pn(4), t(ip) - 0.5);
rn = yy - lgn(pn, t); cn = rn' * rn;
if cn < cost, p = pn; r = rn; dc = cost - cn; cost = cn; lam = lam / 3; stepOK = true; else, lam = lam * 4; end
end
if ~stepOK || dc < 1e-10 * cost, break; end
end
if mode == 1, est(trial, :) = [mtt(p), ttp(p), max(lgn(p, t)) - p(5), p(1)]; else, estLog(trial, :) = [mtt(p), ttp(p)]; end
if trial == 1, keepFit{mode} = struct('y', yy, 'p', p); end
end
end
fprintf('Linear-intensity fits (%d trials, Ns = %d): MTT %.2f +- %.2f s, TTP %.2f +- %.2f s, peak %.1f +- %.1f, AUC %.0f +- %.0f\n', ...
nTrial, Ns, mean(est(:,1)), std(est(:,1)), mean(est(:,2)), std(est(:,2)), mean(est(:,3)), std(est(:,3)), mean(est(:,4)), std(est(:,4)));
fprintf('Log-compressed fits (DR = %d dB): MTT %.2f +- %.2f s (x%.2f), TTP %.2f +- %.2f s\n', ...
DR, mean(estLog(:,1)), std(estLog(:,1)), mean(estLog(:,1)) / mtt(th), mean(estLog(:,2)), std(estLog(:,2)));
%% Flash-replenishment (Wei et al. 1998 model)
tf = (0:1/fr:15)'; thF = [50, 0.8, 2]; % A, beta (1/s), bg
mdlF = @(p, t) p(1) * (1 - exp(-p(2) * t)) + p(3);
estF = zeros(nTrial, 3);
for trial = 1:nTrial
y = noisy(mdlF(thF, tf));
p = [max(y) - y(1), 1, y(1)]; lam = 1e-2; r = y - mdlF(p, tf); cost = r' * r;
for it = 1:200
J = zeros(numel(tf), 3);
for k = 1:3
dp = zeros(1, 3); dp(k) = 1e-6 * max(1, abs(p(k)));
J(:, k) = (mdlF(p + dp, tf) - mdlF(p - dp, tf)) / (2 * dp(k));
end
Hm = J' * J; g = J' * r; stepOK = false;
while ~stepOK && lam < 1e12
d = (Hm + lam * diag(diag(Hm) + 1e-12)) \ g; pn = p + d'; pn(2) = max(pn(2), 1e-3);
rn = y - mdlF(pn, tf); cn = rn' * rn;
if cn < cost, p = pn; r = rn; dc = cost - cn; cost = cn; lam = lam / 3; stepOK = true; else, lam = lam * 4; end
end
if ~stepOK || dc < 1e-10 * cost, break; end
end
estF(trial, :) = p;
if trial == 1, keepF = struct('y', y, 'p', p); end
end
fprintf('Flash-replenishment fits: A %.1f +- %.1f (true 50), beta %.3f +- %.3f 1/s (true 0.8), A*beta %.1f +- %.1f (true 40)\n', ...
mean(estF(:,1)), std(estF(:,1)), mean(estF(:,2)), std(estF(:,2)), mean(estF(:,1).*estF(:,2)), std(estF(:,1).*estF(:,2)));
%% Figures
figure('Position', [50 50 960 330]);
subplot(1, 3, 1); plot(t, keepFit{1}.y, '.', 'Color', [0.6 0.6 0.6]); hold on; plot(t, lgn(th, t), 'k--', t, lgn(keepFit{1}.p, t), 'b'); grid on;
xlabel('time (s)'); ylabel('linear intensity (a.u.)'); title('bolus: lognormal fit'); legend('data', 'truth', 'fit');
subplot(1, 3, 2); plot(t, keepFit{2}.y, '.', 'Color', [0.6 0.6 0.6]); hold on; plot(t, lgn(keepFit{2}.p, t), 'r'); grid on;
xlabel('time (s)'); ylabel('gray level (40 dB map)'); title('same data, log-compressed');
subplot(1, 3, 3); plot(tf, keepF.y, '.', 'Color', [0.6 0.6 0.6]); hold on; plot(tf, mdlF(thF, tf), 'k--', tf, mdlF(keepF.p, tf), 'b'); grid on;
xlabel('time after flash (s)'); ylabel('linear intensity (a.u.)'); title('flash-replenishment');
What to look for: over 100 noise realizations the linear fit gives MTT 11.96 ± 0.33 s (truth 11.97) and TTP 9.94 ± 0.16 s (truth 9.98); the same fit on 40 dB log-compressed data gives MTT ≈ 26.5 s (×2.2) and TTP 10.5 s; the flash fit returns \(A\) = 50.1 ± 1.8, \(\beta = 0.79 \pm 0.10\ \mathrm{s^{-1}}\).
m28_ulm_concentration.mMonte Carlo of PSF isolation versus density, optimal density, conversion from blood concentration, capillary waiting times, Doppler of slow bubbles · change f, PSF, k, C, vessels · runtime ≈ 7 s
%% m28_ulm_concentration.m -- ULM: isolation probability, optimal density, capillary waiting time
% Purpose : bubbles in the image plane are a Poisson process of areal density n (per mm^2).
% A bubble is usable when no other bubble lies inside the ellipse
% (dx/(k Fx))^2 + (dz/(k Fz))^2 < 1 (Fx, Fz = PSF FWHM). Then
% P_iso = exp(-n A), A = pi k^2 Fx Fz, isolated per mm^2 = n exp(-n A),
% maximal at n* = 1/A with P_iso = 1/e. Monte Carlo (periodic box) checks this.
% Second constraint: a capillary of diameter d and speed v receives bubbles at the
% rate C Q = C (pi d^2/4) v, so its first detection waits about 1/(C Q).
% Change : f, PSF widths (in wavelengths), k, blood concentration C, slice thickness, vessel table.
% Expect : Monte Carlo P_iso within about 0.01 of exp(-nA); at 5 MHz (Fx = 1.0, Fz = 0.75 lambda, k = 1)
% A = 0.224 mm^2 and n* = 4.5 per mm^2; at C = 140 per uL an 8 um capillary at
% 1 mm/s waits about 140 s between bubbles.
% Provenance: Poisson Monte Carlo checked against the exact exp(-nA) law.
% Runtime : about 6 s in Octave.
clear; close all; rng(1);
c = 1540; f = 5e6; lam = c / f * 1e3; % mm
Fx = 1.0 * lam; Fz = 0.75 * lam; k = 1; % PSF FWHM (mm) and exclusion factor
A = pi * k^2 * Fx * Fz; % exclusion area (mm^2)
fprintf('lambda = %.3f mm, Fx = %.3f mm, Fz = %.3f mm, exclusion area A = %.4f mm^2, n* = 1/A = %.2f per mm^2\n', lam, Fx, Fz, A, 1/A);
%% Monte Carlo in a periodic 10 x 10 mm box
L = 10; ns = [0.5 1 2 4.47 8 16]; % per mm^2
fprintf('\n n (/mm^2) | nA | P_iso MC | exp(-nA) | isolated per mm^2 (MC)\n');
nTrial = 60; Pmc = zeros(size(ns)); Iso = zeros(size(ns));
for in = 1:numel(ns)
cnt = 0; tot = 0;
for tr = 1:nTrial
N = round(ns(in) * L^2); % fixed count: P_iso = (1 - A/L^2)^(N-1), which is exp(-nA) to < 0.1 %
xy = rand(N, 2) * L;
dx = abs(xy(:, 1) - xy(:, 1)'); dx = min(dx, L - dx); % periodic distances
dz = abs(xy(:, 2) - xy(:, 2)'); dz = min(dz, L - dz);
inside = (dx / (k*Fx)).^2 + (dz / (k*Fz)).^2 < 1;
inside(1:N+1:end) = false;
cnt = cnt + sum(~any(inside, 2)); tot = tot + N;
end
Pmc(in) = cnt / tot; Iso(in) = cnt / (nTrial * L^2);
fprintf(' %8.2f | %5.2f | %.4f | %.4f | %.3f\n', ns(in), ns(in)*A, Pmc(in), exp(-ns(in)*A), Iso(in));
end
fprintf('Theoretical maximum: %.3f isolated per mm^2 at n* (= 1/(e A))\n', 1/(exp(1)*A));
%% From blood concentration to image density, and the capillary waiting time
C = 140; % bubbles per uL (= per mm^3); about 2.4 mL x 3e8/mL diluted in 5 L
slice = 0.5; % mm, elevation slice thickness
vf = 0.05; % vascular volume fraction
n = C * slice * vf; % per mm^2 of image
fprintf('\nC = %d per uL, slice %.1f mm, vascular fraction %.0f %%: n = %.2f per mm^2, nA = %.2f, P_iso = %.2f\n', C, slice, 100*vf, n, n*A, exp(-n*A));
ves = [8 1; 20 3; 50 10; 200 50]; % diameter (um), speed (mm/s)
fprintf('\n vessel d (um) | v (mm/s) | Q (nL/s) | bubbles per s at C | mean wait (s)\n');
for iv = 1:size(ves, 1)
d = ves(iv, 1) * 1e-3; v = ves(iv, 2); % mm, mm/s
Q = pi * d^2 / 4 * v; % mm^3/s
fprintf(' %6.0f | %5.1f | %8.4f | %9.4g | %9.3g\n', ves(iv, 1), v, Q * 1e3, C * Q, 1 / (C * Q));
end
%% Doppler shift of a slow bubble: why SVD cannot separate capillary flow from tissue motion
for v = [0.5 1 10] * 1e-3
fprintf('Axial velocity %.1f mm/s at 5 MHz: Doppler %.1f Hz; displacement per frame at 1 kHz = %.1f um\n', v*1e3, 2*v*f/c, v/1000*1e6);
end
%% Figure
nn = logspace(-1, 2, 200);
figure('Position', [50 50 900 330]);
subplot(1, 2, 1); semilogx(nn, exp(-nn*A), 'b', ns, Pmc, 'ko'); hold on;
semilogx(nn, nn .* exp(-nn*A) * A * exp(1), 'r'); grid on;
xlabel('bubbles per mm^2 of image'); ylabel('fraction');
legend('P_{iso} theory', 'P_{iso} Monte Carlo', 'isolated per area / max', 'Location', 'west'); title('5 MHz, k = 1');
subplot(1, 2, 2); Cs = logspace(0, 4, 200); d = 8e-3; v = 1; Q = pi*d^2/4*v;
loglog(Cs, 1 ./ (Cs * Q), 'b'); grid on; xlabel('blood concentration (bubbles per \muL)'); ylabel('mean wait (s)');
title('8 \mum capillary at 1 mm/s');
What to look for: \(A\) = 0.2235 mm² and \(n^{*}\) = 4.47 per mm² at 5 MHz; the Monte Carlo \(P_{\mathrm{iso}}\) follows \(e^{-nA}\) within about 0.01 and the isolated count peaks near 1.65 per mm²; \(C\) = 140 per µL gives \(n\) = 3.5 per mm² and \(P_{\mathrm{iso}}\) = 0.46; the 8 µm capillary waits 142 s, a 50 µm arteriole 0.36 s; a 1 mm/s bubble has a 6.5 Hz Doppler at 5 MHz.
28.14Self-check
Q1What is the resonance of a free 3 µm air bubble by Minnaert and with surface tension?
3.29/3 = 1.10 MHz by Minnaert; with \(\sigma\) = 0.072 N/m, \(\rho\omega_0^{2}R_0^{2} = 3(1.4)(101.3\ \mathrm{kPa})\) \(+\, 2(3.2)(0.072)/(3\ \mathrm{\mu m})\) = 425.6 + 153.6 kPa, so \(f_0\) = 1.28 MHz (+17 %).
Q2A lipid bubble with \(R_0\) = 1.5 µm, \(\kappa\) = 1.07, \(\chi\) = 0.5 N/m and \(\sigma_0\) = 0.02 N/m. Estimate \(f_0\).
\(3\kappa P_0\) = 0.325 MPa; \(2(3\kappa + 1)\sigma_0/R_0 = 2(4.21)(0.02)/1.5\ \mathrm{\mu m}\) = 0.112 MPa; \(4\chi/R_0\) = 1.333 MPa. Sum 1.77 MPa; \(f_0 = \sqrt{1.77\times 10^{6}/1000}\,/\,(2\pi\times 1.5\ \mathrm{\mu m})\) = 4.46 MHz (28.3).
Q3For a free 2 µm air bubble at resonance, give \(Q\) and \(\sigma_s/\pi R_0^{2}\).
\(\delta_{\mathrm{rad}} = \omega R_0/c = 0.017\), \(\delta_{\mathrm{vis}} = 4\mu/(\rho\omega R_0^{2}) = 0.078\) at 2.04 MHz: \(\delta\) = 0.095, \(Q\) = 10.5, \(\sigma_s/\pi R_0^{2} = 4/\delta^{2} = 440\) (thermal damping would lower both).
Q4Which sequence keeps the cubic term, with what weight, and in which band does it appear?
AM (\(e(1) - 2e(\tfrac12)\)) and AM + PI (\(e(1) + 2e(-\tfrac12)\)) both keep \(\tfrac34 a_3p^{3}\); \(p^{3}\) has components at \(f_0\) and \(3f_0\), and the \(f_0\) part is the "nonlinear fundamental". PI cancels all odd orders.
Q5Tissue moves at 5 mm/s, PRF 2 kHz, f = 3 MHz. What is the PI tissue residue?
\(\tau = 2v/(c\cdot\mathrm{PRF}) = 0.01/(1540\times 2000)\) = 3.25 ns; \(2\lvert\sin(\pi f\tau)\rvert = 2\sin(0.0306) = 0.061\), i.e. −24.3 dB.
Q6A bolus fits a lognormal with \(\mu = \ln 8\) (s) and \(s\) = 0.5, \(t_0\) = 4 s. Give MTT and TTP.
\(\mathrm{MTT} = 8e^{0.125}\) = 9.07 s (from \(t_0\)); \(\mathrm{TTP} = t_0 + 8e^{-0.25} = 4 + 6.23\) = 10.23 s (28.11).
Q7At 7.5 MHz with \(F_x = \lambda\), \(F_z = 0.75\lambda\) and \(k\) = 1, what bubble density maximizes isolated localizations, and how many per mm² per frame?
\(\lambda\) = 0.205 mm; \(A = \pi(0.205)(0.154)\) = 0.0993 mm²; \(n^{*}\) = 10.1 per mm²; maximum \(1/(eA)\) = 3.7 isolated bubbles per mm² per frame.
Q8Why does a 1 µm free bubble need more pressure than a 2 µm one to emit a subharmonic at \(2f_0\)?
Viscous damping scales as \(1/(\omega R_0^{2})\): \(\delta_{\mathrm{vis}}\) = 0.134 at 1 µm versus 0.078 at 2 µm (each at its own resonance), so the parametric growth needs a larger drive. In the model at 200 kPa the 1 µm bubble's \(f/2\) is −26 dB, the 2 µm bubble's −4 dB.
28.15Figures from the literature
Measured data behind the models of this module. Figure numbers checked against the publisher or author copies; links only, no images reproduced.
The buckling–rupture shell, and a real bubble that follows it
P. Marmottant, S. van der Meer, M. Emmer, M. Versluis, N. de Jong, S. Hilgenfeldt and D. Lohse, J. Acoust. Soc. Am. 118, 3499–3505 (2005). doi:10.1121/1.2109427
Look at Fig. 2 (the \(\sigma(R)\) law of 28.3, (28.8): zero when buckled, elastic slope, water tension after rupture), Fig. 3 (equilibrium radius versus ambient pressure for \(R_0\) = 2 µm, \(R_b\) = 1.98 µm, \(\chi\) = 1 N/m, i.e. \(\sigma_0\) = 0.020 N/m, the value used here) and Fig. 5 (Brandaris high-speed radius–time curve of a SonoVue bubble at 2.9 MHz and 130 kPa against the model). Why: it is the measurement that makes the model worth using; the SonoVue fit used \(\chi\) = 1 N/m and \(\kappa_s = 1.5\times 10^{-8}\) kg/s, so this module's \(\chi\) = 0.5 N/m and \(\kappa_s = 5\times 10^{-9}\) kg/s are deliberately softer and less damped.
Compression-only behaviour seen optically
N. de Jong, M. Emmer, C. T. Chin, A. Bouakaz, F. Mastik, D. Lohse and M. Versluis, Ultrasound Med. Biol. 33, 653–656 (2007). doi:10.1016/j.ultrasmedbio.2006.09.016
Look at Fig. 2 (diameter–time curves and spectra of a regular bubble and a compression-only bubble) and Fig. 3 (expansion/compression ratio versus diameter at 1, 1.8 and 4 MHz and 50–200 kPa). Why: about 40 % of the 101 bubbles recorded compressed far more than they expanded, from 50 kPa: the asymmetry of Figure 28.B and Rule "Lipid shells buckle" is a measured property of real agents, not a model artefact.
Shell nonlinearity at a few kPa: skewed resonance curves
M. Overvelde, V. Garbin, J. Sijl, B. Dollet, N. de Jong, D. Lohse and M. Versluis, "Nonlinear shell behavior of phospholipid-coated microbubbles," Ultrasound Med. Biol. 36, 2080–2092 (2010). doi:10.1016/j.ultrasmedbio.2010.08.015
Look at Fig. 6 (measured resonance curves of a 3.2 µm bubble at 7.5, 12.5 and 25 kPa against uncoated, linear viscoelastic and buckling models) and Fig. 7 (frequency of maximum response versus pressure). Why: the resonance of a lipid bubble moves with pressure already at tens of kPa, so the linear formula of 28.2 (and Demo 28.1) is a small-amplitude limit; only the buckling model reproduces the skew. Their Fig. 2 parameters (\(\chi\) = 2.5 N/m, \(\kappa_s = 6\times 10^{-9}\) kg/s) show again how much fitted shells differ.
Fitted shell parameters depend on bubble size
A. A. Doinikov, J. F. Haac and P. A. Dayton, "Modeling of nonlinear viscous stress in encapsulating shells of lipid-coated contrast agent microbubbles," Ultrasonics 49, 269–275 (2009). doi:10.1016/j.ultras.2008.09.007 (open copy: PMC2826140)
Look at Fig. 1 (best-fit shell viscosity \(\kappa_s\) and elasticity \(\chi\) versus initial radius) and Fig. 4 (fit of a measured compression-only radius–time curve). Why: a constant \(\kappa_s\) cannot describe all sizes, which is the limit stated under the low-Q rule; a shear-thinning shell viscosity is one proposed fix.
The first in vivo ultrafast ULM
C. Errico, J. Pierre, S. Pezet, Y. Desailly, Z. Lenkei, O. Couture and M. Tanter, Nature 527, 499–502 (2015). doi:10.1038/nature16066
Look at Fig. 1 (principle: ultrafast frames, bubble detection, localization, accumulation, the pipeline of Figure 28.E) and Fig. 2 (spatial resolution and quantification in the rat cortex through a thinned skull). Why: vessels below 10 µm at more than 10 mm depth from more than 500 frames per second and 75 000 images, about \(10^{6}\) events per coronal plane: the numbers that make the density and time trade-offs of 28.8 concrete.
Acquisition time is set by microvascular flow
V. Hingot, C. Errico, B. Heiles, L. Rahal, M. Tanter and O. Couture, Sci. Rep. 9, 2456 (2019). doi:10.1038/s41598-018-38349-x (open access)
Look at Fig. 1 (one 50 µm vessel reconstructed with 5 µm pixels after 240 s) and Fig. 5 (vessel reconstruction versus acquisition time). Why: it measures the saturation that Rule "Capillaries need minutes" and the last visual-summary card (fraction visited) compute from \(1 - e^{-CQt}\): small, slow vessels fill last.
Benchmarking the localization step (PALA)
B. Heiles, A. Chavignon, V. Hingot, P. Lopez, E. Teston and O. Couture, "Performance benchmarking of microbubble-localization algorithms for ultrasound localization microscopy," Nat. Biomed. Eng. 6, 605–616 (2022). doi:10.1038/s41551-021-00824-8
Look at Fig. 2 (simulation framework and metrics), Fig. 3 (global score per algorithm) and Fig. 5 (PSF non-uniformity in direction and its effect on localization). Why: it turns "which localizer?" into numbers on shared simulated and in vivo data, and shows that PSF shape, not only width, biases sub-pixel estimates: the step 5 of Figure 28.E that innovation most often targets.
28.16Go deeper
- T. G. Leighton, The Acoustic Bubble, Academic Press (1994). The reference for bubble dynamics, damping and cross-sections.
- M. Minnaert, "On musical air-bubbles and the sounds of running water," Philosophical Magazine 16, 235–248 (1933).
- P. S. Epstein and M. S. Plesset, "On the stability of gas bubbles in liquid-gas solutions," J. Chem. Phys. 18, 1505–1509 (1950).
- J. B. Keller and M. Miksis, "Bubble oscillations of large amplitude," J. Acoust. Soc. Am. 68, 628–633 (1980).
- C. C. Church, "The effects of an elastic solid surface layer on the radial pulsations of gas bubbles," J. Acoust. Soc. Am. 97, 1510–1521 (1995).
- P. Marmottant, S. van der Meer, M. Emmer, M. Versluis, N. de Jong, S. Hilgenfeldt and D. Lohse, "A model for large amplitude oscillations of coated bubbles accounting for buckling and rupture," J. Acoust. Soc. Am. 118, 3499–3505 (2005).
- N. de Jong, M. Emmer, C. T. Chin, A. Bouakaz, F. Mastik, D. Lohse and M. Versluis, ""Compression-only" behavior of phospholipid-coated contrast bubbles," Ultrasound Med. Biol. 33, 653–656 (2007).
- J. Sijl, B. Dollet, M. Overvelde, V. Garbin, T. Rozendal, N. de Jong, D. Lohse and M. Versluis, "Subharmonic behavior of phospholipid-coated ultrasound contrast agent microbubbles," J. Acoust. Soc. Am. 128, 3239–3252 (2010).
- M. Versluis, E. Stride, G. Lajoinie, B. Dollet and T. Segers, "Ultrasound contrast agent modeling: a review," Ultrasound Med. Biol. 46, 2117–2144 (2020).
- D. H. Simpson, C. T. Chin and P. N. Burns, "Pulse inversion Doppler: a new method for detecting nonlinear echoes from microbubble contrast agents," IEEE Trans. Ultrason. Ferroelectr. Freq. Control 46, 372–382 (1999).
- K. Wei, A. R. Jayaweera, S. Firoozan, A. Linka, D. M. Skyba and S. Kaul, "Quantification of myocardial blood flow with ultrasound-induced destruction of microbubbles administered as a constant venous infusion," Circulation 97, 473–483 (1998).
- M.-X. Tang, H. Mulvana, T. Gauthier, A. K. P. Lim, D. O. Cosgrove, R. J. Eckersley and E. Stride, "Quantitative contrast-enhanced ultrasound imaging: a review of sources of variability," Interface Focus 1, 520–539 (2011).
- C. Demené et al., "Spatiotemporal clutter filtering of ultrafast ultrasound data highly increases Doppler and fUltrasound sensitivity," IEEE Trans. Med. Imaging 34, 2271–2285 (2015).
- C. Errico, J. Pierre, S. Pezet, Y. Desailly, Z. Lenkei, O. Couture and M. Tanter, "Ultrafast ultrasound localization microscopy for deep super-resolution vascular imaging," Nature 527, 499–502 (2015).
- V. Hingot, C. Errico, B. Heiles, L. Rahal, M. Tanter and O. Couture, "Microvascular flow dictates the compromise between spatial resolution and acquisition time in ultrasound localization microscopy," Sci. Rep. 9, 2456 (2019).
- M. Overvelde et al., "Nonlinear shell behavior of phospholipid-coated microbubbles," Ultrasound Med. Biol. 36, 2080–2092 (2010).
- B. Heiles et al., "Performance benchmarking of microbubble-localization algorithms for ultrasound localization microscopy," Nat. Biomed. Eng. 6, 605–616 (2022). Open code and data (PALA).
- K. Christensen-Jeffries et al., "Super-resolution ultrasound imaging," Ultrasound Med. Biol. 46, 865–891 (2020).
- US prescribing information for Lumason, Definity and Optison (composition, size and concentration figures in 28.1); EFSUMB CEUS guidelines (non-hepatic update 2017) for low-MI practice.