Low-intensity focused ultrasound (LIFU) is an increasingly applied method for achieving non-invasive brain stimulation. However, transmission of ultrasound through the human skull can substantially affect focal point characteristics of LIFU, including dramatic attenuation in intensity and refraction of focal point location. These effects depend on a high-dimensional parameter space, making these effects difficult to estimate from previous work. Instead, focal point properties of LIFU experiments are often estimated using numerical simulation of LIFU sonication through skull. However, this procedure presents many entry barriers to even computationally savvy investigators and often requires expensive computational hardware, impeding LIFU research. We present a novel MATLAB toolbox for rapidly estimating beam properties of LIFU transmitted through bone. Users provide specific values for frequency of LIFU, bone thickness, angle at which LIFU is applied, depth of the LIFU focal point, and diameter of the transducer used and receive an estimation of the degree of refraction/attenuation expected for the given parameters.

We defined a five-dimensional parameter space, comprising carrier frequency, bone thickness, trajectory, transducer size, and focal depth (detailed in Figure 2A), in which we simulated acoustic waves through bone. Some portions of this parameter space placed the theoretical focus either inside bone or outside the head entirely; inputting these values into SMART_FUS will return an error. Similarly, inputting values outside the parameter space studied will also return an error. In all, 12,096 simulations were run using an NVIDIA RTX 3080 GPU with each simulation taking roughly 10 seconds to complete, making for approximately fourteen days of simulation time.

We applied the MATLAB-based pseudospectral solver toolbox k-Wave⁷ to perform the simulations. A simulation space of 256 × 256 × 256 voxels with isotropic dimensions of 0.5mm each was created. Using dimensions that were powers of two maximized the speed of the Fast Fourier Transform used in k-Wave’s Fourier collocation method. A perfectly matched layer (PML) was created to absorb energy at the edge of the simulation; default values of 10 voxels and a PML_alpha of 2 were used. A Courant-Friedrichs-Lewy (CFL) of 0.3 (the default value for k-Wave) was chosen. A trade-off between runtime and simulation fidelity exists for CFL values. A value of 0.3 is appropriate for simple simulations, such as these, and maintains feasible simulation times. When ignoring the PLM layer, our simulations represent a space of 11.8 cm in each direction, capable of subsuming the necessary components of each simulation: a transducer up to 80mm in diameter with a focal depth up to 80mm, and a bone thickness up to 12mm. These dimensions, given the frequencies sampled, result in a points-per-wavelength (PPW) ranging from 11.86 (at the lowest frequency) to 2.96 (at the highest frequency). Homogenous 3D simulations have been observed to converge (in terms of recorded intensity and focal location) at a PPW of 3.5, even with more complex skull models¹² with values lower than this introducing some error. As our medium contains non-aliased skull-models (discussed more below) and a lower CFL of 0.3 (resulting in higher fidelity), compared to the cited study’s CFL of 0.5, a slightly lower PPW of 2.96 (compared to 3.5) at our high-range remains valid.

For each simulation, a single-element transducer was modeled and assigned as the source of the acoustic waveform. The modeled transducer has dimensions that corresponded to a certain radius of curvature and diameter. The radius of curvature was determined using the desired focal depth of the transducer with the simple formula:

Radius of Curvature= (Transducer Focal Depth)²+ (Transducer Radius)²

Transducers were generated in voxel space using k-Wave’s makeBowl function and placed such that they met the upper-left corner of the PML. Bone was assigned to a flat wall of voxels in front of and meeting the transducer face. While transducers assigned to a perpendicular trajectory met bone along its entire face, those assigned a non-perpendicular trajectory were tilted away from the PML so their foci fully avoid the perimeter of the space. Regardless of trajectory, the voxel structure of bone remained flat to avoid an aliased surface (see Figure 1). Since aliasing has been characterized as the most significant source of error in simulations of this type¹², this likely enhances the fidelity of our simulations considerably.

Each voxel in the space not corresponding to bone was assigned a density of 1000 kg/m³ and speed of sound of 1482 m/s, corresponding to that of water. A simulation time of 8.0375e^-5s  was chosen as this allowed the simulated acoustic wave to reach the opposing end of the simulation space. The intensity emitted from the transducer was set to 1.0558MPa. Because all outputs from this function reflect the relative attenuation/refraction between water and bone simulations, this value is arbitrary but was selected from measurements taken from a previous empirical study³ (see Discussion). Because the absorption coefficient for water changes in relation to the frequency of sound being simulated, we calculated it using k-Wave’s water absorption function by inputting the frequency used and a temperature of 37c (body temperature). Bone attenuation was set to 85Np Mhz^-1m^-1 , while speed of sound in bone was set to 2850m/s, and Bone density was set to 1732 kg/m³ as has been defined previously for Bulk Bone^5,13 and appropriate for homogeneous simulations.

SMART_FUS includes two functions: SMART_FUS and SMART_FUS_vis2d.

SMART_FUS takes as input values for each dimension studied and returns an estimate of attenuation and refraction based on linear interpolation of two datasets that sample the attenuation and refraction in the parameter space (see below). Additionally, the nearest simulation point is visualized, providing the user an estimate of focal shape for the given parameters (see Figure 1A).

SMART_FUS_vis2d takes as input values for any three of the five dimensions studied and returns two 2D images depicting the refraction and attenuation expected at points across the two unspecified dimensions, which represent 2D slices through the full 5D parameter space. These 2D arrays are extracted and linearly interpolated to upscale the output matrix and provide a smooth space for visualization (see Figure 1B).