The temperature distribution was calculated through the Pennes bioheat conduction equation written as (Pennes, 1998; Hallaj and Cleveland, 1999):

(9) ρ ⁢ C r ⁢ ∂ ⁡ T / ∂ ⁡ t = r ⁢ ∇ 2 ⁡ T + Q - W B ⁢ C B ⁢ ( T - T 0 )

where C_r [J⋅(kg°C)^–1] and r [W⋅(m°C)^–1] are the specific heat and the thermal conductivity of the medium, respectively, T is the transient temperature of the acoustic medium, T₀ is the initial temperature and set as 37°C in the simulation, Q is the volumetric energy loss which is equal to 2αI, where I=1tp⁢∫0tpp22⁢ρ⁢c⁢dt and t_p is the acoustic wave period, W_B is the blood perfusion rate and C_B is the heat capacity of the blood. Other constant parameters used in the simulation are shown in Table 3.

The Montalibet theoretical equation is (Montalibet et al., 2001):

(10) J = σ ⁢ v z ⁢ B x ⁢ s ⁢ i ⁢ n ⁢ ( ω ⁢ t - ϕ ) / ( 1 + t ⁢ a ⁢ n 2 ⁢ ϕ )

where J is the current density; σ is the conductivity of the water, scalp and skull, and the σ of cerebrospinal fluid and brain tissue are 1.0, 0.33, 0.042, 1.0, and 0.33 S⋅m^–1, respectively; v_z is the proton vibration velocity; ϕ is the time constant and B_x is the intensity of the static magnetic field perpendicular to the direction of ultrasonic propagation. Because v_z=p/ρc, the cell level tanφ and φ are femtoseconds and can be ignored. In summary, the Montalibet theoretical equation can be transformed into:

(11) J ≈ σ ⁢ B ⁢ p / ( ρ ⁢ c )

Establish single neuron models of the STN, GPe, GPi and Th based on the H-H model. The BG-Th neural network model is constructed as shown in Figure 3 based on the structured sparse connection method, where I_{app_GPe}, I_{app_STN,} and I_{app_GPi} are the input currents of other nuclei of the brain to the GPe, STN and GPi, respectively, and I_{app_SM} is the irregular pulse input current of the sensorimotor cortex to the thalamus. The black, yellow, light blue, and green dots are the Gpe, STN, GPi, and Th neurons, respectively. The red line arrow is the excitatory input current, the blue line point is the inhibitory input current, and the dashed box is a nerve nucleus. The BG-Th neural network model structure adopted from prior study (Rubin and Terman, 2004; So et al., 2012). The synaptic connection of the BG-Th neural network is that each STN neuron excites two GPe and two GPi neurons through excitatory synaptic connections. Each GPe neuron inhibits two GPe, two STN and two GPi neurons through inhibitory synaptic connections. Each GPi neuron inhibits one Th neuron, and each Th neuron receives excitability information from the sensorimotor cortex. Each neuron is assumed to be a spherical cell existing within an isotropic medium. Coupling of ultrasound field and static magnetic field on charged ions in nerve tissues jointly generates current I_TMAS. This current has been used in studies of neural tissue stimulation (Yuan et al., 2016) and also in simulation studies of H-H neuron models (Lu et al., 2020). The mathematical model of each neuron in the BG-Th neural network model is shown in follow:

(12) C m ⁢ d ⁢ v G ⁢ P ⁢ e d ⁢ t = - I L - I N ⁢ a - I K - I T - I C ⁢ a - I A ⁢ H ⁢ P - I G ⁢ P ⁢ e → S ⁢ T ⁢ N + I a ⁢ p ⁢ p ⁢ _ ⁢ G ⁢ P ⁢ e

(13) C m ⁢ d ⁢ v G ⁢ P ⁢ i d ⁢ t = - I L - I N ⁢ a - I K - I T - I C ⁢ a - I A ⁢ H ⁢ P + I S ⁢ T ⁢ N → G ⁢ P ⁢ i - I G ⁢ P ⁢ e → G ⁢ P ⁢ i + I a ⁢ p ⁢ p ⁢ _ ⁢ G ⁢ P ⁢ i

