Ex vivo optimisation of a heterogeneous speed of sound model of the human skull for non-invasive transcranial focused ultrasound at 1 MHz

Abstract

Transcranial brain therapy has recently emerged as a non-invasive strategy for the treatment of various neurological diseases, such as essential tremor or neurogenic pain. However, treatments require millimetre-scale accuracy. The use of high frequencies (typically ≥1 MHz) decreases the ultrasonic wavelength to the millimetre scale, thereby increasing the clinical accuracy and lowering the probability of cavitation, which improves the safety of the technique compared with the use of low-frequency devices that operate at 220 kHz. Nevertheless, the skull produces greater distortions of high-frequency waves relative to low-frequency waves. High-frequency waves require high-performance adaptive focusing techniques, based on modelling the wave propagation through the skull.

This study sought to optimise the acoustical modelling of the skull based on computed tomography (CT) for a 1 MHz clinical brain therapy system.

The best model tested in this article corresponded to a maximum speed of sound of 4000 m.s⁻¹ in the skull bone, and it restored 86% of the optimal pressure amplitude on average in a collection of six human skulls. Compared with uncorrected focusing, the optimised non-invasive correction led to an average increase of 99% in the maximum pressure amplitude around the target and an average decrease of 48% in the distance between the peak pressure and the selected target. The attenuation through the skulls was also assessed within the bandwidth of the transducers, and it was found to vary in the range of 10 ± 3 dB at 800 kHz and 16 ± 3 dB at 1.3 MHz.

Keywords:

• High intensity focused ultrasound
• ultrasound
• thermal ablation
• modelling

Introduction

Transcranial focusing of therapeutic ultrasonic waves remains a challenge in the field of medical ultrasound because heterogeneities in the speed of sound through biological tissue and skull bone induce distortions within the ultrasonic wave field [1]. Such distortions result in the partial destruction of the focusing pattern. To restore focusing quality, aberration corrections are performed via ultrasonic arrays, with different signals estimated and emitted on each element of the array [2].

Several aberration correction techniques have been developed over the past two decades, and they are reviewed and discussed in Ref [3]. Time reversal correction (or phase conjugation at one frequency) [4] has been particularly adapted to high intensity focused ultrasound (HIFU) aberration correction because it maximises the pressure at the focus for a specific transmitted acoustic energy [5]. In certain situations, when the ultrasonic array incorporates both transmitting and receiving channels, the echoes of a bright reflector [6] or a point-like transducer located inside the biological tissue may be used [7] because the signals received on the array can be phase conjugated to focus back to the initial position of the source. To preserve the non-invasiveness of the focused ultrasound, a variant of the technique was implemented in vitro for transcranial aberration correction using a cavitation bubble generated either by low negative pressure [8] or by injected liquid droplet vaporisation [9]. However, the in vivo cavitation threshold in the brain is too high for the clinical implementation of bubble-generating techniques [10,11].

The only non-invasive option consists of simulating the acoustic wave propagation through the skull. Several simulation schemes have been developed [3] and can be ranked by increasing numerical complexity as follows: ray tracing, angular spectrum methods and full wave models. All of these models rely on image-based estimations of the skull density and the speed of sound.

As an initial approximation, the skull can be modelled as a homogeneous medium with varying thicknesses. This method was initially based on magnetic resonance imaging [12]; however, computed tomography (CT) provides more information on the bone structure, and 59% of the optimal amplitude (compared with a hydrophone-based correction) was restored at the focus by calculating the thickness-averaged speed of sound from the CT skull density [13] calculated at 750 kHz. The model was improved by considering three homogeneous layers in the skull (2900 m/s in the inner and outer tables and 2500 m/s in the central layer), allowing for 79% of the optimal amplitude to be restored when investigated at 500 kHz [14]. Further models of the speed of sound in the skull were estimated using transcranial simulations based on CT data coupled with a genetic algorithm to fit phase measurements [15], and these simulations have been used in recent studies based on the normal incidence between the ultrasound beam and the skull surface [16].

