DOI: 10.5068/d1qd60
title: SMART FUS: Surrogate model of attenuation and refraction in transcranial focused ultrasound

description: 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.
description: 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<sup>7</sup> 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<sub>alpha </sub>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<sup>12</sup> 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: <i>Radius of Curvature=
 </i><i>(Transducer Focal
 Depth)</i><sup><i>2</i></sup><i>+
 </i><i>(Transducer
 Radius)</i><sup><i>2</i></sup> 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<sup>12</sup>, 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<sup>3</sup> and speed of sound of 1482 m/s, corresponding to that of water. A simulation time of 8.0375e<sup><i>-5</i></sup><i>s</i>  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.<sub> </sub>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<sup>3</sup> (see <b>Discussion</b>). 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 85<i>Np </i><i>Mhz</i><sup><i>-1</i></sup><i>m</i><sup><i>-1</i></sup> , while speed of sound in bone was set to 2850m/s, and Bone density was set to 1732 kg/m<sup>3</sup> as has been defined previously for Bulk Bone<sup>5,13</sup> and appropriate for homogeneous simulations.
description: 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).