Firing activities and magnetic stimulation effects in a Cortico-basal ganglia-thalamus neural network

1.   Introduction

Parkinson's disease (PD) is the second most common neurodegenerative disease in the elderly (after Alzheimer's disease (AD)). According to the research, the number of new cases per year of PD has increased over a 30-year period ^[1]. The main clinical manifestations of PD include motor symptoms such as resting tremor and muscle rigidity and non-motor symptoms such as cognitive disorders and anxiety ^[2]. The significant pathological feature of PD is the degeneration of dopaminergic neurons in the substantia nigra compacta (SNc) ^[3,4], which leads to abnormal $\beta$-band oscillations (13–30 Hz), and excessive synchrony in the basal ganglia (BG) ^[5]. The BG is composed of striatum, SNc, subthalamic nucleus (STN), globus pallidus externa (GPe), and globus pallidus interna (GPi). According to a previous study, STN is currently thought to play a prominent role in the pathophysiology of PD. $\beta$-band oscillations and excessive synchrony in the STN are related to bradykinesia and rigidity in PD ^[6]. Additionally, the firing pattern of other nuclei also varies in the pathological state. In fact, electrophysiologic recordings in monkeys treated with 1-methyl-4-phenyl-1, 2, 3, 6-tetrahydropyridine (MPTP) showed firing rates in the STN increased and GPi but firing rates in the GPe and thalamus (TH) decreased ^[6,7]. Current studies related to the control of abnormal oscillations in PD mainly focus on deep brain stimulation ^[8], optogenetics ^[9] and magnetic stimulaiton ^[10]. With the success of non-invasive magnetic stimulation in the treatment of PD in clinic, its treatment mechanism remains to be investigated.

The motor cortex, which mediates most motor symptoms, is the main output target of the BG. In previous studies, growing activity and bursting also occurred in the motor cortex in the MPTP primate model of PD ^[11]. Meanwhile, the motor cortex plays an input role to STN through the hyperdirect pathway ^[8], that is, the connection from pyramidal neuron (PY) to STN. Leblois et al. ^[12] found that the abnormal oscillation in PD is dominated by the hyperdirect pathway. Pavlides ^[13] further suggested that both resonance and feedback models in the hyperdirect pathway can reproduce anomalous frequency in the animal model of PD. Moreover, Lu ^[14] found that the cortex and hyperdirect pathway contribute to abnormal oscillation in PD.

Clinical and animal models have confirmed that electromagnetic therapy such as transcranial magnetic stimulation (TMS) is effective in the treatment of many neurologic diseases such as PD, AD, and major depressive disorder. On the one hand, TMS can change cortical excitability. The temporary suppression of $\beta$-band activity and desynchronization in the human STN following cortical TMS was reported ^[10]. On the other hand, it can also induce dopamine release in the caudate nucleus ^[15]. As a result, movement disorder in PD can be partly improved after magnetic stimulation treatment. From a micro perspective, magnetic stimulation can change ion concentrations inside and outside of the membrane, which further influences the electrical activity of the nervous system. In addition, electromagnetic induction should be considered in the neuron as the fluctuation in membrane potential can induce changes of the electromagnetic field or the distribution of ion concentrations in the neuron. Ma et al. ^[16] established a magnetic induction neuron model by modeling the magnetic effect as memristor's effect to transfer the magnetic flow into the transmembrane current. Their model could represent the magnetic induction effect induced by the nervous system and external magnetic stimulation. To investigate the underlying mechanism of magnetic stimulation in the treatment of PD, we improve a cortico-basal ganglia-thalamus (CBGT) network considering the magnetic induction. This paper is organized as follows: In Section 2, the model of CBGT network coupled with magnetic induction model is proposed. Section 3 discusses the firing activity and properties in each nucleus under physiological and pathological states. We further explore how the increased synaptic connection in the hyperdirect pathway affects the firing activities in the STN. In Section 4, we investigate the influence of magnetic stimulation on PY. We analyse the statistical results under different parameters of magnetic stimulation. Then, we divide the parameter region into two regions and fit the boundary curves. Section 5 summarizes the results and discusses the perspective of this work.

2.   Model of CBGT neural network

The structure of the CBGT neural network is discussed in Section 2.1. Then, in Section 2.2, we improve the model by introducing the magnetic induction model to consider the magnetic effect on the basis of the CBGT model. The parameters are given, which can represent features in the physiological and pathological states.

