The Effects of the Structural and Acoustic Parameters of the Skull Model on Transcranial Focused Ultrasound

Transcranial focused ultrasound (tFUS) has great potential in brain imaging and therapy. However, the structural and acoustic differences of the skull will cause a large number of technical problems in the application of tFUS, such as low focus energy, focal shift, and defocusing. To have a comprehensive understanding of the skull effect on tFUS, this study investigated the effects of the structural parameters (thickness, radius of curvature, and distance from the transducer) and acoustic parameters (density, acoustic speed, and absorption coefficient) of the skull model on tFUS based on acrylic plates and two simulation methods (self-programming and COMSOL). For structural parameters, our research shows that as the three factors increase the unit distance, the attenuation caused from large to small is the thickness (0.357 dB/mm), the distance to transducer (0.048 dB/mm), and the radius of curvature (0.027 dB/mm). For acoustic parameters, the attenuation caused by density (0.024 dB/30 kg/m3) and acoustic speed (0.021 dB/30 m/s) are basically the same. Additionally, as the absorption coefficient increases, the focus acoustic pressure decays exponentially. The thickness of the structural parameters and the absorption coefficient of the acoustic parameters are the most important factors leading to the attenuation of tFUS. The experimental and simulation trends are highly consistent. This work contributes to the comprehensive and quantitative understanding of how the skull influences tFUS, which further enhances the application of tFUS in neuromodulation research and treatment.


Introduction
Transcranial focused ultrasound (tFUS) is a non-invasive, localized, non-ionizing technique with applications in particular in oncology and neurosurgery. It has been investigated for over 60 years, with seminal studies by Lynn and Fry establishing the potential of applying ultrasound to cerebral tissue [1][2][3]. Depending on stimulation parameters, tFUS is able to determine a wide spectrum of effects, ranging from suppression or facilitation of neural activity to tissue ablation [4]. According to the acoustic intensity of the focal area, tFUS can be divided into low-intensity focused ultrasound (LIFU) and high-intensity focused ultrasound (HIFU).
Numerous experimental studies have demonstrated that LIFU is an emerging tool for non-invasive neuromodulation. It can transmit low-intensity ultrasound through the Sensors 2021, 21, 5962 2 of 14 skull to temporarily and safely modulate regional brain activity [5][6][7][8]. In 2010, Tufail used targeted LIFU to stimulate neuronal activity in the intact hippocampus. The experimental results show that ultrasound-stimulated neuronal activity was sufficient to evoke mice's motor behaviors [9]. In 2013, Deffieux demonstrated the feasibility of using LIFU stimulation to causally modulate behavior in the awake nonhuman primate brain for the first time [10]. In 2014, Legon probed the influence of targeted LIFU to the human primary somatosensory cortex on sensory-evoked brain activity and sensory discrimination abilities. Electroencephalographic recordings showed that LIFU can be used to focally modulate human cortical function [11]. LIFU can not only achieve neuromodulation in normal animals and humans, but also has shown positive and effective effects in the prevention and treatment of neurological diseases such as epilepsy [12], Parkinson's disease [13], and depression [14].
As a non-invasive tumor treatment technology, HIFU has been widely used in the treatment of solid tumors and has broad application prospects [15]. In initial trials in three glioblastoma patients in 2010, multiple HIFU exposures were applied up to the maximum acoustic power available. However, the skull makes ultrasound delivery difficult or even impossible and thus seemed to not achieve thermal coagulation [16]. In 2014, Coluccia applied HIFU under MR imaging guidance, and the partial tumor ablation could be achieved. This work is the first to demonstrate the feasibility of using tFUS to safely ablate substantial volumes of brain tumor tissue. Unfortunately, the applications of tFUS in the brain have been inhibited by the skull, as well as the inability to target and monitor treatments. Over the years, with the goal of understanding the impact of the skull on tFUS, scholars have done a lot of work, including mapping the geometry of the skull [17], exploring the acoustic parameters of the skull [18,19], trying skull phantom [20], and exploring the effect of skull thickness and density on tFUS [21][22][23]. Although previous studies have achieved certain results, these works have not been systematically and comprehensively studied on the influence of the skull on tFUS.
Therefore, the present study aims to investigate this aspect. This article mainly consists of two parts: the first part focuses on the quantitative description of the influence of structural parameters (thickness, curvature radius, and distance to transducer) on tFUS, while the second part contains an accuracy assessment of the impact of acoustic parameters (density, acoustic speed, and absorption coefficient).