Refined models of the speed of sound based on the 3D external and internal structures of the skull as measured by clinical CT have been introduced [17,18]. The corresponding non-invasive correction was shown to restore 90% of the optimal amplitude (compared with 32% without any correction) based on previous testing with a 300-element array on one human skull at 1 MHz with a transcranial speed of sound between 1500 and 3100 m/s [19]. Pichardo et al. estimated the speed of sound as a function of the skull density for a discrete set of frequencies (270, 836, 1402, 1965 and 2526 kHz) [16]. These authors used a two-step optimisation process based on a total of 33 measurement locations on seven skulls (3–5 locations on each skull) that included the coarse optimisation of a large set of random functions using a multi-layer propagation model followed by the fine optimisation using a finite-difference time-domain solution of the Westervelt equation.

In this article, a speed of sound model for a 3D finite difference time domain fluid simulation was optimised and tested at 1 MHz based on six different human skulls, which included more than 3000 measurement points. Additionally, the attenuation through the skulls was assessed between 800 kHz and 1.3 MHz within the bandwidth of the transducers. The 1 MHz frequency was selected based on the objective of optimising a high-frequency non-invasive clinical device for brain therapy [20].

Methods

Six human skulls were provided by the Institut d’Anatomie UFR Biomédicale des Saints-Pères Université Paris Descartes, Paris. This study was approved by the ethics committee of the Centre du Don des Corps, Université Paris Descartes, Paris. Donors provided informed consent before death. The average age of the donors was 88 years.

To assess the accuracy of the numerical model of the skull, the phase aberration of the skull was numerically estimated, the corresponding aberration correction was experimentally applied, and the pressure amplitude around the geometrical centre was measured. Five 3D models were tested for the speed of sound in the skull. The corresponding phase aberration corrections were compared using a unique experimental propagation matrix per skull.

The sequential steps of the method used for each skull are as follows:

1. Measure the propagation matrix around the geometrical centre of the multi-element array in water only without the skull;

2. Place the skull on a stereotactic frame (Elekta, Stockholm, and Sweden);

3. Degas the skull on the stereotactic frame in the water;

4. Place the skull on the stereotactic frame in the water tank in front of the multi-element array (Figure 1);

5. Measure the propagation matrix around the geometrical centre of the multi-element array through the skull;

6. Perform CT scan of the skull on the stereotactic frame using the Open CT indicator (Elekta, Stockholm, Sweden);

7. Numerically simulate the propagation through the skull with a varying maximum speed of sound (c_max) for the model used to compute the 3D speed of sound map in the skull bone;

8. Determine the pressure field through the skull for each tested c_max using propagation matrices and phases from each simulation; and

9. Calculate the maximum pressure, the focusing volume and the distance to the geometrical centre of the multi-element array for each tested c_max.

Figure 1. (A) Diagram of the experimental setup used to measure each propagation matrix. The circles around the hydrophone tip are examples of pressure measurement locations. (B) Photograph of the experimental setup showing the skull mounted on a stereotactic frame in front of the multi-element array.

Step 1 is the calibration step, and step 2 is performed to ensure that the skull cannot move relative to the stereotactic frame during the entire procedure. The degassing in step 3 is required to ensure wave propagation and prevent the trapping of air bubbles inside the skull bone. In step 4, the stereotactic frame is placed at a unique and accurate mechanical position relative to the multi-element array, and then in step 5, the propagation matrix through the skull is measured. In step 6, the geometrical model of the skull and its position relative to the multi-element array are determined. The aberration correction based on numerical simulations is calculated in step 7 using the geometrical model of step 6 as input. Step 8 uses the simulation output of step 7 with the measurements acquired in steps 1 and 5 to determine the pressure field around the geometrical centre of the multi-element array. Finally, step 9 utilises the generated pressure fields to estimate the focusing quality of each aberration correction simulation.

Measurement of propagation matrices

Prior to all of the measurements, each skull is degassed in water with a near-vacuum pressure of 60 mmHg for 24 h.