2.1.   Structure of CBGT neural network

Computational modeling of the CBGT network can be used to investigate the mechanisms in the magnetic treatment of PD. Plenty models such as neural population model ^[17] and conductance-based neural network model ^[18] have been established to investigate brain activities. Many conductance-based neural models, such as Hodgkin–Huxley model and Hindmarsh–Rose model, are employed in revealing the underlying mechanisms of brain function ^[8]. However, these models are computationally costly in large-scale network modeling. Here, we use a simplified neural model that not only can represent the firing properties, but also is computationally efficient in a large-scale network. In this study, we choose the Izhikevich neural model, which is a simplified two-dimensional model that can reproduce the dynamic properties of electrical activity for classical Hodgkin–Huxley model and is suitable for the large spiking neural network simulation.

Lu ^[14] established a CBGT network model based on sparse network of CBGT to represent PD features. The feed-forward CBGT network consists of six nuclei, namely, inhibitory interneurons (IN), PY, STN, GPe, GPi, and TH. Each nucleus has 100 single-compartment Izhikevich neurons, as shown in Figure 1. Arrows represent synapse connection between nuclei: red arrows indicate excitability synapse connection, whereas yellow arrows represent inhibitory synapse connection. In the cortex area, each IN neuron receives excitability inputs from four PY neurons. Each PY neuron receives inhibitory synaptic inputs from one IN neuron and excitability neuron from two TH neurons. In the basal ganglia, each STN neuron receives inhibitory inputs from two GPe neurons and excitability inputs from two PY neurons through the hyperdirect pathway. Each GPe neuron receives inhibitory inputs from two other GPe neurons and excitability inputs from three STN neurons. Each GPi neuron receives excitability inputs from one STN neuron and inhibitory inputs from one GPe neuron. Each TH neuron receives inhibitory inputs from eight GPi neurons.All input neurons are randomly chosen. The effect of the striatum in direct pathway and indirect pathway is simplifed as a constant current applied on GPe and GPi respectively ^[14]. The hyperdirect pathway from cortex to STN is represented by the connection from PY to IN.

2.2.   Izhikevich neural model under magnetic induction in CBGT network

We use the Izhikevich model as the basic node to simulate neurons in the CBGT network. Electromagnetic radiation could affect the electrical activity of neurons. According to Ma ^[16], electromagnetic induction is caused by changes of intercellular and extracellular ion concentrations. Due to Faraday's law, magnetic induction is a result of fluctuations in membrane potential. Electromagnetic induction can be introduced to the Izhikevich model to adjust the membrane potential via a memristor effect. We introduce Ma's magnetic induction model ^[16] to describe the influence of electromagnetic induction in this study. Thus, the Izhikevich model under magnetic effect can be written as Eq 2.1.

Here, another magnetic flux variable $\phi_i$ is introduced to denote the effects of magnetic induction. Details of modeling theory is presented in ^[16]. ${I_i^{mf}}$ is the induction current induced by the magnetic induction. $k$ is the induction coefficient which is determined by media. To represent the continuous variety of media environment, $k$ is set as $k = 0.1$ for half of the neurons and $k = 0.09$ for the rest in each nucleus. $k_1 = 0.0001$, and $k_2 = 0.01$ are parameters that denote feedback gain of magnetic flux induced by membrane potential and the leakage of magnetic flux respectively. Function $\rho$($\phi_i$) represents the memory conductance value related to the magnetic flux. $\alpha_\phi = 0.1$ and $\beta_\phi = 0.02$ are constants in the magnetic induction model. ${\phi_{ext}}$ is the external magnetic induction composed of the magnetic induction coupling flow and the environment magnetic induction. Here, we do not consider the magnetic coupling within each nucleus. Thus, ${\phi_{ext}}$ only includes the environment magnetic induction if it exists, which is normally regarded as 0 if there is no external electromagnetic field.

the after-spike resetting equation satisfies:

