Abstract

Low-intensity focused ultrasound (FUS), combined with microbubbles, is able to locally, and noninvasively, open the blood-brain barrier (BBB), allowing nanoparticles to enter the brain. We present here a study on the diffusion process of gadolinium-based MRI contrast agents within the brain extracellular space after ultrasound-induced BBB permeabilization. Three compounds were tested (MultiHance, Gadovist, and Dotarem). We characterized their diffusion through in vivo experimental tests supported by theoretical models. Specifically, by estimation of the free diffusion coefficients from in vitro studies and of apparent diffusion coefficients from in vivo experiments, we have assessed tortuosity in the right striatum of 9 Sprague Dawley rats through a model correctly describing both vascular permeability as a function of time and diffusion processes occurring in the brain tissue. This model takes into account acoustic pressure, particle size, blood pharmacokinetics, and diffusion rates. Our model is able to fully predict the result of a FUS-induced BBB opening experiment at long space and time scales. Recovered values of tortuosity are in agreement with the literature and demonstrate that our improved model allows us to assess that the chosen permeabilization protocol preserves the integrity of the brain tissue.

1. Introduction

The in vivo characterization of gadolinium-(Gd-) based MRI contrast agent (MR-CA) diffusion within the brain tissue is of great interest for the understanding of drug transport mechanisms in the brain parenchyma, in the framework of the recent pharmaceutical developments targeting entral nervous system (CNS) diseases. Despite increasing efforts and encouraging results, drug delivery to the CNS remains a challenging task. Indeed, the blood-brain barrier (BBB) not only prevents neurotoxic substances from entering the brain but also limits the passage of therapeutic products to the CNS [1, 2]. Many strategies have been studied to overcome this obstacle, including direct injections [3, 4], transient BBB disruption using chemical agents [5, 6], or molecular engineering [7]. More recently, a promising technique has been proposed, allowing the delivery of various compounds to the brain using low-intensity focused ultrasound combined with circulating microbubbles [8].

However, once the molecules have crossed the barrier, they have to diffuse in a highly constrained media, the extracellular space (ECS), to reach their targets [9]. Moreover, since the ECS architecture can change in case of pathologies [10, 11], the characterization of the hindrance experienced by molecules within the brain tissue is essential when designing new therapeutic compounds or diagnostic molecules for brain diseases. Diffusion constraints can be studied by estimating the ECS tortuosity (λ). This parameter compares the apparent diffusion coefficient (ADC) of a molecule within the complex architecture of the ECS to its diffusion coefficient in a free medium Dfree [12]. Different strategies have been proposed to measure the ADC. The most widely used method is real-time iontophoresis [9, 13], using tetramethylammonium (TMA+) as a probe. This technique not only permits the in vivo characterization of the ADC but also, thanks to the small size of both detection electrodes and injection micropipette, proves minimally invasive, with consequent preservation of the integrity of the tissues. Its main drawback consists in the measurement relying on just one spatial point. More recently, diffusion-weighted magnetic resonance imaging (DW-MRI) has been proposed to noninvasively measure the ADC of water molecules in the brain [14, 15]. In comparison to the previous techniques, DW-MRI allows ADC measurements in deeper areas of the brain with a high (typically 2 mm isotropic) spatial resolution [16]. However, contrary to TMA+ and other techniques using labelled molecules that diffuse only across the ECS, DW-MRI detects water, which is also present in the intracellular compartment. To benefit from the advantages offered by MR in acquiring deep volumes of the brain, a new method has been recently introduced by our team, which allows us to detect molecular diffusion only in the ECS structure [17]. To do so, MR-CAs are directly injected into the brain tissue, and their diffusion is followed by acquisition of several longitudinal relaxation-time (T1) parametric maps. MR-CA concentration maps at different diffusion times are then calculated, and from these, the ADC is estimated. When compared to the typical diffusion-based MRI techniques, our method investigates larger areas of the brain with a higher spatial resolution (about 0.2 × 0.2 mm2 in plane and 1 mm in thickness). However, a major issue raised by this procedure consists in intracerebral injections inducing edema, which modifies the diffusion properties of brain tissue.

