PFOB sonosensitive microdroplets: determining their interaction radii with focused ultrasound using MR thermometry and a Gaussian convolution kernel computation

Abstract

Purpose: Micron-sized perfluorocarbon droplet adjuvants to focused ultrasound therapies allow lower applied power, circumvent unwanted prefocal heating, and enhance thermal dose in highly perfused tissues. The heat enhancement has been shown to saturate at increasing concentrations. Experiments were performed to empirically model the saturating heating effects during focused ultrasound.

Materials and methods: The measurements were made at varying concentrations using magnetic resonance thermometry and focused ultrasound by circulating droplets of mean diameter 1.9 to 2.3 µm through a perfused phantom. A simulation was performed to estimate the interaction radius size, empirically.

Results: The interaction radius, representing the radius of a sphere encompassing 90% of the probability for the transformation of acoustic energy into heat deposition around a single droplet, was determined experimentally from ultrasonic absorption coefficient measurements The simulations suggest the interaction radius was approximately 12.5-fold larger than the geometrical radius of droplets, corresponding to an interaction volume on the order of 2000 larger than the geometrical volume.

Conclusions: The results provide information regarding the dose–response relationship from the droplets, a measure with 15% precision of their interaction radii with focused ultrasound, and subsequent insights into the underlying physical heating mechanism.

Keywords:

• focused ultrasound
• MR thermometry
• thermal ablation
• perfluorocarbon emulsions
• sonosensitizers
• microdroplets

1. Introduction

Perfluorocarbon emulsions have gained much interest in recent years for focused ultrasound therapies. The microdroplets have the potential to reduce applied energy and side effects. The potential focused ultrasound applications of these emulsions are relevant to ultrasound guided treatment and magnetic resonance guided treatment, and include: drug delivery, [1–3] lithotripsy, [4] sonothrombolysis, [5] and ablation [6,7]. These pharmaceuticals have not yet entered into clinical trials for ablative therapy and hyperthermia. Further characterization and testing might help translate the technique into a pilot study.

Focused ultrasound applications are complicated by high-power requirements, respiratory motion, unwanted beam path heating, and fluid perfusion difficulties [8–15]. Among the drawbacks to clinical hepatobiliary focused ultrasound studies have been skin burns, the potential for rib resection, or the need to administer intrapleural infusion of warm saline to circumvent unwanted heating [14,16]. By enhancing the acoustic absorption with sonosensitive droplets, the applied power, associated side effects, and treatment times might be reduced.

Ultrasound contrast agents in routine clinical use are gas-filled perfluorocarbon microbubbles with either sulfur hexafluoride, perflutren, or perfluorobutane core and a phospholipid or albumin shell. Applications include echocardiography for left ventricular opacification and endocardial border delineation; along with imaging the liver, breast, and blood perfusion. The four common agents are Optison (GE Healthcare), Sonazoid (GE Healthcare), Sonovue (Bracco Imaging, S.p.A.), and Definity (Lantheus Medical Imaging, Inc.). These microbubbles, also referred to as microspheres, are about 1.2 to 3.1 µm in diameter and are safely administered intravenously at concentrations between 0.65% v:v and 4.4% v:v [17].

The ultrasound contrast agents are nontoxic to kidney and thyroid, with only rare adverse effects including anaphylaxis and hypersensitive immune response [18]. Other research has suggested the possibility of adverse effects at high mechanical index: possible microvascular damage and arrythmia in echocardiography [18]. Sonovue microspheres containing sulfur hexafluoride gas of boiling point −51 °C administered at a 0.8% v:v dose (0.3 ml/kg), exhibit a mean diameter of 2.5 µm, with 90% of bubbles being less than 6 µm, and 99% of microbubbles being less than 11 µm [19]. Sonovue microbubbles adjuvant to ultrasound-guided focused ultrasound ablation have been shown to improve clinical outcomes for uterine fibroids, [20–24] adenomyosis, [25] and hepatocellular carcinoma [26]. The major benefit from the application of the microbubbles was the reduction of sonication times and treatment times without significant adverse effects such as skin burns.

Preclinical investigations of gas-filled perfluorocarbon microbubbles showed that thermal enhancement in biological tissues was created by nonlinear bubble oscillations when subjected to focused ultrasound. The use of focused ultrasound with gas-filled microbubbles demonstrated encouraging results to reduce the cavitation threshold and increase the thermal effect. However, the in vivo half-life of a few minutes is too short for many tumor ablation procedures. The near-field interactions with microbubbles and long treatment times has also shown adverse effects in tissue mimicking gels [27]. Research groups have reported many preclinical study results with the phase-shift liquid perfluorocarbon droplets for use with focused ultrasound ablation [28–36]. To date, no pilot studies or randomized controlled trials have utilized perfluorocarbon phase-shift droplets for focused ultrasound ablation. A number of theoretical models have been developed for acoustic droplet vaporization of phase-shift droplets [37–42]. Most of the models assume that nucleation and vaporization occur within the liquid perfluorocarbon and result in an increase in the bubble volume, causing an oscillatory behavior of expansion and compression that contributes to heat generation.