where $i = 1, 2, 3, ..., N(N = 600)$ for each nucleus composed of 100 neurons. $v_i$ is a fast membrane voltage variable, and $u_i$ represents a slow recovery variable of $i$th neuron. $a_i$ and $b_i$ have connection with the time scale of $u_i$ and affect the sensitivity of $v_i$ in subthreshold state. $c_i$ and $d_i$ denote the after-spike reset value of $v_i$ and $u_i$. That is, $a_i$, $b_i$, $c_i$ and $d_i$ determine the neuron type together. $I_i^{bias}$ is the basic current of the nucleus received from other areas in the brain. $\sum\limits_j I_{ij}^{syn}$ is the synapse current from multiple neurons to $i$th neuron in the CBGT network.

In the CBGT network, neurons receive excitability ($\alpha$-Amino-3-hydroxy-5-methyl-4-isoxazolepropionic acid (AMPA) -mediated glutamatergic) and inhibitory (Gamma aminobutyric acid (GABA)-mediated GABAergic) synaptic currents from other neurons. The details of the synaptic connection are shown in Figure 1. The synapse current can be modeled as Eq 2.2, where

The summation is taken over all presynaptic neurons input from neuron $j$ to $i$, and ${s_j}$ is the presynaptic variable. $g_{ij}$ is the strength of synaptic connection (as shown in Figure 1 as yellow and red arrows) from neuron $j$ in nucleus $J$ to neuron $i$ in nucleus $I$ ($I, J\in \left\{PY, IN, STN, GPe, GPi, TH\right\}$), and $g_{ij} = G_{J\to I}$ as listed in Table 2. $E_{type}$ is reversal potential of synapse. $E_{ex} = 0$ is excitability synapse for the synapses from neurons $j$ in PY to neurons $i$ in IN, from neurons $j$ in PY to neurons $i$ in STN, from neurons $j$ in STN to neurons $i$ in GPe, from neurons $j$ in STN to neurons $i$ in GPi, from neurons $j$ in TH to neurons $i$ in PY (as shown in Figure 1 as red arrows), while $E_{in} = -80$ is inhibitory synapse for the synapses from neuron $j$ in IN to neurons $i$ in PY, from neuron $j$ in GPe to neurons $i$ in STN, from neuron $j$ in GPe to neurons $i$ in GPe, from neuron $j$ in GPe to neurons $i$ in GPi, and from neuron $j$ in GPi to neurons $i$ in TH (as shown in Figure 1 as yellow arrows). The values of $\alpha_{j}$ and $\beta_{j}$ listed in Table 1 are the same in both physiological and pathological states for each nucleus.

Referencing the value of model parameters in ^[14], the physiological parameter values in Eq 2.1 are listed in Table 1, and the basic physiological firing patterns of six nuclei within the CBGT network are shown in Figure 2(a). The parameters of GPi and GPe are set as the same value as the similar properties of neurons in globus pallidus, except for the constant current that $I_{GPi}^{bias} = 15$ set greater than $I_{GPe}^{bias} = 11$ (as shown in Table 1) to follow the stronger inputs in direct pathway from striatum. On the one hand, in the pathological state, the loss of dopaminergic neurons in the SNc induces changes of synapse connection strength between each nucleus ^[7,14]. Specifically, the decrease of dopamine in PD induces the increase of excitability projections from STN to GPe and GPi, the increase of excitability projections from PY to STN and excitability projections from TH to PY ^[7,14]. And It causes the increase of inhibitory projections from GPi to TH as well. Besides, the decrease of dopamine reduces intrapallidal inhibition from GPe to GPi and lateral inhibition in the GPe, which mimics a higher level of enkephalin after dopamine decrease, and induces the decrease of inhibitory projections from GPe to STN ^[7,14]. On the other hand, the loss of dopaminergic neurons in the SNc induces firing property changes of neurons in the CBGT network ^[14]. We increase after-spike reset the value $c$ of $v_{GPe}$, $v_{GPi}$ and $v_{TH}$, and increase after-spike reset value $d$ of $u_{GPe}$ and $u_{TH}$, which ensure to reproduce the burst firing appearing in GPi and TH for pathological state as it was analysed by Liu ^[5]. In addition, the decreasing activity of the cortex in PD is mimicked by reducing $I_{PY}^{bias}$. The parameters in the pathological state are listed in Tables 2 and 3. The basic physiological firing patterns of six nuclei in the CBGT network are shown in Figure 2(b). The membrane potentials of individual neurons for the closed-loop CBGT network are shown in Figure 2. In Figure 2(a), the STN, GPe, GPi, TH, and cortex exhibit regular spiking at a nearly constant frequency in the normal physiological condition. Conversely, in the pathological state (Figure 2(b)), bursting occurs in the GPe, GPi, TH. The inter-burst intervals of GPe, GPi and TH are much longer. Furthermore, the STN displays a more irregular spiking of occasional bursting and a higher firing frequency, whereas the PY exhibits a slower bursting firing pattern in the pathological state. To interpret the regularity of spiking, the inter-spike interval of neuron $i$ is defined by $T_{ISI_i}$, and the statistical index of inter-spike interval series of neuron $i$ can be defined as: $CV_i = \sqrt{ < T_{ISI_i}^2 > - < T_{ISI_i} > ^2}/ < T_{ISI_i} >$, which reflects the regularity of neuronal potentials. And the regularity of nucleus can be defined as: $CV = < CV_i >$ for neuron $i$ in STN and PY. The smaller value of CV corresponds to more regular of neuronal spiking. $CV$ of STN and PY in physiological state are 0.050 and 0.022, whereas in pathological state, and $CV$ of STN and PY is larger: 0.354 and 0.090, which reflect more irregular spiking of STN and PY in pathological state. The firing properties of neuron in six nucleus consist with the results in ^[5,14]. Compared with ^[14], inside magnetic electromagnetic radiation of neuron could induce an additional current $I^{mf}$, as shown in Figure 2(c) and (d), causing less spiking of STN in physiological state and a little bursting of STN in pathological state.