In the present study, we have used two noninvasive methods for the in vivo estimation of the ADC of different Gd chelates diffusing in the ECS after a FUS-induced BBB opening experiment. In both cases, contrast agent diffusion is recorded through dynamic acquisitions of MRI concentration maps. In the first method, the ADC evaluation is performed as in [17], e.g., by fitting a 2D Gaussian curve to the image intensity at different time points. However, diffusion of molecules delivered to the brain with the aid of FUS-induced BBB permeabilization depends on many factors, such as tissue and particle properties, as well as acoustic parameters. For this reason, as a second approach to estimate contrast agent diffusion, we introduce here the first diffusion model able to fully describe and predict at long space and time scales the result of a FUS-induced BBB opening experiment. This model takes into account acoustic pressure, particle size, blood pharmacokinetics, vascular permeability as a function of time, and diffusion process occurring in brain tissue.

Starting from ADC estimation performed with the help of both methods and the evaluation of Dfree for all the compounds by means of in vitro experiments, it is possible to calculate tortuosities in the target region of rats’ brains, to evaluate the effect of the selected BBB permeabilization protocol on the properties of brain tissue.

2. Materials and Methods

2.1. Experimental Procedures

All magnetic resonance acquisitions were performed by using a 7 T/90 mm Pharmascan scanner (Bruker, Ettlingen, Germany). The in vitro acquisitions have been performed by using a 1H transmit-receiver volume coil (Bruker). The in vivo experiments have been conducted by using a dedicated ultrasound single-loop radiofrequency coil [18], whose diameter was wide enough for the ultrasound beam to pass through it and for extensive displacement of the transducer above the rat’s head. A heater device was used to keep temperature at the physiological value (37°C), as monitored by a temperature probe that was inserted inside the magnet (see Figure 1).

Three different gadolinium (Gd) chelates were studied: Dotarem® (Gd-DOTA, Guerbet, France), Gadovist® (Gd-DO3A-butrol, Bayer, Germany), and MultiHance® (Gd-BOPTA, Bracco, Italy). First, we assessed their longitudinal relaxivities (r1) at 7 T and 37°C, using phantoms made of bundles of tubes containing different contrast agent (CA) concentrations in 0.3% w/w agarose gel for each compound. For these phantoms, T1 values were measured by means of an inversion-recovery fast gradient echo (IR-FGE) sequence [17, 19] (echo time (TE)/repetition time (TR1) = 2.5/5 ms, 6 segments, 90 inversion times (TI) varying from 75 ms to 8975 ms, flip angle (FA) = 5°, matrix = 120 × 120 × 5 with resolution = 0.250 × 0.250 × 1.25 mm3, delay between the acquisitions of two segments (TR2) = 15000 ms, and number of averages (NA = 6)). Resulting relaxivities are summarized in Table 1. This table also includes the hydrodynamic diameter (dH) of each CA, measured by dynamic light scattering (DLS). DLS experiments were performed using a NanoZS equipement (Malvern, France) operating at an angle of 173°. For each Gd chelate, the DLS acquisitions were performed at 25°C by using concentrations of 0.5 M for MultiHance and Dotarem and of 1.0 M for Gadovist, e.g., without diluting the samples. We performed five different DLS measurements for each sample. The mean dH and the standard deviation evaluated over the five measures are reported in Table 1.

Evaluations of the Dfree of each compound were done with a method already presented in a previous work [17]. For each product, 10 μL of a 5 mM solution was injected with a Hamilton syringe (diameter = 1 mm) into a tube filled with 0.3% w/w agarose gel. A stereotactic system was used to make the injection central and vertical with respect to the tube. The free diffusion of the CA was then dynamically followed by acquisition of five T1 parametric maps after injection (IR-FGE sequence with the following parameters: TE/TR1 = 2.5/5 ms, 6 segments, 60 TI from 88 ms to 5100 ms, FA = 5°, matrix = 128 × 104 × 14 with res = 0.225 × 0.225 × 1 mm3, TR2 = 9000 ms, NA = 1, and total duration = 12.5 min). A T1 map acquired before the injection was used as a reference.

The number of TI values has been chosen to ensure an accurate estimation of T1 values for a large range of T1. In particular, thanks to this sequence, we are able to detect Gd concentrations with a sensitivity threshold estimated around 2.5 μM [17]. The spatial and temporal resolutions of this mapping sequence were set in order to ensure a sufficient space and time sampling of CA diffusion process. Further details about the optimization of this MRI sequence can be found in [17]. In all T1-parametric maps, all voxels with a T1 value larger than 5000 ms, that is, much larger than both the T1 of gray and white matter at 7 T [20], have been masked and considered as Not-a-Number.