Providing experimental evidence of the efficiency of the optimisation of a simulation parameter is a delicate and time-consuming process. In particular, the propagation medium (in our case, the CT-based co-registered human skulls) cannot move or undergo significant changes (temperature, dryness, etc.) during the optimisation process. We exploited the recording of the propagation matrix that was previously introduced for ultrasonic imaging [21,22] and for therapy [23,24]. For clarity, the salient steps are detailed in this subsection.

For each skull, we recorded the linear operator relating the K = 512 transducers of the multi-element array to a set of M = 16 184 control points embedded around the geometrical centre of the array.

An impulse response, h_mk(t), is defined for each pair (m, k) that constitutes a control point and a transducer element such that h_mk(t) is the signal received on the mth control point after a temporal delta function is applied on the kth transducer of the array (Figure 1). This response includes all the propagation effects through the medium under investigation and the acousto-electric responses of the array element and the receiver. Because the transformations are assumed to be linear and invariant under a time shift, the M × K temporal functions h_mk(t) for 1 ≤  m ≤ M and 1 ≤ k ≤ K can describe any transmitter-receiver operation for the same arrangement. This set of M×K temporal functions characterises the propagation operator.

Let e_k(t), 1 ≤ k ≤ K, be the K input signals to the array. Then, the output signals f_m(t), 1≤ m ≤ M, received in the control domain are given as follows: (1) where * is a temporal convolution operator.

A temporal Fourier transform leads to the following relation: (2)

Equation 2(2) can be written in matrix form: (3) where E(ω ) = (E_k (ω ))₁, ≤, k, ≤, K is the column vector of the Fourier transform of the transmitted signals and F(ω) = (F_m (ω))₁, ≤, m, ≤, M is the column vector of the Fourier transform of the received signals. The transfer matrix H(ω) = H_mk (ω)₁, ≤ m ≤ M, 1 ≤ k ≤ K describes the propagation in the medium from the multi-element array to the set of control points for a chosen frequency component; thus, it is referred to as the monochromatic propagation matrix [21,22].

The control points are located on a 3D grid centred on the geometrical centre of the array in a 6 × 6×12 mm³ volume of interest. As the focal spot is axially elongated (and not isotropic), a 1 mm spatial step was chosen in the axial X direction, larger than the 0.188 mm spatial step chosen in the lateral Y and Z directions. The lateral step was set to less than half the active diameter of the hydrophone. The impulse responses h_mk(t) are acquired for each skull. The electric emission is a 1 μs total duration bipolar square waveform, and reception is performed using a needle hydrophone (HCN400 model, Onda Corporation, Sunnyvale, CA; active diameter of 0.4 mm) mounted on a motorised 3-axis orthogonal translations axes gantry. The transfer matrix H(ω ) is then calculated using the Fourier transform.

The 1 MHz HIFU brain therapy array (Imasonic, Voray sur l’Ognon, France) is composed of 512 elements with an active diameter of 6 mm and a frequency bandwidth of 800 kHz to 1.3 MHz. The elements are arranged on a portion of a sphere. The multi-element array is mounted on a degassed water tank using a mechanical system on which each skull is positioned by a stereotactic frame (Figure 1(B)).

The multi-element array is controlled by a brain therapy electronic system (SuperSonic Imagine, Aix en Provence, France) and programmed with a MATLAB (MathWorks, Natick, MA) interface. One channel of the same electronic system acquires the signal measured by the hydrophone so that the volume measurements can be automated.

To calibrate each measurement, a propagation matrix in water without a skull is measured prior to each transcranial propagation matrix measurement.

Speed of sound and density model

When the propagation matrix is acquired, the stereotactic frame with the skull in place is detached from the water tank and scanned using a clinical X-ray CT scanner (Siemens Sensation Cardiac 64 scanner, Siemens, Munich, Germany), using a voltage of 80 kV.