Liquid-core perfluorocarbon microdroplets are expected to increase the stability and half-life for potential ablation applications. Perfluorooctyl bromide (PFOB) has high oxygen solubility, low toxicity, a short biological half-life, a high boiling point (142 °C), low volatility, was temporarily FDA approved as an oxygen-carrying blood substitute, [43] and was recently granted orphan designation by the EMA for treatment of respiratory distress syndrome [44] and congenital pulmonary hypoplasia [45]. Thus, emulsified PFOB has been considered a good adjuvant to focused ultrasound therapies.

The mean diameter of emulsified droplets can be adjusted by modifying the surfactant nature and concentration and the manufacturing technique. The size was chosen such that the emulsions are transpulmonary, to target the tumor microvasculature and remain in continuous systemic circulation [43, 46–48]. The kidneys tend to excrete particles between approximately 1 to 6 nm, the liver tends to filter particles between approximately 50 to 100 nm, and the spleen filter particles between about 200 to 500 nm [49]. As these droplets are about 2 µm in diameter, they might have a tendency to mechanically aggregate in and be retained in small capillaries [50]. The liver filters substances through a capillary network of sinusoids expressing fenestrae pores that have been measured at 107 nm with transmission electron microscopy of healthy human liver biopsies [51]. These sinusoids are a common location of metastatic tumor sites. Blood from the hepatic arterial network [52] arrives at the liver sinusoids and metastatic tumor cells accumulate in the sinusoids [53]. Droplets less than about 100 nm can also undergo pinocytosis and enter into endothelial cells [54]. Previous microbubble studies showed pulmonary capillary networks are about 8 µm in diameter and can create a sieve effect on microbubbles larger than this size, impeding circulation, and increasing dissolution rate [55]. These droplets are large enough to prevent uptake by the liver sinusoids and spleen, but small enough to cross the pulmonary capillaries. Hence, the droplets should remain in circulation until being degraded by the mononuclear phagocyte system (MPS). Lecithin-PFOB droplets have been shown to accumulate in regions with high amounts of macrophages, such as tumor sites [56]. Methods to slow the uptake by the MPS include unreactive emulsifiers, surface coating with polyethylene glycol (PEG) moieties, and reducing droplet size [54].

The perfluorocarbon core is immiscible in aqueous solution and requires emulsification to stabilize the substance in blood plasma. The toxicity of PFOB emulsion has been reported to have a LD₅₀ between 14.72–45 g/kg in rodents [57–59]. Intravenous administration of 0.1–0.2 µm diameter lecithin-PFOB emulsions in rodents showed a half-life of 6 h after receiving doses of 2.5 g/kg, with 99.8% of the PFOB eliminated through the lungs, the remainder eliminated through feces, and none eliminated through urine [60, 61]. The half-life of PFOB in humans is about 3–4 days [62].

Unemulsified oral PFOB doses of 2–12 ml/kg showed no toxic effects within three days in 60 human subjects [60] and about 95% of PFOB is passed through the gastrointestinal tract within 24 h [63]. Lecithin-PFOB emulsions with 0.1–0.2 µm diameter showed no detectable toxicity in humans after 1 g/kg intravenous dose [60]. The applications of PFOB in human subjects include: a substitute for allogeneic blood transfusion during surgery, [64–66] a bronchodilator for mild atopic asthma, [67, 68] treatment of respiratory distress in premature infants, [69, 70] a CT contrast agent, [60, 71] an ultrasound contrast agent, [72] and an oral gastrointestinal MRI contrast agent [57, 73, 74].

The size of the droplets is very important not only for heating effects, but also for pharmacokinetics, biodistribution, and clearance. Pharmacokinetic studies of lecithin-PFOB emulsions indicate that the elimination route from the blood occurs through macrophage phagocytosis and subsequent elimination by pulmonary expiration [75]. Opsonization of the lecithin-PFOB emulsions is characterized by a coronal layer of absorbing plasma proteins expressing receptors for macrophages that eventually engulf the emulsions via phagocytosis [76]. Phagocytosis for these emulsions has been observed to occur primarily via Kupffer cells and splenocytes that results in perfluorocarbon deposits in the liver and spleen, where the molecules can remain at small amounts, relative to the intravenous dose, far longer than their circulatory half-life [77]. The macrophage immune response to the lecithin-PFOB emulsions results in relatively large deposits of PFOB in the liver and spleen, soon after administration [77]. The emulsions are then degraded, and unmetabolized perfluorocarbon molecules diffuse back into the blood stream where the molecules bind to blood plasma lipids [77]. The PFOB-bound lipoproteins are transported to the lungs, where the molecules are expelled through respiration [62]. It has been shown that in lecithin-PFOB droplets, essentially all PFOB is cleared from the lungs within a few days [60, 61].