3.   Basic state of CBGT neural network

In this section, we demonstrate the firing activities of the CBGT network (Eq 2.1) in the physiological and pathological states and further study the effect of the hyperdirect pathway. The parameters are listed in Section 2. We use Brina2 ^[19] to simulate the network with fourth-order Runge–Kutta algorithm (time step $h = 0.1ms$) in Python 3.7. The initial value of $v_i$ is set as a random number taken from a uniform distribution of [-5, 5], whereas the initial values of $u_i$ and $\phi_i$ are chosen from uniform distribution of [0, 1]. To describe the excessive synchrony in the CBGT neural network, we use Kuramoto order parameter to define the synchronous coefficient (Eq 3.1).

where the magnitude $r(t)\in [0, 1]$ is the order parameter. $r(t)$ measures the coherence of oscillator population, and $i$ is the imaginary unit. $t_{j, k}$ is the $k$th spiking time of $j$th oscillator of the network, and $\theta_j$ is the phase. Thus, the time average order $r$ of the network is $r = {\langle r(t) \rangle}$, which could define the phase synchronization of each nucleus in the CBGT network. That is, $r = 0$ indicates asynchronous, whereas $r = 1$ is fully synchronized.

3.1.   Physiological and pathological states in CBGT network

To further illustrate the dynamic characteristics of each nucleus, the firing activities in the physiological and pathological states are shown in Figures 3, 4 and 5. The parameters in the two states are listed in Tables 1, 2, and 3. Figure 3(a) and (b) describe the spatiotemporal patterns of each neuron in the CBGT network in two states. The nucleus exhibits more coherent oscillations in the pathological state than in the physiological state. The firing patterns are basically determined by individual neuron parameters and slightly different as a result of synapse current injection. Experiments in MPTP-treated monkeys have found that the firing rates increased in the STN and GPi but decreased in the GPe, TH, and cortex. In the physiological state (Figure 5(a)), the average firing rates are 8, 70, 77, 17, 59, and 71 Hz for the STN, GPe, GPi, TH, PY, and IN, respectively. Conversely, in the pathological state (Figure 5(b)), the firing rates are 23, 36,101, 13, 20, and 69 Hz. As shown in Figure 4 and Figure 5(a), the dominant frequency of STN transfers from the $\alpha$-band (8–13 Hz) to the $\beta$-band (13–30 Hz) as a result of PD. The dominant frequency of PY decreases to the $\beta$-band, whereas the peak frequencies of GPe and GPi appear in the $\beta$-band. The model reproduces the typical $\beta$-band oscillation and average firing rates in basal ganglia. This finding indicates that our model in the physiological and pathological states is in agreement with the statistical results of firing frequency in the previous work ^[14].

3.2.   Effect of hyperdirect pathway