During the acoustic measurements, transportation and CT scans, each skull is kept mounted at the same position inside its stereotactic frame to secure its mechanical position relative to the multi-element array. The CT scan is performed after placement of the clinical Open CT Indicator (Elekta, Stockholm, Sweden) [25]. The software detects the Open CT indicator so that the skull and multi-element array geometrical positions can be determined based on the Elekta stereotactic frame references.

CT images are obtained at a resolution of 0.45 × 0.45 × 0.4 mm³, which corresponds to λ/3 in water and λ/9 in the cortical bone of the skull at 1 MHz. Note that in our case, this resolution is achieved by setting the volume of interest as low as possible because the number of points reconstructed in the CT image is set to 512 points.

The 3D geometry is then calculated by the brain system planning software (SuperSonic Imagine, Aix en Provence, France), which automatically segments the skull and the Open CT Indicator. The position of the multi-element array relative to the Open CT Indicator is determined to deduce the coordinates of the elements in the numerical skull model.

The numerical model is based on the density and the speed of sound at each point of a Cartesian grid. The medium around the skull is modelled as water: ρ = 1000 kg/m³ and c = 1500 m/s. The density is calculated from the CT images (in Hounsfield units (HU)) using Equation 4(4) , where HU_min is set to −1024, HU_max is set to 2400, ρ_min is set to 1000 kg/m³, and ρ_max is set to 2700 kg/m³. (4)

The HU throughout the images are thresholded to the minimum and maximum values (HU_min and HU_max, respectively) prior to applying Equation 4(4) .

The speed of sound is calculated as a linear function of the apparent density of the skull using Equation 5(5) . The speed of sound is considered to be minimal in the most porous areas where voxels contain only water; therefore, c_min was set to 1500 m/s and c_max was adjusted. To optimise the speed of sound model, c_max was varied between 3000 and 5000 m/s as described in detail in the following subsection. (5)

Examples of density map and the corresponding speed of sound through the skull are shown in Figure 2, as well as density plots across the bone.

Figure 2. Coronal views of skull density and speed of sound used in the simulation for skulls #1, #2 and #3. Bottom right: density path point along the vertical axis (distance =140 mm on the coronal views): ‘−’ skull 1, ‘.’ skull 2, ‘−•’ skull 3, ‘- -‘skull 4,•‘ ^’ skull 5 and ‘×’ skull 6. For plot, the traversed distance is normalised by the thickness of the skull bone measured at each location.

Numerical simulation

A 3D numerical simulation of a spherical acoustic wave originating from the geometrical centre of the multi-element array is simulated through the skull. The central frequency of the pulse is set to 1 MHz. The fluid linear wave in Equation 6(6) is discretised by a second-order finite difference method in the spatial domain and a first-order method in the temporal domain. (6)

Equation 6(6) . Heterogeneous wave equation used for the simulation.

The 3D CT digital maps are interpolated to dx = λ/10 = 0.15 mm. The time step of the simulation is set to dt = dx/c_max/√3 to comply with the digital stability condition [26], where c_max is the maximum speed of sound through the skull. The surface of each of the 512 transducers is discretised into 317 virtual point sources in the simulation with a pitch of λ/5 (0.3 mm), and the signal is integrated over the surface of the transducer to take into account the directionality of the elements.

Two different acoustical media are considered: homogeneous soft tissue (brain and skin) and heterogeneous tissue (skull bones). The 3D simulations are only performed from the outer surface of the brain and not from the original point source. Specifically, a set of synchronised sources are located at the closest grid points to a sphere located close to the skull assuming a homogeneous medium in the brain. In this study, both types of tissue are considered as fluids, whereas in Pulkkinen et al. [27], the skull was modelled as a solid by coupling acoustic and visco-elastic wave equations.

The computer code was developed by SuperSonic Imagine, and approximately 3 h are required to perform the full 3D propagation through the skull on a dual Intel Xeon© computer running at 2.6 GHz with 24 GB of memory.