Emulsions composed of F₈TAC₁₈-PFOB were used by Desgranges et al. [10] as sonosensitizers for enhancing the energy deposition. The microdroplets ranged from 0.67 to 4.07 µm mean diameter and were designed to circulate in the endovascular system, rather than accumulate in tumors via the enhanced permeability and retention effect. Liquid-core PFOB microdroplets were embedded in tissue mimicking gels at different concentrations and demonstrated significant thermal enhancement, with reproducible and reversible heating effects. Lorton et al. confirmed these results after administration of 0.1 to 0.2% v:v PFOB-filled microdroplets in in vitro and ex vivo perfused tissue models. The study reported stable microdroplet size for at least 30 days in appropriate storage conditions. An increase in the energy deposition by a factor of 2.5 to 3 was observed after microdroplet injection in the perfusion fluid of the tissue-mimicking model, freshly excised pig liver, and excised pig kidneys. The study demonstrated reproducible thermal effects when repeating sonications. This indicated that emulsions were stable for at least 90 min under those conditions. At this moment, there are no experimental data on the biological half-life of these micron-size droplets. Their diameter is comparable to standard ultrasound contrast agents suggesting no risk of embolism.

Emulsions composed of 1.61% v:v lecithin-PFOB were shown to increase the attenuation coefficient by 7.36 m⁻¹ in static human blood when measured with 4.7 MHz diagnostic ultrasound [78]. This is equivalent to 1.61 m⁻¹ at 1.031 MHz. This value includes both scattering and absorption contributions for low energy diagnostic ultrasound. Lorton et al showed a maximum increase of HIFU absorption coefficient from baseline of 1.2 m⁻¹ with 1.42-2.16 µm F₈TAC₁₈-PFOB droplets at 0.235% v:v in ex vivo pig kidneys [43].

Presented here is a technique to measure the ultrasonic absorption enhancement experimentally, as a function of sonosensitizer concentration. The study compares the model to experimental magnetic resonance thermometry data acquired from the previously described perfused tissue models. The primary interest is to determine the effective interaction radius of such a microdroplet with focused ultrasound using macroscopic quantities. Secondary, this technique might generate predictive in vivo behavior. Such an application could be the estimation of the administered concentration needed to achieve a desired thermal dose.

2. Materials and methods

The microdroplets were prepared using a fluorinated surfactant (F₈TAC₁₈) and a PFOB core as described in previous publications [10, 43, 79, 80]. The F₈TAC₁₈ fluorinated surfactant features a hydrophobic fluorinated tail group and a hydrophilic poly-tris(hydroxymethyl)aminomethane (or polyTRIS) head group. The F₈TAC₁₈ surfactant was synthesized by a one-step free radical polymerization [81] of the acryloyl monomer tris(hydroxymethyl)acrylamidomethane (THAM) using 1H,1H,2H,2H-perfluorodecanethiol (telogen) as transfer reagent and AIBN (α,α-azobisisobutyronitrile) as radical initiator [10]. As the surfactant molecular formula is tunable and influences the droplet size, it was chosen to yield the suitable droplet size for intra-capillary circulation (i.e., below 3–4 µm). An emulsion was prepared by sonicating a mixture of water, F₈TAC₁₈, and PFOB with a PT 3100 homogenizer (Kinematica, Luzern, Switzerland), with surfactant amount above its critical micelle concentration of 0.0027 mmol/L. Size distribution was characterized with a Mastersizer 2000 laser diffraction particle size measurement system (Malvern Instruments, Orsay, France) via Mie light scattering theory. The particles had mean diameters of 1.9 and 2.3 µm for the samples used in the experiments.

A bi-compartmental perfused phantom was prepared of 1.5% w:w agarose for the matrix, bringing rigidity and reducing permeability. The intra-capillary circulation compartment was prepared of 3% w:w agarose, 11% w:w glycerol and 15% w:w condensed milk by casting gel strands in a biopsy sheath. The perfused model mimicked the tissues represented by the rigid gel and intra-capillary circulation. The ratio between the circulating fluid and the tissue mimicking gel was 1:9. More details are described by Lorton et al. [43].

These were placed inside a custom container allowing stabilization and attached to a perfusion chamber to regulate fluid flow. The gel was cast within this phantom and a region in the center was embedded with gel strands for circulating the droplets. It is necessary to remove phantom air bubbles to prevent MRI artifacts and HIFU beam scattering. The presence of air bubbles is one of the main considerations for proper response. Air bubbles in the phantom are shown in the left image of Figure 1. The image on the right also shows the phantom after air bubbles were removed by continuous fluid exchange with a vacuum pump degassing compartment.

Figure 1. MR T1-weighted gradient recalled echo 3 D images. Left: Phantom with air bubbles. Right: Phantom without air bubbles. 1. Agarose-based matrix gel. 2. Agarose and condensed milk strands for fluid perfusion. 3. Ultrasound coupling gel. 4. Water. 5. HIFU transducer.