Mathematical Equation for Acoustic Field
Based on the self-programming simulation, we have already accurately reproduced the tFUS propagation through the skull [24,25]. Throughout the paper, the nonlinear Westervelt equation is used to describe acoustic fields in the simulation [26,27]. It incorporates the effects of nonlinearity, attenuation, reflection, refraction, and scattering of waves, and can be written as In this equation, p is acoustic pressure, c is acoustic speed, t is time, ρ is medium density, δ is the thermoviscous diffusivity, β is the nonlinearity coefficient. The equation is calculated by the finite-difference time-domain method (FDTD) in 3D. Since the center frequency of the ultrasound is 0.9 MHz, the simulation space step is set to one-eighth of the ultrasonic wavelength in water (0.2 mm), and the time step is 10 ns.
In the COMSOL Multiphysics 5.4 (COMSOL Inc., Stockholm, Sweden) simulation part, this paper uses the Helmholtz equation to calculate the acoustic field In this equation, r is radial coordinate, z. is axial coordinate, p is acoustic pressure, ω is angular frequency, ρ c is density, and c is acoustic speed. The remaining settings of the two simulations are exactly the same. Numerical simulations are performed on the OMNISKY MultiGPU CSC4450 server platform with an internal configuration of dual Xeon E5-2697 V3 processor (2.60 GHz). Figure 1 schematically shows the numerical simulation model, including the transducer, skull, and water, which size of the numerical domain is 100 × 100 × 100 mm. The Mur first-order absorbing boundary condition was used at the boundary to avoid the reflection of ultrasound [28]. Among them, we first focus on the thickness of the skull (T: 2:4:22 mm), the radius of curvature (R: 60:10:100 mm), and the distance to transducer (D: 11:1:20 mm) according to different positions (such as temporal bone, occipital bone, parietal bone, etc.) of the skull [17]. According to previous studies, skull density, acoustic speed, and absorption coefficient trended toward a reduction with increasing age [29,30]. According to the suggestion in the literature, the density of the skull (ρ: 1700:30:2000 kg/m 3 ) [19], the acoustic speed (c: 2000:300:5000 m/s) [18], and the absorption coefficient (α: 200:480:5000 dB/m) [31] are also considered in the present study [32,33]. The physical properties of the objects are shown in Table 1. In this equation, is radial coordinate, is axial coordinate, is acoustic pressure, is angular frequency, is density, and is acoustic speed. The remaining settings of the two simulations are exactly the same. Numerical simulations are performed on the OMNISKY MultiGPU CSC4450 server platform with an internal configuration of dual Xeon E5-2697 V3 processor (2.60 GHz). Figure 1 schematically shows the numerical simulation model, including the transducer, skull, and water, which size of the numerical domain is 100 × 100 × 100 mm. The Mur first-order absorbing boundary condition was used at the boundary to avoid the reflection of ultrasound [28]. Among them, we first focus on the thickness of the skull (T: 2: 4: 22 mm), the radius of curvature (R: 60: 10: 100 mm), and the distance to transducer (D: 11: 1: 20 mm) according to different positions (such as temporal bone, occipital bone, parietal bone, etc.) of the skull [17]. According to previous studies, skull density, acoustic speed, and absorption coefficient trended toward a reduction with increasing age [29,30]. According to the suggestion in the literature, the density of the skull (ρ: 1700: 30: 2000 kg m ⁄ ) [19], the acoustic speed (c: 2000: 300: 5000 m s ⁄ ) [18], and the absorption coefficient (α: 200: 480: 5000 dB m ⁄ ) [31] are also considered in the present study [32,33]. The physical properties of the objects are shown in Table 1.   [18] 1482 [18] 0.02 [19] 3.5 [24] Acrylic 1203 [18] 2600 [18] 18 [34] 4.4 [24]

