Non-spherical oscillations drive the ultrasound-mediated release from targeted microbubbles

Ultrasound-driven microbubbles are attractive for a variety of applications in medicine, including real-time organ perfusion imaging and targeted molecular imaging. In ultrasound-mediated drug delivery, bubbles decorated with a functional payload become convenient transport vehicles and offer highly localized release. How to ef ﬁ ciently release and transport these nanomedicines to the target site remains unclear owing to the microscopic length scales and nanoseconds timescales of the process. Here, we show theoretically how nonspherical bubble oscillations lead ﬁ rst to local oversaturation, thereby inducing payload release, and then to microstreaming generation that initiates transport. Experimental vali-dation is achieved through ultra-high-speed imaging in an unconventional side-view at tens of nanoseconds timescales combined with high-speed ﬂ uorescence imaging to track the release of the payload. Transport distance and intrinsic bubble behavior are quanti ﬁ ed and agree well with the model. These results will allow for optimizing the therapeutic use of targeted microbubbles for precision medicine.

C linical ultrasound probably meets all essential requirements to be the ideal medical imaging modality: an excellent safety record, good imaging resolution and penetration depth, and a real-time, or even faster, acquisition speed. However, its specificity remains limited. Microbubble ultrasound contrast agents can alleviate this shortcoming. Clinical use of microbubbles has been made possible by their core physical property, namely resonance. Upon exposure to the oscillatory pressure variations of an ultrasound wave, compressibility of the gas core combined with inertia of the surrounding liquid leads to a mass-spring-like resonance behavior whose Minnaert eigenfrequency obeys an inverse relation with the bubble radius 1 . The unique acoustic properties of microbubbles that present a typical radius of 1−3 μm have also been used for two decades to boost perfusion imaging in the clinic, making them unrivaled blood pool agents 2 .
In vivo, microbubbles have to circulate in the blood stream until reaching the target area. Extensive efforts have therefore been put into enhancing the stability of microbubbles by coating them with polymers, proteins, or lipids. Among these coatings, phospholipid shells have been the most widely used owing to their high biocompatibility, flexibility, and enhanced non-linear acoustic properties 3 . Furthermore, a lipid coating can be easily functionalized, e.g. by including stealth capability by adding polyethylene glycol (PEG), which prevents interactions with the immune system and prolongs the circulation time 4 . Moreover, targeting ligands can be attached to the periphery of the shell to specifically recognize and adhere to diseased cells and tissues 5,6 , thereby bringing a molecular imaging dimension to the clinical use of microbubbles. The same ligands can also be used to load the microbubbles with various nanoconstructs, turning them into efficient microcarriers. The number of formulations for drugloaded microbubbles has vastly expanded. These include a broad range of interfacial structures, from simple drugs loaded onto the bubble shell 7,8 to more complex nanoconstructs that can entrap genes or chemotherapeutics 9,10 .
Driven near its resonance frequency, a bubble displays maximal radial response and generates secondary effects, such as harmonics and subharmonics 11,12 , streaming, acoustic radiation forces, shape instabilities 13 , and non-spherical oscillations 14,15 . These effects are of prime importance in applications such as cleaning, (bacterial) biofilms removal 16 , mechanical destruction of thrombus 17 or tumors 18,19 or inducing vessel wall permeation, e.g for blood brain barrier opening 20 . The combination of the mechanical action of targeted microbubbles adherent to a surface with the release of a drug payload is of great interest for therapeutic applications as medicine is becoming increasingly focused towards personalized 21 and localized therapy 22,23 .
Although there is convincing experimental evidence for controlled release of the payload, the physical mechanisms behind it remain unclear. Some attempts were made to explain these phenomena: Borden et al. 24 proposed that the surface area reduction leads to the expulsion of excessive shell material, which is, however, difficult to assess owing to the nanoseconds timescales and nanometer length scales governing the problem. O'Brien et al. 25 attempted to explain the shedding in relation to bubble stability by looking at molecular viscosity, while Kwan et al. 26 concluded that the release could originate from a cyclic nucleation and aggregation of lipid folding events. This, however, is difficult to combine with the observation that shedding can occur within just a few ultrasound cycles 27 .
Here, we take a major step forward, both in the experimental observation of the microbubble shedding and in its physical description within the relevant pressure regime 28 . Time-resolved observation of the release was performed in an unconventional side view, which allowed the visualization of the release in the relevant plane of sight. The setup combines two high-speed cameras in order to simultaneously record the bubble oscillations at a tens of nanoseconds timescale, together with high-speed fluorescence imaging of the release and transport of the bubble payload at a tens of microseconds timescale. To this end, a fluorescent dye (DiI) was incorporated into the bubble shell to serve as a model drug to unravel the release and transport mechanisms of the payload at a range of acoustic parameters. We also study the controlled release theoretically, where we show good agreement with the experimental observations. These observations shed light on the physical mechanisms involved in drug delivery with microbubbles, and demonstrate the importance of the physical environment for the efficacy of the controlled release.