For each measured propagation matrix, several simulations reproduce the specific geometry of the experiment with different parameters for the speed of sound through the skull. The minimum speed is always set to 1500 m/s because it is assumed that the lowest skull densities contain porous zones filled with water. Five values for the maximum speed of sound are tested: 3000, 3500, 4000, 4500 and 5000 m/s.

Focusing using propagation matrices and aberration correction

The monochromatic propagation matrices H(ω) is derived from the corresponding time signals using fast Fourier transform. To match a typical HIFU treatment using monochromatic signal, only the 1 MHz frequency is selected and used in this study.

For each skull, the simulation signals are then normalised and multiplied by the propagation matrix to yield the corresponding pressure field in the focal area according to Equation 2(2) .

The results are compared with an invasive experimental time reversal operation, which is an actual physical experiment considered the gold-standard reference correction. The reference experiment is also based on the previously measured propagation matrix [28].

The therapeutic array is used as a time reversal mirror: a signal s(t) is emitted at the geometric centre of the array, which corresponds to the emission vector s(t) ^t[0,… 0,1,0,…., 0], where m = M/2 is the only value set to one in the vector. We denote S(ω) as the Fourier transform of this emission vector. H(ω) is the propagator from the therapeutic array to the target area, and ^tH(ω) is the propagator from the target to the array because of spatial reciprocity: h_mk^2→1(t) = h_km^1→2(t). According to the previous section, the signals recorded on the array are determined as follows: (7)

The received signals are then time reversed. This time-reversal operation is mathematically described by the variable change (t) →(−t), which is equivalent to a phase conjugation operation in the frequency domain. The re-emitted signals are then determined as follows: (8)

The set of signals E(ω) are then propagated towards the focus area, and the signals obtained around the geometrical centre are finally determined as follows: (9)

For each 3D field obtained using aberration correction, the maximum pressure and its distance to the geometrical centre of the multi-element array are determined. In each case, the maximum amplitude ratio is calculated relative to the maximum amplitude obtained with the corresponding optimal correction. In addition, the focusing volume is calculated for each correction as the volume for which pressure is higher than half the maximum.

For each skull the pressure amplitude at the acoustical focus was measured with the hydrophone for each of the 512 elements. The average pressure was calculated and compared to the value without the presence of the skull bone: attenuation was estimated for each skull as the ratio of the average pressure with the skull to the average pressure without the skull provided.

Results

The transcranial pressure field obtained at 1 MHz for the six skulls using the propagation matrix were calculated. The typical impact of focusing on the skull can be observed in Figure 3 (left) when transmission is performed without correction. In this case, the maximum pressure is 41% of the optimal pressure (Figure 3, right), which is located more than 2 mm from the geometrical centre of the multi-element array, and side lobes are present. The average maximum pressure without correction is 45% of the optimal correction amplitude, and the average shift is 1.05 mm. Detailed results are provided in Tables 1 and 2.

Figure 3. Slices in the three orientations when focusing through skull #3. Left: without correction; centre: with simulation correction; right: with the measured correction (invasive time reversal). Scales are normalised to the maximal amplitude of the optimal focusing.

Table 1. Maximal transcranial pressure with and without correction normalised by the maximal amplitude of the optimal focusing for each skull.

Skull #	1	2	3	4	5	6	Average relative pressure
Relative maximal pressure without correction	0.456	0.429	0.407	0.506	0.469	0.458	0.454
Relative maximal pressure with simulation correction	0.820	0.845	0.873	0.885	0.880	0.870	0.862

Table 2. Distance of the maximum amplitude relative to the geometrical centre of the multi-element array in the three directions without aberration correction.

Skull #	1	2	3	4	5	6	Average absolute distance (mm)
Distance in axial direction X (mm)	2.00	0.00	2.00	0.00	0.00	0.00	0.667
Distance in lateral direction Y (mm)	−0.188	0.00	0.563	0.375	−0.188	0.563	0.313
Distance in lateral direction Z (mm)	−0.188	0.375	−0.563	−0.750	0.188	−0.375	0.469
Distance (mm)	2.02	0.375	2.15	0.836	0.266	0.677	1.05

