Shell viscosity estimation of lipid-coated microbubbles

Abstract

Understanding the shell rheology of ultrasound contrast agent microbubbles is vital for anticipating their bioeffects in clinical practice. Past studies using sophisticated acoustic and optical techniques have made enormous progress in this direction, enabling the development of shell models that adequately reproduce the nonlinear behaviour of the coated microbubble under acoustic excitation. However, there have also been puzzling discrepancies and missing physical explanations for the dependency of shell viscosity on the equilibrium bubble radius, which demands further experimental investigations. In this study, we aim to unravel the cause of such behaviour by performing a refined characterisation of the shell viscosity. We use ultra-high-speed microscopy imaging, optical trapping and wide-field fluorescence to accurately record the individual microbubble response upon ultrasound driving across a range of bubble sizes. An advanced model of bubble dynamics is validated and employed to infer the shell viscosity of single bubbles from their radial time evolution. The resulting values reveal a prominent variability of the shell viscosity of about an order of magnitude and no dependency on the bubble size, which is contrary to previous studies. We find that the method called bubble spectroscopy, which has been used extensively in the past to determine the shell viscosity, is highly sensitive to methodology inaccuracies, and we demonstrate through analytical arguments that the previously reported unphysical trends are an artifact of these biases. We also show the importance of correct bubble sizing, as errors in this aspect can also lead to unphysical trends in shell viscosity, when estimated through a nonlinear fitting from the time response of the bubble.

---

1. Introduction

During the past decades, the development of stabilised gas microbubbles smaller than a red blood cell (<10 μm), and thus able to navigate the thinnest capillaries, has opened unique opportunities to utilise the effects of ultrasonic cavitation in medicine. Upon high-amplitude ultrasound driving, the microbubbles are set into a periodic succession of growth and collapse phases. The microbubbles thus pump energy out of the acoustic wave and release it locally in the form of sound (even at frequencies other than the driving frequency) and mechanical action. This extraordinary energy-focusing ability makes a simple bubble an extremely capable device for a large number of biomedical applications. Coated microbubbles are currently used as ultrasound contrast agents to improve the contrast in perfusion imaging of the myocardium,^1,2 liver,³ brain,⁴ kidneys⁵ and other organs, providing unique information on blood volume and velocity owing to their large echogenicity. The imaging potential of microbubbles can be augmented by equipping the coating with specific ligands. Site-targeted contrast agents can be used to non-invasively image molecular events in vivo such as inflammation,^6,7 angiogenesis^8,9 and early tumour formation.^10,11 Moreover, microbubbles are presently investigated for therapeutic applications. The high-energy mechanical effects of bubble cavitation (e.g. acoustic microstreaming¹² and microjetting¹³) can be harnessed to promote clot lysis in acute ischemic stroke,^14,15 deliver plasmid DNA and drugs in the immediate perivascular region for myocardial infarction,¹⁶ atherosclerosis¹⁷ and tumours treatment^18,19 or to locally permeabilise the blood–brain barrier to enhance the delivery of therapeutic agents for brain cancer,^20–22 Alzheimer's disease²³ and amyotrophic lateral sclerosis²⁴ treatment.

A prolonged intravascular stability against dissolution and coalescence is achieved by encapsulating the microbubbles with a protein, biopolymer or lipid shell and by using a high-molecular-mass gas content such as perfluorocarbons and sulphur hexafluoride, which are less prone to outward loss compared to air.²⁵ To date, most clinically approved ultrasound contrast agents are coated with a phospholipid shell which is more flexible than other types of coatings, yielding less stable but more echogenic microbubbles. Lipids naturally adsorb on gas–liquid interfaces and self-assemble into a monolayer which lowers the gas–liquid surface tension and thus the large Laplace pressure, arresting the gas efflux. Moreover, it imparts rheological properties to the surface.²⁶ Monolayers of DPPC (1,2-distearoyl-sn-glycero-3-phosphocholine), a frequent main lipid constituent of the bubble shell, show at low harmonic driving frequencies a dominantly viscous behaviour similar to three-dimensional liquids. On the other hand, at higher frequencies their microstructure dominates the interfacial rheology and induces an additional elastic contribution.²⁷ This elasticity is of intrinsic rheological origin and is added to the Gibbs (apparent) monolayer elasticity, which arises when interface deformations yield variations of surfactant concentration.²⁸ The properties of the shell have a fundamental impact on the dynamical behaviour of ultrasound contrast agent microbubbles, and consequently, on their biophysical effects. Compared to uncoated bubbles, coated bubbles show a higher resonance frequency due to a stiffer interface, a lower amplitude of response owing to viscous dissipation mechanisms between lipid molecules and a plethora of nonlinear effects.²⁹

Unfortunately, characterising the shell of a microbubble under clinically-relevant working conditions poses significant challenges owing to the micrometric size and the extremely high driving frequencies (∼10⁶ Hz) to which these bubbles are subjected. Direct experimental rheological techniques such as Langmuir trough and capillary pressure tensiometry methods are limited to significantly lower driving frequencies (∼1 Hz and ∼10² Hz, respectively)³⁰ and at least one order of magnitude smaller film curvatures.³¹ Therefore, as current direct approaches are not adequate for characterising the shell properties, only indirect methods can be used. The seminal works of de Jong et al.,³² Frinking & de Jong,³³ Hoff et al.³⁴ and Gorce et al.³⁵ provided the first estimates of the elastic and viscous contribution of the shell by performing acoustic attenuation measurements of polydisperse suspensions of coated microbubbles. This approach was later refined by Parrales et al.³⁶ and Segers et al.³⁷ by employing monodisperse bubble suspensions produced by microfluidic flow-focusing techniques. Moreover, Helfield et al.³⁸ estimated the shell properties from scattered pressure measurements of single bubbles. Morgan et al.³⁹ provided the first optical characterisation by fitting a shell model to the radial evolution, captured on a streak camera image, of single ultrasound-driven microbubbles. The same system was later also used by Doinikov et al.⁴⁰ Only with the development of the groundbreaking Brandaris camera⁴¹ was the curtain on the wide-field time-resolved dynamics of the microbubble finally lifted. van der Meer et al.⁴² and Overvelde et al.⁴³ were the first to leverage this new potential to characterise the shell of single DSPC-based microbubbles, followed by Luan et al.⁴⁴ to investigate the effect of liposomes loaded on the shell and van Rooij et al.⁴⁵ to evaluate the influence of the carbon chain length of the phospholipids on the shell properties. A side laser scattering system was used by Tu et al.^46,47 to characterise the shell properties of commercial phospholipid-coated microbubbles by measuring their response to an ultrasound pulse. The sensitivity of this technique has been further improved by employing a forward laser-scattering system that enabled Lum et al.⁴⁸ to measure the decaying free oscillations and Lum et al.,⁴⁹ Supponen et al.⁵⁰ and Daeichin et al.⁵¹ to retrieve the frequency response of single bubbles set into small-amplitude oscillations via photoacoustic forcing. Such small oscillations (few tens of nanometers) simplify the characterisation, as the absence of nonlinear effects allows the use of linear models.

In contrast to the dilatational modulus, for which there is broad empirical evidence regarding its independence from the bubble size, the reported dilatational viscosity values have been shown to increase with the bubble equilibrium radius and with a growing rate that strongly varies from study to study.^{36–40,42,44–47,50,51} This observed dependency is widely regarded as unphysical, as it is illogical for a material property to be dependent on the quantity of the material. van der Meer et al.,⁴² Doinikov et al.⁴⁰ and Tu et al.⁴⁷ also showed the shell viscosity to decrease with the (maximum) strain rate of the bubble, suggesting a strain-thinning behaviour to possibly explain its questionable dependency on the bubble size. However, this conclusion contradicts the observations from Helfield et al.³⁸ and van Rooij et al.⁴⁵ who found no significant relationship between shell viscosity and strain rate. Moreover, the results reported by Tu et al.⁴⁷ are particularly puzzling: as the parameters constituting the maximum strain rate (Ṙ/R)_max ≈ 2πfΔR_max/R₀—bubble resting radius R₀, maximum radial expansion (through pressure) ΔR_max and driving frequency f—vary, the shell viscosity gives rise to different families of curves instead of collapsing on the same universal curve. In our view, this means that the dependence of shell viscosity on the initial radius cannot be explained by a dependence on the strain rate. Furthermore, the assumption that the shell viscosity is dependent on the strain rate throughout the entire radial excursion of the bubble pursuant to the Cross model, as proposed by Doinikov et al.⁴⁰ following the result from van der Meer et al., is questionable. In fact, during the expansion phase of the bubble, the phospholipids separate from each other, leading to a loss of the microstructure responsible for the non-Newtonian characteristics of the shell beyond a certain level of interfacial strain. Conversely, during the compression phase of the bubble, the in-plane strain of the microstructure may have a minor role compared to the out-of-plane buckling. Most studies reporting a dependency of shell viscosity on the bubble size have estimated the shell properties by analysing the bubble response in the frequency domain, adopting the “bubble spectroscopy” technique introduced by van der Meer et al.⁴² or a variation of it. The remaining works instead inferred the shell properties by fitting the experimental bubble radius–time curves.

In this work, we aim to clarify the apparent dependency of the shell viscosity on the microbubble size and the unexplained discrepancies in the values reported in different studies for similar bubbles. We use ultra-high-speed microscopic imaging and optical trapping to directly observe the response of single ultrasonically driven coated microbubbles in an unbound fluid environment across a range of bubble radii. We characterise the shell dilatational viscosity examining the bubble response in the time domain, instead of adopting the common bubble spectroscopy approach. Contrary to other studies which fitted the bubble time response, we first confronted our experimental data against the theoretical predictions to verify that the chosen bubble model is capable of describing the salient features of the bubble response, in particular the bubble resonance size. The shell dilatational viscosity is then inferred by fitting the experimental response with the theoretical one. Finally, the results are compared with previous studies and an explanation for the origin of the unphysical dependence of shell viscosity on the microbubble size is provided.

2. Theory

The radial motion of a spherical bubble of radius R(t) is described by a modified form of the Rayleigh–Plesset equation for mildly compressible Newtonian media,^52–55 which includes a pressure term Σ(R, Ṙ) to account for the generalised interfacial stresses σ_s(R, Ṙ) and readswhere over-dots denote time differentiation, ρ_l the liquid density, c_l the speed of sound in the medium, p_g the gas pressure inside the bubble, p_∞ the undisturbed ambient pressure, p_d(t) the acoustic driving pressure and μ_l the dynamic viscosity of the medium.