Results
Microbubble dynamics theory and modeling. The dynamics of free bubbles is described by the Rayleigh−Plesset equation and has been investigated for a century. It can theoretically account for large oscillation amplitudes and water compressibility 29 , nonadiabatic compression 30 and complex viscoelastic behavior of the bubble shell 31,32 . In practice, however, microbubbles are never in free space and further modifications were made to the primary models to account for the presence of a substrate. The so-called method of images is such an example in potential flow theory that considers the substrate as an acoustic reflector and offers an easy way to account for its effect on the bubble dynamics 33 . More extended spherical models were developed 34,35 and tested 36,37 to investigate the effect of the distance between the bubble and the substrate. However, the present physical problem requires going beyond the spherical models by considering the non-spherical dynamics of microbubble oscillations and the subsequent release and transport of the bubble payload.
Successfully explaining the physical aspects involved in shedding from microbubbles implies an understanding and description of the microbubble oscillations without the constraint of spherical symmetry. In this section we therefore build a potential flow model based on a simple idealization of the flow kinematics for the microbubble dynamics near a substrate, where the bubble is allowed a non-spherical axisymmetric motion (see Fig. 1).
First, we evaluate the expected displacement of the substrate as a result of the arrival of ultrasound, which can be estimated from a no-slip boundary condition and from liquid compressibility. For the studied range of pressures (150−400 kPa) at a driving frequency of 1 MHz, the amplitude of motion of a water fluid particle induced by the ultrasound wave is not more than 50 nm 38 . Such a displacement is negligible as compared to the observed bubble oscillation amplitude that typically reaches a few Fig. 1 Schematic description of the model. a Description of a bubble targeted to a substrate including the geometrical parameters. b The bubble undergoes volumetric oscillations and oscillatory translations under ultrasound exposure, thereby generating streaming micrometers. Consequently, the displacement of a substrate with larger acoustic impedance is therefore also negligible. In addition, an oscillating bubble will preferentially displace water rather than inducing a local deformation of the substrate owing to the substrate's elasticity (E = 3 GPa for polystyrene). Thus, the substrate can then be identified as a rigid plane where we neglect its deformation along the contact area with the bubble. The center of this contact area is then taken as the origin to build our axisymmetric model. We now describe the radius R not only as a function of time but also of the elevation angle ϑ, Fig. 1a. Note that in the following the geometry used for the mathematical derivation remains spherical, since that condition is required to make use of the Rayleigh−Plesset dynamics. We consider in first approximation the velocity potential as radial in this coordinate system. The error associated with this assumption is discussed in Supplementary Note 1. This now allows us to rewrite the Rayleigh −Plesset equation in spherical coordinates (r, ϑ): with ρ ' the liquid density, and R(ϑ, t) the time-dependent and angle-dependent radius with the overdots representing its time derivatives. P i is the internal gas pressure: Rðϑ; tÞ 2 À ΔPðϑ; tÞ þ P 0 and P e is the external gas pressure: Here, μ ' is the liquid dynamic viscosity, P 0 is the ambient pressure, P A is the applied ultrasound pressure, V is the bubble volume with V 0 the initial bubble volume, γ is the polytropic exponent of the gas, and κ R 0 is the shell viscosity. ΔP is the local pressure drop across the interface and the Young−Laplace equation is used to generalize the pressure drop due to surface tension: where σ w is the surface tension of water, r loc is the timedependent local curvature on the bubble surface in the observation plane and R c is the time-dependent radius of curvature in the orthogonal direction. The radius of curvature R c can be estimated from the bubble volume V ¼ ð4πR 3 c Þ=3. The exact derivation of the curvature of the bubble and a discussion on this estimation can be found in Supplementary Note 2.
For numerical evaluation each bubble can be discretized and divided into angular segments of length R = R(ϑ, t) and width dϑ that link the surface of the bubble to the origin at an elevation angle ϑ (see Fig. 1b). In general, each segment of the microbubble is considered to obey Eq. (1) and all segments together are coupled through the gas volume V, that in turn determines the internal gas pressure Eq. (2), and through r loc , that is calculated using a finite difference method from the two closest neighbors. In many cases, the use of functionalized bubbles includes targeting to specific cell receptors. As a consequence, bubbles adhere to the functionalized supports upon contact and this adherence translates into the pinning of the triple line separating liquid, gas, and the substrate. The proposed way of modeling the non-spherical bubble dynamics gives a convenient and straightforward way of handling this pinning by setting the velocity of the first segment to zero. The resulting set of differential equations (one equation per segment) is then solved using the ODE113 solver in Matlab.
Describing the microbubble in this segmented picture has implications that we discuss in more detail now. First, the bubble volume being a summation of all segment volumes can be expressed as the sum of the volume fraction belonging to the segment with index i and the remaining volume V r : with dϑ the angle discretization step chosen for the model. When introducing Eq. (5) into Eq. (2), it becomes clear from Eq. (1) that (i) each segment will be a second-order oscillator whose eigenfrequency and phase behavior depends entirely on its initial length, identical to the Minnaert eigenfrequency of a spherical bubble. And (ii), that there is now an additional driving term V r in Eq. (1) that is governed by the dynamics of all the other segments and that interferes (physically) with the original acoustic driving term. As a direct consequence of (i), the phase of the oscillations will increase for increasing initial segment length and therefore with increasing ϑ (see Supplementary Note 3). The segments with the smallest angles expand first, leading to an oblate shape during microbubble growth. The segments with the smallest angles also retract first, leading to a prolate shape upon collapse. This phase difference along the surface of the microbubble will prove crucial for both the release of the shell material and its subsequent transport. The numerical model treats coupled oscillating segments and therefore allows for the existence of surface waves at the bubble interface. These surface waves are known to attenuate over a typical time that scales with the wavelength as To account for this, we implement a damping correction for each segment in the numerical model with a damping term of the form κ _ R: with κ = 2.3×10 −5 m 2 s −1 the only empirical parameter of the model.
Non-spherical bubble oscillations. Figure 2a shows a typical radius−time curve R c (t) of a 3.9 μm radius targeted microbubble driven by a 100-cycle ultrasound pulse at a pressure of 166 kPa. The recording was taken at a frame rate near 7.5 million frames per second, and are shown here together with the simulated response at this pressure and for this particular bubble size. Figure 2b shows the corresponding simulated and measured contours R(ϑ, t). Figure 2c, d also shows simulated and experimental contours of microbubbles exposed to pressures of 249 and 331 kPa, respectively. Supplementary Movie 4 shows the simulated bubble dynamics for a 3.1-μm bubble exposed to 84, 166, and 210 kPa and Supplementary Movie 5 shows the simulated bubble dynamics for a 2.4, 3.1, and 4-μm bubble exposed to 210 kPa. From these contours, one can appreciate how the bubbles take an oblate shape while growing and a prolate shape while collapsing, as predicted. Such observations were made earlier 14,15 and the physical explanation is now provided here.
Lipid oversaturation. An important implication of considering the phase as a function of the elevation angle is that the collapse of the bubble is now regarded as non-spherical. This becomes obvious when taking a closer look at the oscillation dynamics near the substrate, where the contact line is pinned. Until now, ultrasound-mediated release mechanisms of lipids or nanoparticles from microbubbles have been described in the context of spherical oscillations alone [24][25][26] and the most commonly considered mechanism to explain detachment of the shell material from the surface of a microbubble involves oversaturation of the lipid layer due to the surface area reduction upon compression, as it is the most intuitive. As mentioned before, in practice, both in vitro and in vivo microbubbles are never in free space. In the free field, the total collapse of a spherical bubble is prevented by a rapid increase of pressure and temperature in the gas core. In the presence of highly non-spherical oscillations, however, the earlier collapse of the bubble segment at small elevation angles is not prevented by such a pressure increase, owing to the phase lag of the bubble segments at a larger elevation angle. Thus, given a sufficiently high acoustic driving pressure, one can expect that nothing prevents this particular segment of the bubble from fully collapsing. This cusp formation process is demonstrated in Fig. 3, showing a simulated lipid shell collapse for such a pinned microbubble for increasing pressures. From Fig. 3a it also becomes clear that the aspherical collapse leads to an oversaturation of phospholipids at a localized position between the cusp and the substrate. Plotted in an angle-resolved view, Fig. 3b shows the change in surface area of each segment in the collapse phase for increasing acoustic pressure. The corresponding increase in the local phospholipid concentration exceeds a factor of ten at pressures in excess of 166 kPa. Figure 3c also shows that the cusp formation becomes less pronounced for larger bubble sizes, away from resonance, as will be detailed in the following section.
Although it is well-known that the presence of a substrate can change bubble resonance behavior 34,39 , little investigations were performed on its effect on the sphericity of the oscillations in this pressure range and in this plane of observation. Unlike our preliminary study 27 where we studied microbubbles floating up against a substrate by buoyancy, here targeted microbubbles were used, resulting in the pinning of the area in contact with the Opticell™ membrane. It was found experimentally, see Supplementary Note 4, that the radius of the pinned contact line was on average half the bubble radius with no significant dependency on the bubble size or other quantifiable parameter. r contact = R 0 /2 was therefore used for the simulation. It should be kept in mind that there is a fair amount of variability and the contact radius can vary from 25% up to 70% of the bubble radius for the extreme cases. On the other hand, the size of the contact area has limited impact on the simulation output, see Supplementary Note 5. Even removing the pinning altogether (while the bubble kept contact with the substrate) in the simulation leads to the very same conclusions. With decreasing contact area the shape of the bubble at small ϑ angle turns from a cusp into a tip for which all the following arguments also apply. We note that this tip-like shape was observed before 15 for the bubbles presenting the smallest contact areas. Also, the pinning behavior of the targeted bubble did not change the previously observed release behavior 27 , as confirmed by our preliminary experimental verifications and by the simulations.
Microstreaming and transport. It was proposed before that microstreaming surrounding the microbubble during ultrasound exposure is an important transport mechanism following the release events 27 . The axisymmetric velocity field resulting from simplified non-spherical bubble oscillations was calculated before by Marmottant and Hilgenfeldt 40 by assuming the bubble motion to be a superposition of a volumetric oscillation and an oscillatory translation perpendicular to the substrate and with a crucial phase difference Δφ between the two. The resulting flow field, u(r, z), see Fig. 1, is given by a combination of three Stokes singularities: a stokeslet (stk), a dipole (dip), and a hexadecapole (hexdp) 41 . The streamlines calculated for a 2.5 μm radius bubble exposed to a burst of 100 cycles at 1 MHz are shown in Fig. 4 and explain the trajectory followed by the released material: the phospholipids detach due to the local oversaturation near the substrate (see subsection Lipid oversaturation) and are dragged along the microbubble surface with the local flow field. Once at the top of the bubble, the material is transported along the centerline (z-axis in our nomenclature) and away from the bubble as illustrated in Fig. 4a. The experimental evidence of this unintuitive behavior is captured in a side view recording using high-speed fluorescence imaging at a frame rate of 50,000 frames per second. Figure 4b shows a typical example of fluorescent material transported during several tens of microseconds. In about one third of the cases, a secondary pinched-off bubble, itself carrying some fluorescent material, was traveling along the central streamline together with the shed material without, however, having a measurable effect on the transport. The velocities observed were The velocity of the fluorescently labeled molecules drops to zero as soon as ultrasound stops. This is a consequence of the highly viscous behavior of the surrounding water at these ultrashort length scales.
To describe the centerline transport along the z-axis (r = 0) the general flow field description u(r, z) can be simplified to give: ð9Þ with z as described in Fig. 1, d the initial distance of the center of the bubble to the substrate, thus slightly smaller than the bubble radius due to the pinning, and where the overbars are used to indicate that the space variable is non-dimensionalized with the initial bubble radius R 0 . ϵ z and ϵ R are the relative volumetric oscillation amplitudes and the relative oscillatory translation amplitudes, respectively, as defined in ref. 40 . The prefactor in Eq. (7) gives the typical velocity of the streaming and is proportional to the sine of the phase difference Δφ between volumetric oscillations and oscillatory translations. Thus, in this description, streaming can indeed only arise from non-spherical oscillations. The amplitude of the volumetric oscillations and the oscillatory translations, as well as their phases, and their phase difference, follow directly from the simulation. Figure 5a displays the volumetric oscillation and oscillatory translation amplitudes as a function of the bubble radius at a driving pressure of 100 kPa and a frequency of 1 MHz. It is evident that it displays strong resonance behavior. Similarly, the corresponding phase of the simulated volumetric oscillations and oscillatory translations, Fig. 5b, displays the characteristics of a resonant system. In addition, Fig. 5b plots the calculated phase difference for different pressures. Note that the bubble size corresponding to the maximum phase difference is found to increase with increasing pressure.
Bubble streaming-induced transport. The transport distance s(t) of a fluid element along the centerline (z-axis) is calculated by integrating the Lagrangian velocity of the particle, Eq. (7), from time zero to t: The transport distance is measured directly from the high-speed fluorescence recordings. Figure 6a plots the experimental transport distance vs. time. The integration in Eq. (11) is carried out over the duration of the ultrasound burst, in this case 100 μs. The calculated transport distance is valid for bubbles in a stationary oscillation regime and in a stationary streaming flow field, i.e. in the absence of transient effects, and with the bubble kept at a fixed location. Sometimes, and even more so at elevated pressures, the bubbles break-up, jet, pinch-off daughter bubbles or move rapidly along the substrate while disturbing the streaming field. For the targeted microbubbles that do fulfill the conditions of the model, the measured transport distance is plotted in Fig. 6c-e together with the predicted transport distance (see Supplementary Note 6 for details), where we find overall very good agreement. Note that the experimental results show some variability in the observed transport distance, owing to the large variations in the nonlinear bubble response and subsequent second-order effects, which includes bubble streaming. Finally, coming back to the time evolution of the bubble streaming-induced transport, Fig. 6b shows the calculated transport distances at resonance, showing that most efficient transport is achieved within the first 200 μs and that the effect of longer pulse duration, and notably increased pressure, is only marginal, since a doubling of the pressure only leads to a 25% increase of the transport distance.