Experimental Setup
Since the acoustic properties of the acrylic plate (such as density, sound speed, and attenuation coefficient, etc.) are close to those of the skull, the acrylic plate phantom mimics the skull [18]. Based on computer numerical control (CNC) technology, acrylic plates with different thicknesses and curvature radius are designed ( Figure 2). The experimental schematic and device diagram are shown in Figure 3. The ultrasound transducer is a custom designed single-element piezoelectric ceramic focused transducer (Fluids Engineering Laboratory, Department of Mechanical Engineering, the University of Tokyo, Japan), which has a center frequency of 0.9 MHz, an opening diameter of 100 mm, a focal length of 80 mm, and a central hole diameter of 50 mm. The power amplifier (AG1021 T&C, US) is used to amplify the signal induced by the signal generator (33522A Agilent, US), and further drive the transducer using the continuous wave (CW) mode to generate a 10s continuous sine wave. The ultrasound field around the geometrical focus is measured by polyvinylidene fluoride (PVDF) membrane needle hydrophone (NCS-1, China) mounted on

Experimental Setup
Since the acoustic properties of the acrylic plate (such as density, sound speed, and attenuation coefficient, etc.) are close to those of the skull, the acrylic plate phantom mimics the skull [18]. Based on computer numerical control (CNC) technology, acrylic plates with different thicknesses and curvature radius are designed ( Figure 2). The experimental schematic and device diagram are shown in Figure 3. The ultrasound transducer is a custom designed single-element piezoelectric ceramic focused transducer (Fluids Engineering Laboratory, Department of Mechanical Engineering, the University of Tokyo, Japan), which has a center frequency of 0.9 MHz, an opening diameter of 100 mm, a focal length of 80 mm, and a central hole diameter of 50 mm. The power amplifier (AG1021 T&C, USA) is used to amplify the signal induced by the signal generator (33522A Agilent, Santa Clara, CA, USA), and further drive the transducer using the continuous wave (CW) mode to generate a 10s continuous sine wave. The ultrasound field around the geometrical focus is measured by polyvinylidene fluoride (PVDF) membrane needle hydrophone (NCS-1, Beijing, China) mounted on a three-axis stage (Sigma Koki, Tokyo, Japan). All experiments are carried out in a degassed water tank.

Results
The results section is divided into three parts. First, the consistency between the simulation and experiment is explored. The second and third parts are the influence of skull model structural parameters and acoustic parameters on tFUS.

Transducer Surface Pressure
To ensure consistency of the initial conditions of the experiment and simulation, we first use a hydrophone close to the transducer in pure water to measure its surface sound pressure. Then, the acoustic pressures at 12 points distributed on the surface of the transducer according to the clock dial are measured under different excitation conditions (Figure 4a), and its average value is used as the initial acoustic pressure for the simulations. Figure 4b,c show the acoustic pressure value at each point and overall trends under different excitation power. It can be seen that the transducer surface pressure increased with increasing excitation power, and fluctuated slightly. a three-axis stage (Sigma Koki, Japan). All experiments are carried out in a degassed wa tank.

Results
The results section is divided into three parts. First, the consistency between the si ulation and experiment is explored. The second and third parts are the influence of sk model structural parameters and acoustic parameters on tFUS.

Transducer Surface Pressure
To ensure consistency of the initial conditions of the experiment and simulation, first use a hydrophone close to the transducer in pure water to measure its surface sou pressure. Then, the acoustic pressures at 12 points distributed on the surface of the tra ducer according to the clock dial are measured under different excitation conditions (F ure 4a), and its average value is used as the initial acoustic pressure for the simulatio Figure 4b,c show the acoustic pressure value at each point and overall trends under d ferent excitation power. It can be seen that the transducer surface pressure increased w increasing excitation power, and fluctuated slightly.

Results
The results section is divided into three parts. First, the consistency between the simulation and experiment is explored. The second and third parts are the influence of skull model structural parameters and acoustic parameters on tFUS.

Transducer Surface Pressure
To ensure consistency of the initial conditions of the experiment and simulation, we first use a hydrophone close to the transducer in pure water to measure its surface sound pressure. Then, the acoustic pressures at 12 points distributed on the surface of the transducer according to the clock dial are measured under different excitation conditions (Figure 4a), and its average value is used as the initial acoustic pressure for the simulations. Figure 4b,c show the acoustic pressure value at each point and overall trends under different excitation power. It can be seen that the transducer surface pressure increased with increasing excitation power, and fluctuated slightly.