The acoustic field was modulated at the specified duty cycle by a pseudo-random phased array transducer (Imasonic, Besançon, France, 131 mm focal length, 140 mm aperture, 1031 kHz) and controlled by a beam former unit (Image Guided Therapy, Pessac, France). An impedance matching unit optimized the energy transfer between the beam former and the ultrasound transducer. The transducer was placed below the phantom apparatus within the scanner bore (3 T, Prisma Fit, Siemens, Erlangen, Germany) and an 11-cm loop coil was used as the receiver antenna.

The perfusate was degassed and assessed with MRI to confirm bubble removal. Temperature maps were generated from the phase maps using proton resonance frequency shift thermometry [43, 82]. Real-time control and positioning were monitored with Thermoguide software (Image Guided Therapy, Pessac, France). Maximum focal spot temperature elevation and the cumulative temperature field, referred to as integral temperature, were measured and equated to the acoustic absorption coefficient. The temperature maps parallel or perpendicular to the HIFU beam were measured at various particle concentrations ranging from 0.0-0.3% v:v in the perfusate. The acquisition parameters are provided in Lorton et al. [43].

The spatial integral of the temperature increase between the beginning and end of sonication was measured from the MRT temperature measurements. The ultrasonic absorption coefficient was formulated, relative to the baseline control condition, at each concentration, from the integral temperature according to Equation (1)(1) $$\alpha =\frac{\sum _1^N{\mathrm{T}}_i}{E\Delta Z}{c}_m\rho V_p$$(1) . This coefficient is useful to compare varying sonication conditions. The global (or, volumetric) temperature increase in the region of interest encompassing the sonication pattern, at the end of a sonication period was considered, aiming to overcome the thermal diffusion effect. A pixel-wise map of the absorption coefficient is not achievable with this method, however, a spatial variation of the absorption coefficient in our perfused model is not expected. Each sonication resulted in a measurement of the absorption coefficient, with about six sonications at each concentration. The averages and standard deviations were assessed. The values were in good agreement with an alternative method from Cornelis et al. [83] (1) $$\alpha =\frac{\sum _1^N{\mathrm{T}}_i}{E\Delta Z}{c}_m\rho V_p$$(1) $\alpha$=acoustic absorption coefficient, m⁻¹, $T_i$=temperature in voxel above threshold of 0.5 °C, N=number of voxels in MR image, $E$=FUS transducer energy emitted during sonication, J, $c_m$=specific heat capacity of tissue or phantom, J·kg⁻¹·°C⁻¹, $\rho$ =average voxel density of tissue or phantom, kg·m³${, V}_p$=voxel volume, m³$, \Delta Z$=voxel and slice thickness, m.

The interaction radius is conceptually the region beyond droplet surface an increased probability of temperature enhancement is expected. The interaction radius is defined as the radius of a sphere encompassing 90% of the probability for the transformation of acoustic energy into heat deposition around a single droplet (See Appendix). It illustrates the relationship between the physical droplet radius and the extended regime where heating effects are present. This model forms the basis of the following analysis used to generate the simulation. Each individual droplet was modeled as a sphere with radius, R_o. The radial distance to the edge of the effective heating regime is denoted as the interaction radius, R_int. The actual model uses a three-dimensional spherical Gaussian normal probability density function to relate the probability of acoustic energy transformation into heat near the particle.

A simplified one-dimensional description of the 3 D interaction radius model is shown in Figure 2 for illustration purpose.

Figure 2. (A) Depiction of relationship between droplet physical surface and theoretical interaction radius. (B) Interaction regime represented as a one-dimensional normal distribution of probability (P) for acoustic energy transformation into heat.

Figure 3 illustrates the outcome from overlapping interaction regimes of neighboring models. Figure 3(A) shows how two interaction regimes can overlap when the physical droplet surface is very close to adjacent droplets. The cumulative result of the convolution function is given in Figure 3(B). The simulation was first compared with multiple approaches using additive superposition without thresholding, subadditive cumulative effects, and with the thresholding effect shown in Figure 3(C). The actual simulation was performed in a three-dimensional Cartesian coordinate grid.

Figure 3. Depiction of theories for overlapping interaction regimes. Red lines represent individual heating regimes. Green lines represent effective regimes. P is the probability for acoustic energy transformation into heat.

The ‘critical concentration’ approximates the transition point from linear to nonlinear effects and is denoted as ‘C_o’. Linear effects would be that applying a double dose results in twice the amount of heating. Nonlinear effects describe that applying a double dose results in noticeably less than twice the amount of heating. This is the approximate concentration when interaction regimes begin to overlap from adjacent droplets. The acoustic parameters, for example, frequency, pressure at the source surface, F-number of the beam were identical in the peer-to-peer comparative experiments with/without sonosenzitizers.