The recent study ^[14] suggested that over-activity of the hyperdirect pathway between cortex and STN has played a key role in generating pathological oscillation in PD and following dopamine degeneration, and abnormal oscillations of STN can be driven by cortex. Increased coupling between cortex and STN has also been observed by measuring the functional connectivity in animal models of PD ^[14]. To investigate the effects of the hyperdirect pathway to synchronous oscillations, synchronization, dominant frequency and $\beta$- band power for each nucleus are observed in Figure 6 with the growth of $G_{PY \to STN}$.

In Figure 6, under each set of parameter, we simulate over 10 times as $G_{PY \to STN}$ varying from 0.05 to 1 with step 0.05 to investigate the effect of hyperdirect pathway. Under the parameters of physiological state and with the increasment of $G_{PY \to STN}$, it shows that the increasing synaptic strength of hyperdirect pathway has relatively small effect on the synchronization (Figure 6(a)), dominant frequency (Figure 6(b)) and $\beta$-band power (Figure 6(c)) to the CBGT network. The fluctuation of $r$ within ten simulations is caused by the random factors, ie., random initial value and random sparse connection of the CBGT network. The synchronization of nucleus in BG is less than 0.4, while the dominant frequencies of STN, GPe, GPi are not within the $\beta$-band and corresponding $\beta$-band power is small, which indicate the change of hyperdirect pathway doesn't induce the pathological firing in the CBGT network. However, in the pathological state, the increasing synaptic strength of hyperdirect pathway has a relatively small effect on the dominant frequency and $\beta$-band power of nucleus except for STN in basal ganglia. As shown in Figure 6(d), the synchronization coefficient of STN rising with $G_{PY \to STN}$. And when $G_{PY \to STN} > 0.7$, $r_{STN}$ is around 0.7 and STN represents high synchronization. In Figure 6(e), The dominant frequency of STN rapidly increase from $\alpha$ band to $\beta$-band when $0.1\leq G_{PY \to STN}\leq 0.2$. As shown in Figure 6(f), the $\beta$-band power of STN gradually increase along with $G_{PY \to STN}$. That is, hyperdirect pathway contributes to the abnormal oscillation of STN in PD. And it has relatively small effect on the synchrony of TH, PY, and IN. Differing from the synchrony in the previous work ^[14], as a result of inter-magnetic induction, the synchronous coefficient $r$ of STN has a fluctuation change with the increasement of $G_{PY \to STN}$.

4.   Control of magnetic stimulation in CBGT neural network

Magnetic stimulation such as TMS has been widely used in the treatment of neurological disorders. In this section, we mainly focus on the effect of magnetic stimulation on the CBGT neural network and provide a control strategy.

4.1.   Effects of magnetic stimulation

Here, we adopt the periodic exponentially damped sine function as our magnetic shape, which was used in previous studies due to the fact that TMS induces I-wave and D-wave in the cortex ^[20] (Figure 7).

Applying Ma's model, external magnetic stimulation can be expressed by external magnetic induction $\phi_{ext}$. We model the pulse $\phi_{ext} = H(t)$ as waveform of Eq 4.1 in each period (Figure 7). In Figure 7, the waveform has a stimulation pulse width of 5 ms initially and then rests for another 4 ms in a period.

As magnetic stimulation mainly affects PY with a long axon structure ^[20], we apply $\phi_{ext}$ on PY. Amplitude $A$ in Eq 4.1 is set to $2.5$ and period $T = 25$, whereas $\omega = 0.5$ and $\tau = 1.0$. Figure 8 shows the simulation results of the CBGT network under magnetic stimulation effects in the pathological state. Compared with the corresponding firing activities in PD, external magnetic induction mainly induces the dominant frequency of STN and PY rising from $\beta$-band to $\gamma$-band (40.0 Hz) and successfully reduces the $\beta$-band power of STN and PY green bars in Figure 9(b). Under the pathological state with external magnetic induction (Eq 4.1), the synchronization of STN, GPi and PY could be reduced (Figure 9(a)). And the dominant frequency of STN is no longer in $\beta$-band (Figure 5(a)). Magnetic stimulation nearly has no positive effect on the activity of GPe. Figure 9(a) plots the synchronous coefficient $r$ in the physiological state (blue circles), pathological state (red stars), and magnetic stimulation of the pathological state (green triangles).