An example of the pressure field obtained using the numerical 3D simulation correction is plotted for skull #3 in Figure 3 (centre). The restored focusing pattern is similar to the optimal case (Figure 3, right), and the maximum position is less than 0.2 mm from the geometrical centre.

The maximal pressures obtained using the different simulation correction methods are plotted in Figure 4. The effectiveness of the optimal simulation corrections varies from 82 to 88% and appears to be consistent among the skulls. The maximum pressures exhibit an average amplitude of 86% ± 2.7% and a 0.54 mm shift compared with the optimal correction (details provided in Tables 1 and 3). The highest amplitude is obtained for a maximum speed of either 4000 or 4500 m/s. In addition, the maximum amplitude is located an average of 0.125, 0.188 and 0.333 mm from the geometrical centre in the Y and Z lateral directions and X axial direction, respectively, as detailed in Table 3. This shift was two-fold smaller than the shift without correction (Table 2).

Figure 4. Normalised maximum amplitudes obtained with and without the 3D numerical simulation correction. Normalisation is relative to the maximal amplitude obtained by optimal focusing (invasive time reversal).

Table 3. Distance of the maximum amplitude relative to the geometrical centre of the multi-element array in the three directions, with aberration correction by simulation.

Skull #	1	2	3	4	5	6	Average absolute distance (mm)
Distance in axial direction X (mm)	1.00	0.00	0.00	0.00	−1.00	0.00	0.333
Distance in lateral direction Y (mm)	0.00	0.00	0.188	0.188	−0.188	0.188	0.125
Distance in lateral direction Z (mm)	−0.188	0.188	0.00	−0.375	0.00	−0.375	0.188
Distance (mm)	1.02	0.188	0.188	0.42	1.02	0.42	0.543

For each skull, the phase differences between the optimal simulation (maximum 4000 m.s⁻¹ speed of sound in the skull) and the experimental data were calculated (Table 4). The standard deviation of the absolute phase differences varied from 0.51 to 0.62, illustrating the robustness of the simulation among the six human skulls. The spatial distribution of the transcranial phase estimation over the transducer surface obtained by the optimal simulation is plotted for skull #4 in Figure 5. Similar phase patterns were observed for the experiments and the simulations (left and centre). In addition, the corresponding phase differences plotted in Figure 5 (right) show areas with a nearly constant phase offset (top and bottom of the multi-element array) and few elements with higher phase errors.

Figure 5. Spatial distribution of the phase on the array surface for skull #4 with a maximal speed of sound of 4000 m/s. From left to right: experimental phases, simulation phases and their differences.

Table 4. Phase differences and absolute phase differences between the simulation and the experimental data for each skull.

Skull #	1	2	3	4	5	6
Mean phase differences (rad)	0.0004	−0.0097	0.0075	0.0369	0.1063	0.0108
Standard deviation of the phase differences (rad)	0.788	0.973	0.734	0.766	0.747	0.873
Mean absolute phase differences (rad)	0.599	0.754	0.527	0.599	0.551	0.686
Standard deviation of the absolute phase differences (rad)	0.511	0.615	0.512	0.479	0.514	0.542

The calculated focusing volume for each correction is plotted in Figure 6. For all skulls, the focusing volume without correction is many times larger than the focusing volume calculated with the optimal correction, thus demonstrating the defocusing impact of the skull. However, the simulation-based aberration correction leads to a focusing quality near the optimal correction.

Figure 6. Normalised focusing volume obtained with and without the 3D numerical simulation correction. The normalisation is relative to the focusing volume obtained by the optimal focusing (invasive time reversal).

The mean attenuation measured for each skull is presented in Table 5 and Figure 7. The frequency dependence of the attenuation can be clearly observed in Figure 6 and corresponds to highly similar slopes among the skulls. The average transmission coefficient at 836 kHz was 0.34, which is 6% higher than the value of 0.32 found in Ref [16].

Figure 7. Mean attenuation for each skull within the bandwidth from 800 kHz to 1.3 MHz.

