^{1}School of Biomedical Engineering and Technology, Tianjin Medical University, Tianjin, China^{2}Department of Radiotherapy, Cancer Institute and Hospital of Tianjin Medical University, Tianjin, China

**Objective:** Parkinson’s disease (PD) is a degenerative disease of the nervous system that frequently occurs in the aged. Transcranial magnetoacoustic stimulation (TMAS) is a neuronal adjustment method that combines sound fields and magnetic fields. It has the characteristics of high spatial resolution and noninvasive deep brain focusing.

**Methods:** This paper constructed a simulation model of TMAS based on volunteer’s skull computer tomography, phased controlled transducer and permanent magnet. It simulates a transcranial focused sound pressure field with the Westervelt equation and builds a basal ganglia and thalamus neural network model in the PD state based on the Hodgkin-Huxley model.

**Results:** A biased sinusoidal pulsed ultrasonic TMAS induced current with 0.3 T static magnetic field induction and 0.2 W⋅cm^{–2} sound intensity can effectively modulate PD states with RI ≥ 0.633. The magnitude of magnetic induction strength was changed to 0.2 and 0.4 T. The induced current was the same when the sound intensity was 0.4 and 0.1 W⋅cm^{–2}. And the sound pressure level is in the range of −1 dB (the induced current difference is less than or equal to 0.019 μA⋅cm^{–2}). TMAS with a duty cycle of approximately 50% can effectively modulates the error firings in the PD neural network with a relay reliability not less than 0.633.

**Conclusion:** TMAS can modulates the state of PD.

## Introduction

Parkinson’s disease (PD) is a common degenerative neurological disease of the nervous system in middle-aged and elderly people. Its incidence is the second highest among neurodegenerative diseases affecting the elderly. Its prevalence among people over 65 is approximately 1.7%. Its incidence and prevalence both increase with age (Zhang et al., 2005). At present, the main therapies for the treatment of early PD are a combination of multiple anti-PD drugs, nerve nucleus damage surgery, and deep brain stimulation therapy. Drug combination therapy needs to gradually increase the drug dose as the course of the disease progresses, and the corresponding side effects such as dyskinesia are more frequent (Angeli and Merkel, 2008). Surgical treatment of nerve nucleus damage is effective, but it is invasive and irreversible. As the disease progresses, some patients require further treatment due to recurrence or even worsening of symptoms (Knovich et al., 2009). In 2002, the United States Food and Drug Administration approved implantable deep brain stimulation (DBS) therapy. The main stimulation targets of this therapy are the internal globus pallidus (GPi) and the subthalamic nucleus (STN) (Hirsch et al., 1988; Aubry et al., 2003). During treatment, it is necessary to drill a hole in the skull, and puncture the brain to implant electrodes to stimulate the GPi or STN. A slight deviation of the electrode stimulation point may induce depression, and there are risks of bleeding, infection, wire breakage. Currently only 1.6–4.5% of severe PD patients are treated with this modality (Morgante et al., 2007; Strutt et al., 2015). Nonimplantable transcranial magnetic stimulation and direct current stimulation have low spatial resolution and are suitable for stimulation of the motor cortex or superficial nerve tissue (Chris et al., 2019).

In 2003, Norton et al. proposed the transcranial magnetoacoustic stimulation (TMAS) method and conducted a theoretical derivation. This technology combines magnetic and focused ultrasound based on the Hall effect, and uses the induced current generated by the coupling of transcranial focused ultrasound and static magnetic field to excite or inhibit neurons in the target area (Norton, 2003). It has the advantages of noninvasiveness, deep brain stimulation, and high spatial resolution. Yuan et al. (2017) studied the influence of TMAS parameters on neuron desynchronization based on Hodgkin-Huxley (H-H) and Hindmarsh-Rose neuron model simulations. The simulation results show that the cycle is similar to the neuron firing cycle and that the duty cycle is 40–70%. Zhang et al. (2018) studied the effect of TMAS current density on the action potential of excitatory or inhibitory cortical neurons based on the Izhikevich model. The results showed that with increasing current density, the firing interval of cortical neurons decreases and the firing frequency increases (Zhang et al., 2018). Wang H. et al. (2019) used TMAS with different spatial peak time average sound intensities under the condition of a static magnetic field intensity of 0.3 T to stimulate the motor cortex of mice. The results showed that when the sound intensity was 25–144 mW⋅cm^{–2}, the TMAS induced electric field can be reduced. The EMG signal of the mice is induced under the intensity of the ultrasound (Wang H. et al., 2019). In the same year, Wang Y. et al. (2019) conducted experiments in the hippocampus of TMAS PD mice, and the results showed that TMAS could improve neuroplasticity through postsynaptic regulation. Liu et al. (2020) studied the effect of DBS on the PD state basal ganglia-thalamus (BG-Th) neural network based on the Izhikevich model and used the relay reliability index (RI) to evaluate the stimulation results. The larger the RI, the better the stimulation effect. TMAS can effectively stimulate the activity of brain neurons (Liu et al., 2020). At present, only TMAS numerical simulations and experiments on rats and mice without considering the distribution of induced currents have been carried out. They cannot reflect the state of each neuron in the BG-Th neural circuit, and cannot reflect the effective sound field and its generation under static magnetic field conditions. In mouse or rat animal experiments, because the mouse skull is very thin, it has little effect on the ultrasound transcranial focused sound pressure field, whereas the human skull is approximately 2.33–19.08 mm thick and has heterogeneity, which will distort the ultrasound transcranial focused sound pressure field. Unfavorable phenomena such as phase distortion and defocusing occur, which cause distortion of the focus position and the shape of the TMAS induced current (Ding et al., 2015). In 2006, Yand et al. recorded the frequency of EEG in adult mice at 24–42°C for 60 min, and found that brain tissue damage is reversible when the temperature is less than 40°C (Yang and Qi, 2006).