The measurements of the ADC were performed in vivo on 9 Sprague Dawley male rats (3 rats/compound, 120–140 g, Janvier, Le Genest-Saint-Isle, France). Animal testing complied with the recommendations of the European Community (86/609/EEC) and French legislation (decree no. 87/848). The experimental setup is shown in Figure 1. The rats were anesthetized by means of 1.5–2% isoflurane in a mixture of air and oxygen, and their heads were chemically shaved to ensure a proper coupling with the ultrasound transducer. They were then placed in prone position in a cradle, integrating a stereotactic frame and a dedicated radiofrequency coil (Figure 1(b)). A custom build catheter (25G) was inserted into the caudal vein to perform injections from outside the MRI scanner. Temperature monitoring and breathing monitoring were performed using a rectal temperature probe and a respiration probe (Figure 1(c)), respectively.

A MR-compatible focalized transducer with 1.5 MHz central frequency (diameter: 25 mm, focal depth: 20 mm, focal spot dimensions: 1.1 mm in-plane, 6 mm thickness, Imasonic, France) was coupled to any animal’s head via a balloon filled with degassed water. The transducer was mounted on a mobile stage, and its position could be tuned from outside the magnet by using MR-compatible motors (see Figure 1(b)). The movement of the motors and ultrasound parameters were controlled by a dedicated software (Thermoguide®, Image Guided Therapy, France) (Figure 1(a)). All acoustic pressures were estimated from previous calibration of the transducer, taking a skull transmission factor varying with animals’ weight [21].

In Figure 2, the experimental protocol is shown. After rat installation, an acoustic radiation force imaging (ARFI) sequence [22, 23] was performed to localize the ultrasound focal point in rats’ brains, consisting in a standard multislice multiecho sequence (MSME; TE/TR = 28/1080 ms, matrix = 64 × 64 × 5, and res = 0.5 × 0.5 × 2 mm3) modified by the addition of two motion-sensitizing gradients (MSGs; duration of one MSG = 8 ms and duration of the ultrasound bursts = 4 ms). Knowing the current position of the focal spot, the transducer was moved using the motors so as to focalize ultrasound in the left striatum of the rats. This location has been chosen to ensure a high acoustic transmission through the skull as detailed in a recent work published by our team [21]. A second ARFI image was acquired to assess the good positioning of the ultrasound focal spot. T1-weighted (T1w) anatomical images were acquired, before the BBB opening, by using an MSME (TE/TR = 8.3/300 ms, matrix dimension = 256 × 256 × 10, resolution of 0.125 × 0.125 × 1 mm3, and 3 averages). This was followed by a bolus injection of Sonovue® microbubbles (Bracco, Milan, Italy; 1.5 × 108 bubbles/mL, 1.6 mL/kg, 3 s) via tail vein catheter, approximately 5 s before transcranial sonication (3 ms burst every 100 ms over a period of one minute; estimated focal acoustic pressure in the brain = 0.8 MPa). 30 seconds after the end of the ultrasound session, Gd chelates were intravenously injected via bolus (5 seconds, 0.5 M and 1.6 mL/kg for MultiHance and Dotarem; 1 M and 0.8 mL/kg for Gadovist). T1-weighted (T1w) images were acquired 30 seconds after the CA injection to verify the BBB disruption. Using the same IR-FGE sequence as the one used for in vitro diffusion, T1 parametric maps were acquired before and after sonication in order to dynamically follow the diffusion of the Gd chelates in the brain. At the end of each experimental session, a T2-weighted (T2w) image was acquired to verify the absence of any hemorrhage or edema due to ultrasound-induced BBB disruption. A rapid acquisition with relaxation enhancement (RARE) sequence was used with the following parameters: TE/TR = 10/3800 ms, RARE factor = 8, and matrix = 128 × 128 × 32 with resolution = 0.225 × 0.225 × 0.5 mm3.

2.2. Data Analysis

From T1 maps, the corresponding concentration maps were calculated using the following relationship between the longitudinal relaxation rates, 1/T1, and the Gd-chelate concentrations, [CA] [24]:where 1/T1,0 is the relaxation rate of the sample without CA, i.e., before the injection. From this equation, CA concentration maps were then obtained. All voxel values in the T1 or T1,0 maps larger than 5000 ms were considered as Not-a-Number in the CA maps. These voxels were not considered in the CA-diffusion analysis.

In all cases (both for in vivo and in vitro acquisitions), we have assigned to each CA map the time elapsed between the CA injection (in agarose gel or in the caudal vein for the in vitro and the in vivo acquisitions, respectively) and the beginning of the CA-map acquisition sequence.