Acoustic Pressure Field Comparison
Under the same initial conditions (excitation power 1w), we have explored the influence of the presence of skull model (T8 mm, R60 mm, D10 mm) on tFUS. As shown in Figure 5a-d, both simulation and experiment accurately focus at the geometric focal point 80 mm under without skull model condition. The influence of the acrylic plate as an elastic body on the distortion of focus is shown in Figure 5e-h; while Figure 5e is a self-programming acoustic pressure field, the Figure 5f is COMSOL acoustic pressure field. In the presence of the bone plate, the experimental and simulated focal points have shifted about 2-3 mm from the geometric focal point to the transducer side, which is consistent with the results of previous studies [35]. The attenuation caused by the acrylic plate was very strong, with 56.22%, 57.65%, and 56.87% for self-programming, experimental, and COM-SOL, respectively (Table 2), which is very similar to the results in the literature [19,36]. In other words, the existence of the skull model will significantly change the focus position and the focusing region energy.  In the presence of the bone plate, the experimental and simulated focal points have shifted about 2-3 mm from the geometric focal point to the transducer side, which is consistent with the results of previous studies [35]. The attenuation caused by the acrylic plate was very strong, with 56.22%, 57.65%, and 56.87% for self-programming, experimental, and COMSOL, respectively (Table 2), which is very similar to the results in the literature [19,36]. In other words, the existence of the skull model will significantly change the focus position and the focusing region energy.   Figure 5f is COMSOL acoustic pressure field. In the presence of the bone plate, the experimental and simulated focal points have shifted about 2-3 mm from the geometric focal point to the transducer side, which is consistent with the results of previous studies [35]. The attenuation caused by the acrylic plate was very strong, with 56.22%, 57.65%, and 56.87% for self-programming, experimental, and COM-SOL, respectively (Table 2), which is very similar to the results in the literature [19,36]. In other words, the existence of the skull model will significantly change the focus position and the focusing region energy.

The Influence of the Skull Model Structural Parameter
Under the condition of ensuring the rest of the parameters remained unchanged, we have changed the thickness of the skull model, the radius of curvature, and the distance to transducer one by one.

Thickness
The influence of acrylic plate thickness on ultrasound propagation is examined next. The fixed distance of the acrylic plate (Excitation power: 1 w, R: 60 mm) from the open end of the transducer (D) is 10 mm and its thickness grew from 2 mm to 22 mm. Figure 6a displays the linear fitting results where different colors represent different numerical simulations and experiments. The self-programming curve is shown in blue; the experimental and COMSOL curves are in orange and green, respectively. The numerical simulation and experimental results show that the focus acoustic pressure decreased with increasing acrylic plate thickness. The focus acoustic pressure in self-programming fell by 81.13% compared with 78.84% in experimental, and by 83.04% in COMSOL. According to the fitting formula, we calculated the focus acoustic attenuation value for a 1mm increase in thickness. Compared with the initial acoustic pressure, the focus sound pressure will be reduced by 0.357 dB, 0.378 dB, and 0.353 dB, respectively.

Curvature Radius
We then studied the influence of an acrylic plate curvature radius on the propagation of ultrasonic waves. The distance of the acrylic plate (excitation power: 1 w, T: 8 mm) from the open end of the transducer (D) is 10 mm and its curvature radius grew from 60 mm to 100 mm. Figure 6b shows the linear fitting curves of the three results. The numerical simulation and experimental results show that the focus acoustic pressure decreased with increasing acrylic plate curvature radius. The focus acoustic pressure in self-programming fell by 12.41% compared with 14.83% in experimental, and by 5.24% in COMSOL. From the fitting results, the focus acoustic pressure will decrease by 0.027 dB, 0.032 dB, and 0.011 dB for every 1 mm increase in curvature radius, respectively.