Transcranial magnetic acoustic stimulation is a method of neuromodulation based on Hall effect coupling TMS and FUS to generate induced currents, which may be accompanied by the stimulating effect of ultrasound during TMAS treatment. Verhagen et al. (2019) stimulated motor association cortex and granular prefrontal cortex of six healthy macaques under pulsed ultrasound conditions with a fundamental frequency of 1.4 MHz, a repetition frequency of 10 Hz, a duration of 20 s, and a time-averaged spatial peak acoustic intensity (I_{SPTA}) of 7.2 and 9.5 W cm^{–2}, respectively, and showed that changes in neural activity could be induced. Folloni et al. (2019) stimulated the amygdala and anterior cingulate cortex of 11 healthy macaques with pulsed ultrasound with ISPTA of 19.5 and 5.63 W cm^{–2}, respectively, at a fundamental frequency of 250 KHz, pulse duration of 30 ms, duty cycle of 30%, and duration of 40 s. The results showed an effect on neural activity in the target area, while there was no effect on non-target (Folloni et al., 2019).

Based on CT images of volunteer’s skull, this paper establishes a TMAS numerical simulation model comprised of human skull, 128-element phase-controlled transducer and permanent magnet. Numerical simulations transcranial sound pressure field, and then coupled with the static magnetic field to obtain the TMAS induction electric field distribution. It stimulates the STN in the BG-Th neural network in the PD state by using the induced currents at the focal point and the acoustic axis, and explored the parameters of the sound pressure field by changing the waveform, duty cycle, and repetition frequency of the transcranial focused ultrasound when fundamental frequency is fixed at 500 kHz. The effect on various neurons in the BG-Th neural network was evaluated, and screening of effective parameters that can make Th respond normally to cortical control motor behavior signals was conducted.

## Models and Methods

The numerical simulations were all performed on a Lenovo Think Station D30 workstation (Lenovo Group Ltd., Beijing, China) with an Nvidia TitanX GPU (NVIDIA Corporation, Santa Clara, CA, United States). Simulations were performed based on computer programming using CUDA C on the platform of Visual Studio Community 2013 (Microsoft Corporation, Redmond, Washington, United States). The simulation of the BG-Th model was carried out in MATLAB (The Mathworks Inc., Natick, MA, United States).

### Sound Pressure Field

#### Numerical Simulation Model

Volunteer’s head CT data (49-year-old male, scanning parameters were 120 kV and 100 mA, scanning thickness was 3 mm) was used to establish a human head, 128-element concave spherical phase-controlled transducer. The numerical simulation model of TMAS transcranial focusing comprised of water and permanent magnets is shown in Figure 1. Among them, the regular distribution array transducer has an opening diameter of 112 mm, a radius of curvature of 86 mm, an element radius of 4 mm, an operating frequency of 0.7 MHz, and a focal depth of 66 mm, which can reach the positions of the STN (Starr, 1999). The numerical simulation area is 112 mm × 112 mm × 100 mm, and the acoustic axis is the z-axis. The spatial step size of the numerical simulation model is dx = dy = dz = 0.25 mm, and the time step size is dt = 10 ns. The boundary of the model is processed by the Mur first-order boundary absorption condition.