To evaluate the Dfree value of injected molecules, the following bidimensional Gaussian function was fitted to concentration-map data for each time point after the CA injection [17, 25, 26]:where A is the Gaussian amplitude and (x0, y0) are the coordinates of its center along the absolute axes (x, y). a, b, and c are functions depending on the Gaussian widths (σX and σY) along its main axes (X and Y) and on the angle θ between (X, Y) and (x, y):

The regression algorithm used to fit the data with Gaussian functions is the Levenberg–Marquardt algorithm [27], available in the GSL, GNU Scientific Library (https://www.gnu.org/software/gsl/doc/html/nls.html). In particular, we used the version of this algorithm implemented in the scaled LMDER routine in MINPACK, written by Jorge J. More, Burton S. Garbow, and Kenneth E. Hillstrom (https://people.sc.fsu.edu/∼jburkardt/f_src/minpack/minpack.html).

Defining and as the molecular mean square displacements along X and Y, the diffusion coefficients along these axes, Dfree,X and Dfree,Y, are given by Fick’s law:where t is the instant time after injection, i.e., the diffusion time. Dfree values were then calculated as the mean value of Dfree,X and Dfree,Y components:

The first method used to evaluate the ADC consisted in placing a mask surrounding the disruption site in CA maps, to which the same Gaussian fitting procedure was applied. For any compound, the ADC was estimated in any rat’s striatum as the average:

The second ADC estimation took into account how the BBB permeabilization changes after the ultrasound application, together with CA pharmacokinetics after injection. A homemade MATLAB code was used to simulate CA diffusion within the ECS after the BBB opening.

The code comprises the following components:(i)A source function S (x, y, z, t), describing the contrast agents that move from the blood to the brain, was modeled aswhere α is a proportionality constant, requiring a first guess on its value, QCA(x, y, z, t) is the amount of CA crossing the BBB [28], and CAblood(t) describes CA pharmacokinetics. For a Gd chelate of hydrodynamic diameter (dH), QCA(x, y, z, t) is defined as [28, 29]where σ0 is the standard deviation of the distribution of the gap sizes generated in the BBB by ultrasound and k is the BBB closure rate (k = 1.54e−5·s−1). Since it has been demonstrated that blood-brain barrier disruption is characterized by a mechanical index (MI), which is linearly dependent on the effective acoustic pressure (Pex) [30], we considered the same dependence for σ0. In particular, according to the work published by Marty et al. in 2013 [28], we applied the relationship σ0 = 2.1·Pex. Starting from the simulated acoustic pressure map, we obtained the σ0(x, y, z) distribution.The kinetic term in equation (7) can be expressed bysince our time resolution in the acquisition of CA maps (12.5 min) allows for just considering the wash out of CAs, in obedience to Tofts’ two-compartment kinetic model [31]. CAblood(t) depends on the injected CA concentration (CAinj) and on its clearance rate from the blood, b. CAinj was estimated for each animal by taking into account its weight and an average blood volume of 6.86 mL/100 g [32], while b was fixed at 25 minutes [33].(ii)Introducing the source term S(x, y, z, t) into Fick’s second law, the evolution of CA-concentration long time was found by simulating the equationfor a temporal and spatial resolution higher than those characterizing [CA] maps. Specifically, an isotropic spatial resolution (dx, dy, dz) equal to 0.125 μm was selected, while the temporal step dt was set at 5 s. While α has been guessed, the initial ADC values used for simulations were chosen from the equationstarting from tortuosity values of the target region of the brain recovered from the literature [34] and molecular Dfree values retrieved from our experiments.(iii)The simulated CA volume was downsampled in space and time to the MRI acquisition resolution and then coregistered to the experimental three-dimensional [CA] distribution in the CA-concentration maps.

Due the large focal-spot length (∼6 mm), CA concentration can be considered as constant along this direction (called z) for all the slices taken in account. This makes the CA gradient negligible along z, as well as the related diffusional process (see Figures S1 and S2 in the Supplementary Materials). For this reason, the previous equation can be considered as reasonably describing the following bidimensional dynamics:

This last equation was integrated in order to estimate [CA](x, y, t). Through a cumulative fit including the experimental CA maps for the central slice, the ADC components along x and y and the proportionality constant α were found. Different ADCs around the value suggested by equation (11) were simulated until the fit algorithm converged.

To evaluate the quality of the experimental approach chosen to mimic molecular free diffusion (i.e., the injection of the compound in 0.3% w/w of agarose gel), it is worth estimating the hydrodynamic diameter of the molecules, using the Stokes–Einstein equation:where k = 1.38·10−23 Pa·m3·K−1 is the Boltzmann constant, T is the temperature in Kelvin degrees, and η is the viscosity of the agar gel (6.92·10−4 Pa·s).

From the mean ADC recovered through the two aforementioned methods, the tortuosity values were estimated with the help of equation (11).

3. Results

Figure 3 shows an example of in vitro diffusion data and their analysis. Concentration maps (Figure 3(a)) were acquired 4 to 56 minutes after the injection of MultiHance. These data were fitted by means of the bidimensional Gaussian function reported in equation (2). The simulated Gaussian distributions resulting from the fit are shown in Figure 3(b). Taking into account the voxel values in the central row of the Gaussian spots pictured in Figures 3(a) and 3(b), it is possible to assess the quality of the fit, as illustrated in Figure 3(c), where the black dots represent the data and the red curve their Gaussian fit.

Fick’s law (equation (4)) was used to fit the squares of the fitted Gaussian widths (σx and σy) as a function of the diffusion time, in order to obtain an estimation of Dfree,X and Dfree,Y (Figure 3(d)). The Dfree values found for each compound are given by the average of the two components and are summarized in Table 1.

The ADCs were estimated by analyzing in vivo concentration maps, as the ones shown in the upper panel of Figure 4. Specifically, these maps were acquired 2 to 84 minutes after bolus injection of Dotarem. Prior to compute Gaussian fits on concentration maps, a mask including only the BBB disruption site was applied (Figure 4(b)). The first method for ADC evaluation consists in fitting 2D Gaussian functions to such maps. The resulting distributions are shown in Figure 4(c).

As for the in vitro measurements, the overlapping between data and fit curve is shown (see Figure 4(d)). By comparing, through a two-sample Kolmogorov–Smirnov test, the data shown in Figure 4(c) with the respective Gaussian profiles at each time point, we obtained values equal to 5.6e − 4; 0.258, 0.258, 0.440, and 0.2581, meaning that only at the first time point the Gaussian fit results to be different from the data. We also evaluated the ADC values without taking into account the first time point. However, since the values obtained with and without the first time point varied less than the error estimated by the respective linear fits and less than the variations inside the n = 3 rat pools, we also considered the first time point to estimate the ADCs.

The temporal evolution of the squared Gaussian widths is shown in Figure 4(e), together with their fits by Fick’s Law. Starting from ADCX and ADCY values, the ADC in each rat’s striatum was found. By average over the entire set of rats, the mean ADCs, reported in Table 2, were estimated, as well as brain tortuosity, λI.

The second method proposed to evaluate brain diffusional properties is based on a model taking in account both the temporal changes in BBB permeabilization after ultrasound application and CA blood pharmacokinetics.

Figure 5 shows an example of CA distributions inside the brain, obtained by fitting this model to experimental concentration maps obtained by diffusion measurements on Multihance.

Once again, Figure 5(a) reports the masked concentration maps used to evaluate brain tortuosity, while the maps in Figure 5(b) are obtained via model. The ADCs estimated by average of model results obtained for each compound are shown in Table 2 (ADCII). In the same table, the values obtained for the proportionality constant α of the source term are included. Entering Dfree and ADCII values found by this second approach in equation (11), brain tortuosity is once again retrieved (λII in Table 2).

For the sake of comparison, in Figure 6, the distribution profiles extracted from the centers of [CA] maps are shown, as previously done in Figure 4(d). This dataset refers to an experiment on Gadovist with black dots representing experimental data and [CA] red and blue profiles representing theoretical data obtained from the first and second method, respectively. By comparing, through a two-sample Kolmogorov–Smirnov test, the data with the simulated and the Gaussian profiles, we obtained, at different time points, values equal to 0.985, 0.374, 0.147, 0.047, 0.147, and 0.047 for method I and equal to 0.675, 0.675, 0.736, 0.736, 0.736, and 0.736 for method II.

These results show that method II allows for obtaining distribution shapes that are more similar to data at all the time points than Gaussian fits in method I.

4. Discussion

This work introduces two new methods suitable for the in vivo characterization of molecular diffusion processes taking place in the ECS after transient BBB permeabilization with low-intensity focused ultrasound in order to deliver MR-contrast agents to the brain. We used MRI to record MR-CA diffusion. By measuring DFree (free-medium diffusion) and ADC values within the ECS, brain tissue tortuosity was calculated in order to have information on brain architecture.

To assess the quality of the experimental approach chosen to evaluate molecular free diffusion, it is worth comparing the hydrodynamic diameter of the molecules, dH(S-E), obtained through equation (13), to the ones found by using DLS. As can be noticed from Table 1, the hydrodynamic diameter found through these two methods agree, which means that the diffusion of the compounds in 0.3% w/w of agarose gel can be considered as free. In addition, Dfree values in Table 1 can be compared to the analogous ones already published in the literature. Specifically, Marty et al. [17] have found the same Dfree for Dotarem, whereas Thorne and Nicholson [35] have estimated a free diffusion coefficient equal to (2.22 ± 0.16)·10−10 m2/s for a molecule with hydrodynamic diameter of 2.95 ± 0.02 nm, which is comparable to one that was found for a slightly smaller molecule of MultiHance (dH = 2.3 ± 0.1 nm and Dfree = (2.8 ± 0.2) · 10−10 m2/s).

Table 2 shows that, irrespective of the applied method, ADC values scale correctly with molecular size, decreasing at increasing dH (ADCDotarem > ADCGadovist > ADCMultiHance), as expected from comparison to the literature [35]. Furthermore, all ADC values are smaller than their associated Dfree, which confirms the hindrance experienced by diffusion across the ECS.

Tortuosities obtained by method I and II (λI and λII) are compared to those appearing in the literature in order to assess the goodness of ADC estimation.

λI and λII obtained for the different molecules turn out constant, which agrees with the literature. Indeed, all of our test molecules have a hydrodynamic diameter ten times smaller than the intracellular gap d, which is typically comprised between 20 and 64 nm in healthy rats’ brains [35, 36].

In this case, the stationary wall-drag effect, expected for larger molecules by virtue of viscosity theory, affects neither molecular diffusion [36] nor tortuosity, whose value only depends on the ECS structure and not on the size of the diffusion probes.

4.1. Limitations and Future Perspectives

In the present work, both the methods used to estimate the molecular apparent diffusion coefficients are based on a protocol validated by our team in 2013 [17], e.g., the dynamic acquisitions of CA-concentration maps through an IR-FGE MRI sequence. Although this sequence has been accurately tuned to be sensitive to a large range of CA concentrations and to have a sufficiently high temporal and spatial resolution to record molecular diffusion, further work is needed to improve such resolutions. For example, a suitable way to increase the speed of the MRI sequence currently used is by using compressed sensing MRI techniques [37]. Doing so, we expect to reduce the acquisition time and therefore to get access to diffusion data of MRI contrast agents at high temporal resolution.

The second limitation of our experimental approach is related to the possibility to evaluate CA diffusion only in two dimensions. Indeed, our method allows us to estimate the transversal components (x and y) of the ADC but not to evaluate CA diffusion processes along z-axis. This is due to the gradient concentration and to the relatively low spatial resolution in this direction. In order to improve our delivery method and to be more sensitive to Gd concentration gradients along z-axis, future experiments can be performed by using multielement transducer to produce a controlled steering of the ultrasound beam in the z direction (see Figure S3 in the Supplementary Materials). With this steering approach, it will be possible to permeabilize the BBB in a smaller region of the brain. In addition, by improving the spatial resolution in z of the concentration maps, it will be then possible to characterize the particle diffusion also along this direction.

Another limitation of our work concerns the capability of method II to fully predict the amount of particles getting in the brain after a FUS-induced BBB opening experiment. Indeed, from a qualitative point of view, one can expect the inclusion of the source term to provide a better data description when the blood-to-ECS flux is larger, i.e., for CAs of smaller size, since the QCA expression is a monotonically decreasing function with the molecular hydrodynamic diameter, dH. However, the amount of particles getting in the brain after a FUS-induced BBB permeabilization is dependent from many factors, some of them being difficult to precisely control. For example, if the coupling of the water balloon between the transducer and the head, or if the position of the transducer, slightly changes between two experiments, the transmitted acoustic power could vary inside the brain and, consequently, the amounts of particles delivered to brain tissue [21, 38]. McDannolds et al. [39] have recently shown that even the level of oxygen used as a carrier gas for anesthesia during the experiments can change microbubble activity and BBB disruption. All these aspects, varying among experiments, change the value of the constant of proportionality α. For this reason, in order to use our model to simulate an experimental outcome, the simulations need to be performed by varying α between 0 (e.g., the worst-case scenario corresponding to a failure of the experiment) and 0.07 (e.g., the maximum value of α found in this work).

4.2. Comparison between the Two ADC Estimation Methods

The first method consists in fitting Gaussian distributions to CA-map data in the brain region where diffusion occurs. From this fit, the molecular square displacements, and so their ADC, can be evaluated. This kind of postprocessing is already accepted in the literature [16], although originally applied to CA diffusion patterns acquired after intracerebral injection of compounds. However, this method presents some limitations. The first one concerns the application of this fit to CA maps with low signal-to-noise ratio (SNR).

In particular, we define the SNR in each slice of the CA maps, as the ratio between the maximum CA delivered in the slice and the standard deviation in a region (20 voxels × 20 voxels) located in the contralateral hemisphere. The Gaussian fit overestimates the distribution widths for SNR smaller than 10. This is the case, for example, of the acquisition shown in Figure 6. The errors committed by method I on the estimation of the distributions widths are confirmed by the values obtained when comparing the Gaussian profiles to the respective data points, through a two-sample Kolmogorov–Smirnov test. Indeed, the values resulted to be smaller than 0.05 at two time points. The same issue does not affect ADC estimations when the compounds are intracerebral injected, as in [17]. Indeed, in this latter case, the SNR is higher than the one obtained through BBB opening since the CA concentration diffusing within the ECS is 100 times larger than the CA delivered through BBB-opening.

On the other hand, when method II is applied to analyze the same dataset, it is possible to obtain particle distributions more similar to the experimental ones, as confirmed by the values larger than 0.05, resulting from the same kind of statistical test.

In addition, to fit the data through the first method, we use the version of Levenberg–Marquardt algorithm implemented in the scaled LMDER routine in MINPACK [27]. This scaled LMDER routine makes use of both the function and its derivative, so it could explain why in some cases, as the one shown in Figure 6, the main differences between the data and the respective Gaussian fit can be found near the peak.

With respect to the first method, the second ADC estimation method presented in this work is based on a diffusion model that includes a source term. The source term describes the flux from the blood to the ECS only, which is appropriate if the two pools have a large concentration difference. This approximation can be quantitatively justified. Indeed, the CA concentration injected in the blood system is around 3 mM, while, as can be noticed from Figures 46, the maximum CA delivered in the brain is estimated to be approximately 100 times smaller. In addition, the CA concentration in blood is much higher than the ECS concentration, during the duration of whole of the experiments (about 1 hour) (see Figure 4 in Supplementary Materials).

Another possible way to compare the two methods is to compare the different tortuosity values, λI and λII, shown in Table 2. It has been recently shown with histology that low-intensity pulsed ultrasound could be used to transiently enlarge the ECS width [40]. In particular, by estimating the overall volume of distribution of different nanoparticles, Frenkel et al. found an enhanced volume of 36% in average. The volume where particles diffuse in ECS is characterized by the volume fraction υ = VECS/VT, defined as the ratio between the volume of ECS (VECS) and the volume of the whole tissue measured in a small region of the brain (VT) (Sykova, Physiol Rev. 2008). In healthy brain tissue, the ECS volume fraction υ is estimated around 0.20. However, by considering the study proposed by Frenkel et al. [40], the volume fraction enlarges of 36% after FUS application, leading to a volume fraction of υ = 0.27. Since the relationship between the tortuosity value, λ, and υ is the following as given by [41]:and the expected value of brain tortuosity after a FUS-induced BBB permeabilization experiment is equal to 1.32, e.g., more similar to the values obtained through method II than the ones estimated through the Gaussian fit.

5. Conclusions

In this study, we used two methods to characterize the contrast agent bidimensional diffusion within the brain after ultrasound-induced BBB opening. These techniques allow to investigate macromolecules biodistribution within the ECS with a slow time scale, suitable for the study of cellular uptake and transport, as well as of the potential clearance processes related to bulk flow or glymphatic pathway. Although it is well known that focused ultrasound combined with microbubbles permits to transiently and noninvasively break tight junctions, locally increasing the BBB permeabilization and so promoting drug delivery into the brain [8, 28, 4244], so far no study has been performed to fully characterize, on a macroscopic space and time scale, the distribution of a compound when it enters the brain.

By using a motorized and MR-compatible ultrasound system, we were able to target the right striatum of 9 rats in a very precise and reproducible manner, in order to study diffusion processes in a specific area of the brain. Three commercially available MR-CAs were tested (Dotarem®, Gd-DOTA; Gadovist®, Gd-DO3A-butrol; MultiHance®, Gd-BOPTA). Their diffusion from the BBB-disruption site was followed by acquisition of several CA maps within 1 hour from application of ultrasound. The tested compounds are characterized by a similar hydrodynamic diameter (about 1–2 nm), which resulted in a similar hindering of diffusion in the ECS. Since the CA distribution depends on the diffusion properties of brain tissue, we have evaluated its tortuosity, a parameter comparing molecular ADC inside the tissue to its free-diffusion counterpart in a media without obstacles. The methods proposed here to estimate λ are both based on data processing of MR-CA maps. The first approach does not describe the dependence of molecular diffusion neither on fundamental biological aspects nor on the specific protocol used to permeabilize the BBB.

For this reason, we have presented a mathematical model able to fully predict time evolution of CA distributions within the brain after BBB permeabilization induced by FUS. Our model takes into account different biological features concerning the BBB-opening mechanism, such as the gap distribution between endothelial cells, in turn depending on the effective acoustic pressure transmitted through the skull and the shape of the focal spot, the BBB closure rate, and the CA concentration in blood after bolus injection and its physiological rate of clearance. The match with the experimental data allows us to introduce this approach as a new tool to successfully predict and plan drug distribution after a BBB-opening experiment, for any particle size and acoustic parameter, in all brain regions.

Abbreviations

BBB:Blood-brain barrier
US:Ultrasound
FUS:Focused ultrasound
CA:Contrast agents
CA map:CA concentration map
ADC:Apparent diffusion coefficient
ECS:Extracellular space.

Data Availability

The MRI data used to support the findings of this study are available from the corresponding author upon request.

Disclosure

Earlier results of the present work have been presented at IEEE International Ultrasonics Symposium (IUS) in 2016 [26] and at the conference NeWS in 2017.

Conflicts of Interest

The authors declare no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Authors’ Contributions

AC, SM, and BL planned the experiments. AC and RM performed MRI acquisitions and FUS experiments. AC performed data analysis. AC wrote the simulation code. SM provided MRI data analysis pipelines. NT performed DLS experiments. SM, BL, and SDP managed the overall project.

Acknowledgments

AC received an Enhanced Eurotalents postdoctoral fellowship (Grant Agreement no. 600382), part of the Marie Sklodowska-Curie Actions Program, cofunded by the European Commission and managed by the French Atomic Energy and Alternative Energies Commission (CEA).

Supplementary Materials

Supplementary Figure 1: (A) acoustic pressure field at 1.5 MHz. The pressure field is normalized by the maximum pressure (obtained at the focal spot). The acoustic pressure field shown on the right has been rescaled to the spatial resolution of the concentration maps. (B) Axial and sagittal views of the Gd-concentration maps: it can be noticed that the focal spot dimension along z-axis is comparable to the thickness of the rat brain in the area where the BBB was permeabilized. Moreover, the spatial resolution in this direction is much lower than the in-plane resolution (around 4.4 times). This significantly lower resolution along z-axis does not allow to precisely quantify the variations of CA concentration in this direction. Supplementary Figure 2: left: sagittal views of the Gd-concentration maps; right: concentration profiles extrapolated from the center of the BBB opening site. These profiles show that the Gd concentration changes only slightly with the z position. This is due to the low resolution of the Gd-concentration map along z and also to the dimensions of the focal spot along this direction. Supplementary Figure 3: acoustic pressure fields at 1.5 MHz for the concave transducer (F/D = 0.8) without steering (left) and with a 2.5 mm steering toward the transducer (right). Both pressure fields are normalized by the maximum pressure obtained at the focal spot in case of the absence of steering. With a 2.5 mm steering toward the transducer, the volume of the focal spot is decreased by 20% and the maximum pressure is increased by 10% compared to the experiment without steering. Supplementary Figure 4: comparison between the CA concentration in blood (picture in blue) and the maximum CA delivered during a BBB opening experiment (in red). In particular, experimental points refer to the data shown in Figure 4, while the trend of CAblood along time, t, has been derived through the equation CAblood(t) = CAinj·exp(−t/b), with b = 25 minutes (Aime and Caravan, JMRI, 2009). (Supplementary Materials)