A relationship between the critical concentration, interaction radius, and actual droplet radius can be estimated analytically in a first approximation. If the droplets are given a body-centered-cubic ideal distribution in space, and the critical concentration corresponds to the concentration where one interaction radius is just beginning to touch neighboring interaction radii, Equation (2)(2) $$\frac{{\mathit{R}}_{\mathit{int}}}{R_o}=\frac{\sqrt{3}}{4}\mathrm{*}\sqrt[3]{\frac{8\pi }{3}}\mathrm{*}\frac{1}{\sqrt[3]{C_o}}$$(2) can be derived: (2) $$\frac{{\mathit{R}}_{\mathit{int}}}{R_o}=\frac{\sqrt{3}}{4}\mathrm{*}\sqrt[3]{\frac{8\pi }{3}}\mathrm{*}\frac{1}{\sqrt[3]{C_o}}$$(2) ${\mathit{R}}_{\mathit{int}}$ = interaction radius, m$, R_o$ = initial radius, m$, C_o$ = critical concentration, v:v

Numerical computing was performed with MATLAB (The MathWorks, Inc., Natick, Massachusetts, USA) using three-dimensional convolution. This approach proved most useful as the function is compatible with GPU acceleration and tall arrays. Performing the simulations on the high-performance computing (HPC) cluster reduced the simulation time from four days to one hour. The illustration of the method is shown in Figure 4. The matrix was cubic with a side length of 1024 elements. Each matrix element represents a 250 nm voxel.

Figure 4. 3D Matrix Convolution. A zeros matrix was generated. Then values of 1 were inserted at random coordinates. Next, the binary matrix was convolved with a Gaussian kernel. The output matrix has Gaussian elements placed at each location the values of 1 were located. Inherent to convolution, overlapping Gaussians were superimposed, additively. The final matrix was set to a threshold value equal to the maximum of the Gaussian distribution.

All matrix voxels were initialized to zero. Particle positions corresponding to concentrations between 0.0 and 0.2% v:v, at 0.02% v:v intervals, were assigned with the MATLAB randi function If the distance to the center of the nearest droplet was less than twice the geometrical droplet radius, valued at 2 µm, the position was regenerated. Each generated particle was marked by setting its corresponding matrix voxel to 1, forming a binary matrix. The actual droplets used in the experiments had a mean diameter of 2.3 µm or 1.9 µm depending on the batch.

Numerical integration of the three-dimensional Gaussian probability distribution in spherical coordinates reveals that 90% of the distribution volume is contained within 2.5 standard deviations. This convention is analogous to hydrogen atomic radius to help standardize the interaction radius to represent the probability of temperature increase surrounding the droplets. Obviously, the numerical value of the interaction radius depends on the underlying convention. The simulated Gaussian kernel standard deviations ranged from 2-9 µm, which is equivalent to 90% cumulated probability interaction radii between 5-22.5 µm. The range of the standard deviations used for the simulations were initially estimated from Equation (2)(2) $$\frac{{\mathit{R}}_{\mathit{int}}}{R_o}=\frac{\sqrt{3}}{4}\mathrm{*}\sqrt[3]{\frac{8\pi }{3}}\mathrm{*}\frac{1}{\sqrt[3]{C_o}}$$(2) .

The spherical Gaussian kernel was generated with the normal distribution centered within the kernel matrix. The kernel was cubic with the side length of four standard deviations and rescaled to a maximum amplitude set conventionally to 255. Every droplet was assigned the same interaction radius at a given time point of the simulation, as illustrated in Figure 4, that is, same convolution kernel. Then the concentration (or, the number) of droplets was adjusted iteratively.

Then, the binary matrix was convoluted with the Gaussian kernel. This has the effect of placing the kernel at each particle position and summating all particle kernels. After summation, voxels with values greater than 255 were set to the threshold value of 255, to represent the saturating effect from the interaction of neighboring interaction radii. The illustration is shown in Figure 4. The simulation was iterated over multiple concentrations and droplet interaction radii. The final value of all-droplets cumulated probability was normalized to a maximum of 1, meaning a division by 255 times the size of the 3d matrix, 1024³. The experimental data were normalized to the horizontal asymptote versus the volume concentration and overlaid on the simulations.

Further, a new independent step of simulation was performed using a modified interaction radius.

The calculation required about 24 gigabytes of memory and jobs were submitted in a batch job array on the HPC cluster, with 8 CPUs per task and 4 GB per CPU. This allowed the calculations to be performed in parallel, substantially reducing the simulation time. A portrayal of simulation slices at constant concentration and increasing interaction radii is illustrated in Figure 5.

Figure 5. 2D median section through the 3D volume of simulation used to create a data point for the simulation results in Figure 7, with 250 nm resolution and 0.04% v:v droplets. Shown FOV is 250 µm. The interaction radius increases at (a) 7.5 µm, (b) 12.5 µm, (c) 17.5 µm, and (d) 22.5 µm, respectively. Nonlinear threshold-based 3D summation of overlapping probabilities was applied, see Figure 3(C). Note the saturating effect from adjacent droplets in particular for larger interaction radii.