Distance to Transducer
In the next step, the influence of an acrylic plate distance to transducer has been regarded as a research target in this part. We use a three-axis stage to change the acrylic plate (excitation power: 1 w, T: 8 mm, R: 60 mm) distance from the open end of the transducer (D: 11 mm-20 mm). Figure 6c shows the linear fitting curves of the three results. The numerical simulation and experimental results show that the focus acoustic pressure decreased with increasing distance to the transducer. The focus acoustic pressure in self-programming fell by 5.53% compared with 10.23% in experimental, and by 2.15% in COMSOL. According to the fitting formula, the focus acoustic pressure will decrease by 0.048 dB, 0.089 dB, and 0.019 dB for every increase of 1 mm between the acrylic plate and transducer.  Table 3, which mainly include the coefficient of determination (R-square), the fitting formula, the correlation coefficient between fitting curves, and the attenuation of unit distance. X is the independent variable, representing the variation of thickness, radius of curvature, and distance from the transducer, respectively. Y is the focus acoustic pressure. In addition, the higher the Pearson correlation coefficient, the higher the degree of agreement between the curve and the self-programming fitting curve. It can be seen from the slope of the fitted curve and the attenuation results per unit distance that the importance of the impact is:

The Influence of the Skull Model Acoustic Parameter
We then explore the influence of the skull model acoustic parameters on tFUS based on self-programming and COMSOL. Additionally, corresponding structural parameters are set as excitation power1 w, T8 mm, R60 mm, D10 mm. First, we change all three  Table 3, which mainly include the coefficient of determination (R-square), the fitting formula, the correlation coefficient between fitting curves, and the attenuation of unit distance. X is the independent variable, representing the variation of thickness, radius of curvature, and distance from the transducer, respectively. Y is the focus acoustic pressure. In addition, the higher the Pearson correlation coefficient, the higher the degree of agreement between the curve and the self-programming fitting curve. It can be seen from the slope of the fitted curve and the attenuation results per unit distance that the importance of the impact is: (1) thickness; (2) distance to the transducer; (3) radius of curvature. The self-programming results show that the attenuation of the thickness of the skull model (0.357 dB/mm) is about 7.44 times the distance to the transducer (0.048 dB/mm) and about 13.22 times the radius of curvature (0.027 dB/mm).

The Influence of the Skull Model Acoustic Parameter
We then explore the influence of the skull model acoustic parameters on tFUS based on self-programming and COMSOL. Additionally, corresponding structural parameters are set as excitation power 1 w, T8 mm, R60 mm, D10 mm. First, we change all three parameters of the linear interval simultaneously. Then, to evaluate the effect of varied acoustic parameters, we vary only one parameter at a time, while keeping others constant.

Comprehensive Influence of Three Parameters
In the first part, we take the values of density ρ (1700:30:2000 kg/m 3 ), acoustic speed c (2000:300:5000 m/s), and absorption coefficient α (200:480:5000 dB/m) in linear intervals, and change these three acoustic parameters simultaneously. Figure 7a,c are selfprogramming and COMSOL acoustic pressure fields and the variation trend of focus acoustic pressure under different parameters are displayed in Figure 7b,d. Figure 7e shows the histogram of the change in the focus acoustic pressure for the two simulation methods under different acoustic parameters. Our results indicate that the declining trend of the focus acoustic pressure gradually becomes smooth as three parameters are increased at the same time. It needs to be noted here that it is difficult to form a distinguishable focal region after ρ 1880 kg/m 3 , c 3800 m/s, and α 3080 dB/m. At this point, the self-programming acoustic pressure decrease by 94.41%, and COMSOL decline by 94.68%. To determine which acoustic parameter is most influential for the determination of the attenuation of the acoustic pressure, we have varied parameters one by one. parameters of the linear interval simultaneously. Then, to evaluate the effect of varied acoustic parameters, we vary only one parameter at a time, while keeping others constant.

Comprehensive Influence of Three Parameters
In the first part, we take the values of density ρ (1700:30:2000 kg/m 3 ), acoustic speed c (2000:300:5000 m/s), and absorption coefficient α (200:480:5000 dB/m) in linear intervals, and change these three acoustic parameters simultaneously. Figure 7a,c are self-programming and COMSOL acoustic pressure fields and the variation trend of focus acoustic pressure under different parameters are displayed in Figure 7b,d. Figure 7e shows the histogram of the change in the focus acoustic pressure for the two simulation methods under different acoustic parameters. Our results indicate that the declining trend of the focus acoustic pressure gradually becomes smooth as three parameters are increased at the same time. It needs to be noted here that it is difficult to form a distinguishable focal region after ρ 1880 kg/m 3 , c 3800 m/s, and α 3080 dB/m. At this point, the self-programming acoustic pressure decrease by 94.41%, and COMSOL decline by 94.68%. To determine which acoustic parameter is most influential for the determination of the attenuation of the acoustic pressure, we have varied parameters one by one.