The pressure contribution of the interface Σ(R, Ṙ) is given by the radial component of the surface divergence of the interfacial stress^26,29

Σ(R, Ṙ) = [∇_s·σ_s(R, Ṙ)]_r, (2)

where is the surface gradient operator, is the surface unit tensor and n is the unit vector normal to the interface. As opposed to earlier models that treated the shell as a three-dimensional solid,^34,39,56 Chatterjee and Sarkar^57,58 were the first to suggest modelling the bubble coating as a two-dimensional continuum. This approach is considered the most appropriate for molecular membranes (i.e. phospholipid⁵⁹ and protein⁶⁰ shells) for which the assumption of continuum in the normal direction is questionable.⁶¹ They proposed to model the microbubble shell as a viscoelastic solid interface, expressing the surface stress tensor, owing to the spherical symmetry of the problem, as

σ_s(R, Ṙ) = [σ₀ + E_s(J − 1)]I_s + κ_s(∇_s·Ṙn)I_s. (3)

The first term represents the isotropic part of the Evans–Skalak model for elastic interfaces⁶² without differentiating the thermodynamic contribution from the rheological one^63,64 and the second the isotropic part of the Boussinesq–Scriven model for viscous interfaces²⁶ which is the surface equivalent of the Newtonian model for compressible bulk fluids. σ₀ is the interfacial surface tension at the equilibrium bubble radius R₀, E_s is the interfacial dilatational modulus, κ_s is the interfacial dilatational viscosity and J = R²/R₀² is the relative area deformation. Upon substitution of eqn (3) into eqn (2) and by making use of the following differential identities, valid for a spherically symmetric problem, ∇_s·n = 2/R and ∇_s·I_s = −(2/R)n, the pressure contribution given by the interface readsTo account for the strong nonlinear effects that coated microbubbles exhibit, such as subharmonic oscillations,^65–69 compression-only behaviour,⁷⁰ thresholding behaviour⁷¹ and resonance frequency shifting,⁴³ which a quasi-linear shell model as eqn (4) cannot represent, Marmottant et al.⁷⁰ proposed a phenomenological extension inspired by the behaviour of liquid–gas interfaces covered by insoluble surfactants,²⁷ yieldingwhich accounts for: (i) the buckling of the coating when the bubble is compressed below its shell compression limit radius R_buckling resulting in a tensionless interface, (ii) the rupturing of the coating when the bubble is expanded beyond R_rupture where it experiences a surface tension equal to that of a pure gas–water interface σ_water = 72.8 mN m⁻¹, and (iii) the elastic regime in between. The Marmottant model has proven to be a compelling engineering representation of the microbubble shell behaviour, capable of explaining a broad spectrum of nonlinearities,^43,69 and to be sophisticated enough to accurately predict the behaviour of bubbles in the present study, without the need for additional free parameters. However, this model is not entirely free from physical inconsistencies, an important one being the assumption of a constant shell viscosity for varying the free molecular area. Nevertheless, this can be justified as being an average viscosity over the radial excursion of the bubble.

In order to close the problem, the gas pressure contained in the bubble must be determined, in principle, by solving the full Navier–Stokes equation inside the bubble. This can be done analytically only for small-amplitude oscillations for which the governing equations can be linearised^72,73 or numerically at a high computational cost. For these reasons, extensive use is made of the crude approximation of a polytropic process:

with p_g,0 the bubble internal pressure at rest and n the polytropic index. In spite of its appealing simplicity, this approximation suffers from serious limitations.⁷⁴ First, its range of validity is restricted to low-amplitude oscillations and to global Péclet numbers Pe = R₀²ω/D_g (where ω = 2πf is the driving angular frequency and D_g is the gas thermal diffusivity) much smaller or much bigger than one, and therefore to the particular cases where the gas follows an isothermal process (n = 1) or an adiabatic one (n = γ, where γ is the gas specific heat ratio), respectively, for the whole radial motion of the bubble. And second, its use results in no energy loss associated with the change of temperature of the gas. These restrictions are unfortunate in the context of ultrasound contrast agent microbubbles because under clinical conditions they often undergo large nonlinear oscillations, characterised by a thermal motion defined by a global Pe ∼ 1 with both isothermal and adiabatic behaviours occurring at different phases of oscillation, and subjected to a thermal loss of the same order of magnitude of the viscous damping. To address the latter point, it is customary in the literature to artificially increase the liquid viscosity doubling its value.⁷⁵ These considerations motivate us to take a step beyond the polytropic approximation. In an attempt to simplify the problem, Nigmatulin et al.⁷⁶ and Prosperetti et al.⁷⁴ independently found the spatial variation of the pressure inside the bubble to be negligible in most cases (Δp_g/p_g ∼ O(RṘ/λc_g, Ṙ²/c_g²), where λ and c_g are representative values of the sound wavelength and speed in the bubble), allowing, for a perfect gas, to exactly express the gas radial velocity as:and from this result recover an exact relationship for the gas pressure:where r is the radial coordinate, T_g is the gas temperature, and K_g is the thermal conductivity of the gas. These two relationships allow to reduce the number of governing partial differential equations (PDEs) from three to one, the one necessary to find the temperature gradient at the bubble surface, i.e. ∂_rT_g|_R, which is the term that governs the heat exchange and thus closes the problem. One can solve either the energy equation or the continuity equation coupled with the equation of state for perfect gases. The first alternative was chosen by Kamath & Prosperetti,⁷⁷ who solved the energy equation by using a Galerkin–Chebyshev spectral method. Recently, the second alternative was pursued by Zhou,⁷⁸ who found the temperature profile inside the bubble to be divided into three regions: (1) an internal layer characterised by uniform temperature, (2) a buffer layer, and (3) a layer adjacent to the bubble surface characterised by a linear temperature distribution, which in turn suggests dividing the computational domain into as many cells. In this way, the continuity equation

[small rho, Greek, dot above]_g + ∇·(ρ_gu_g) = 0, (9)

where ρ_g denotes the gas density, can be converted into two ordinary differential equations (ODEs)where m_g,i is the gas mass in cell i and f_j is the mass flux across the cell interface j. ρ^uw_g,j is the density of the neighboring cell of interface j in the upwind direction, u^rel_g,j is the convective velocity and S_j is the interface j surface area. The velocity u_g can be obtained from eqn (7), while the temperature T_g computed through the equation of state for ideal gaswhere is the the specific gas constant. The change in bubble surface temperature is only a small fraction of that of the gas and can, therefore, be neglected, i.e. T_g|_R ≈ T_l.⁷⁴ At this point, the motion of the bubble can be effortlessly solved by integrating in time the system of ODEs consisting of eqn (1), (8) and (10). The reader can refer to the article by Zhou⁷⁸ for full details of the method. The simplicity and accuracy of the model make it suitable for this work.

3. Methods

3.1 Microbubble synthesis

Lipid-coated microbubbles with a C4F10 (perfluorobutane, Fluoromed) gas core are prepared in-house by probe sonication. The lipid coating consists of 90 mol% DPPC (1,2-distearoyl-sn-glycero-3-phosphocholine, NOF EUROPE) and 10 mol% of DPPE-PEG5K (1,2-dipalmitoyl-sn-glycero-3-phosphoethanolamine-N-[methoxy(polyethylene glycol)-5000]), Avanti Polar Lipids). The lipids are first dissolved in chloroform. The solution is placed under vacuum at 35 °C overnight to remove the solvent and produce a dry lipid film. The dry lipid film is rehydrated with 1× PBS (phosphate buffered saline, Boston BioProducts) to yield a total lipid concentration of 2 mg mL⁻¹. The solution is then sonicated (SFX550, Branson; 20 kHz, 550 W) at low power (30%) for five minutes to convert the multilamellar vesicles to unilamellar liposomes. Microbubbles are formed by probe-sonicating the surface of the lipid solution at full power for ten seconds while simultaneously flowing C4F10 gas over it. The microbubble suspension is then cooled down to room temperature and washed using centrifugation. Differential centrifugation is used to size-isolate the microbubbles into a specific size distribution (1–4 μm-radius), based on the protocol developed by Feshitan et al.⁷⁹ Fluorescent-labeled microbubbles are fabricated by adding the lipid dye DiI (1,1′-dioctadecyl-3,3,3′3′-tetramethylindocarbocyanine perchlorate, Sigma-Aldrich) to the lipid solution at a concentration of 5 μg mL⁻¹ before probe sonication.

3.2 Experimental setup

To study the dynamics of microbubbles, we developed a multimodal microscope (“Mscope”, Fig. 1) which allows for wide-field ultra-high-speed imaging, holographic optical trapping and wide-field fluorescence.

Fig. 1 Schematics of the experimental setup. (M) Microscope subsystem, (F) fluorescence subsystem, (S) acousto-optic shutter subsystem, (HOT) holographic optical tweezers subsystem. (AOTF) Acousto-optic tunable filter, (BD) beam dump, (C) camera, (DB) dichroic beamsplitter, (EF) emission filter, (FL1–FL2) flash lamp, (HI) halogen illuminator, (HW) half-wave plate, (I) iris, (L) laser, (L1–L8) lens, (OL) objective lens, (SF) spatial filter, (SLM) spatial light modulator, (TC) test chamber, (TL) tube lens, (US) ultrasound transducer, (WB) water bath, (XF) excitation filter. (0*) zeroth-order laser beam image plane, (1*) first-order laser beam image plane. The zoomed-in inset depicts the pre-test conditions, with a single microbubble being optically trapped >50 μm from the top window of the test chamber.