An exponential function was fit to the experimental absorption coefficient averages according to Equation (3)(3) $$\alpha \left(x\right)=b\left(1-e^{-\frac{x}{C_o}}\right)$$(3) . (3) $$\alpha \left(x\right)=b\left(1-e^{-\frac{x}{C_o}}\right)$$(3) $x$ = measured emulsion concentration, v:v, $b$ = horizontal asymptote, $C_o$ = critical emulsion concentration, v:v.

Equation (3)(3) $$\alpha \left(x\right)=b\left(1-e^{-\frac{x}{C_o}}\right)$$(3) enables the computation of two important quantities. First, the asymptote, designated by the variable ‘b’, indicates the maximum absorption coefficient. Also obtained is the fit constant ‘C_o’. The variable is denoted as the ‘critical concentration’ and occurs at 68% of the maximum value. The fit constants were determined by nonlinear least squares curve fitting with the MATLAB lsqcurvefit function

3. Results

The thresholding approach, shown in Figure 3(C), best fits the experimental data. The other simulation results using subadditivity and superposition were omitted (data not shown). Figure 5 shows an example of the convolution kernel simulation for solution containing 0.04% v:v droplets. It is intended to illustrate how the interaction radius and ultrasound absorption surrounding the droplets might appear. The second image corresponds to the best fit to the experimental data and is likely most representative of the actual phenomenon. The summation of the matrices shown in Figure 5 were used as the data points of the solid lines in Figure 7.

Figure 6 shows the experimental acoustic absorption coefficients fit to Equation (3)(3) $$\alpha \left(x\right)=b\left(1-e^{-\frac{x}{C_o}}\right)$$(3) . The data are relative to the control baseline and normalized to the overall maximum. Each data point represents an average and standard deviation of about 6 sonications performed with the same sonication conditions and droplet concentration. The data in Figure 6 was derived using Equation (1)(1) $$\alpha =\frac{\sum _1^N{\mathrm{T}}_i}{E\Delta Z}{c}_m\rho V_p$$(1) , which includes the total energy emitted from the transducer.

Figure 6. Acoustic absorption coefficient measurements for a series of perfusate emulsion concentrations and focused ultrasound sonication conditions. The values illustrated relative to baseline for the control condition without droplets present and normalized to overall maximum.

Figure 7 illustrates the Gaussian kernel simulation data along with experimental data measurements. Condition 1, Condition 2, and Condition 3 refer to the experimental conditions listed in Table 1. It is evident that a threshold was reached at higher concentrations, where no substantial improvement was obtained from adding more emulsion to solution. Table 1 lists the sonication parameters, with the fit constants from Equation (3)(3) $$\alpha \left(x\right)=b\left(1-e^{-\frac{x}{C_o}}\right)$$(3) , and the analytical interaction radii determined from Equation (2)(2) $$\frac{{\mathit{R}}_{\mathit{int}}}{R_o}=\frac{\sqrt{3}}{4}\mathrm{*}\sqrt[3]{\frac{8\pi }{3}}\mathrm{*}\frac{1}{\sqrt[3]{C_o}}$$(2) . The analytical solution using a simplified binary model of body centered cubic network predicts the particle interaction radii between about 7-9 µm estimating the order of magnitude. The convolution kernel simulation results in Figure 7 matched very closely the experimental data within a radial distance of 12.5–13.75 µm for the actual convention of 90% cumulated probability, that is a factor of about 12.5 larger than the geometrical droplet radius.

Figure 7. Gaussian Kernel Simulations for a 1024x1024x1024 Matrix. Solid lines represent an interaction radius that is equal to 2.5 standard deviations of a spherical Gaussian distribution. These were normalized to the maximum possible attenuation in the simulation volume. The points represent the first order exponential fit line through the experimental data that has been normalized to the horizontal asymptote.

Table 1. Sonication parameters, empirical fit constants, and model interaction radius comparison.

	Condition 1	Condition 2	Condition 3
Formula	Emulsion 1	Emulsion 1	Emulsion 2
Range, % v:v	0.0-0.29	0.0-0.2	0.0-0.29
Time on, ms	990.0	950.0	950.0
Time off, ms	10.0	50.0	50.0
Executions	33.0	33.0	33.0
Total Time, sec	33.0	33.0	33.0
Duty Cycle, %	99.0	95.0	95.0
Applied Power, W	94.5	137.0	137.0
Emitted Energy, kJ	3.09	4.29	4.29
Flow Rate, mL/s	0.10	0.10	0.10
Slice Orientation	Coronal	Coronal	Transverse
Total Data points	36	38	50
Fit Constant, b, m^-1	0.47	0.32	0.90
Fit Constant, Co , v:v	0.0013	0.0018	0.0017
Actual Droplet Radius, Ro, µm	1.15	1.15	0.96
Body centered cubic model, Rint/Ro	8.01	7.21	7.33

4. Discussion