Influence of a Single Coefficient: Density
In the second part, skull density ρ is increased from 1700 to 2000 kg/m 3 at 30 kg/m 3 intervals, while maintaining the fixed acoustic speed c 3500 m/s and absorption coefficient α 2600 dB/m. Figure 8a,c are self-programming and COMSOL acoustic pressure fields and the variation trend of focus acoustic pressure under different densities are displayed in Figure 8b,d. Figure 8e,f show the fitted curves and histograms of the change in focus acoustic pressure for the two simulation methods at different density conditions. The results show that with the increase of density, the focus acoustic pressure shows a slight downward trend. The self-programming acoustic pressure decrease by 3.12% in total, and 13.02% in the COMSOL. From the linear fitting diagram, we calculated the focus acoustic attenuation value for every 30 kg/m 3 increase in density. Compared with the initial acoustic pressure, the focus sound pressure will be reduced by 0.024 dB, and 0.115 dB, respectively.

Influence of a Single Coefficient: Density
In the second part, skull density ρ is increased from 1700 to 2000 kg/m 3 at 30 kg/m 3 intervals, while maintaining the fixed acoustic speed c 3500 m/s and absorption coefficient α 2600 dB/m. Figure 8a,c are self-programming and COMSOL acoustic pressure fields and the variation trend of focus acoustic pressure under different densities are displayed in Figure 8b,d. Figure 8e,f show the fitted curves and histograms of the change in focus acoustic pressure for the two simulation methods at different density conditions. The results show that with the increase of density, the focus acoustic pressure shows a slight downward trend. The self-programming acoustic pressure decrease by 3.12% in total, and 13.02% in the COMSOL. From the linear fitting diagram, we calculated the focus acoustic attenuation value for every 30 kg/m 3 increase in density. Compared with the initial acoustic pressure, the focus sound pressure will be reduced by 0.024 dB, and 0.115 dB, respectively.

Influence of a Single Coefficient: Acoustic Speed
The third part is the effect of acoustic speed, skull acoustic speed c is increased from 2000 to 5000 m/s at 300 m/s intervals, while maintaining the fixed density ρ 1850 kg/m 3 and absorption coefficient α 2600 dB/m. Figure 9a,c are self-programming and COMSOL acoustic pressure fields and the variation trend of focus acoustic pressure under different acoustic speed are displayed in Figure 9b,d. Figure 9e,f show the fitted curves and histograms of the change in focus acoustic pressure for the two simulation methods at different density conditions. The results showed that with the increase of acoustic speed, the focus acoustic pressure shows a significant downward trend. The self-programming sound pressure decrease by 25.01% in total, and 50.14% in the COMSOL. From the linear fitting diagram, we calculated the focus acoustic attenuation value for every 30 m/s increase in acoustic speed. Compared with the initial acoustic pressure, the focus sound pressure will be reduced by 0.021 dB, and 0.055 dB, respectively.

Influence of a Single Coefficient: Acoustic Speed
The third part is the effect of acoustic speed, skull acoustic speed c is increased from 2000 to 5000 m/s at 300 m/s intervals, while maintaining the fixed density ρ 1850 kg/m 3 and absorption coefficient α 2600 dB/m. Figure 9a,c are self-programming and COMSOL acoustic pressure fields and the variation trend of focus acoustic pressure under different acoustic speed are displayed in Figure 9b,d. Figure 9e,f show the fitted curves and histograms of the change in focus acoustic pressure for the two simulation methods at different density conditions. The results showed that with the increase of acoustic speed, the focus acoustic pressure shows a significant downward trend. The self-programming sound pressure decrease by 25.01% in total, and 50.14% in the COMSOL. From the linear fitting diagram, we calculated the focus acoustic attenuation value for every 30 m/s increase in acoustic speed. Compared with the initial acoustic pressure, the focus sound pressure will be reduced by 0.021 dB, and 0.055 dB, respectively.