**Figure 1.** Numerical simulation model of transcranial focusing of a concave spherical phased transducer with 128 elements.

#### Sound Wave Equation

The Westervelt acoustic wave nonlinear propagation equation is (Westervelt, 1963; Kang et al., 2008),

where ∇^{2} is the Laplacian, *p* (Pa) is the acoustic pressure, *ρ* (kg⋅m^{–3}) and c (m⋅s^{–1}) are density and sound velocity of the acoustic medium, respectively, and δ = 2*c*^{3}α/*ω*^{2} is the acoustic diffusion coefficient where α (dB⋅mm^{–1}) is the acoustic attenuation coefficient, β is the nonlinear coefficient, ω = 2πf (rad⋅s^{–1}) is the angular frequency, f is the drive frequency of the transducer and t is the irradiation time.

In this paper, the parameters of the skull and brain tissue such as density (ρ), sound speed (c) and attenuation coefficient (α) were obtained from the bone porosity (φ) converted from the Hounsfield unit (H) of the CT images and the calculation method was as follows (Aubry et al., 2003):

where the *ρ _{bone}*,

*c*

_{bone}and

*α*are the density, speed of sound and the attenuation of cortical skull bone, respectively, and

_{bone}*ρ*,

_{water}*c*

_{water}and

*α*are density, speed of sound, and the attenuation of water, respectively. Other constant parameters used in the simulation are shown in Table 1.

_{water}#### Element Driving Signals

Based on the time reversal (TR) method (Ding et al., 2015), place the pulse wave virtual point sound source *S*_{0}(*t*) at the target focus F, as shown in Figure 2,

where *p(t)* = *p*_{0}sin(2π*ft*) is the instantaneous sound pressure, *p*_{0} is the sound pressure amplitude, *T*_{1} is the sine pulse width, *T*_{s} is the repetition period, *T*_{1}/*T*_{s} is the duty cycle, *T*_{A} is the irradiation duration, *b* is the bias parameter, and the pulse function *y*(*t*) is:

The bias parameter b and the pulse signal *y(t)* jointly determine the pulse ultrasonic wave shape, see Table 2 for details.

Figures 2A,B show the processes of driving signal acquisition and focusing. Point F was set as the focal targets. Sequentially record the sound pressure signal, *p*_{i}(*t*-*T*_{p}), of the sound wave emitted by the point sound source propagating to element i. Reverse each element according to the time series T_{p} at different times and obtain the TR signal *p*_{i}(*T*_{P}-*t*) of each element in the transducer. T_{p} is a time delay sequence containing head information. The relative initial phase delay Δt_{i} of p_{i} (T_{P}-t) over a period of time is calculated using the least squares function fitting method, and then the sinusoidal signal amplitude is modulated with the same input sound intensity. The excitation signal of each element after phase adjustment using the numerical fitting method of TR is:

### Temperature Field

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

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_{0} 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=\frac{1}{{t}_{p}}{\int}_{0}^{{t}_{p}}\frac{{p}^{2}}{2\mathrm{\rho}c}\mathit{d}t$ 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.

### Electric Field

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

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:

### Basal Ganglia-Thalamus Neural Network

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:

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:

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*:

where *g _{m}_{→}_{n}* is the maximum synaptic conductance,

*E*is the synaptic reversal voltage, ${\sum}_{j}{S}_{\mathrm{\alpha}}^{j}$ represents the overall synaptic conductance of all presynaptic neurons,

_{m}_{→}_{n}*S*

_{α}is synaptic variables, and for

*I*,

_{GPe→GPe}*I*and

_{GPe→GPi}*I*are set as:

_{GPe→STN}where *H*_{∞} is a step function. For *I _{STN→GPe}*,

*I*and

_{STN→GPi}*I*,

_{GPi→Th}*S*

_{α}is set as:

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 [*C**a*]′ are the derivatives of the gating variables, inducing the gating variables to switch between 0 and 1 (So et al., 2012).

### Evaluation of the Stimulation Effects

#### Firing Rates

By counting the number of spikes for each neuron in a nucleus, we can obtain an average rate, *RF*, for every nucleus in the BG-Th neural network,

where*fires* is the total number of firings of the neuron during *T*_{A}.

#### Thalamic Relay Fidelity

The reliability index (*RI*) is another primary index to quantify the Parkinsonian state by measuring the fidelity of Th throughput. *RI* is defined as follows (Rubin and Terman, 2004):

where *n*_{SM} denote the total number of *I*_{SM} pulses, and the *n*_{errors} is the inaccurate responses of thalamic neurons to sensorimotor cortex impulses including missed, delay and burst firing corresponding to no excitation, excitation 5 ms after input, and multiple excitations within 25 ms.

In healthy BG-Th, optimum performance of the Th neural network can be achieved, and the *RI* is equal to 1, when each input pulse from the sensorimotor cortex *I*_{SM} results in a single action potential in each Th neuron. In addition, in the PD state, 0 < *RI* < 1. Due to the randomness of neuron connections and the value of the initial membrane potential in this paper, all evaluation indicators are averaged five times.

## Results

### Transcranial Magnetic Acoustic Stimulation Induction Electric Field

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.

**Figure 4.** The influence of input sound intensity on TMAS, **(A)** sound pressure field, **(B)** temperature field, **(C)** TMAS-induced field, **(D)** sound pressure curve at the sound axis, **(E)** is the temperature rise at the focus and the skull, **(F)** is the maximum density of the TMAS-induced current at the focal point (*B* = 0.3 T, *t* = 10 s).

### Stimulate Single Subthalamic Nucleus

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_{1} 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.

**Figure 5.** The influence of TMAS input sound intensity on a single STN membrane firing. Columns 1–3 ultrasonic input sound intensity *I*_{0} = 0.1, 1.0 and 2.0 W⋅cm^{–}^{2}, rows 1–4 are the upper, lower half sine wave, sine wave and offset sine wave stimulation. The blue and black curves are STN membrane potential and the TMAS-induced current, **(I _{1})** is the curve of STN neuron firing rate versus TMAS input sound intensity (

*f*= 0.5 MHz,

*f*

_{s}= 10 Hz,

*T*

_{1}/

*T*

_{s}= 50%,

*B*= 0.3 T).

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_{1} 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.

**Figure 6.** The effect of TMAS repetition frequencies on a single STN discharge. Columns 1–3 ultrasonic pulse repetition frequencies *f*_{s} = 5, 10 and 20 Hz, rows 1–4 are the upper, lower half sine wave, sine wave and offset sine wave stimulation, the blue and black curves are STN membrane potential and TMAS-induced current, **(I _{1})** is the curve of STN neuron firing rate versus TMAS repetition frequency (

*f*= 0.5 MHz,

*I*

_{0}= 0.1 W⋅cm

^{– 2},

*T*= 50%,

_{1}/T_{s}*B*= 0.3 T).

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_{1} 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.

**Figure 7.** The effect of TMAS duty cycle on the discharge of a single STN. Columns 1–3 Ultrasonic Duty Cycle, *T*_{1}/*T*_{s} = 5%, 50 and 95%, rows 1–4 are the upper, lower half sine wave, sine wave and offset sine wave stimulation, the blue and black curves are STN membrane potential and the TMAS-induced current, **(I _{1})** is the curve of STN neuron firing rate versus TMAS duty cycle (

*f*= 0.5 MHz,

*f*

_{s}= 10 Hz,

*I*

_{0}= 0.1 W⋅cm

^{– 2},

*B*= 0.3 T).

### Stimulate the Subthalamic Nucleus in the Basal Ganglia-Thalamus Neural Network

#### Health and Parkinson’s Disease

The BG-Th neural network is constructed based on the H-H model and the structured sparse connection method. Each nucleus contains 10 neurons and the duration is 1,000 ms. The membrane potential of the 1st neuron in the STN, GPe and GPi nuclei and ten neurons of Th membrane potentials in the BG-Th neural network in healthy and PD states changes with time as shown in Figure 8. Figures 8A–C are the membrane potentials of the STN, GPe, and GPi with time in health and Figures 8D–F in PD states. Figures 8D,H are the curves of the membrane potential of 10 neurons in the Th nucleus in the healthy and PD states. The red, black, and green arrows indicate empty firing, delayed firing, and burst firing, respectively. From Figure 8, in the PD state, the firing times of STN and GPi neurons in the BG-Th neural network increase, while the firing times of GPe neurons decrease. In the PD state, Th neurons cannot respond to *I*_{SM} one by one.

**Figure 8.** The curve of neuron membrane potential over time, single neuron of STN, GPe, GPi and ten neurons of Th membrane potentials under **(A–D)** health and **(E–H)** PD states.

#### Influence of the Ultrasound Parameters

Under the conditions of a pulse repetition frequency of 10 Hz and an input sound intensity of 0.1 W⋅cm^{–2}, the random function rand is used to determine the initial value of the neuron membrane potential and the position of the synaptic connection, and the *RI* value of the Th nucleus after the STN nucleus is stimulated by the TMAS uniform electric field with different ultrasonic duty cycles is shown in Table 4. In Table 4, the first column is the *RI* value corresponding to the original PD state 10 times, and columns 2 to 11 are the *RI* value after TMAS stimulation with different duty cycles. The bold data are larger than the *RI* value in the first column corresponding to the PD state, which means that the stimulus is effective. From Table 4, when the PD *RI* ≥ 0.567, the BG-Th neural network has a corrective effect when using a TMAS stimulating current with a certain duty cycle. When the PD *RI* ≥ 0.667, the TMAS duty cycle is approximately 50%, which can be used for correction in most cases. Follow-up ultrasound parameters were studied for PD status with PD *RI* ≥ 0.667.

**Table 4.** *RI* value after TMAS stimulation with different duty cycles (*f* = 0.5 MHz, *f*_{s} = 10 Hz, *I*_{0} = 0.1 W⋅cm^{–2}, and *B* = 0.3 T).

Under the condition that the duty cycle of the bias pulse sinusoidal ultrasound is 50% and the repetition frequency is 10 Hz, uniform TMAS with different ultrasonic input sound intensities at the geometric focus stimulates the STN nucleus in the BG-Th neural network, and the Th discharge result is shown in Figure 9. Figures 9A–C show the delay, burst and missed firing times of Th nerve nucleus after stimulation of TMAS-induced current with an input sound intensity of 0.1–2.0 W⋅cm^{–2}, and Figure 9D shows the change curve of *RI* with ultrasonic input sound intensity. The point is a single simulation result, the curve is the fitting curve after five averages, and the point line is the average number of various abnormal responses and Th relay reliability *RI* under the five PD states. Figure 9 shows that as the input sound intensity increases, the abnormal response oscillation increases. When the input sound intensity is between 0.1–0.4 W⋅cm^{–2}, the *RI* value after stimulation is greater than the corresponding original PD *RI* value. When the input sound intensity is 0.2 W⋅cm^{–2} with a corresponding maximum induced current density of 22.3 μA⋅cm^{–2}, the *RI* value after stimulation is the largest, and TMAS stimulation has the best effect in improving the abnormal response of PD.

**Figure 9.** The influence of the input sound intensity of the biased ultrasound TMAS on the abnormal discharge of the Th nerve nucleus in the BG-Th neural network. **(A–C)** Number of delay, burst and missed firing, **(D)** *RI* from Th. The points and curves are the average of single and five times after stimulation, the dashed line is the average of five times in the PD state, (*f* = 0.5 MHz, *f*_{s} = 10 Hz, *T _{1}/T_{s}* = 50%,

*B*= 0.3 T).

When the bias pulse sinusoidal ultrasonic input sound intensity is 0.2 W⋅cm^{–2}, the duty cycle is 50%, and the repetition frequency is 10 Hz, the acoustic axis sound pressure or induced current density distribution area is as shown in Figure 10A, and then sampling is performed 10 times around the focus point, which is the highest stimulation current. The embedded diagram in the Figure 10A shows the sound pressure and nonuniform TMAS stimulation current received by ten STN neurons, and the delay firing, burst firing and miss firing of Th firing after stimulation are shown in Figures 10B–D. Figure 10E shows the variation curve of *RI* with ultrasonic input sound intensity. From Figure 10, when the sampling interval is 0.25, 0.50, 0.75, 1.00, 1.25, and 1.50 mm, the maximum current difference corresponding to the induced current of TMAS is 0.003, 0.019, 0.044, 0.077, 0.119, and 0.159 μA⋅cm^{–2}, respectively. As the TMAS-induced current difference increases, the average number of Th nuclei delayed firing, empty firing, and burst firing gradually increases, and the *RI* gradually decreases. When the maximum current difference of the TMAS-induced current is less than or equal to 0.019 μA⋅cm^{–2}, whereas the sampling interval is 0.50 mm, the *RI* value after stimulation is greater than the corresponding original PD *RI* value, and TMAS stimulation can effectively regulate the abnormal response of PD status. When the maximum current difference is less than or equal to 0.019 μA⋅cm^{–2}, the corresponding sound pressure level is −1 dB, and the effective nerve stimulation range is less than the −6 dB sound focal range.

**Figure 10.** The influence of the biased ultrasound distribution on the abnormal discharge of the Th nerve nucleus in the BG-Th neural network. **(A)** is the acoustic axis sound pressure and the TMAS-induced current density after sampling (block diagram), **(B–D)** Number of delay, burst and missed firing, **(E)** *RI* from Th, (*f* = 0.5 MHz, *f*_{s} = 10 Hz, *T _{1}/T_{s}* = 50%,

*I*

_{0}= 0.2 W cm

^{– 2},

*B*= 0.3 T).

#### Influence of the Static Magnetic Field Strength

When the fundamental frequency of the bias pulse sinusoidal ultrasound is 0.5 MHz, the duty cycle is 50%, the input sound intensity is 0.2 W⋅cm^{–2}, the repetition frequency is 10 Hz, the maximum density of the TMAS-induced current is 22.3 μA⋅cm^{–2}, and the current density distributions with different static magnetic field strengths are as shown in Figure 11. Figure 11A shows the change curve of the induced current density with ultrasonic input sound intensity under different magnetic field strength conditions. Figure 11B shows the spatial peak time-averaged sound intensity variation curve with ultrasound output sound intensity at the focus point. Figures 11C–F shows the induced current density distribution under static magnetic field strengths of 0.1, 0.2, 0.3 and 0.4 T. The white curve is the skull, the white cross point is a straight line to locate the focal point, and the black curve is the ultrasonic −1 dB focal range. From Figure 11, to ensure that the maximum induced current density is 22.3 μA⋅cm^{–2}, when the static magnetic field intensity is 0.1, 0.2, 0.3, and 0.4 T, the corresponding ultrasonic input sound intensity is 1.5, 0.4, 0.2, and 0.1 W⋅cm^{–2}, corresponding to an ultrasonic −1 dB focal length The field size is basically the same. With increasing magnetic field strength, the side lobes of the induced current density gradually increase; however, the TMAS-induced current density distribution is basically the same. When the ultrasound input sound intensity is less than 0.5 W⋅cm^{–2}, the focal I_{SPTA} is less than 5.63 W⋅cm^{–2}, and the stimulation effect of ultrasound in TMAS is negligible.

**Figure 11.** The distribution of induced current density under different static magnetic field strengths, **(A)** is the change curve of induced current density with the input sound intensity, **(B)** is the spatial peak time-averaged sound intensity variation curve with output sound intensity, **(C–F)** are the distribution of induced current density (*J* = 22.3 μA⋅cm^{– 2}, *t* = 10 s).

## Discussion

### Ultrasound Mechanical Effects in Transcranial Magnetic Acoustic Stimulation

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.

### The Feasibility of Biased/Rectified Ultrasonic Waveforms

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 Limitations of the Used Basal Ganglia-Thalamus Network

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.

### The Limitations of Acoustic Modeling

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.

## Conclusion

Based on the Westervelt acoustic wave nonlinear equation, this paper constructs the ultrasonic transcranial focused sound pressure field, superimposes the uniform static magnetic field, and obtains the TMAS-induced current distribution according to the Montalibet theory. The construction of the PD state BG-Th neural network model is based on the H-H neuron model. The TMAS-induced current is used to stimulate a single STN neuron and the STN nucleus in a neural network, and parameters such as the ultrasonic waveform, sound strength, duty cycle, and repetition frequency are changed according to the firing state and relay reliability judgment stimulus effect. The results show that when the duty cycle of the ultrasonic wave is approximately 50%, the sound pressure field is −1 dB, the static magnetic field strength is 0.3 T and the ultrasonic input sound strength is 0.2 W⋅cm^{–2}, The magnitude of magnetic induction strength was changed to 0.2 and 0.4 T. The induced current was the same when the sound intensity was 0.4 and 0.1 W⋅cm^{–2}. The TMAS-induced current in the focal range can modulate the PD state with an *RI* not less than 0.633. TMAS could be effective for PD stimulation therapy. These results provide theoretical reference data for the clinically safe and effective treatment of PD by TMAS.

This paper refers to an existing TMAS experiment in rats and mice, and simulates the sound pressure field based on sinusoidal ultrasound before and after modulation and the establishment of the TMAS induction electric field. We propose that when the input sound intensity is between 0.1–0.3 W⋅cm^{–2}, TMAS has a modulation effect on PD, which is consistent with the TMAS mouse experiment results of Wang H. et al. (2019). The stimulation effect of other waves needs to be further studied. This paper is only based on the CT data from one volunteer’s skull to conduct TMAS modeling and simulation. The next step will be to explore the influence of skull parameters on TMAS PD based on the CT data of multiple volunteers’ skulls.

## Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

## Ethics Statement

This study was reviewed and approved by the Ethics Committee of Tianjin Medical University, China. Written informed consent was obtained from the individual(s) for the publication of any potentially identifiable images or data included in this article.

## Author Contributions

YZ designed the study and wrote the manuscript. MZ and ZL collected the relevant literatures. PW provided the CT data. XJ reviewed and edited the manuscript. All authors read and approved the submitted manuscript.

## Funding

This work was partially funded by National Natural Science Foundation of China (Grant No. 81272495) and the Natural Science Foundation of Tianjin (Grant No. 16JC2DJC32200).

## Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

## Acknowledgments

Numerical simulations were all performed on a Lenovo Think Station D30 workstation (Lenovo Group Ltd., Beijing, China) with an Nvidia TitanX GPU (NVIDIA Corporation, Santa Clara, CA, United States).

## Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnins.2021.761720/full#supplementary-material

## References

Angeli, P., and Merkel, C. (2008). Pathogenesis and management of hepatorenal syndrome in patients with cirrhosis. *J. Hepatol.* 48, S93–S103. doi: 10.1016/j.jhep.2008.01.010

Aubry, J. F., Tanter, M., Pernot, M., Thomas, J. L., and Fink, M. (2003). Experimental demonstration of noninvasive transskull adaptive focusing based on prior computed tomography scans. *J. Acoust. Soc. Am.* 113, 84–93. doi: 10.1121/1.1529663

Chris, B., Brem, A. K., Arns, M., Brunoni, A. R., and Bennabi, D. (2019). Repetitive transcranial magnetic stimulation treatment for depressive disorders: current knowledge and future directions. *Curr. Opin. Psychiatry* 32, 1–7. doi: 10.1097/YCO.0000000000000533

Ding, X., Wang, Y., Zhang, Q., Zhou, W., Wang, P., Luo, M., et al. (2015). Modulation of transcranial focusing thermal deposition in nonlinear HIFU brain surgery by numerical simulation. *Phys. Med. Biol.* 60, 3975–3998. doi: 10.1088/0031-9155/60/10/3975

Folloni, D., Verhagen, L., Mars, R. B., Fouragnan, E., Constans, C., Aubry, J. F., et al. (2019). Manipulation of subcortical and deep cortical activity in the primate brain using ranscranial focused ultrasound stimulation. *Neuron* 101, 1109–1116. doi: 10.1016/j.neuron.2019.01.019

Hallaj, I. M., and Cleveland, R. O. (1999). FDTD simulation of finite-amplitude pressure and temperature fields for biomedical ultrasound. *J. Acoust. Soc. Am.* 105, L7–L12. doi: 10.1121/1.426776

Hirsch, E., Graybiel, A. M., and Agid, Y. A. (1988). Melanized dopaminergic neurons are differentially susceptible to degeneration in Parkinson’s disease. *Nature* 334, 345–348. doi: 10.1038/334345a0

Kang, I. L., Sim, I., Kang, G. S., and Min, J. C. (2008). Numerical simulation of temperature elevation in soft tissue by high intensity focused ultrasound. *Mod. Phys. Lett. B* 22, 803–807. doi: 10.1142/S0217984908015413

Knovich, M. A., Storey, J. A., Lan, G. C., Torti, S. V., and Torti, F. M. (2009). Ferritin for the clinician. *Blood Rev.* 113, 95–104. doi: 10.1016/j.blre.2008.08.001

Liu, C., Zhao, G., Wang, J., Wu, H., and Loparo, K. A. (2020). Neural network-based closed-loop deep brain stimulation for modulation of pathological oscillation in Parkinson’s disease. *IEEE Access* 99:1. doi: 10.1109/ACCESS.2020.3020429

Lu, M., Wei, X., Che, Y., Wang, J., and Loparo, K. A. (2020). Application of reinforcement learning to deep brain stimulation in a computational model of Parkinson’s disease. *IEEE Trans. Neural Syst. Rehabil. Eng.* 28, 339–349. doi: 10.1109/TNSRE.2019.2952637

Montalibet, A., Jossinet, J., Matias, A., and Cathignol, D. (2001). Electric current generated by ultrasonically induced Lorentz force in biological media. *Med. Biol. Eng. Comput.* 39, 15–20. doi: 10.1007/BF02345261

Morgante, L., Morgante, F., Moro, E., Epifanio, A., Girlanda, P., Ragonese, P., et al. (2007). How many parkinsonian patients are suitable candidates for deep brain stimulation of subthalamic nucleus? Results of a questionnaire. *Park. Relat. Disord.* 13, 528–531. doi: 10.1016/j.parkreldis.2006.12.013

Norton, S. J. (2003). Can ultrasound be used to stimulate nerve tissue? *Biomed. Eng. Online* 2:6. doi: 10.1186/1475-925X-2-6

Pennes, H. H. (1998). Analysis of tissue and arterial blood temperatures in the resting human forearm. *J. Appl. Physiol.* 85, 5–34. doi: 10.1152/jappl.1998.85.1.5

Rubin, J. E., and Terman, D. (2004). High frequency stimulation of the subthalamic nucleus eliminates pathological thalamic rhythmicity in a computational model. *J. Comput. Neurosci.* 16, 211–235. doi: 10.1023/B:JCNS.0000025686.47117.67

So, R. Q., Kent, A. R., and Grill, W. M. (2012). Relative contributions of local cell and passing fiber activation and silencing to changes in thalamic fidelity during deep brain stimulation and lesioning: a computational modeling study. *J. Comput. Neurosci.* 32, 499–519. doi: 10.1007/s10827-011-0366-4

Starr, P. (1999). Magnetic resonance imaging-based stereotactic localization of the globus pallidus and subthalamic nucleus. *Neurosurgery* 44, 303–313.

Strutt, A. M., Simpson, R., Jankovic, J., and York, M. K. (2015). Changes in cognitive−emotional and physiological symptoms of depression following STN−DBS for the treatment of Parkinson’s disease. *Eur. J. Neurol.* 19, 121–127. doi: 10.1111/j.1468-1331.2011.03447.x

Tyler, W. J., Tufail, Y., Finsterwald, M., Tauchmann, M. L., Olson, E. J., and Majestic, C. (2008). Remote excitation of neuronal circuits using low-intensity, low-frequency ultrasound. *PLoS One* 3:e3511. doi: 10.1371/journal.pone.0003511

Verhagen, L., Gallea, C., Folloni, D., Constans, C., Jensen, D. E., Ahnine, H., et al. (2019). Offline impact of transcranial focused ultrasound on cortical activation in primates. *Elife* 8:e40541. doi: 10.7554/eLife.40541

Wang, H., Zhou, X., Cui, D., Liu, R., Tan, R., Wang, X., et al. (2019). Comparative study of transcranial magneto-acoustic stimulation and transcranial ultrasound stimulation of motor cortex. *Front. Behav. Neurosci.* 13:241. doi: 10.3389/fnbeh.2019.00241

Wang, Y., Feng, L., Liu, S., Zhou, X., Yin, T., Liu, Z., et al. (2019). Transcranial magneto-acoustic stimulation improves neuroplasticity in hippocampus of Parkinson’s disease model mice. *Neurotherapeutics* 16, 1210–1224. doi: 10.1007/s13311-019-00732-5

Yang, F., and Qi, W. (2006). Electron microscopic study of brain tissue and EEG of mouse after acute heat stress Medical. *J. Chin. Peoples Liberation Army* 4, 165–168. doi: 10.1016/S1872-2040(06)60020-0

Yuan, Y., Chen, Y. D., and Li, X. L. (2016). A new brain stimulation method: noninvasive transcranial magneto-acoustical stimulation. *Chin. Phys. B* 25, 226–231. doi: 10.1088/1674-1056/25/8/084301

Yuan, Y., Pang, N., Chen, Y., Wang, Y., and Li, X. (2017). A phase-locking analysis of neuronal firing rhythms with transcranial magneto-acoustical stimulation based on the hodgkin-huxley neuron model. *Front. Comput. Neurosci.* 11:1. doi: 10.3389/fncom.2017.00001

Zhang, S., Cui, K., Zhang, X., Shi, X., Ge, M., Zhao, M., et al. (2018). Effect of transcranial ultrasonic-magnetic stimulation on two types of neural firing behaviors in modified izhikevich model. *IEEE Trans. Magn.* 54, 1–4. doi: 10.1109/TMAG.2017.2773086

Keywords: transcranial magnetoacoustic stimulation, Parkinson’s disease, neural network, Hodgkin-Huxley model, numerical simulation

Citation: Zhang Y, Zhang M, Ling Z, Wang P and Jian X (2021) The Influence of Transcranial Magnetoacoustic Stimulation Parameters on the Basal Ganglia-Thalamus Neural Network in Parkinson’s Disease. *Front. Neurosci.* 15:761720. doi: 10.3389/fnins.2021.761720

Received: 20 August 2021; Accepted: 28 September 2021;

Published: 18 October 2021.

Edited by:

Xiaoli Li, Beijing Normal University, ChinaCopyright © 2021 Zhang, Zhang, Ling, Wang and Jian. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Xiqi Jian, jianxiqi@tmu.edu.cn