In essence, these findings indicate that in the treatment of PD, magnetic stimulation could affect the firing activities directly on the cortex and further control the abnormal oscillations in STN through the hyperdirect pathway.

4.2.   Statistical results of magnetic stimulation

On the basis of analysis in Section 4.1, magnetic stimulation could mainly further influence the firing properties of STN in the BG through the hyperdirect pathway. To explore the relationship between various parameters of external magnetic stimulation and the firing properties of STN, we simulate the pathological CBGT neural network under external magnetic stimulation with the waveform in Figure 7. Amplitude $A$ of external magnetic induction $H(t)$ changes from 0.1 to 5 with a step of 0.1, whereas period $T$ varies from 10 to 110 with intervals of 2. As the randomness of the network structure, simulation under each set of parameters repeats 10 times. In order to further understand the effects of the pulse period and amplitude of magnetic stimulation, the synchronous coefficient $r$, dominant frequency, $\beta$-band power and max power of STN are adopted to identify the controllability of magnetic stimulation.

4.2.1.   $\beta$-band oscillation and synchronization of STN under various magnetic stimulation

Figure 10(a) plots the statistical results of the synchronous coefficient $r$ of STN under various magnetic stimulation. Once the synchronous coefficient is greater than 0.5, we regard that the STN is in an excessive synchronous state related to PD. The yellow and red area surrounded by black contour denotes that $r_{STN}$ is greater than 0.5, which represents the excessive synchronization of STN corresponding to pathological state. To control the excessive synchronization of STN, the parameters of magnetic stimulation could be chosen from the blue areas of Figure 10(a). In addition to the synchronization of STN, external magnetic stimulation also affects the power spectrum density of STN. In Figure 10(b), the sky blue area surrounded by white curves denote dominant frequency in $\beta$-band, indicating that the dominant frequency is within 13–30 Hz. Accompanied with the increase of amplitude $A$ and decrease of period $T$, dominant frequency and $\beta$-band power shows a gradual decline trend in Figure 10(b) and (c). Figure 10(d) is the power of STN in dominant frequency under various magnetic stimulation. It shows that though the dominant frequency is in $\beta$-band, the $\beta$-band power and power of dominant frequency in STN decrease when $A$ increases and $T$ decreases. To control the abnormal firing in $\beta$-band, the parameters of magnetic stimulation could be chosen from the dark blue areas of Figure 10(b–d).

4.2.2.   Control area of STN under various magnetic stimulation

To obtain the control strategy, the synchronous coefficient $r$, dominant frequency, and $\beta$-band power of STN in Figure 10 are considered together. Under each set of parameters, we consider that the controllability of STN is equal to 1 once dominant frequency of STN is within 13–30 Hz, the $r_{STN}$ is greater than 0.5, and $\beta$-band power is greater than 16. Accordingly, the controllability is 0 if the dominant frequency is less than 13 Hz or greater than 30Hz, $r_{STN}$ is smaller than 0.5, and $\beta$-band power less than 16, corresponding to successful control of STN. In Figure 11, black represents a value of 1 to denote the failed control of magnetic stimulation, whereas white represents a value of 0 when the activity of STN is well controlled.

The region of parameters in Figure 11 can be divided into two parts: controllable region I and uncontrollable region II. Region I in the bottom right corner with color blue is defined as the totally controllable area, whereas the rest is uncontrollable region II. In region I, the abnormal activity of STN is well controlled by magnetic stimulation, that is, a totally controlled region. In region II, magnetic stimulation is ineffective in most areas, and few parameter points can be controlled, that is, uncontrolled region. The controllable points in region II are various in 10 numerical simulations. It shows that a lower $T$ and higher $A$ can better control the abnormal activity of STN. The boundary curve between region I and region II (lower white curve in Figure 11) has the same trends in 10 simulations. The polynomial function

is chosen to fit the lower boundary curve in Figure 11, which has a R-square close to 1. The parameters of the fitted curve are listed in Table 4. In the fitted boundary curve over 10 simulations, one $\sigma$ ranges for $p_1$, $p_2$, $p_3$ and $p_4$ are $(-0.5597, -0.3689)$, $(3.893, 4.604)$, $(-4.864, -1.881)$, and $(9.57, 11, 24)$, respectively. The fitted parameters from No.1 to No.10 fall within one $\sigma$ range of No.All (white curve of Figure 12) except for the $p3$ of No.3. It indicates that the boundary curve between two region is stable from a statistical perspective. In summary, regions I and II are different in various simulations. The parameters in region II cannot be totally controlled mainly come from dominant frequency. To well control the activity of STN related to PD, magnetic stimulation parameters in region I should be used.