Discussion
The presence of a shell has two effects on bubble dynamics: an increased damping due to shell viscosity and an increased stiffness due to shell elasticity. Several models were proposed to account for bubble shells 42,43 . The present model includes a shell viscosity (see Eq. (2)), which mainly increases the required pressure to reach the same oscillation amplitude and has otherwise negligible impact on the response. Rigorously considering the shell elasticity for non-spherical oscillation is complex as the response of a lipid monolayer to bending and the phospholipids dynamics on the bubble surface at MHz rates and on the microscale is unknown. Here, the coating is almost entirely shed within the first few microseconds, and therefore shell elasticity is neglected to compute the streaming. Adding shell elasticity to our computation, according to the Marmottant model 31 , induces a shift in the resonant radius but resulted in comparable surface reduction plots at resonance (Fig. 3a) and led to the same physical description of the shedding process. A contribution of shell elasticity was therefore not included. The release of material from the bubble coating is often accompanied by the pinch-off of smaller submicron-sized bubbles. A direct implication is the immediate size reduction of the mother microbubble. This pinch-off process was also reported in earlier work 44,45 . In the present experiments, these events were very difficult to control and highly irreproducible. It appears that the secondary bubbles carry some fluorescent shell material with them, while moving away from the bubble. Pinch-off occurred in about one third of the cases, with no obvious dependency on bubble size or oscillation amplitude.
For larger bubbles (≥4 μm) and higher pressures (≥350 kPa), jet formation was often observed. Such jets were reported before 46 : following asymmetrical collapse, the jet penetrates the bubble core and impinges onto the substrate. The role of jetting on cell poration is often discussed, but its resulting effects on the microbubble shell or on the delivery of a payload has received little attention. During collapse, the fluorescent coating appears to be following the jet, leading to its deposition onto the supporting substrate. Although jetting is not our present focus, this effect could certainly be of further interest for drug delivery with microbubbles.
The present study was performed on a model drug inserted between the phospholipids that coat the microbubble. A previous study 27 has demonstrated that microbubbles loaded with liposomes or polymeric nanoparticles quantitatively display identical behavior when compared in terms of microbubble oscillation amplitude. Changing the type of payload is therefore not expected to change the conclusions presented here. The ability of the targeted microbubbles to deliver various payloads over a distance several times larger than their own size presents great potential for targeted drug delivery and theranostic applications.
The behavior of phospholipid molecules on the microbubble surface is highly complex and the subject of much attention from researchers in the field of interfacial chemistry. Quasi-static Langmuir trough measurements state that the shell should collapse when reducing the bubble size by only a few percent 26,47 . It is however known not to hold for ultrasound-driven microbubbles, in particular due to the fast dynamical aspects 48 . On short timescales, the residence time of molecules at the interface will be related to the desorption potential barrier and to the viscosity associated with molecular motion 49 . Also, the compression of the phospholipids in a localized area of the surface can be expected to result in molecular transport along the shell, motion that is also subjected to molecular friction. Motion of individual molecules and the associated molecular processes remain to be explored. Interesting insights here may be provided, e.g. by molecular dynamic simulations.
Finally, this work highlights the importance of the substrate that supports the microbubble. Note that in the proposed model, the substrate is considered rigid in a fluid dynamical sense, which is justified by the very short length scales involved, as well as what  is observed in experiment. The results presented here improve our understanding of the delivery of a drug from targeted microbubbles and thereby our control over microbubble shedding. The present configuration of payload release along the symmetry axis of the adherent bubble can be highly efficient in small vessels and capillaries. However, this configuration may be suboptimal for sonoporation or sonothrombolysis, where the transport should preferentially be directed toward the target cells.
In order to optimize drug delivery, next steps should extend the present work towards more compliant substrates, as found in vivo, and define more advanced strategies to control the transport direction.