We are reporting an experimental method for estimating the interaction radius of high boiling point liquid core sonosensitizers with focused ultrasound, based on macroscopic measurements and a numerical model of overlapping effects. The results obtained on three independent batches are suggesting an accuracy of approximately 15%, based on the deviation of the critical concentrations C_o provided in Table 1, the error propagation via Equation (2)(2) $$\frac{{\mathit{R}}_{\mathit{int}}}{R_o}=\frac{\sqrt{3}}{4}\mathrm{*}\sqrt[3]{\frac{8\pi }{3}}\mathrm{*}\frac{1}{\sqrt[3]{C_o}}$$(2) using the cubic root and the intrinsic variation of the average geometrical radius of droplets R_o in different batches, approximately 10%. The interaction radius for a 90% heating volume was found approximately 12.5 times larger than the geometrical radius of the particle, corresponding to an interaction volume on the order of 2000 larger than the geometrical volume. The saturation effect occurs therefore at low concentration. This is a fundamental result which explains the high efficacy of these sonosensitizers. A previous report mentioned that the apparent acoustic absorption in the volume of one droplet is at least 10,000 higher than biological tissue. If this energy was converted into heat inside the droplets, this would yield very rapid evaporation. However, previous reports demonstrated that these droplets are not evaporated during their interaction with HIFU [10]. Therefore, the only explanation is that the heating effect is distributed in a larger volume outside the particle, three orders of magnitude larger than the geometrical volume of the particle itself, which is the main highlight of the current study.

Investigating the intrinsic physical mechanism of enhanced absorption of ultrasound around the microdroplets was outside the purpose of this report, that relies on macroscopic measurement of temperature elevation at the scale of millimeter and a numerical model of non-linear summation of probabilities. We recall that the concerned microdroplets trapped in a semi-rigid gel demonstrated a stable response to repeated sonications [10] and their enhancement of ultrasound absorption was a factor of 6 higher in a liquid emulsion versus the semi-rigid gel [43]. A combination between thermoviscous heating [84] and nucleation of stable cavitation in water [85] is compatible with this behavior, but further observational studies at microscopic scale are required to determine the balance between the two effects and the range of acoustic parameters that applies.

Measuring the thermal enhancement at conditions below the threshold for both stable and inertial cavitation can characterize the contributions from thermoviscous energy losses to attenuation. This might be performed with a method for determining the cavitation threshold using passive cavitation detection, Doppler ultrasound, or B-mode imaging [5, 34, 86–88]. Traditionally, this was done by monitoring the ultraharmonic and broadband emissions. Broadband emissions correspond to the inertial (also called transient) cavitation signal and sub-/ultraharmonic emissions correspond to the stable cavitation signal [5, 85, 89–91]. Methods have been developed to detect and quantify the contribution of cavitation to the ultrasonic absorption during MR-HIFU [87, 88]. A hydrophone can be used to monitor the pressures at which the emissions begin to occur. By recording the absorption above and below, these thresholds, the contributions can be deduced. This technique would elucidate the individual contributions to overall acoustic absorption at expected in vivo pressures and is the subject of future work with the droplets.

The concentrations used in the study are less than lecithin-PFOB emulsion concentrations administered during previous clinical trials. Assuming an average human body weight is 62 kg and the average human total blood volume is 4400 ml, [92] the concentrations were formulated in terms of g/kg. The density of PFOB is 1920 mg/mL. A 0.1% v:v concentration of F₈TAC₁₈-PFOB emulsion in the total blood volume is equivalent to 4.4 ml or 8.5 g of PFOB in circulation. This results in a 0.137 g/kg dose of PFOB emulsion. This dose is a factor of 14 less than previous studies with lecithin-PFOB emulsions administered at 1.2–1.8 g/kg [64].

According to the discussion detailed in Lorton et al. [43] the corresponding doses of PFOB and F₈TAC₁₈ surfactant would be several hundred and several thousand, respectively, lower the extrapolated LD₅₀ from rats. The droplet shell in this study is composed of biocompatible fluorinated surfactant (F₈TAC₁₈) and the core consists of perfluorooctyl bromide. The beneficial features of this F₈TAC₁₈ surfactant include biocompatibility in the absence of hemolytic activity, LD₅₀ less than 4.5 g/kg in rats, and a half-life of about 40 h in rats [93, 94]. The molecule was mainly eliminated unmetabolized in urine, did not cross the blood–brain barrier, and tended to distribute in tissues at amounts near the blood plasma concentration, without affinity for specific tissues [94].

Further studies should address the biological half-life of the F₈TAC₁₈-PFOB emulsions and would help to optimize the surfactant formulation and the strategy for concomitant administration of sonosensitizer and HIFU sonication. The droplets might be cleared by MPS and/or deposited in the liver and spleen after administration, in vivo. This would prevent the long circulation times seen in ex vivo studies. The surfactant used here, F₈TAC₁₈, has been shown to be quite inert in rodents (unpublished data) and may ensure sufficient duration of the droplet circulation into the bloodstream. Techniques such as pegylation might be further considered to prevent rapid MPS clearance and prolong circulation times [95].