The main body of the system (Fig. 1(M)) is a custom-built upright microscope realised using modular optomechanics components (Thorlabs, cage system), installed on an optical table (T1220C, Thorlabs) with active isolators (PTS603, Thorlabs) to minimise environmental vibrations. The microscope is equipped with a water-dipping objective lens of focal length (CFI Plan 100XC W, Nikon) and a tube lens (TL400-A, Thorlabs) for a total magnification of 200×. An ultra-high-speed camera (HPV-X2, Shimadzu) allows for recordings at 10 million frames per second over 25.6 μs of continuous visualisation of a 64 μm × 40 μm field of view with a 160 nm pixel resolution. Backlight illumination is provided by a continuous halogen illuminator (OSL2, Thorlabs) for live imaging and two Xenon flash lamps used sequentially (MVS-7010, EG&G) for video recording, all combined into a single optical fiber output and focused on the sample through a custom-built condenser (L1 and L2, , AC127-025-A, Thorlabs). A water bath (T_l ≈ 22 °C) filled with deionised water is accommodated at the base of the microscope. A custom-designed 3D-printed test chamber with optically and acoustically transparent windows on top and bottom (low-density polyethylene, acoustic impedance z = 1.79 × 10⁶ Pa s m⁻¹) is suspended in the water bath and its position controlled by a three-axis motorised microtranslation stage (PT3/M-Z8, Thorlabs). A broadband focused ultrasound transducer (PA1612, Precision Acoustics; 1.5 MHz of center frequency, 85 mm of focal length, ≈4 mm of beam width at −6 dB) is positioned in the water bath at an angle of 30 with respect to the horizontal plane to minimise acoustic reflections and manoeuvred using a manual three-axis microtranslation stage (#12694, #66-511, Edmund Optics). A function generator (LW 420B, Teledyne LeCroy) is used to produce the driving pulse, which is then amplified by a radiofrequency power amplifier (1020L, E&I). A calibrated needle hydrophone (0.2 mm, NH0200, Precision Acoustics) is employed to measure the driving acoustic pressure inside the test chamber and to align the acoustic focal point with the optical field.

The holographic optical tweezers and fluorescent microscopy are realised using a 532 nm optically pumped semiconductor continuous-wave laser (Verdi G10, Coherent). The laser beam enters an acousto-optic tunable filter (AOTF.NC-VIS/TN, AA Opto Electronic) that functions as a microsecond-fast electronic shutter (Fig. 1(S)). The transmitted beam (0th order) is ceased with a beam dump while the diffracted beam (1st order), whose intensity is dependent upon the radiofrequency drive power provided by a radiofrequency driver (MOD.8C.10.b.VIS, AA Opto Electronic), is selectively sent to the optical tweezers or the fluorescence setup by means of a movable mirror.

Three-dimensional optical trapping of single microbubbles has been demonstrated in the pioneering works of Prentice et al.⁸⁰ and Garbin et al.⁸¹ by focusing an optical vortex laser beam obtained by converting the laser transverse mode from Gaussian to Laguerre–Gaussian by means of a diffractive optical element. Compared to the infrared laser employed in those studies, using a green laser reduces the local heating of water by three orders of magnitude.^82,83 In our implementation (Fig. 1(HOT)), the optical path of the optical tweezers setup is kept as short as possible to achieve the highest pointing stability. The diffractive optical element (DOE) is provided by a reflective spatial light modulator (SLM) (PLUTO-2.1 VIS-096, HOLOEYE Photonics). The SLM is placed after a 7.5× galileian telescope, made of a pair of two lenses (L3, , ACN127-020-A, Thorlabs; L4, , AC254-150-A, Thorlabs), and an iris (SM1D12, Thorlabs) which expand and spatially filter the laser beam to completely fill the active area of the SLM with an approximately uniform-intensity distribution which helps to generate a strong optical trap. The linear polarisation of the laser is rotated by using a zero-order half-wave plate (WPHSM05-532, Thorlabs) to match the alignment of the liquid crystalline molecules of the SLM microdisplay. The laser beam mode is converted into a Laguerre–Gaussian one by implementing an optical vortex phase mask on the SLM. In order to minimise the diffraction effects, a system is used to relay the SLM image plane to the objective back focal plane (L5, , AC254-400-A, Thorlabs; L6, , AC254-200-A, Thorlabs). The magnification of the system is 0.5× in order to slightly overfill the objective lens back aperture and thus obtain the best optical trapping performance. The 0th-order component of the beam, unavoidably present due to the pixelated structure of the SLM, is suppressed by performing spatial filtering after having separated the 1st-order-component image plane from the 0th one by superimposing a phase mask for a lens onto the SLM. Upon reflection on a reflective-band dichroic beamsplitter (ZT532dcrb, Chroma), the beam is focused on the sample by the objective lens. The optical trap can be freely positioned within the sample image plane by superimposing a phase mask for a prism onto the SLM.

With the use of movable mirrors the optical tweezers setup can be bypassed and make use of the laser for performing epifluorescence microscopy (Fig. 1(F)). The size of the field of illumination and consequently the laser-excitation power density can be adjusted varying the spacing of a lens relay system (L7 and L8, , AC254-100-A, Thorlabs). Undesirable light wavelengths are prevented from illuminating the specimen by means of a narrow passband excitation filter (ZET532/10x, Chroma). The dichroic beamsplitter reflects the excitation light into the objective lens. The laser line is removed from the specimen image with an emission filter (ET590/50m, Chroma).

4. Results

4.1 Microbubble radius measurement

Measuring accurately the radius of a microbubble is essential to equitably compare its experimental behaviour with theoretical prediction. However, making this measurement with bright-field microscopy is not trivial as the size of a microbubble is comparable to the wavelength of light. The light waves passing around the microbubble generate a Fresnel diffraction pattern that appears as concentric bright and dark circles around the bubble and thus hinders the accurate estimation of its dimension. Fortunately, fluorescence microscopy does not suffer from the interference between the light rays and the specimen because the specimen itself is self-emitting. By adding a fluorophore to the shell composition, it is possible to leverage fluorescence microscopy to image the microbubble. However, the brightness intensity of fluorophores is not sufficient for ultra-high speed recordings. Fluorescence microscopy is therefore only employed to evaluate the systematic error present in the bubble radius measurements obtained with bright-field microscopy. The error is assessed by comparing the bubble radius as measured with fluorescence and bright-field microscopy. Fig. 2 compares the two, showing that bright-field microscopy overestimates the bubble radius by R_err = (140 ± 40) nm across all bubble sizes. It is important to mention that this error value is characteristic to this specific optical setup and follows from the particular combination of the light source and objective lens numerical apertures used. It should not be surprising that a discrepancy smaller than the pixel resolution can be recovered, owing to the sub-pixel accuracy of the edge detection algorithm employed for detecting the bubble contour.

Fig. 2 Comparison between bubble radii measured from a fluorescence (R_fluoro) and bright-field (R_bright) microscopy image (orange dots). The green solid line represents the zero-offset line, while the orange dashed line represents a fit with an offset of 140 nm. Snapshots of the same microbubble taken with bright-field and fluorescence microscopy are shown in the orange and green boxes, respectively.

4.2 Microbubble response upon ultrasonic driving

The lipid-coated microbubbles are injected at a very low dilution ratio in the test chamber in order to avoid inter-bubble radiative forces (secondary Bjerknes forces) upon acoustic excitation. Because of buoyancy, the bubbles float up and adhere to the top window of the test chamber. Individual bubbles are then optically trapped and moved away from the top window by >50 μm using optical tweezers in order to avoid wall interferences. The order of the optical vortex and the laser intensity are adjusted manually for each bubble size. The bubble under investigation is acoustically driven by a 20-cycle sinusoidal pulse at 1.5 MHz and with 40 kPa of peak negative pressure produced by the ultrasound transducer. A video recording at 10 million frames per second of the bubble response is performed using the ultra-high-speed videomicroscopy apparatus. The optical trap is deactivated during the recording by using the microsecond-fast electronic shutter. Fig. 3(a) summarises the temporal chain of events occurring during a single experiment, including the measured pressure driving pulse (black), the measured laser intensity (i.e., the optical trap strength) (green) and the recording time window (orange).

Fig. 3 (a) Sequence and time of events in the experiment. The black line represents the driving pressure pulse produced by the ultrasound transducer and recorded with a hydrophone (1.5 MHz, 40 kPa, 20 cycles). The green line shows the optical trap intensity around the microbubble in arbitrary units measured with a photodetector. The orange-shaded area depicts the video recording time window. (b) Normalised maximum radial expansion for different equilibrium bubble radii R₀ extracted from time-resolved video recordings of single bubble responses to the driving pressure pulse (orange dots). The orange-shaded area represents the envelope of maximum radial expansion simulated numerically with the model detailed in Section 2 using the range of shell modulus and viscosity values given in the legend. The grey and black solid lines depict the numerically computed maximum radial expansion for a low and high shell viscosity, respectively, and a fixed shell modulus (values given in the legend). (c) Radial time evolution of two bubbles having approximately the same size but largely different radial excursions (orange line). Selected snapshots of the dynamics of the two bubbles, extracted from Movies S1 and S2 (ESI†), are shown in the orange boxes. The grey and black solid lines represent the numerically computed bubble radial time evolution for a low and high shell viscosity, respectively, and a fixed shell modulus (values given in the legend). The remaining parameters used in the simulations are: p₀ = 102.2 kPa, σ_w = 72.8 mN m⁻¹, c_l = 1481 m s⁻¹, μ_l = 9.54 × 10⁻⁴, ρ_l = 997.8 kg m⁻³, σ₀ = 0 N m⁻¹, γ = 1.4, K_g = 0.026 W m⁻¹ K⁻¹, , T_g|_R = 295 K.

Radius–time curves are extracted from the recordings using a feature-extraction image-processing algorithm. The normalised maximum bubble radial expansion R_max/R₀ − 1 is taken as the exemplary feature of the bubble response and as the target for the inverse problem of inferring the shell viscosity from the bubble response. This choice is more robust than considering the average maximum expansion over several oscillation cycles because it prevents the bubble dissolution, which may occur in the final cycles to a varying extent, or shape modes, which develop for large bubbles again in the final cycles, from influencing the bubble response evaluation and the consequent shell viscosity estimation. In fact, bubble shrinkage leads to a reduction in radial excursion while shape modes cause the loss of the bubble's spherical symmetry and thus the inability to correctly relate the bubble's cross-sectional area with its volume and hence to an equivalent radius. Nevertheless, these issues always occur after the maximum bubble expansion is reached, which we can easily verify visually from individual videos. By selecting the maximum bubble expansion as the target for the inverse problem of inferring the shell viscosity from the bubble response, we effectively prevent these concerns from compromising the estimation of shell viscosity. Fig. 3(b) depicts with orange dots the distribution of the normalised maximum bubble radial expansion across a range of bubble equilibrium radii R₀. The bubble radius is adjusted by considering the systematic sizing error coming from the Fresnel interference and deemed equal to R_err = 140 nm.

As expected, the curve exhibits a pronounced peak at a specific bubble radius that represents the resonance bubble radius. However, it also unveils a large variability in the response for microbubbles of similar size. A representative example of this arbitrariness is detailed in Fig. 3(c), which shows the large difference in the radial time evolution of two bubbles of nearly identical size (R₀ = 2.5 μm), reporting their radius–time curves in solid orange lines and selected snapshots of their dynamics, extracted from Movies S1 and S2 (ESI†), in the orange boxes. This irregular bubble behaviour can be explained with a dispersion in the shell parameters: dilatational modulus but most prominently dilatational viscosity. The orange-shaded area in Fig. 3(b), computed using the bubble dynamics model reported in Section 2, shows how a range of dilatational modulus E_s = 0.3 − 2 N m⁻¹ and dilatational viscosity κ_s = 1.5 × 10⁻⁹–12 × 10⁻⁹ kg s⁻¹ allows us to fully describe the variability in the experimental results. This figure also suggests a useful qualitative classification of the bubble response into two main regimes, i.e., shell-dominated and inertia-dominated, based on the bubble size.

In the shell-dominated regime, the shell properties have a large influence on the bubble dynamics. Therefore, minor changes in the shell rheology may lead to profoundly different bubble response outcomes. This leverage effect is amplified in the resonance sub-regime where a bubble experiences the largest radial excursion. It should be noted that high radial excursion values can possibly result in the elastic shell regime to play a marginal role in the overall bubble response compared to the prevailing buckling and rupture regimes. This is also the case for the combination of shell composition and negative peak pressure employed in this study. In fact, as shown in Fig. 3(b), the bubble response in the resonance regime is not influenced by the choice of the dilatational modulus. The grey-shaded area, representing the dispersion of the bubble response for the same range in shell viscosity reported before but now for a specific dilatational modulus value E_s = 0.6 N m⁻¹, matches finely the orange-shaded one. Therefore, for clinical applications, where the bubble suspension is subjected to relatively strong ultrasound forcing and the resonance is leveraged to get the maximum result at the minimum mechanical index, characterising the shell viscosity is certainly more important than doing so for the dilatational modulus. Nonetheless, for sub-resonant bubbles, the response is less sensitive to shell viscosity variations and therefore these alone—even with a larger range of values than the one used for Fig. 3(b)—cannot fully explain the scatter found in the experimental results, which must therefore also be attributed to the dilatational modulus.

On the other hand, in the inertia-dominated regime, the shell plays a marginal role and the radial excursion of the bubble is mainly dependent on its equilibrium radius, and therefore on its inertia, as evident from the different curves representing different values of shell properties collapsing on a single one. As a consequence, the elevation of this curve can only be controlled by the driving pressure. Therefore, finding a good match between experiments and theory in this regime means that the measured and theoretical pressure are coinciding, propping our acoustic pressure measurement. Also note how the variability of the bubble response drops in this regime which correlates with the diminished influence of shell properties on the response itself. This observation further reinforces the idea that the observed variability in the bubble response stems from the scattering of shell properties from one bubble to another.

The measured bubble radial time evolutions agree well with the theoretical predictions as shown by the representative radius–time curves in Fig. 3(c). The bubble with the largest excursion is best modelled with a shell viscosity κ_s = 1.8 × 10⁻⁹ kg s⁻¹. The discrepancy in the final cycles and the new equilibrium radius are caused by the net gas loss of the bubble, which is likely facilitated by the partial shedding of the shell. These two phenomena are not accounted for in the theoretical model. Consequently, as the ultrasound burst progresses, the bubble experiences changes in its resonance frequency and shell properties compared to the initial conditions, leading to deviations in its response from the theoretical predictions. Conversely, the bubble with the smallest radial excursion is simulated with a much higher shell viscosity κ_s = 10 × 10⁻⁹ kg s⁻¹. In this case, the bubble does not undergo shrinkage, resulting in excellent agreement between the theoretical and experimental curves, even in the tail region. This result further strengthens our hypothesis regarding shell viscosity being the prime actor in explaining the variability in bubble response. We emphasise that bubble dissolution only occurs during intense radial oscillations. Otherwise, the bubbles remain stable for extended periods, tens of minutes, within the aqueous solution.

The theoretical bubble response is calculated by setting the initial surface tension σ₀ to zero because no finite value can encompass all data points, whereas a null value successfully does so, as illustrated in Fig. 3(b). It is worth noting that a tensionless or nearly tensionless (null Laplace pressure) bubble is also consistent with its long-term stability against dissolution in a saturated medium.²⁵ Although the microbubbles are originally fabricated with a perfluorobutane gas core, it is likely that this has been replaced by air due to the extended period of time (few tens of minutes) that the bubbles reside in an air-saturated aqueous medium. A study by Kwan and Borden,⁸⁴ supports this notion, revealing that the heavy-molecular-weight gas initially present in the gas core of microbubbles is substituted by the air dissolved in the surrounding aqueous medium within a few minutes. Accordingly, we set γ = 1.4, K_g = 0.026 W m⁻¹ K⁻¹ and . The driving acoustic signal used in the simulations is the actual experimentally recorded pressure pulse (Fig. 3(a)). The remaining parameters used are: p₀ = 102.2 kPa, σ_w = 72.8 mN m⁻¹, c_l = 1481 m s⁻¹, μ_l = 9.54 × 10⁻⁴ Pa s, ρ_l = 997.8 kg m⁻³, T_g|_R = 295 K.

4.3 Shell dilatational viscosity estimation

The shell dilatational viscosity of each tested microbubble is inferred from the observation of its radial time evolution by solving an inverse problem. A grid search optimisation algorithm is used to find the theoretical bubble radial time evolution that results in the minimum deviation between empirical and theoretical maximum bubble expansions. The free parameters are the dilatational shell viscosity κ_s and modulus E_s. The other parameters are held fixed at the values presented in the previous section. We would like to emphasise that due to the typical shape of pressure signals generated by ultrasound transducers, inferring shell viscosity from the maximum expansion of the bubble is a more robust approach than inferring it from the decaying rate of the bubble response tail. The response tail, in fact, does not truly represent the free dissipative relaxation of the bubble but rather the forced response to the tail of the pressure driving signal (see Fig. 3(a)). On the contrary, the (maximum) radial expansion of a bubble correlates robustly with its shell viscosity as shown in Fig. 3(b).

Fig. 4 reports with the orange dots the values of shell viscosity resulting from solving the inverse problem, which show a large variability of one order of magnitude and an average value of . In contrast to previous studies (e.g. in red van der Meer et al.,⁴² in green Daeichin et al.⁵¹), no dependence on the bubble size is found.

Fig. 4 Shell dilatational viscosity κ_s for different equilibrium bubble radii R₀. The orange dots represent the values gathered in the present study by solving the inverse problem with an optimisation algorithm. The objective function is to minimise the deviation in the maximum normalised radial expansion between the experimental and synthetic radial time evolutions. The green dots represent the dataset collected by van der Meer et al.⁴² by extracting the shell viscosity from the width of the frequency response of single microbubbles reconstructed by scanning through a range of driving frequencies. The red dots represent the dataset collected in Daeichin et al.⁵¹ by extracting the shell viscosity from the width of the frequency response of single microbubbles directly acquired using impulsive driving. The coloured solid lines represent linear fits to the corresponding datasets.

We feel confident about the reliability of our shell viscosity estimates, i.e. the absence of biases potentially being introduced by the theoretical model used to solve the inverse problem, because we have proven with Fig. 3(c) that the model employed is capable of describing the (nonlinear) dynamics of single microbubbles and, most importantly, with Fig. 3(b), that it is also able to faithfully represent the overall microbubble behaviour across the whole size range with a single set of parameters, and in particular, to match the position of the experimental resonance peak, which is crucial for avoiding systematic errors. From this perspective, one can finally understand the importance of taking into account even the small systematic error R_err in the sizing of the bubbles, which causes a minimal but still significant shift in the resonance peak position.

It should be noted that these values of shell viscosity do not necessarily represent the ground truth, but instead, a single-number portrayal of the complex dissipative processes occurring in the shell through the lens of the Marmottant model. Nothing prevents new rheological developments concerning the effect of dynamical parameters (e.g. the strain of the surface area) on the shell viscosity, from changing this quantitative picture.

We do not provide the result of the inverse problem regarding the dilatational modulus because it would be inconclusive. As explained in the previous section, at these acoustic driving conditions, for a wide range of bubble sizes, the elastic regime of the shell plays a marginal role and therefore the modulus can take on an arbitrary value. To correctly infer the modulus, it is necessary to probe the linear response of the bubble. Unfortunately, resolving such small oscillations exceeds the limits of optical microscopy. One must therefore resort to laser scattering techniques, which, however, compromise on the side of versatility. Notwithstanding, the range of values found in this study, as shown in Fig. 3(b), matches well the measurements carried out with the aforementioned laser scattering techniques (see e.g. Daeichin et al.⁵¹).

5. Discussion

5.1 The shell dilatational viscosity shows large variability and no dependency on the bubble size

We believe that the variability in the bubble response is genuine and not a byproduct of the experimental methodology in light of the following arguments: (i) absence of significant arbitrary interference on the bubble dynamics from the pressure waves that reflect off from the top surface of the test chamber owing to the comparable acoustic impedances of water and the membrane material (low-density polyethylene, z = 1.79 × 10⁶ Pa s m⁻¹) and as confirmed by the pressure pulse recorded inside the test chamber with a needle hydrophone inserted through a side opening and shown in Fig. 3(a); (ii) absence of significant arbitrary interference on the bubble dynamics from pressure waves being scattered by neighbouring bubbles by reason of the extremely low dilution ratio; (iii) absence of arbitrary thermal effects affecting the shell rheological properties caused by the optical trap because of the very low localised temperature increase, estimated in the order of a few tens of millikelvin based on the work by Peterman et al.⁸⁵ and the very low electromagnetic absorption coefficient of water at the wavelength of 532 nm.^82,83

Most probably, the variability in bubble response can be rooted in the arbitrariness of the rheological properties of the shell. In support of this theory, it is noteworthy that the variability of the bubble response changes in accordance with the influence exerted by the shell properties on it, as clearly illustrated in Fig. 3(b). The rheological properties of the shell, in turn, are dependent on its microstructure. At surface tension values close to zero, a phospholipid monolayer self-assembles in a two-dimensional polydomain liquid crystal with different tailgroup orientations. High variation in the size and distribution of such phospholipid domains on the bubble interface has been observed by Borden et al.⁸⁶ by using fluorescence and electron microscopy and by Kooiman et al.⁸⁷ with confocal fluorescence microscopy. The edge of the domains, where phospholipids with different tailgroups interact with each other, can be considered as defects of the bubble shell. A larger amount of defects results in more frequent interactions between phospholipids with different tailgroup orientations which, in turn, could yield a higher shell viscosity as experimentally proven by Hermans & Vermant²⁷ by testing the interfacial shear rheology of DPPC at different preshear values (i.e. at different levels of crystalline microstructure refinement). Based on these findings, we can conclude that the rheological properties of the shell and consequently the acoustic behaviour of the bubble are uniquely determined by the characteristic processing history of the shell. It follows that we cannot rule out the possibility that the repeated driving pressure pulses to which the bubbles in the test chamber are subjected during the test may change the original morphology of the shell (through clustering, partitioning or even shedding off of shell material). Nevertheless, the time interval between each ultrasound pulse in our experiment (a few minutes) is designed to allow for the completion of any dissolution processes and the restoration of full lipid coverage at the bubble interface (i.e. null surface stresses). Therefore, we believe that the effect of shell aging, which inevitably occurs during the operation of microbubbles in any case, does not significantly contribute to the observed variability in shell properties that is already naturally present.

The independence of shell viscosity on the bubble size found in our measurements is consistent with the unchanging shell thickness (consisting, most probably, of only a monolayer of phospholipids) and therefore shell properties with the size of the bubble. However, our shell viscosity estimates contradict previous studies that reported a smaller scatter and an upward trend in shell viscosity as the bubble equilibrium radius increases, albeit with varying rates across studies.^{36–40,42,45–47,50,51} We believe that these divergent results are an artifact caused by the measurement techniques adopted. In the next sections, we present arguments supporting our hypothesis.

5.2 Bubble spectroscopy: overestimating the frequency response bandwidth affects the shell viscosity estimation

Most of studies which found the upward trend with the bubble equilibrium radius have estimated the shell viscosity by analysing the bubble response in the frequency domain,^{36–38,42,45,50,51} adopting the “bubble spectroscopy” technique introduced by van der Meer et al.⁴² or an equivalent of it for acoustic attenuation measurements.^36,37 The bubble spectroscopy method propounds to estimate the shell viscosity from the half-energy bandwidth of the microbubble frequency response treating the microbubble as a linear harmonic oscillator. We now show how the upward trends are caused by the overestimation of bandwidth due to experimental biases. To align ourselves with these previous works and for ease of exposition from here on, we adopt the polytropic process approximation for the gas contained in the bubble (i.e., eqn (6)).

If one assumes that the bubble undergoes small radial oscillations, eqn (1) can be linearised and reduced to:

ẍ + 2δω₀ẋ + ω₀²x = −p_d(t)/ρ_lR₀², (12)

where x(t) = (R(t) − R₀)/R₀ is the normalised radial excursion, ω₀ = 2πf₀ is the natural angular frequency of the system expressed asand δ is the damping coefficient consisting of three contributions: the contribution from sound reradiationthe contribution from the liquid viscosityand the contribution from the shell dilatational viscosityAssuming a harmonic acoustic driving p_d(t) = p_a sin(ωt) and applying the Fourier transform to eqn (12), one obtains the frequency response of the bubble-oscillator:where ω = 2πf is the ultrasound angular frequency. If the damping δ is small enough (δ < 0.1), this can be derived from the frequency response energy aswhere is the full width of the frequency response energy at half peak response (half-energy bandwidth) and f_res = f₀(1 − 2δ²)^1/2 is the peak response frequency or resonance frequency. Finally, the shell viscosity κ_s can be derived from the damping δ asExemplary frequency response energies for microbubbles obeying the linearised model (eqn (17)), with equilibrium radii R₀ = 2–5 μm and shell viscosities κ_s = 1 × 10⁻⁹ kg s⁻¹ and 1 × 10⁻⁸ kg s⁻¹, are given in Fig. 5(a). These two values are chosen because they represent approximately the limit values of shell viscosity found in the present study. It is good to be reminded, however, that this linear approximation only holds for markedly small acoustic driving. Based on the Marmottant shell model, as early as p_a = 5 kPa, the nonlinear frequency response deviates from the linear one, as shown in Fig. 5(b) for a microbubble with radius R₀ = 3.5 μm and shell viscosity κ_s = 5 × 10⁻⁹ kg s⁻¹. Furthermore, this threshold in p_a decreases the closer the microbubble at rest is to the buckled state (i.e., in Marmottant model terms, the closer σ₀ is to zero). This suggests that the statement made in the study of van der Meer et al.⁴² regarding the absence of substantial differences in the frequency response of a microbubble driven at p_a = 1 kPa and one at p_a = 40 kPa is inaccurate. This oversight was caused by the use of a linearised model to simulate the frequency response of the microbubble driven at p_a = 40 kPa.

Fig. 5 (a) Frequency response energies computed with the linearised model (eqn (17)) for microbubbles with equilibrium radii R₀ = 2–5 μm and shell viscosities κ_s = 1 × 10⁻⁹ kg s⁻¹ (grey line) and 1 × 10⁻⁸ kg s⁻¹ (black line). (b) Effect of the acoustic driving pressure amplitude p_a on the frequency response energy of a bubble with equilibrium radius R₀ = 3.5 μm and shell viscosity κ_s = 5 × 10⁻⁹ kg s⁻¹. The grey dashed line is the result of the linearised model (eqn (17)). The pink solid lines are the results of the modified Rayleigh–Plesset equation (eqn (1)) with the Marmottant shell model (eqn (5)) and the polytropic process approximation for the gas contained in the bubble (eqn (6)) for increasing acoustic driving pressure amplitude p_a = 1 kPa, 5 kPa and 40 kPa. The remaining physical parameters used, so as to align with earlier works such as van der Meer et al.,⁴² are: p₀ = 100 kPa, σ_w = 72 mN m⁻¹, c_l = 1500 m s⁻¹, μ_l = 0.002 Pa s, ρ_l = 1000 kg m⁻³, E_s = 0.54 N m⁻¹, σ₀ = 0.02 N m⁻¹, n = 1.07.

In an experimental investigation, the frequency response of a microbubble can be reconstructed by scanning through a range of driving frequencies^{38,42,44,45,49,50} or directly acquired using impulsive driving.⁵¹ If the response thus obtained follows the linearised model described by eqn (17) and depicted in Fig. 6(a) in the solid line, the physical parameters are known exactly and the response is free of experimental biases, then the bubble spectroscopy method, through the use of eqn (18) and (19), allows the shell viscosity to be inferred correctly, as shown for κ_s = 1 × 10⁻⁹ kg s⁻¹ and 1 × 10⁻⁸ kg s⁻¹ in Fig. 6(b) in solid grey and black lines, respectively. However, if the experimental parameters used to survey the microbubble frequency response are not chosen carefully, the spectroscopy method can lead to an erroneous evaluation of the shell viscosity. Large frequency steps when scanning through frequencies, short observation times of the bubble dynamics leading to spectral leakage and high driving pressure amplitudes that would negate the linear oscillator hypothesis can cause an overestimation of the half-energy bandwidth and, as a consequence, of the shell viscosity. If we assume that these inaccuracies predominantly affect the half-energy bandwidth by a systematic deviation Δf_err, the estimated shell viscosities will be impacted by an amount κ_s,err = (π/2)ρ_lR₀³Δf_err(1 − 2δ²)⁻¹. For small damping values, which is usually the case for phospholipid-coated microbubbles, this expression can be approximated as:

It is clear that the cubic dependency on the equilibrium bubble radius implies that even small systematic errors on the half-energy bandwidth yield significant variations in the shell viscosity across a range of bubble sizes. The impact on the shell viscosity estimate of an experimental systematic error Δf_err = 100 kHz in the half-energy bandwidth of frequency responses, as illustrated in Fig. 6(a) in the dashed line, is given in Fig. 6(b) as dashed lines. This result can explain the unphysical increase of shell viscosity with the bubble size and the lower variability compared to the present work that many studies have observed.

Fig. 6 Effectiveness of the bubble spectroscopy method in inferring the shell viscosity of microbubbles with frequency responses affected by errors. (a) The solid line represents the error-free frequency response energy of a microbubble that obeys the linear model (eqn (17)). The dashed line represents the frequency response energy of the same bubble but affected by a bandwidth overestimation Δf_err = 100 kHz. The dotted-dashed line represents the error-free frequency response magnitude of the same bubble. All the lines pertain to a bubble with equilibrium radius R₀ = 3.5 μm and shell viscosity κ_s = 1 × 10⁻⁸ kg s⁻¹. (b) The thin solid lines represent the shell viscosity values to be inferred, in black κ_s = 1 × 10⁻⁸ kg s⁻¹ and in grey κ_s = 1 × 10⁻⁹ kg s⁻¹. The thick solid lines are the result of applying the bubble spectroscopy method (eqn (18) and (19)) to error-free . The dashed lines depict the result of bubble spectroscopy applied to affected by Δf_err = 100 kHz. The dotted-dashed lines depict the result of bubble spectroscopy applied to instead of . The physical parameters used are: p₀ = 100 kPa, σ_w = 72 mN m⁻¹, c_l = 1500 m s⁻¹, μ_l = 0.002 Pa s, ρ_l = 1000 kg m⁻³, E_s = 0.54 N m⁻¹, σ₀ = 0.02 N m⁻¹, n = 1.07.

We would also like to clarify, in view of some confusion in previous studies, that to correctly estimate the damping δ from the frequency response full width at half maximum (FWHW) the frequency response energy has to be used. To estimate the damping from the frequency response magnitude , showcased in Fig. 6(a) in the dotted-dashed line, the full width at -maximum must be employed instead. The consequences of estimating the damping from the half-magnitude bandwidth instead of the half-energy bandwidth can be evaluated computing from eqn (17) the frequency at half magnitude , which for small values of damping can be approximated as , as opposed to the frequency at half energy f_HE ≈ f_res(1 ± δ). This translates into an error in the bandwidth and a departure from the true value of the shell viscosity equal to:

The effect of this error on the shell viscosity estimate is given in Fig. 6(b) as dotted-dashed lines. Again, the error causes an upward shift and an upward trend in the shell viscosity estimate as the equilibrium radius of the bubble increases. We also emphasise that an inappropriate choice of the bubble dynamics model or inaccuracies in the physical parameter measurements can also lead to unphysical trends in the shell viscosity.

We now examine two previous studies, in which the bubble spectroscopy method was employed to infer the shell viscosity, in order to verify whether the values found are correct or are, instead, possibly a byproduct of the inaccuracies listed in the previous section. The two studies chosen are those of van der Meer et al.⁴² and Daeichin et al.⁵¹ because they stand as exemplary models for a large class of other studies which employed the bubble spectroscopy method. The values of shell viscosity they found are shown in Fig. 4, along with the values from the present study. The work of van der Meer et al. pioneered the use of bubble spectroscopy and represents the linchpin for other closely related studies (e.g. Luan et al.,⁴⁴ Helfield et al.³⁸ and van Rooij et al.⁴⁵). The authors reconstructed the frequency response of single microbubbles point-by-point by scanning the frequency of the driving acoustic pulse produced by a broadband ultrasound transducer. At each input frequency, the bubble dynamics was recorded by using an ultra-high-speed videomicroscopy facility, the radial time evolution extracted from it and its Fourier transform computed to finally obtain the response at the input frequency. The work of Daeichin et al. represents the last and most developed iteration of the laser-scattering measurement technique, pioneered by Dove et al.⁸⁸ and already employed by Lum et al.⁴⁹ and Supponen et al.⁵⁰ The authors here have developed a brilliant technique to characterise in one fell swoop the whole microbubble frequency response. The microbubble was probed with a broadband photoacoustic wave generated by a nanosecond-pulsed laser illuminating an optical absorber and the resulting oscillations were detected by light scattering of a continuous laser.

As a first step, we study the scaling laws of the possible deviations of the shell viscosity values reported in the two studies from an assumed true value in order to identify the nature of the potential errors made. In view of our results reported in Section 4, we assume a value κ_s,true = 5 × 10⁻⁹ kg s⁻¹, constant across bubble sizes, as the true shell viscosity value. Fig. 7(a) reveals that the values found by van der Meer et al. deviate from κ_s,true following a trend (dotted green line) that varies proportionally to R₀³ (solid green line). This tendency, in view of eqn (20), advocates that the reported values for shell viscosity are predominantly affected by an approximately constant systematic overestimation of the . The elevation of the fitting line gives an estimate of such overestimation, which yields Δf_err ≈ 270 kHz. The left part of Fig. 7(b) compares the shell viscosity values reported by van der Meer et al. (green dots) with the shell viscosity estimates (dotted lines) inferred from synthetic microbubble frequency responses (eqn (17)) affected by a bandwidth error Δf_err = 270 kHz. The black lines correspond to a shell viscosity value κ_s = 1 × 10⁻⁸ kg s⁻¹, while the grey lines correspond to a value κ_s = 1 × 10⁻⁹ kg s⁻¹. The excellent agreement found suggests that a systematic overestimation of the , approximately constant with bubble size, may be the dominant cause for the unphysical trends in the shell viscosity found in several studies. Fig. 7(a) also reveals that the deviation of the shell viscosities found by Daeichin et al. follows a quadratic trend in logarithmic scale (dotted red line) which results from the sum of a first contribution that scales with R₀³ω₀δ and a second one that scales with R₀³ (sum in solid red line, single contributions in grey lines), unveiling, in view of eqn (21) and (20), that the FWHM is measured from the magnitude instead of the energy of the frequency responses and the bandwidth itself is possibly overestimated, on average, by an amount Δf_err ≈ 190 kHz. The right part of Fig. 7(b) compares the shell viscosity values reported by Daeichin et al. (red dots) with the shell viscosity estimates (grey lines) inferred from the magnitude of synthetic microbubble frequency responses (eqn (17)) affected by a bandwidth error Δf_err = 190 kHz, yielding again a remarkably good agreement.

Fig. 7 (a) Scaling laws of the deviation of the shell viscosity values reported in two previous studies from the assumed true value κ_s,true = 5 × 10⁻⁹ kg s⁻¹. The experimental values are given with the dots (van der Meer et al. in green and Daeichin et al. in red). The dotted lines represent a fit of the experimental values (linear in green, quadratic in red). The green solid line scales with R₀³, representing the deviance caused by an error in the (eqn (20)). The red solid line results from the sum of a first contribution that scales with R₀³ω₀δ (dotted-dashed grey line) and a second one that scales with R₀³ (dashed grey line). The latter one represents the deviance caused by measuring the FWHM from the frequency response magnitude instead of energy (eqn (21)). (b) (left) Comparison between the shell viscosity values reported in van der Meer et al. (green dots) and the values inferred from linear frequency response energies affected by an error Δf_err = 270 kHz in the (dashed lines). The black lines correspond to a shell viscosity κ_s = 1 × 10⁻⁸ kg s⁻¹, while the grey lines to a shell viscosity κ_s = 1 × 10⁻⁹ kg s⁻¹. (right) Comparison between the shell viscosity values reported in Daeichin et al. (red dots) and the values inferred from linear frequency response magnitudes affected by an error Δf_err = 190 kHz in the (dashed lines). The physical parameters used are: p₀ = 100 kPa, σ_w = 72 mN m⁻¹, c_l = 1500 m s⁻¹, μ_l = 0.002 Pa s, ρ_l = 1000 kg m⁻³, E_s = 0.54 N m⁻¹, σ₀ = 0.02 N m⁻¹, n = 1.07.

The comprehensive specifications reported in the work of van der Meer et al. regarding their shell characterisation procedure allow us to directly examine its accuracy by replicating and applying it on synthetic responses of microbubbles of known shell viscosity. The values of the physical parameters and the specifics of the acoustic driving and of the spectroscopy analysis are taken directly from the information reported in their study. The model used for the bubble dynamics is nonlinear and based on the modified Rayleigh–Plesset equation (eqn (1)) with the Marmottant shell model (eqn (5)) and the polytropic process approximation for the gas contained in the bubble (eqn (6)). The latter allows a smooth carryover of the physical parameters used in the study in the bubble model. Note that employing the full model presented in Section 2 would still lead to the same results that we are going to present. The shell viscosity values chosen are again κ_s = 1 × 10⁻⁹ kg s⁻¹ and κ_s = 1 × 10⁻⁸ kg s⁻¹.

The results of this endeavour are presented in Fig. 8(a) which reveals with the dotted-dashed lines that their methodology leads to an inaccurate estimate of the shell viscosity: its value gets inflated, the scatter for a given bubble size is reduced and a strong dependency on the initial radius appears. The estimated shell viscosities compare remarkably well with the values reported by the authors (green dots), endorsing the idea that the apparent dependency of the shell viscosity on the bubble equilibrium radius is indeed only the result of an experimental artefact. The methodological sources of bias occurred in the experimental procedure are diverse: (i) E₁: the linear model was applied to bubbles behaving nonlinearly (p_a = 40 kPa); (ii) E₂: the frequency response was sampled with too coarse a resolution (Δf_sam = 50 kHz); (iii) E₃: the bubble dynamics was observed for too short a period of time (T_obs = 4.2 μs) leading to spectral leakage in the frequency spectrum of its radial time evolution; (iv) E₄: the resonance peak area of the Fourier transform of the radial time evolution of the microbubble under investigation was inappropriately considered as the measure of its response at the specified input frequency. The aftereffect of each source of bias listed above is presented in Fig. 8(b) for a bubble with shell viscosity κ_s = 1 × 10⁻⁹ kg s⁻¹ and two different equilibrium radii R₀ = 2 μm and R₀ = 5 μm. Note that the arithmetic sum of these values does not exactly return the results shown in Fig. 8(a) because the errors combine nonlinearly. To further confirm that the cause of the departure from the prescribed shell viscosity lies in the experimental parameters adopted, the solid thick lines in Fig. 8(a) show that reducing the methodological inaccuracies, detailed above, to negligible values (i.e. using p_a ≈ 1 kPa, Δf_sam < 5 kHz, T_obs > 100 μs) the spectroscopy method would be capable of recovering the correct shell viscosity. Note, however, that the small oscillations resulting from such a low driving pressure would be beyond the resolving capabilities of any existing optical microscopy apparatus.

Fig. 8 (a) Comparison between the shell viscosity values reported by van der Meer et al.,⁴² in green dots, and the values estimated applying the spectroscopy methodology detailed in their work to synthetic nonlinear frequency responses of microbubbles of known shell viscosities κ_s = 1 × 10⁻⁹ kg s⁻¹ and κ_s = 1 × 10⁻⁸ kg s⁻¹, in dashed black and grey lines, respectively. The specifications used, taken from their study, are: p_a = 40 kPa, Δf_sam = 50 kHz, T_obs = 4.2 μs. The solid lines are the result of the same procedure but having reduced the methodological sources of bias using the following specifications: p_a = 1 kPa, Δf_sam = 5 kHz, T_obs = 100 μs. (b) Error breakdown in the estimation of a shell viscosity κ_s = 1 × 10⁻⁹ kg s⁻¹ for two different bubble sizes, R₀ = 2 μm and R₀ = 5 μm on the basis of the methodological inaccuracies E₁–E₄ present in the work of van der Meer et al. Please refer to the text for their description. The model of the bubble dynamics used for calculating the frequency responses is based on the modified Rayleigh–Plesset equation (eqn (1)) with the Marmottant shell model (eqn (5)) and the polytropic process approximation for the gas contained in the bubble (eqn (6)). The physical parameters used are: p₀ = 100 kPa, σ_w = 72 mN m⁻¹, c_l = 1500 m s⁻¹, μ_l = 0.002 Pa s, ρ_l = 1000 kg m⁻³, E_s = 0.54 N m⁻¹, σ₀ = 0.02 N m⁻¹, n = 1.07.

In this regard, we acknowledge the great potential of the shell investigation technique introduced in the work of Daeichin et al. The photoacoustic impulsive driving and the detection method based on laser scattering allow one to surmount the limitations due to the use of optical microscopy that the work of van der Meer et al. ran up against and, in principle, to obtain very accurate shell viscosity and modulus estimates from linearly oscillating bubbles. However, we point out that, in the work in question, in addition to having improperly extracted the damping from the half-amplitude instead of half-energy bandwidth, the observation time used is potentially too short and the driving pressure signal might be deviating from an exact Dirac impulse. These are all factors that can compromise the shell viscosity estimation and possibly explain the reported upward trend with the equilibrium radius.

5.3 Radius–time curve fitting: model and bubble sizing bias affect the shell viscosity estimation

Finally, we turn our attention to the works conducted by Morgan et al.,³⁹ Doinikov et al.⁴⁰ and Tu et al.^46,47 which estimated the shell viscosity of single microbubbles using an approach very similar to ours: by finding the theoretical radius–time curve that results in the minimum mean square error from the empirical curve, recorded with a streak camera in the first two cases or with a side laser scattering technique in the second two. These works also found the shell viscosity increasing with the bubble size. Clearly, the error analysis we carried out for bubble spectroscopy cannot be applied to these studies. However, we believe that their results can also possibly suffer from other methodological oversights causing the upward trend.

Morgan et al.³⁹ and Doinikov et al.⁴⁰ did not report whether the theoretical model used is able to capture the behaviour of the microbubbles across a range of sizes, as we did and shown in Fig. 3(b) and, in particular, whether it can identify the correct resonance radius for the driving frequency employed. Specifically, the microbubbles were driven at a high pressure amplitude (over 100 kPa) and thus into strong nonlinear oscillations, while the theoretical model employed is based on a linear representation of the coating physics. Therefore, the resonance peak from the theoretical model is expected to be at a larger radius than that found in the experiments. This leads the optimisation algorithm to assign a low viscosity to small-diameter bubbles and vice versa, and thus to an upward trend in the shell viscosity similar to that found in their study. Tu et al.^46,47 did employ a nonlinear model, however, they solely conducted indirect measurements of bubble evolution, by converting scattering light signals. This approach raises concerns regarding the reliability of the obtained bubble radius measurements. It is important to note that any inaccuracies in the radius measurement can introduce distortions in the profile of the bubble behaviour across bubble sizes. This, in turn, can result in a mismatch between the experimental data and the theoretical model employed for estimating the shell viscosity, significantly impacting the accuracy and reliability of the inferred shell viscosity. In this regard, Tu et al. also did not compare the theoretical model used with the experimental bubble response across bubble sizes. This omission prevents the resolution of lingering doubts regarding the accuracy of the shell viscosity measurement.

To gain a better understanding of the sensitivity of the shell viscosity estimate on the radius measurement error, let us consider synthetic data of the response of bubbles across different sizes (represented by dashed lines in Fig. 9(a)), generated using a model of bubble behaviour (illustrated by solid lines in Fig. 9(a)), but with a deliberate error R_err = −3% introduced in the radius measurement. If the aforementioned model is subsequently employed to infer the shell viscosity, the error induced by the inaccurate bubble sizing will propagate to the shell viscosity estimate, leading to clear trends correlating with the bubble radius, as shown in Fig. 9(b). This arises from the discrepancy between the positions of the resonance peaks in the model and the synthetic data. It is evident that if the error on the radius has an opposite sign, the trends will exhibit a reversal in their slope. To avoid these undesirable effects in our work, we carefully verified the agreement between the model and the experimental data, and thus the absence of errors on either side, prior to estimating shell viscosity, as portrayed in Fig. 3(b).

Fig. 9 Effectiveness of the radius–time curve fitting in inferring the shell viscosity of microbubbles in the presence of a radius measurement error. (a) The solid lines represent the model of the bubble dynamics (modified Rayleigh–Plesset equation (eqn (1)) with the Marmottant shell model (eqn (5)) and the polytropic process approximation for the gas contained in the bubble (eqn (6))) employed to infer the shell viscosity from radius–time curves of microbubbles. The model is presented for equilibrium radii R₀ = 2–5 μm and shell viscosity values κ_s = 1 × 10⁻⁹ kg s⁻¹ (grey lines) and 1 × 10⁻⁸ kg s⁻¹ (black lines). The dashed lines depict the synthetic dataset generated using the aforementioned model of bubble dynamics but with a deliberate error R_err = −3% in the radius measurement. (b) The solid thin lines represent the shell viscosity to be inferred. The dashed lines are the result of the radius–time curve fitting in inferring the shell viscosity of the aforementioned synthetic dataset. The physical parameters used are: f = 1.5 MHz, p_a = 40 kPa, p₀ = 100 kPa, σ_w = 72 mN m⁻¹, c_l = 1500 m s⁻¹, μ_l = 0.002 Pa s, ρ_l = 1000 kg m⁻³, E_s = 0.54 N m⁻¹, σ₀ = 0.02 N m⁻¹, n = 1.07.

6. Conclusions

In this work, we conducted measurements of the dilatational viscosity of the shell of ultrasound contrast agent microbubbles. By combining ultra-high-speed microscopy imaging, optical trapping and wide-field fluorescence, we obtained recordings of the response of individual microbubbles to an ultrasound driving pulse with unprecedented accuracy. We validated an advanced theoretical model of bubble dynamics based on the compressible Rayleigh–Plesset equation, and augmented by the Marmottant model to account for the shell and the Zhou model to accurately track the pressure variation in the bubble, against the experimental data collected, and found excellent agreement, both with the individual bubble responses and the overall response curve across the entire range of bubble sizes. This model was therefore employed in an optimisation routine to solve the inverse problem of determining the shell viscosity of single bubbles from their radial time evolution. The resulting shell viscosity values show a large, one order of magnitude variability, ranging approximately from κ_s = 1 × 10⁻⁹ kg s⁻¹ to κ_s = 1 × 10⁻⁸ kg s⁻¹, and an average value of . Most interestingly, we found no dependency of shell viscosity on the bubble equilibrium radius. These results contrast with more than two decades of previous studies, which have reported less variability and, very intriguingly, a strong upward trend as the bubble equilibrium radius increases. The latter behaviour is considered to be unphysical, nevertheless, it has thus far never been either explained or proven as such. We believe that the root cause of this issue in the majority of studies lies in the high sensitivity to errors of “bubble spectroscopy”, the method employed to extract shell viscosity from the bubble frequency response bandwidth. We demonstrated with analytical arguments that an overestimation of the bandwidth leads to inflated shell viscosity values, reduced dispersion for a given bubble size, and the appearance of a strong dependence on the equilibrium radius. We pinpointed several methodological inaccuracies in previous studies that lead to an overestimation of bandwidth, including: the use of linear models to describe the nonlinear bubble response, a too short observation time of the bubble response, a too low frequency sampling rate, the use of the half-magnitude bandwidth instead of the half-energy bandwidth, to name a few. We exemplified that resolving these inaccuracies allows the same trend found in our measurements to be recovered. Furthermore, our findings indicate that the unphysical trends observed in the remaining studies, where the shell viscosity was determined by fitting radius–time curves, may be attributed to errors in bubble sizing. This issue leads to a misalignment of resonance peaks between the experimental data and the bubble dynamics model used to infer the shell viscosity. As a consequence, the optimisation procedure assigns incorrect shell viscosity values to the microbubbles, resulting in trends that falsely correlate with the bubble radius. Finally, we suggested that the variability found in the shell viscosity values reflects the arbitrariness in the size and distribution of the phospholipid domains among the bubble shells. The results presented could provide useful guidelines for characterising the shell properties of microbubbles which is critical for predicting their bioeffects in medical applications.

Author contributions

Marco Cattaneo: conceptualisation, methodology, software, formal analysis, investigation, data curation, writing – original draft, visualisation, project administration. Outi Supponen: conceptualisation, methodology, investigation, writing – review & editing, project administration, resources, supervision.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

We thank Dr Gazendra Shakya for the preparation of microbubbles. We also thank Prof. Jan Vermant for insightful discussions about phospholipid monolayer rheology. We acknowledge ETH Zürich for the financial support.

Notes and references

1. K. Wei, A. Jayaweera, S. Firoozan, A. Linka, D. Skyba and S. Kaul, Circulation, 1998, 97, 473–483.
2. R. Vogel, A. Indermühle, J. Reinhardt, P. Meier, P. Siegrist, M. Namdar, P. Kaufmann and C. Seiler, J. Am. Coll. Cardiol., 2005, 45, 754–762.
3. S. Wilson, P. Burns, D. Muradali, J. Wilson and X. Lai, Radiology, 2000, 215, 153–161.
4. S.-J. Rim, H. Leong-Poi, J. Lindner, D. Couture, D. Ellegala, H. Mason, M. Durieux, N. Kassel and S. Kaul, Circulation, 2001, 104, 2582–2587.
5. K. Wei, E. Le, J.-P. Bin, M. Coggins, J. Thorpe and S. Kaul, J. Am. Coll. Cardiol., 2001, 37, 1135–1140.
6. J. Lindner, J. Song, J. Christiansen, A. Klibanov, F. Xu and K. Ley, Circulation, 2001, 104, 2107–2112.
7. B. Kaufmann, J. Sanders, C. Davis, A. Xie, P. Aldred, I. Sarembock and J. Lindner, Circulation, 2007, 116, 276–284.
8. H. Leong-Poi, J. Christiansen, A. Klibanov, S. Kaul and J. Lindner, Circulation, 2003, 107, 455–460.
9. D. Ellegala, H. Leong-Poi, J. Carpenter, A. Klibanov, S. Kaul, M. Shaffrey, J. Sklenar and J. Lindner, Circulation, 2003, 108, 336–341.
10. M. Smeenge, F. Tranquart, C. Mannaerts, T. De Reijke, M. Van De Vijver, M. Laguna, S. Pochon, J. De La Rosette and H. Wijkstra, Invest. Radiol., 2017, 52, 419–427.
11. J. Willmann, L. Bonomo, A. Testa, P. Rinaldi, G. Rindi, K. Valluru, G. Petrone, M. Martini, A. Lutz and S. Gambhir, J. Clin. Oncol., 2017, 35, 2133–2140.
12. P. Marmottant and S. Hilgenfeldt, Nature, 2003, 423, 153–156.
13. P. Prentice, A. Cuschieri, K. Dholakia, M. Prausnitz and P. Campbell, Nat. Phys., 2005, 1, 107–110.
14. A. Alexandrov, A. Demchuk, W. Burgin, D. Robinson and J. Grotta, J. Neuroimaging, 2004, 14, 113–117.
15. C. Molina, M. Ribo, M. Rubiera, J. Montaner, E. Santamarina, R. Delgado-Mederos, J. Arenillas, R. Huertas, F. Purroy, P. Delgado, P. Delgado and J. Alvarez-Sabín, Stroke, 2006, 37, 425–429.
16. L. Qian, B. Thapa, J. Hong, Y. Zhang, M. Zhu, M. Chu, J. Yao and D. Xu, J. Thorac. Dis., 2018, 10, 1099–1111.
17. H. Yuan, H. Hu, J. Sun, M. Shi, H. Yu, C. Li, Y. Sun, Z. Yang and R. Hoffman, In Vivo, 2018, 32, 1025–1032.
18. G. Dimcevski, S. Kotopoulis, T. Bjånes, D. Hoem, J. Schjøt, B. Gjertsen, M. Biermann, A. Molven, H. Sorbye, E. McCormack, M. Postema and O. Gilja, J. Controlled Release, 2016, 243, 172–181.
19. Y. Wang, Y. Li, K. Yan, L. Shen, W. Yang, J. Gong and K. Ding, Chin. J. Cancer Res., 2018, 30, 553–563.
20. A. Carpentier, M. Canney, A. Vignot, V. Reina, K. Beccaria, C. Horodyckid, C. Karachi, D. Leclercq, C. Lafon, J.-Y. Chapelon, L. Capelle, P. Cornu, M. Sanson, K. Hoang-Xuan, J.-Y. Delattre and A. Idbaih, Sci. Translational Med., 2016, 8, 343re2.
21. A. Idbaih, M. Canney, L. Belin, C. Desseaux, A. Vignot, G. Bouchoux, N. Asquier, B. Law-Ye, D. Leclercq, A. Bissery, J.-Y. Delattre and A. Carpentier, Clin. Cancer Res., 2019, 25, 3793–3801.
22. T. Mainprize, N. Lipsman, Y. Huang, Y. Meng, A. Bethune, S. Ironside, C. Heyn, R. Alkins, M. Trudeau, A. Sahgal, J. Perry and K. Hynynen, Sci. Rep., 2019, 9, 321.
23. N. Lipsman, Y. Meng, A. Bethune, Y. Huang, B. Lam, M. Masellis, N. Herrmann, C. Heyn, I. Aubert, A. Boutet, K. Hynynen and S. Black, Nat. Commun., 2018, 9, 2336.
24. A. Abrahao, Y. Meng, M. Llinas, Y. Huang, C. Hamani, T. Mainprize, I. Aubert, C. Heyn, S. Black, K. Hynynen, N. Lipsman and L. Zinman, Nat. Commun., 2019, 10, 4373.
25. K. Ferrara, R. Pollard and M. Borden, Annu. Rev. Biomed. Eng., 2007, 9, 415–447.
26. D. A. Edwards, H. Brenner and D. T. Wasan, Interfacial Transport Processes and Rheology, Elsevier, 1991.
27. E. Hermans and J. Vermant, Soft Matter, 2014, 10, 175–186.
28. H. Manikantan and T. Squires, J. Fluid Mech., 2020, 892, P1.
29. B. Dollet, P. Marmottant and V. Garbin, Annu. Rev. Fluid Mech., 2019, 51, 331–355.
30. S. Derkach, J. Krägel and R. Miller, Colloid J., 2009, 71, 1–17.
31. N. Alvarez, L. Walker and S. Anna, Langmuir, 2010, 26, 13310–13319.
32. N. de Jong, L. Hoff, T. Skotland and N. Bom, Ultrasonics, 1992, 30, 95–103.
33. P. Frinking and N. De Jong, Ultrasound Med. Biol., 1998, 24, 523–533.
34. L. Hoff, P. Sontum and J. Hovem, J. Acoust. Soc. Am., 2000, 107, 2272–2280.
35. J.-M. Gorce, M. Arditi and M. Schneider, Invest. Radiol., 2000, 35, 661–671.
36. M. Parrales, J. Fernandez, M. Perez-Saborid, J. Kopechek and T. Porter, J. Acoust. Soc. Am., 2014, 136, 1077–1084.
37. T. Segers, E. Gaud, M. Versluis and P. Frinking, Soft Matter, 2018, 14, 9550–9561.
38. B. Helfield and D. Goertz, J. Acoust. Soc. Am., 2013, 133, 1158–1168.
39. K. Morgan, J. Allen, P. Dayton, J. Chomas, A. Klibanov and K. Ferrara, IEEE Trans. Ultrason. Eng., 2000, 47, 1494–1509.
40. A. Doinikov, J. Haac and P. Dayton, Ultrasonics, 2009, 49, 269–275.
41. C. Chin, C. Lancée, J. Borsboom, F. Mastik, M. Frijlink, N. De Jong, M. Versluis and D. Lohse, Rev. Sci. Instrum., 2003, 74, 5026–5034.
42. S. Van Der Meer, B. Dollet, M. Voormolen, C. Chin, A. Bouakaz, N. De Jong, M. Versluis and D. Lohse, J. Acoust. Soc. Am., 2007, 121, 648–656.
43. M. Overvelde, V. Garbin, J. Sijl, B. Dollet, N. de Jong, D. Lohse and M. Versluis, Ultrasound Med. Biol., 2010, 36, 2080–2092.
44. Y. Luan, T. Faez, E. Gelderblom, I. Skachkov, B. Geers, I. Lentacker, T. van der Steen, M. Versluis and N. de Jong, Ultrasound Med. Biol., 2012, 38, 2174–2185.
45. T. van Rooij, Y. Luan, G. Renaud, A. van der Steen, M. Versluis, N. de Jong and K. Kooiman, Ultrasound Med. Biol., 2015, 41, 1432–1445.
46. J. Tu, J. Guan, Y. Qiu and T. Matula, J. Acoust. Soc. Am., 2009, 126, 2954–2962.
47. J. Tu, J. Swalwell, D. Giraud, W. Cui, W. Chen and T. Matula, IEEE Trans. Ultrason. Eng., 2011, 58, 955–963.
48. J. Lum, J. Dove, T. Murray and M. Borden, Langmuir, 2016, 32, 9410–9417.
49. J. Lum, D. Stobbe, M. Borden and T. Murray, Appl. Phys. Lett., 2018, 112, 111905.
50. O. Supponen, A. Upadhyay, J. Lum, F. Guidi, T. Murray, H. Vos, P. Tortoli and M. Borden, J. Acoust. Soc. Am., 2020, 147, 3236.
51. V. Daeichin, M. Inzunza-Ibarra, J. Lum, M. Borden and T. Murray, IEEE Trans. Ultrason. Eng., 2021, 68, 2311–2314.
52. A. Prosperetti and A. Lezzi, J. Fluid Mech., 1986, 168, 457–478.
53. A. Lezzi and A. Prosperetti, J. Fluid Mech., 1987, 185, 289–321.
54. M. Brenner, S. Hilgenfeldt and D. Lohse, Rev. Mod. Phys., 2002, 74, 425–484.
55. A. Prosperetti, J. Acoust. Soc. Am., 2013, 133, 3719–3726.
56. C. Church, J. Acoust. Soc. Am., 1995, 97, 1510–1521.
57. D. Chatterjee and K. Sarkar, Ultrasound Med. Biol., 2003, 29, 1749–1757.
58. K. Sarkar, W. Shi, D. Chatterjee and F. Forsberg, J. Acoust. Soc. Am., 2005, 118, 539–550.
59. M. Borden, H. Zhang, R. Gillies, P. Dayton and K. Ferrara, Biomaterials, 2008, 29, 597–606.
60. C. Christiansen, H. Kryvi, P. Sontum and T. Skotland, Biotechnol. Appl. Biochem., 1994, 19, 307–320.
61. C. Pozrikidis, Modeling and Simulation of Capsules and Biological Cells, Chapman and Hall/CRC, 2003.
62. E. A. Evans and R. Skalak, Mechanics and Thermodynamics of Biomembranes, CRC Press, 1980.
63. N. Jaensson and J. Vermant, Curr. Opin. Colloid Interface Sci., 2018, 37, 136–150.
64. N. Jaensson, P. Anderson and J. Vermant, J. Non-Newtonian Fluid Mech., 2021, 290, 104507.
65. O. Lotsberg, J. Hovem and B. Aksum, J. Acoust. Soc. Am., 1996, 99, 1366–1369.
66. P. Shankar, P. Krishna and V. Newhouse, Ultrasound Med. Biol., 1998, 24, 395–399.
67. P. Shankar, P. Krishna and V. Newhouse, J. Acoust. Soc. Am., 1999, 106, 2104–2110.
68. Y. Sun, D. Kruse, P. Dayton and K. Ferrara, IEEE Trans. Ultrason. Eng., 2005, 52, 1981–1991.
69. J. Sijl, B. Dollet, M. Overvelde, V. Garbin, T. Rozendal, N. De Jong, D. Lohse and M. Versluis, J. Acoust. Soc. Am., 2010, 128, 3239–3252.
70. P. Marmottant, S. Van Der Meer, M. Emmer, M. Versluis, N. De Jong, S. Hilgenfeldt and D. Lohse, J. Acoust. Soc. Am., 2005, 118, 3499–3505.
71. M. Emmer, A. van Wamel, D. Goertz and N. de Jong, Ultrasound Med. Biol., 2007, 33, 941–949.
72. B. Chapman and M. Plesset, J. Fluids Eng. Trans. ASME, 1971, 93, 373–376.
73. A. Prosperetti, J. Acoust. Soc. Am., 1977, 61, 17–27.
74. A. Prosperetti, L. Crum and K. Commander, J. Acoust. Soc. Am., 1988, 83, 502–514.
75. M. Versluis, E. Stride, G. Lajoinie, B. Dollet and T. Segers, Ultrasound Med. Biol., 2020, 46, 2117–2144.
76. R. Nigmatulin, N. Khabeev and F. Nagiev, Int. J. Heat Mass Transfer, 1981, 24, 1033–1044.
77. V. Kamath and A. Prosperetti, J. Acoust. Soc. Am., 1989, 85, 1538–1548.
78. G. Zhou, J. Acoust. Soc. Am., 2021, 149, 923–933.
79. J. Feshitan, C. Chen, J. Kwan and M. Borden, J. Colloid Interface Sci., 2009, 329, 316–324.
80. P. Prentice, M. MacDonald, T. Frank, A. Cuschieri, G. Spalding, W. Sibbett, P. Campbell and K. Dholakia, Opt. Express, 2004, 12, 593–600.
81. V. Garbin, D. Cojoc, E. Ferrari, R. Proietti, S. Cabrini and E. Di Fabrizio, Jpn. J. Appl. Phys., Part 1, 2005, 44, 5773–5776.
82. G. Hale and M. Querry, Appl. Opt., 1973, 12, 555–563.
83. R. Pope and E. Fry, Appl. Opt., 1997, 36, 8710–8723.
84. J. Kwan and M. Borden, Soft Matter, 2012, 8, 4756–4766.
85. E. Peterman, F. Gittes and C. Schmidt, Biophys. J., 2003, 84, 1308–1316.
86. M. Borden, G. Martinez, J. Ricker, N. Tsvetkova, M. Longo, R. Gillies, P. Dayton and K. Ferrara, Langmuir, 2006, 22, 4291–4297.
87. K. Kooiman, M. Emmer, T. Kokhuis, J. Bosch, H. De Gruiter, M. Van Royen, W. Van Cappellen, A. Houtsmuller, A. Van Der Steen and N. De Jong, Proceedings – IEEE Ultrasonics Symposium, 2010, pp. 900–903.
88. J. Dove, M. Borden and T. Murray, Opt. Lett., 2014, 39, 3732–3735.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d3sm00871a

---