Table 5. Transmission coefficients and resulting attenuation for each skull at 1 MHz.

Skull #	1	2	3	4	5	6
Amplitude transmission coefficient	0.23	0.27	0.29	0.39	0.35	0.21
Standard deviation of the transmission coefficient	0.11	0.24	0.12	0.18	0.18	0.14
Attenuation (dB)	12.7	11.3	10.6	8.18	9.12	13.6

The most efficient aberration corrections are obtained for skulls four and five, which also present the lowest attenuation. Of note, the highest attenuation values (skulls 1, 2 and 6) are observed in the skulls with the highest average absolute phase difference between the simulation and experimental data.

For each skull, for each simulation, we calculated the speed of the sound distribution inside the skull, and the corresponding histograms are plotted in Figure 8. Although considerable differences are observed for the distributions among skulls, the distributions do not appear to impact the aberration correction method, which is robust for the skull structure.

Figure 8. Speed of sound histograms for each skull, with a maximum speed set to 4000 m/s.

Discussion and conclusions

This study demonstrates that a CT-based aberration correction on a collection of six human skulls can restore 86% of the optimal pressure amplitude at frequencies as high as 1 MHz. The best model for the aberration correction corresponds to a maximum transcranial speed of sound of 4000 m.s⁻¹. The optimal maximal speed of sound of 4000 m.s⁻¹ through the skull bone determined in this study should not be directly compared with previous measurements near 3000 m.s⁻¹ [1]. In our case, the maximal speed of sound was only a parameter of our linear model and should not be considered an actual measurement. In addition, the maximal density values of the external cortical bone were only found in certain areas (like the sides of the skull in Figure 2) and may not correspond to the measurement locations of previous studies.

Compared with the non-corrected case, the non-invasive correction leads to an average increase in the maximum pressure amplitude around the target of 99% and an average decrease in the distance between the peak pressure and the selected target of 48%. In addition, the simulation-based corrections lead to nearly optimal focusing volumes, which confirm the obtained focusing quality. This performance indicates the potential for efficient clinical applications of high-frequency non-invasive focused ultrasound brain therapy.

Comparison with previous studies