5.   Discussion

In order to study the firing activities in CBGT and understand the therapeutic effect of magnetic stimulation in PD, a CBGT network including magnetic induction effect was established by using the Izhikevich neural model as the basic node in this work. The computational model can be used to explore the effects of magnetic stimulation on the firing patterns of biophysical neural activity and control strategy. The role of magnetic stimulation is described by a magnetic induction model, which includes current induced inside and outside the model. The main results of this work are as follows.

Firstly, according to the properties in experiments, we reproduce the basic state in the CBGT neural network. The CBGT neural network in physiological and pathological conditions is proposed and analyzed by time series, spatiotemporal patterns, firing rates, power spectrum, and synchronous coefficient for each nucleus. Secondly, to explore the effect of the hyperdirect pathway, we further study the synchronization, oscillation dominant frequency and $\beta$-band power of the nucleus in various connection strength between PY and STN. The results show that a stronger synaptic connection from PY to STN in the hyperdirect pathway contributes to the excessive synchronization and $\beta$-band firing of STN. Then, to investigate the effect of magnetic stimulation, we apply an external magnetic induction with a periodic exponentially damped sine function on PY. The results show that magnetic induction can suppress $\beta$-band oscillation of STN and PY. And it also reduces synchronous coefficient $r$ of STN, GPi and PY. That is to say, external magnetic induction mainly affects the activity of the cortex directly and further influences the firing properties of STN in the BG through the hyperdirect pathway. To determine the control strategy, we analyze the statistical results under different parameters of magnetic induction including amplitude and period. The statistical synchronous coefficient $r$, dominant frequency, $\beta$-band power and max power of STN show that the increasing amplitude $A$ and decreasing period $T$ could suppress the abnormal firing activity of STN. Considering the $\beta$-band oscillation and synchronization with controllability, the region of magnetic stimulation of parameters can be divided into controllable (region I) and uncontrollable(region II). The fitted boundary curve between regions is given.

The intensity region of magnetic stimulation used in Section 4.2 is feasible in the computational model. However, in clinics, the security of the magnetic parameter should be verified. Previous work has considered various forms of external magnetic stimulation in isolated neurons and found that it can induce multiple patterns and bifurcations ^[16]. We also use the magnetic stimulation wave-shape with the sine function, and the statistical results are similar.

Our preliminary study makes effort to explore the effect of magnetic stimulation under electromagnetic radiation with a neural dynamical network model. However, due to the simplify of our CBGT network model, it fails to control abnormal oscillations with the waveforms of repetitive TMS and theta-burst stimulation due to the limitation of the computational model. This indicates that besides the hyperdirect pathway between cortex and STN, other therapeutic mechanisms of magnetic stimulation in PD should be further take into account. For example, TMS is expected to have an inhibitory effect on PY neurons mediated by activation of interneurons. Its treatment process accross longer time scale may relate to the endogenous dopamine in the striatum induced by magnetic stimulation ^[15], which could contribute to the improvement of synaptic connection and abnormal oscillations in PD. Consequently, how magnetic stimulation influences the ion current of neurons in basal ganglia directly and modulates dopamine level, and how the increase of dopamine contributes to synaptic plasticity quantitatively are valued and remained to be studied for a more specific model.

In this work, we take magnetic induction into account in a CBGT network model to study the therapeutic mechanism of magnetic stimulation in PD by considering internal and external magnetic effects. This work expands the application area of the magnetic induction model and helps to consider firing activities under magnetic effects and the controllability of magnetic stimulation in the CBGT network. In addition, our work provides insights into the effects of magnetic stimulation on PD.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (Nos. 11972292, 12172291, 12072265, and 11872304), 111 Project (No. BP0719007). The authors wish to acknowledge Ms. Zhuan Shen, Ms. Yuzhi Zhao and Mr. Haiwei Hu for their help in interpreting the significance of the results of this study.

Conflict of interest

The authors declare there is no conflicts of interest.