According to Desgranges et al. [10], the same technology allowed synthesizing PFOB-FTAC droplets of smaller size than used in this paper, for example, 620 nm, to be considered if this condition would demonstrate to be advantageous in vivo. Alternatively, a minimally invasive approach of locoregional intra-arterial infusion (e.g., hepatic artery) with fractionated bolus of sonosensitizer prior to each 30 s–60s sonication. This approach is common for contrast enhanced angiography and local chemotherapy delivery [52].

5. Conclusions

The results of the study were able to illustrate the saturating effect to HIFU absorption that perfluorocarbon microdroplets exhibit at increasing concentrations. The experimental results provide a measure of the amount of droplets needed to obtain a particular increase in focused ultrasound absorption. The absorption behavior could be fit with an exponential function to determine a measure of ultrasonic absorption coefficient and allows us to predict the size of the heated volume surrounding the droplet. The model describes an interaction region outside the physical droplet volume where higher probabilities of temperature increase, relative to the bulk medium, were present. Two models were derived to predict the size of this interaction radius; one was determined analytically in first approximation and the other determined semi-empirically from a convolution kernel simulation. The simulation predicts this interaction region to be a factor of 12.5 larger than the geometrical droplet radius. The results provide a measure for a dose-response relationship from the application of the droplets, as well as insight on the underlying physical mechanism for thermal enhancement.

Acknowledgements

The Center for Biomedical Imaging at University of Geneva is acknowledged for magnetic resonance facility access. The computations were performed at University of Geneva on the ‘Baobab’ and ‘Yggdrasil’ HPC clusters.

Disclosure statement

The authors declare that they have no conflict of interest.

Additional information

Funding

This project has received funding from the European Union’s Horizon 2020 research and innovation program under the Skłodowska-Curie grant agreement No 813766, and from the Swiss National Science Foundation grant CR32I3_162722.

Appendix

The kernel was generated by a Gaussian function in three-dimensional Cartesian coordinate space. The binary matrix represents the central coordinates of the droplet probability distributions in the simulations. The positions were generated randomly, ensuring the positions were separated by more than the physical droplet diameters. The output probability distribution matrix was generated by convoluting the binary matrix with the Gaussian kernel.

Three-dimensional Gaussian kernel function represented by a Cartesian coordinate system. $$G(x,y,z)=A·\hspace{0.17em}\text{exp}\hspace{0.17em}\left(-\frac{x^2+y^2+z^2}{2{\sigma }^2}\right)$$

Matrix representing the binary matrix of central coordinates of the droplet locations. $$M(x,y,z)=\mathrm{ }\left[m_{i,j,k}\right]\in \mathbb{E}\left(\left\{0,1\right\}\right)$$

Convolution of the two functions to obtain a matrix representing the probability of converting ultrasound into heat. $$\mathbb{P}(x,y,z)=M⨂G$$

The 90% volume of a three-dimensional Gaussian probability distribution function can be equated to the standard deviation as follows. The volume of the Gaussian kernel can be found by a triple integral. Converting to spherical coordinates simplifies the analysis. The function was integrated at 2.5 standard deviations of the Gaussian distribution and at infinite to show the measure of the 90% interaction radius volume. The integrals were solved numerically with the MATLAB integral3 function with the constants A and σ set to unity.

Triple integral of three-dimensional Gaussian kernel function represented by a Cartesian coordinates system. $$V\left(x,y,z\right)=\iiint A·\hspace{0.17em}\text{exp}\hspace{0.17em}\left(-\frac{x^2+y^2+z^2}{2{\sigma }^2}\right)\mathit{d}\mathit{xdydz}$$

Triple integral of three-dimensional Gaussian kernel function represented by a radial coordinates system. $$V\left(\rho ,\phi ,\theta \right)=\iiint A·\hspace{0.17em}\text{exp}\hspace{0.17em}\left(-\frac{{\rho }^2}{2{\sigma }^2}\right){\rho }^2\mathrm{s}\mathrm{in}\left(\theta \right)\mathit{\text{dρdφdθ}}$$

Integral representing the Gaussian volume at a distance of infinite from the central coordinate. $$V_{\mathrm{\infty }}={\int }_0^{\pi }{\int }_0^{2\pi }{\int }_0^{\mathrm{\infty }}A·\hspace{0.17em}\text{exp}\hspace{0.17em}\left(-\frac{{\rho }^2}{2{\sigma }^2}\right){\rho }^2\text{sin}\left(\theta \right)\mathit{\text{dρdφdθ}}$$

Integral representing the Gaussian volume at a distance of 2.5 standard deviations from the central coordinate. $$V_{2.5\sigma }={\int }_0^{\pi }{\int }_0^{2\pi }{\int }_0^{2.5\sigma }A·\hspace{0.17em}\text{exp}\hspace{0.17em}\left(-\frac{{\rho }^2}{2{\sigma }^2}\right){\rho }^2\text{sin}\left(\theta \right)\mathit{\text{dρdφdθ}}$$

Dividing the two integral to get the proportion of the maximum. $${V_{2.5\sigma }/V}_{\mathrm{\infty }}=0.90$$