In this article, the maximum speed of sound (c_max in Equation 5(5) was adjusted to optimise transskull refocusing. This linear relationship between the speed of sound in the bone and the bone density (Equation (5)(5) has been used in the past [17,19]. c_max was set to 2900 m/s in Ref [17] and to 3100 m/s in Ref [19]. These values appear to be significantly lower than the optimum found in this article (4000 m/s). However, the density of the skulls used in this article is also higher (2800 kg/m³) compared to Ref [17] (2100 kg/m³) and Ref [19] (2200 kg/m³). Figure 9 shows the relationship between the speed of sound in the bone and the bone density for all the parameters tested in this article, and for previous models. It demonstrates that all models are consistent.

Figure 9. Speed of sound as a function of skull density for all the values of the maximum speed of sound (c_max) tested in this article and for previously published parameters (Ref [17], black dotted line and Ref [19], green dotted line).

Limitations of the work

Even if the patterns of the simulated phase are similar to that of experimental patterns, residual errors persist. The origins of the differences between the simulated and experimental phases are not fully understood, although two main limitations could prevent the realisation of an optimal correction. The first limitation is the use of a fluid model that considers the compressional wave only as it propagates through the skull. The angle beam incidences at the outer and inner bone surfaces and the complex geometry of the diploe lead to mode conversions that can influence the phase aberration [29,30]. Three dimensional numerical models have been developed to consider shear mode conversions in heterogeneous media [31–33], although they require more memory and computer time than fluid models [3]. To lower the computational cost, models that couple cost-efficient simulations in fluid (soft) and solid (skull) tissues with acoustic and visco-elastic wave equations have been developed [26]. Optimised computational costs were quantified in Ref [26], in which simulations were performed on an HP CP4000 BL ProLiant cluster Vuori (Finnish IT Center for Science (CSC)). Eight computing nodes that each contained two six-core 2.6 GHz AMD Opteron 64-bit processors were used for each simulation. The computing nodes communicated via an InfiniBand network. The average computation time for a single simulation was 65 ± 3 h for sonications with the skull modelled as a solid, whereas the total memory used in the simulation model was approximately 64 GB. The simulation time required for the use of the simplified model in which the skull was modelled as a fluid was 40 h. For the sake of comparison, the code presented here requires 2–3 h (depending on the skull) to perform the full 3D fluid propagation analysis through the skull on a dual Intel Xeon^© computer running at 2.6 GHz with 24 GB of memory.

The second limitation is the spatial resolution of the clinical CT scanner used to determine the skull density. Based on the speed of sound variations in the skull bone, the resolution corresponds to λ/9 in the cortical area of the skull but only λ/3 in the central porous layer. We carefully set the field of view to its minimum value to fully exploit the 512 × 512 pixel limitation of current clinical scanners. Therefore, the resolution of the CT cannot be improved with current clinical CT scanners, although novel large field micro-CTs [34] have the potential for use in further investigations of the impact of CT image resolution. Pioneering work that considers the bone micro structure was based on synchrotron tomography [30,33].

Additionally, the finite-difference simulation used in this article intrinsically leads to a numerical dispersion under our conditions, in which the medium includes speed of sound heterogeneities [35]. Further investigations should compare the results with frequency domain simulations to determine the impact of these limitations [31,36].

Nonlinearities were not considered in our study but could be implemented [31,37,38]; however, the low-pass filtering effect of the skull bone limits nonlinearities to the focal area [39], contrary to abdominal applications in which nonlinearities develop from the near field to the far field.

Cavitation was also neglected in this study, although it is known to enhance the therapeutic effect of focused ultrasound [40–42]. Cavitation is currently monitored and avoided during clinical transcranial thalamotomy for safety reasons [43].

Our experimental setup was used to estimate the accuracy of a specific numerical model. However, the propagation matrix measurements in this study can be used to test any aberration correction method based on clinical CT scans for all of the frequencies within the bandwidth of the multi-element array. These matrices have recently been used to test a novel random-based focusing technique [44]; however, they could also be used to compare simpler, faster simulation models in addition to more complex, computationally intensive models. The data are too large to be provided as supplementary material; however, the corresponding author will share the data with any group willing to test a numerical approach in an experimental setup.

This study validates our skull aberration correction method and the estimated optimal model parameters and verifies these data for use in vivo during brain therapy treatments on patients. Future work will first focus on estimating the aberration correction for different frequencies and target locations (such as non-central locations, which remain challenging). Ultrasonic attenuation simulations must be included in the model to predict the skull areas that best transmit the acoustic energy. Such predictions may lead to more pertinent selections of clinical patients for brain HIFU therapy and facilitate the optimisation of the position of the multi-element array.

Acknowledgements

This work was supported by the Bettencourt Schueller Foundation and the Agence Nationale de la Recherche under the programme “Future Investments”’ under reference ANR-10-EQPX-15. The Institut Langevin is also supported by LABEX WIFI (Laboratory of Excellence ANR-10-LABX-24) within the French Program “Investments for the Future” under reference ANR-10-IDEX-0001–02 PSL*. The authors wish to thank the anonymous reviewers for their comments and thoughtful suggestions.

Disclosure statement

Mickael Tanter is cofounder and a shareholder of Supersonic Imaging, Aix en Provence, France.

Laurent Marsac and Raphael La Greca are employees of Supersonic Imaging, Aix en Provence, France.

Additional information

Funding

This work was supported by the Bettencourt Schueller Foundation and the Agence Nationale de la Recherche under the programme “Future Investments”’ under reference ANR-10-EQPX-15. The Institut Langevin is also supported by LABEX WIFI (Laboratory of Excellence ANR-10-LABX-24) within the French Program “Investments for the Future” under reference ANR-10-IDEX-0001–02 PSL*.