Influence of a Single Coefficient: Acoustic Speed
The third part is the effect of acoustic speed, skull acoustic speed c is increased from 2000 to 5000 m/s at 300 m/s intervals, while maintaining the fixed density ρ 1850 kg/m 3 and absorption coefficient α 2600 dB/m. Figure 9a,c are self-programming and COMSOL acoustic pressure fields and the variation trend of focus acoustic pressure under different acoustic speed are displayed in Figure 9b,d. Figure 9e,f show the fitted curves and histograms of the change in focus acoustic pressure for the two simulation methods at different density conditions. The results showed that with the increase of acoustic speed, the focus acoustic pressure shows a significant downward trend. The self-programming sound pressure decrease by 25.01% in total, and 50.14% in the COMSOL. From the linear fitting diagram, we calculated the focus acoustic attenuation value for every 30 m/s increase in acoustic speed. Compared with the initial acoustic pressure, the focus sound pressure will be reduced by 0.021 dB, and 0.055 dB, respectively.

Influence of a Single Coefficient: Absorption Coefficient
The final portion of the acoustic parameter is the impact of the absorption coefficient. Skull absorption coefficient α is increased from 200 to 5000 dB/m at 480 dB/m intervals, while maintaining the fixed density ρ 1850 kg/m 3 and acoustic speed c 3500 m/s. Figure 10a,c are self-programming and COMSOL acoustic pressure fields and the variation trend of focus acoustic pressure under different absorption coefficients are displayed in Figure 10b,d. Figure 10e,f show the fitted curves and histograms of the change in focus acoustic pressure for the two simulation methods at different density conditions. The most spectacular downward trend is observed for acoustic pressure as the absorption coefficient increases and it is difficult to form a distinguishable focal region after α 3080 dB/m. At this point, the self-programming acoustic pressure decrease by 93.66%, and COMSOL decline by 89.64%. The manuscript calculates the slope of each segment at an interval of 480 dB/mm based on the fitting formula, as shown in Figure 10e. The results show that after α 3080 dB/m, the slope change becomes significantly smaller and the attenuation of sound pressure tends to be flat.

Influence of a Single Coefficient: Absorption Coefficient
The final portion of the acoustic parameter is the impact of the absorption coefficient. Skull absorption coefficient α is increased from 200 to 5000 dB/m at 480 dB/m intervals, while maintaining the fixed density ρ 1850 kg/m 3 and acoustic speed c 3500 m/s. Figure  10a,c are self-programming and COMSOL acoustic pressure fields and the variation trend of focus acoustic pressure under different absorption coefficients are displayed in Figure  10b,d. Figure 10e,f show the fitted curves and histograms of the change in focus acoustic pressure for the two simulation methods at different density conditions. The most spectacular downward trend is observed for acoustic pressure as the absorption coefficient increases and it is difficult to form a distinguishable focal region after α 3080 dB/m. At this point, the self-programming acoustic pressure decrease by 93.66%, and COMSOL decline by 89.64%. The manuscript calculates the slope of each segment at an interval of 480 dB/mm based on the fitting formula, as shown in Figure 10e. The results show that after α 3080 dB/m, the slope change becomes significantly smaller and the attenuation of sound pressure tends to be flat.  Table 4 shows the fitting results under the three conditions. Where x is the independent variable, representing the variation of density, acoustic speed, and absorption coefficient, respectively. Y is the focal sound pressure. The results show that the effect of the acoustic absorption coefficient on acoustic parameters is the most critical. Among them, the self-programming results show that the attenuation caused by the acoustic speed of the skull model (0.021 dB/mm) is basically the same as the attenuation caused by the density (0.024 dB/mm).  Table 4 shows the fitting results under the three conditions. Where x is the independent variable, representing the variation of density, acoustic speed, and absorption coefficient, respectively. Y is the focal sound pressure. The results show that the effect of the acoustic absorption coefficient on acoustic parameters is the most critical. Among them, the selfprogramming results show that the attenuation caused by the acoustic speed of the skull model (0.021 dB/mm) is basically the same as the attenuation caused by the density (0.024 dB/mm).