(14) C m ⁢ d ⁢ v S ⁢ T ⁢ N d ⁢ t = - I L - I N ⁢ a - I K - I T - I C ⁢ a - I A ⁢ H ⁢ P + I S ⁢ T ⁢ N → G ⁢ P ⁢ e - I G ⁢ P ⁢ e → G ⁢ P ⁢ e - I a ⁢ p ⁢ p ⁢ _ ⁢ S ⁢ T ⁢ N + I T ⁢ M ⁢ A ⁢ S

(15) C m ⁢ d ⁢ v T ⁢ h d ⁢ t = - I L - I N ⁢ a - I K - I T - I G ⁢ P ⁢ i → T ⁢ h + I S ⁢ M

where C_m = 1 μF⋅μm^–2 is the membrane capacitance, v_GPe, v_GPi, v_STN and v_Th are the membrane potentials of the GPe, GPi, STN and Th neurons, respectively. I_L, I_Na, I_K, I_T, I_Ca and I_AHP are the leak current, sodium current, potassium current, low-threshold T-type calcium current, high threshold calcium current, and after hyper polarization K⁺ current. I_SM is the cortical sensorimotor signal received by the thalamus and can be described by:

(16) I S ⁢ M = { A S ⁢ M    a / f S ⁢ M < t < a / f S ⁢ M + T S ⁢ M 0    o ⁢ t ⁢ h ⁢ e ⁢ r ⁢ s , 0 < t < T A , a = 1 , 2 , 3 ⁢ …

where A_SM = 3.5 pA μm^–2 is the amplitude, T_SM = 5 ms is the pulse width, and f_SM is instantaneous frequency. The instantaneous frequencies of the incoming pulses follow a γ distribution with an average rate of 14 Hz and a coefficient of variation of 0.2 (So et al., 2012). I_{app_i} (i∈{GPe, STN, GPi}), the positive constant bias currents, which can be viewed as the net synaptic input to these nuclei from other brain regions, I_{app_STN}, I_{app_GPe} and I_{app_GPi} are set at 33, 21 and 22 pA μm^–2 in the healthy state, and 23, 7 and 15 pA μm^–2 in the PD state, respectively. I_TMAS is the TMAS stimulated current. I_m?n (m, n∈{GPe, STN, GPi}) represents the synaptic currents from presynaptic cell m to postsynaptic cell n:

(17) I m → n = g m → n ⁢ ( v m - E m → n ) ⁢ ∑ j S α j

where g_m→n is the maximum synaptic conductance, E_m→n is the synaptic reversal voltage, ∑jSαj represents the overall synaptic conductance of all presynaptic neurons, S_α is synaptic variables, and for I_GPe→GPe, I_GPe→GPi and I_GPe→STN are set as:

(18) dSα/dt=2(−Sα)H∞(v−20)−0.04Sα

(19) H∞(v)=1/(1+exp((−v−57)/2))

where H_∞ is a step function. For I_STN→GPe, I_STN→GPi and I_GPi→Th, S_α is set as:

(20) dSα/dt=zα

(21) dzα/dt=0.234u(t)−0.4zα−0.04Sα

where z_α is the time derivative of the synaptic variable, and if the presynaptic neuron discharge exceeds the threshold of −10 mV at time t, the value of u(t) is 1; otherwise, it is 0. This paper is based on the random function rand to determine the initial value of the neuron membrane potential and the location of the synaptic connections. The specific parameter settings and variable expressions in the model are shown in Supplementary Appendix. All potentials have units of mV, conductance have the unit of mS⋅cm^–2, currents have units of μA⋅cm^–2, concentration have units of mol⋅m^–3, time constants have units of ms. In addition, n′,h′,r′and [Ca]′ are the derivatives of the gating variables, inducing the gating variables to switch between 0 and 1 (So et al., 2012).

This paper first simulates the condition that the input sound intensity of pulsed sinusoidal ultrasonication is 0.3 W⋅cm^–2, the fundamental frequency is 0.5 MHz, the repetition frequency is 10 Hz, the duty cycle is 50%, and the static magnetic field is 0.3 T. Figure 4 shows the sound pressure field, temperature field and TMAS-induced current density distribution at the geometric focus based on the numerical simulation model of ultrasonic transcranial focusing shown in Figure 1 combined with the TR method. Figures 4A–C show the sound pressure field, temperature field and TMAS-induced current density distribution when the sound intensity is 0.3 W⋅cm^–2, respectively, and Figures 4D–F are the change curves of the sound axis pressure, temperature rise and induced current density with input sound intensity. The skull is inside the white hyperbola, and the cross dashed line locates the focal point. From Figure 4, after being calibrated by the TR method, it can be accurately focused transcranially. After 10 s of irradiation, the temperature only rose by 0.40°C and there were no hot spots on the skull. The position of the maximum induced current is consistent with the focal position of the sound pressure field. As the ultrasonic input sound intensity increases, the sound pressure at the focal point gradually increases, and the TMAS-induced current density gradually increases. When the input sound intensity is less than or equal to 2.16 W⋅cm^–2, the temperature rise at the focal point is less than 3°C after 10 s of irradiation, the temperature of the brain tissue is less than 40°C, and there are no hot spots on the skull.

In this paper, the effect of ultrasound parameters on TMAS treatment was investigated by numerical simulation while keeping the pulsed sinusoidal ultrasound fundamental frequency of 0.5 MHz and the magnetic induction strength of static magnetic field is 0.3 T. Take the TMAS-induced current at the focal point to stimulate a single STN neuron as an example. When the duty cycle is 50%, the repetition frequency is 10 Hz, and the STN neuron membrane potential and firing rate with different TMAS input sound intensities are shown in Figure 5. Columns 1–3 in Figure 5 are the curves of the STN membrane potential and I_TMAS with time under the conditions of an input sound intensity of 0.1, 1.0, and 2.0 W⋅cm^–2, respectively. Lines 1–2, 3–4, 5–6, and 7–8 in Figure 5 are the pulse upper half sine, pulse lower half sine, pulse sine and bias pulse sine ultrasonic TMAS current stimulator results, respectively. Figure 5I₁ shows the change curve of the STN discharge rate with ultrasonic input sound intensity. From Figure 5, with the increase of ultrasonic input sound intensity, the STN discharge rate remains basically unchanged after the lower half sine and sine wave stimulation, the STN discharge rate increases after the upper half sine and offset sine wave stimulation, and the STN discharge rate changes dramatically after the bias sine wave stimulation.

Under the condition that the ultrasound duty cycle is 50% and the input sound intensity is 0.1 W⋅cm^–2, the STN neuron membrane potential changes with the ultrasound pulse repetition frequency of the TMAS stimulation as shown in Figure 6. Columns 1–3 in Figure 6 are the curves of the STN membrane potential and I_TMAS with time under ultrasonic pulse repetition frequencies of 5, 10, and 20 Hz, respectively. Figure 6I₁ shows the curve of the firing rate with the repetition frequency of the ultrasonic pulse. From Figure 6, with increasing ultrasonic pulse repetition frequency, the STN discharge rate is basically unchanged the pulse upper sine wave, lower half sine wave and sine wave stimulation. After bias sine wave stimulation, the STN discharge rate increases, but as the repetition frequency increases, the STN discharges. The rate of change is small.

When the ultrasonic pulse repetition frequency is 10 Hz and the input sound intensity is 0.1 W⋅cm^–2, the STN neuron membrane potential changes with the ultrasonic duty cycle of TMAS stimulation as shown in Figure 7. Columns 1–3 of Figure 7 are the curves of the STN membrane potential and I_TMAS with time under duty cycles of 5, 50, and 95%. Figure 7I₁ shows the discharge rate change curve with the ultrasonic duty cycle. From Figure 7, with the increase of the ultrasound duty cycle, the STN discharge rate is basically unchanged after the lower half sine and sine wave stimulation, the STN discharge rate increases after the upper half sine and bias sine wave stimulation, and the STN discharge rate changes dramatically after the bias sine wave stimulation. Since the upper and lower half sines have opposite effects on neuronal membrane potential excitation and inhibition, sine waves have basically no effect on the neuronal membrane potentials at the considered fundamental frequency of 0.5 MHz. The follow-up study is based on biased sinusoidal pulse ultrasound, which can modulate neurons.

The results of this paper show that the TMAS induced current can effectively suppress the PD state under the conditions of magnetic induction strength of 0.3 T and ultrasonic output sound intensity of 0.2 W⋅cm^–2. The ultrasound focal I_SPTA under this condition is 2.41 W⋅cm^–2, which is much smaller than the 5.63 W⋅cm^–2 with ultrasound stimulation in the literature (Folloni et al., 2019). However, when the magnetic induction strength was less than 0.2 T, the ultrasound output sound intensity was stronger than 0.4 W⋅cm^–2, the focal I_SPTA was greater than 4.17 W⋅cm^–2, and the sound intensity I_SPTA was close to 5.63 W⋅cm^–2. Ultrasound can directly cause action potentials in the absence of a magnetic field by several different mechanisms: acoustic radiation force, mechanosensitivity of ion channels, flexoelectricity, cavitation, etc. (Tyler et al., 2008). In this paper, the stimulation effect of induced current is studied according to Montalibet theory, but the direct mechanical effect of ultrasound is not studied. The ultrasound stimulation effect required further analysis.

In this study, when performing ultrasonic waveform analysis in TMAS, it was found that at a fundamental frequency of 0.5 MHz, the upper half-sine and lower half-sine of sinusoidal ultrasound TMAS had excitatory and inhibitory effects on neurons, respectively, which would result in a weak stimulation effect on neurons. Since half-sine is difficult to implement in the actual transducer, and bias sine can be generated by superimposing a pulsed square wave signal on the pulsed sine signal, only half-sine is used for the principal analysis of the smaller sine wave effect, and bias sine wave ultrasound TMAS is chosen to explore the effect of stimulation parameters on the BG-TH neural network in the PD state.

The BG mainly includes the striatum, STN, GPe, GPi, substantia nigra pars compacta (SNc) and substantia nigra reticulata (SNr) nuclei. In this study, we simulated the STN, GPe, and GPi nuclei in BG based on the H-H neuron model. The sensory-motor cortex input was considered as pulsed electrical input and the signals from other brain regions were considered as direct current input, and the enhancement and weakening of synaptic connectivity signals due to the absence of dopaminergic neurons in the SNc in the PD state were simulated by changing this direct current input. Each nucleus pulposus is composed of 10 H-H neurons. Although it was shown by Lu et al. (2020) that this neural network can be used for the exploration of the effects of PD state stimulation parameters, in order to consider more comprehensively the effects of TMAS stimulation on the BG-TH neural network, neuronal models of the striatum, SNc, SNr, and sensorimotor cortex could be introduced in the next step. In addition, there is an outlier at the sampling interval of 1.50 mm in Figure 10E, which may be caused by the small volume of the neural network, and the next step is to construct a large-scale BG-TH neural network and increase the number of neurons contained in each neural cluster. In addition, in the spherical cell model, all induced Montalibet currents will cross the neuronal cell membrane. This approximation may lead to a decrease in the threshold value of the calculated TMAS. The next step will be to consider the loss of induced current in the flow through the cell membrane by constructing multi-compartment neuron model.

In this study, a non-homogeneous head model is constructed based on real human head CT, and the Westervelt acoustic nonlinear propagation equation be used to simulate the sound pressure field, and the Pennes biological heat conduction equation is sampled to simulate the temperature field. Since vascularity and blood flow are not the focus of attention in this paper, these two factors are not considered in the above model, which may have some effect on the distribution of sound pressure/temperature in the acoustic pressure field/temperature field. In this paper, TMAS modeling simulation was performed based on the cranial CT data of only one volunteer, and the next step will be to investigate the effect of cranial parameters on TMAS PD for the cranial CT data of multiple volunteers.