Methods
Numerical methods. The set of second-order differential equations of Eq. (1)  Sample preparation. The topside of an OptiCell (Thermo Fisher Scientific, Waltham, MA, USA) was coated by adding 1 mL of 1 μg mL −1 solution of Neu-trAvidin (Life Technologies Europe, Bleiswijk, The Netherlands) in phosphatebuffered saline (PBS) (Life Technologies Europe) 24 h prior to the experiment. The surface was then rinsed with PBS to remove all unbound proteins and incubated for 1 h with 1% bovine serum albumin (Sigma-Aldrich) to prevent unspecific binding. Afterwards, the surface was cut into rectangular pieces of 5 mm × 20 mm. Biotinylated microbubbles were injected into the tank and allowed to adhere to the surface by flotation for 5 min.
Experimental methods. The substrate was immersed in a water bath (T ≈ 22°C), filled with demineralized water, under a 45°angle with respect to the vertical, and could be moved with a 3D micropositioning stage. A focused, single-element ultrasound transducer (C302, Panametrics, Waltham, MA, USA) with a center frequency of 1 MHz was fixed on one side of the water tank and used for insonation. A detailed schematic of the setup is given in Supplementary Note 7. For each recording, a single microbubble was exposed to a single ultrasound burst. The acoustic pressures were measured with a calibrated needle hydrophone (0.2 mm, Precision Acoustics, Dorchester, UK). Optical recordings were performed using a microscope equipped with a 20x water-immersion objective (LUMPFL, Olympus, Zoeterwoude, The Netherlands) that was fixed at a 45°angle relative to the vertical (with the optical axis parallel to the substrate). A continuous-wave diode laser (5 W, λ = 532 nm; Cohlibri, Lightline, Osnabrück, Germany) was focused onto the sample through the same objective, using a dichroic mirror. The laser light was gated using an acousto-optic modulator (AOTF.nC-VIS, AA Optoelectronic, Orsay, France) to generate a 500 μs laser pulse. Simultaneously, a KL 2500 LED light source (Schott, Mainz, Germany) was used to superimpose bright-field and fluorescence. High-speed recordings were acquired with a high-speed camera (SA-X, Photron, West Wycombe, UK), operating at 5000 to 50,000 fps. Ultra-highspeed imaging was performed with the Brandaris 128 camera 51 at frame rates ranging from 5 to 10 million fps. Using a beam splitter, 80% of the light was directed into the Brandaris camera, while the remaining 20% was directed towards the Photron camera. A test was performed for 12 targeted microbubbles recorded in a top-view configuration. Identical shedding was observed as compared to our previous study 27 , demonstrating that the targeting did not significantly alter the process.
Analysis. The recordings were contrast-enhanced, and analyzed using a thresholding method in MATLAB to extract the bubble contours R(ϑ, t). The bubble volume was obtained by assuming axial symmetry. Due to the more complex aspect of the released fluorescent material, the transport distance was extracted by manually selecting the center of the fluorescence patch in each frame of the recordings.
Data availability. The authors declare that all the data acquired during this study, as well as the Matlab codes used for simulating microbubbles responses and analyses will remain stored within the servers of the University of Twente and will be made available on request by the authors.