Discussion
In this study, we have analyzed the effects of skull model structural and acoustic parameters on tFUS based on simulation and experiment. These skull parameters will vary with the patient's age, physiological state, and ultrasound incident acoustic window. For variation tendency and the importance ranking, the self-programming results are consistent with the experimental and COMSOL results. To an extent, with the increase of skull model thickness, radius of curvature, distance to transducer, density, acoustic speed, and absorption coefficient, the focus acoustic pressure of the self-programming, experiment, and COMSOL all showed a decreasing tendency of varying degree. Specifically, as the thickness, radius of curvature, and distance to transducer increase by 1mm, the self-programming results show that the focus acoustic pressure drops by 0.357 dB, 0.027 dB, and 0.048 dB, respectively. As the density and acoustic speed increase by 30 kg/m 3 and 30 m/s, the self-programming results show that the focus acoustic pressure drops by 0.024 dB and 0.021 dB, respectively. In addition, as the absorption coefficient increases at an interval of 480 dB/m, the focus acoustic pressure drops from sharp to gentle. Our results indicate that the thickness and absorption coefficient play the most important role in the skull model structural and acoustic parameters.
Nevertheless, we are still seeing some reasonable differences between self-programming simulation, experiment, and COMSOL simulation. After comprehensive consideration, potential major reasons for the observed acoustic pressure field differences between selfprogramming and COMSOL include the following three considerations. First of all, the acoustic field calculation formulas of the two simulations are different, one is Westervelt and the other is Helmholtz. Westervelt equation is used to describe nonlinear propagation in many cases, while Helmholtz equation is mainly used to describe linear propagation. Secondly, the acoustic wave propagation equation of self-programming additionally considers nonlinear factors. The acoustic phenomena are fundamentally nonlinear, and can be described by an approximate linearity when the ultrasound amplitude energy is low, and the nonlinear effect becomes progressively larger as the ultrasound amplitude increases. At last, for the numerical calculation method of the acoustic wave propagation equation, self-programming uses the FDTD method, while COMSOL uses the finite element method (FEM). Both FDTD and FEM are electromagnetic wave simulation algorithms. FDTD is faster in calculation and has average accuracy, FEM is average in calculation and has higher accuracy. About the unavoidable differences between experiments and numerical simulations, these are likely caused by simulation parameters that could be mismatched with respect to the true parameters. Besides, the energy inhomogeneity of the transducer surface is ignored here in the numerical simulation. And due to the spherical shape of the transducer surface. There is a non-right angle between the hydrophone and the transducer, and the pointing of the hydrophone may have an effect on the measurement results.
In addition, this study uses homogeneous acrylic plates to simulate the human skull (mainly cortical bone), ignoring the trabecular bone structure. A smooth and ideal homogeneous acrylic plate may not have the problems of scattering, aberration, and sound field distortion caused by the heterogeneous structure of the real skull, which is different from the real situation. In previous studies, we have used self-programming simulation to build a head model based on real human skull CT data, and further simulated the propagation of tFUS in the skull and brain tissue. The main purpose of this article is to explore the influence of skull model parameters on tFUS one by one, and to screen out the most influential parameters, so it is relatively reasonable to use acrylic board here. With the development of phase-controlled transducer technology, researchers have been able to achieve precise intracranial focusing. However, it generally requires a time-reversal method to obtain the phase information of each element, and the phased array system is generally expensive to manufacture. When exploring the influence of a certain parameter, the advantages of phase-controlled transducers are not so obvious. In this paper, based on a single-element transducer, quantitatively and systematically explored the impact of changes in skull parameters on tFUS due to the patient's age, physiological state, and ultrasound incident acoustic window.
In future research, we will establish the patient's heterogeneous head model based on its head CT (Computed Tomography)/MRI (Magnetic Resonance Imaging) data, and further explore the impact of heterogeneous skull parameter changes on tFUS, especially the problems of scattering and aberration. Related research will help the application of tFUS in neuromodulation research and treatment.

Conclusions
The current study comprehensively explores the effects of the structural and acoustic parameters of the skull model on tFUS. It contains one experimental study and two simulation studies. The results show that the focus acoustic pressure showed different degrees of decreasing trend with the increase of structural parameters or acoustic parameters. Both simulation methods can accurately fit the acoustic pressure variation trend, and the experimental and simulation fitting curves are highly correlated. This study provides a reference for the application of tFUS in neuromodulation and therapy.