Multiple-site deep brain stimulation with delayed rectangular waveforms for Parkinson's disease

1.   Introduction

Parkinson's disease (PD) affects tens of millions of people worldwide, and associated socioeconomic burden of the condition is set to increase as the elderly population grows. PD is a highly debilitating and prevalent neurodegenerative disorder characterized by both motor and nonmotor symptoms [1], with the former mainly including muscle rigidity, 4-7 Hz rest tremor and akinesia [2]. The core, but not exclusive pathology is the degeneration of the dopaminergic neurons in the substantia pars compacta (SNc) of the midbrain that projects to the striatum [3].

The symptoms of PD ultimately result from the death of dopaminergic neurons in a subcortical brain structure called the basal ganglia (BG) [4]. The BG is a collection of nuclei tucked deep in the brain, which is composed of three nuclei of neurons: the internal and external segments of the globus pallidus (GPi & GPe), and the subthalamic nucleus (STN). They have been implicated in a wide variety of tasks - action gating, selection, or automation [5], and in the motor aspects of neurological disorders such as PD [6].

Human patients with advanced PD are often treated by deep brain stimulation (DBS), which can alleviate the disease's motor symptoms [7]. In the 1960s, it was discovered that electrical stimulation to specific areas of the brain can alleviate symptoms in neurological disorders [8]. Today, the treatment technique has grown and DBS has been used in more than 150000 patients worldwide [9]. DBS is an approved treatment for a variety of conditions such as PD, dystonia and obsessive-compulsive disorder [10].

The pulsatile stimulation pattern of high-frequency DBS used recently consists of a first, anodal, activating pulse part and a second, cathodal, charge-balanced pulse part [11]. The second part immediately follows the first part but has a longer duration and smaller amplitude, as shown in Fig. 1. Although DBS is a remarkable therapy for PD today, there is also evidence that chronic high-frequency stimulation may cause long-term tissue damage due to the risky combination of the frequency and amplitude [12].

The traditional DBS waveform is a square wave and always consumes a large amount of energy to operate, which increases the risks of device replacement surgeries. Many investigations indicated that well-designed irregular waveform may have clinical advantages compared with the standard DBS. Hofmann et al. found that the introduction of a gap with a specific and optimized duration in a biphasic pulse could increase the efficiency and therefore reduce the energy consumption of the stimulation [13]. Moreover, the conventional DBS usually performed on a fixed target area of the BG. Guo and Rubin proposed that multi-site stimulation, applied to the STN, regularize GPi outputs and significantly diminish TC relay error [14]. Recent experimental and clinical studies have put forward a closed-loop DBS. Peter A. Tass group showed that permanent pulsatile high-frequency stimulation subjected to an amplitude modulation by linear or nonlinear delayed feedback methods can effectively and robustly desynchronize the STN-GPe network of model neurons and suggest this approach for desynchronizing closed-loop DBS [15]. Furthermore, they find that adaptive pulsatile linear delayed feedback stimulation (apLDF) can be efficient in suppressing abnormal synchronization and also causes a stronger reduction of the LFP beta burst length [16].

Inspired by the multi-site stimulation on the STN proposed by Guo and the work done by Peter A. Tass group, multiple-site stimulation with the pulse-delay-pulse type DBS was discussed in this paper. We introduced a gap after the activating pulse part in the pulsatile stimulation pattern. And according to the different combinations of the stimulus targets, we proposed three different options: STN-GPe, GPe-GPi, and STN-GPi. And it was concluded that DBS at STN-GPe resulted in better performance at lower amplitudes and frequencies. Here, we adopted multiple-site stimulation strategy to get away from the hazard of the high-frequency stimulation. By using distinct stimulation strategies, it helps us to understand how PDP-type DBS can provide therapeutic effects for PD based on a computational model of basal ganglia-thalamus (BG-TH) network.

The content of this paper is organized as follows. The BG-TH network and the DBS waveform are introduced in Section 2. In order to evaluate the stimulation effects, three quantitative measures are proposed in this Section. Analysis of the effects of the PDP-type DBS with multiple-site stimulation are presented in Section 3. Finally, a brief conclusion is drawn in Section 4.

2.   Methods

2.1.   The model of basal ganglia-thalamus (BG-TH) network

The model used in this paper was established by So et al. [17], which is composed of four nuclei of neurons: GPi, GPe, STN and TH as shown in Fig. 2. Each nucleus consists of 10 neurons. Each STN neuron can produce excitatory effects on two GPe neurons and two GPi neurons simultaneously. Each GPe neuron exerts inhibitory effects on two STN neurons, two GPe neurons, and two GPi neurons. Each GPi neuron makes an inhibitory projection to a TH neuron. These projections were chosen in accordance with the topographical organization and the convergence/divergence of synaptic connectivity within the BG system [18].

The dynamics of each cell in each nucleus is represented by a single-compartment conductance-based model. And their membrane potentials were similar to those used by Rubin and Terman [19,20]. Membrane potentials of the individual STN, GPe, GPi, and TH cells were governed by a set of differential equations:

where $C_{m}$ represents the membrane capacitance. Additionally, $I_{L}$, $I_{Na}$, $I_{K}$, $I_{T}$, $I_{Ca}$, $I_{a h p}$ are the leak current, sodium current, potassium current, low-threshold T-type calcium current, high-threshold calcium current, and after hyperpolarization potassium current.

Synaptic current from the presynaptic neuron $\alpha$ to the postsynaptic neuron $\beta$ is modeled by:

where $g_{\alpha \rightarrow \beta}$ is the maximal synaptic conductance between neuron $\alpha$ and neuron $\beta$, and $E_{\alpha \rightarrow \beta}$ is the reversal potential of the synaptic current. The summation is taken over all the presynaptic neurons projected to a given postsynaptic neuron.

For STN to GPi, STN to GPe and GPi to TH, the synaptic variables are governed by a second order alpha synapse to acquire an acute sense of synaptic effects as follows:

Here, $u(t) = 1$ if the presynaptic cell crosses a threshold of $-10$ mV at time $t$, indicating the presence of an action potential in the presynaptic cell. Otherwise, $u(t) = 0$.

For GPe to STN, GPe to GPi, GPe to GPe, the synaptic variables are reserved as the following first order differential equation:

where $H_{\infty}$ is a smooth approximation of the Heaviside step function.

Besides, TH neurons receive excitatory input from the sensory motor cortex (SMC), which is modeled by a series of single-phase current pulses with an amplitude of 3.5 $\mu\mathrm{A} / \mathrm{cm}^{2}$ and a duration of 5 $\mathrm{ms}$. The instantaneous frequencies of incoming pulses were drawn from a gamma distribution with an average rate of 14 Hz and a coefficient of variation of 0.2. The variance is included to simulate the non-regular nature of incoming signals from the cortex [17]. For the parameter values and more details involved in the model, see So et al. [17].

2.2.   DBS waveform

A biphasic, charge-balanced pulse prevents net charge injection into the stimulated tissue. Charge balancing is therefore a necessary requirement to avoid tissue damage [21]. The ratio of pulse part amplitudes also defines the ratio of pulse part widths, which thus guarantees charge balancing [13]. We set the ratio throughout this work is 15:1 to ensure that there is sufficient time to depolarize the membrane potential to maximize DBS efficiency. The amplitudes of the cathodic and anodic phase are 300 $\mu\mathrm{A}$ and -20 $\mu\mathrm{A}$, and the pulse widths are 0.2 $\mathrm{ms}$ and 3.0 $\mathrm{ms}$. The pulse width parameters of the pulse used in this study are within the parameter range used in clinical studies for DBS [22]. In the previous studies, for a gap length of 1.8 $\mathrm{ms}$ an improvement of more than 90% of the maximum possible improvement was obtained [13]. Hence, in this paper, we chose the gap length in our DBS signal as 1.8 $\mathrm{ms}$. With this configuration, we have the ability to restore the thalamic fidelity while having a fixed optimal gap length. The DBS waveform denoted as PDP used in the BG-TH network model is shown in Fig. 3.

2.3.   Evaluation measures of the stimulation

1) Error index (EI)

In healthy state, each input from SMC results in a single action potential in each TH neuron. The performance of the BG network was quantified by measuring how accurately the TH neurons relayed inputs from the SMC as introduced by Terman [19]. The error index (EI) is defined as the total number of errors divided by the total number of excitatory inputs ($I_{SMC}$). Three types of errors in thalamic transmission were considered: misses, bursts and spurious events. A miss occurred when a neuron failed to spike, a burst occurred when a neuron spiked more than once within 25 $\mathrm{ms}$ of a stimulating input, and a spurious event occurred if the thalamic cell fired in the absence of a stimulating input [17]. The network achieves optimal performance when each input pulse from the SMC results in an action potential in each TH neuron, that is EI equals to zero.

2) Power Spectra

Patients with PD exhibit prominent feature of low frequency oscillation in the beta band. And oscillatory power in the beta band correlates with akinesia /bradykinesia [23]. Consequently, by computing the power spectrum, we can detect the pathological oscillations of PD in BG network. Spectral analyses were performed using the "mtspectrumpt" function of the Chronux neural signal analysis package (chronux.org) (sliding 1 $\mathrm{s}$ window, 0.1 $\mathrm{s}$ step size and $[3, 5]$ tapers (3 is the time-bandwidth product and 5 is the number of tapers)) [24]. Oscillatory power in the GPi was calculated by integrating the spectral power of GPi spike times in the $7-35$ Hz frequency band [25].

3) Firing Rates

In the healthy state, the STN, GPe and GPi neurons exhibit regular firing at near constant frequencies. Changes in firing patterns observed in these nucleus upon switching from healthy to Parkinsonian state. That is, the STN neurons exhibit more irregularity with varied interspike interval, and GPe and GPi neurons fired more frequently in burst. So, under the Parkinsonian state, several prominent changes will be reflected in the abnormal firing rates. PD patients would have higher firing rates of their STN neurons [26]. The average firing rate is defined by dividing the number of spikes of neurons by the simulation time. In what follows, we calculate the average firing rate for STN, GPe and GPi by 1000 $\mathrm{ms}$ simulation time.

2.4.   Simulation details

The state of the BG-TH network can be achieved by changing the value of the external current applied to the STN, GPe and GPi. The switch from healthy condition to Parkinsonian condition in the model can be achieved by decreasing the constant bias currents $I_{a p p_{-} i}(i \in\mathrm{S T N, G P e, G P i})$, whose values for two different states are presented in Table 1. And to simulate the role of direct pathway of the basal ganglia, the value of $I_{S t r \rightarrow G P i}$ is set to $-2$ $\mu\mathrm{A} / \mathrm{cm}^{2}$.

Simulations were implemented in Matlab 2019a. The numerical method used was the forward Euler method with a time step of 0.01 $\mathrm{ms}$. Unless otherwise indicated, the results and the plots with error bars are obtained from the average results of 50 independent simulation runs due to the randomness of the initial conditions.

3.   Results

With the multiple-site stimulation, we mainly study the effects of the stimulation frequency of the PDP-type DBS for the purpose of attaining better therapeutic effects at a lower frequency. Through the measures mentioned above, we can find the optimal frequency, which can be used in the future clinic.

3.1.   The effectiveness of the PDP-type DBS

First of all, EI is used as the first indicator to measure the effect of the PDP-type DBS treatment and explore the influence of stimulation frequency. The stimulation frequency of DBS varies from 10 to 200 Hz.

1) STN-GPe DBS

Figure 4 exhibits the effect of stimulation frequency on EI of STN-GPe DBS. It can be seen that when the frequency range is 10 Hz-160 Hz, the value of EI shows a decreasing trend with the increase of the stimulation frequency. PDP-type DBS with frequency being equal to and greater than 30 Hz is sufficient to restore the BG network to the healthy condition. Stimulation on STN and GPe at the same time led to an improvement in thalamic fidelity with lower stimulation frequency.

Actually, STN-GPe DBS increases the tonic excitation to GPi from STN and the tonic inhibition to GPi from GPe simultaneously. The overall effect of inhibition to GPi overcomes the excitatory effect by STN at even low frequencies. Specifically, stimulation of STN and GPe resulted in high frequency regular spiking of neurons in both GPe and GPi due to the excitatory connection from STN to GP. The GPi has inhibitory efferent to the TH cells, and regularization of the output from the GPi resulted in the TH receiving regular tonic inhibition, which led to higher thalamic fidelity.

Therefore, we can choose the lower optimal frequency range of the PDP-type DBS to be 20 Hz-100 Hz.

2) STN-GPi DBS

Figure 5 presents the influence of the stimulation frequency on EI of STN-GPi DBS. It is evident that the trend of EI is more complicated than that in STN-GPe DBS. When the frequency range is 10 Hz-100 Hz, the value of EI first decreases gradually and then increases a little, which indicates that there is an optimal frequency that minimizes the value of EI within this range. The same result occurred when the frequency changes from 100 Hz-200 Hz. When the stimulation frequency is greater than 170 Hz, the value of EI will increase dramatically.

It may suggest that, from the viewpoint of the stimulation effectiveness, adding PDP-type DBS to STN and GPi to modulate the Parkinsonian state is not an applicable technique, which depends greatly on the stimulation parameters. We need to acknowledge the validity of this hypothesis later by other measures. The temporary optimal frequency of the PDP-type DBS can be chosen as 60 Hz-90 Hz.

3) GPi-GPe DBS

Figure 6 shows the influence of GPi-GPe DBS. When the frequency range is 20 Hz-160 Hz, the value of EI is approximately equal to 0. In addition, PDP-type DBS with frequency being greater than 160 Hz makes the BG network in a higher frequency state.

Direct stimulation on GPi caused high frequency regular firing of the GPi neurons, and therefore a regular level of inhibition to the TH would result in higher fidelity. On the other hand, introducing PDP-type DBS to GPe led to inhibition and silencing of neurons in both the STN and GPi, thereby in turn led to a loss of inhibition to the TH neurons, which also resulted in improved fidelity.

This suggests that the low frequency PDP-type DBS on GPi-GPe site may be chosen as an alternative method to modulate the Parkinsonian state. In this situation, we choose the optimal frequency range of the DBS as 20 Hz-100 Hz.

3.2.   The effect of the stimulation frequency on $\beta$ oscillation of GPi

Although EI has been regarded as a significant index to measure the therapeutic effect of DBS in many studies, it is clear that we cannot only rely on the results of screening based on EI to justify that the DBS at these frequencies can restore the entire network to the healthy state. The $\beta$ power of GPi nuclei was used as the second indicator to measure the effect of the DBS. The range of the stimulation frequency was selected according to the results we obtained in the above section.

1) STN-GPe DBS

Figure 7 displays the influence of the stimulation frequency on $\beta$ power of GPi on STN-GPe stimulation. When the frequency is 20 Hz and 30 Hz, the $\beta$ band oscillation in GPi is significantly enhanced, even though the value of EI remains at a lower level. The reason is that $\beta$ oscillation is defined as in the range of $17-35$ Hz, the same frequency added to the network will cause resonance, which makes the power surge.

The reciprocally connected STN-GPe network may be a possible source of pathological low frequency oscillatory activity in PD. STN-GPe stimulation suppressed pathological low-frequency oscillations by exciting some GPi neurons through the STN-GPi pathway and inhibiting other GPi neurons through the STN-GPe-GPi pathway. Excited GPi neurons showed a decrease in pathological burst activity and exhibited a more regularized firing, while inhibited GPi neurons did not transmit the pathological activity to the TH neurons. Therefore, pathological low-frequency oscillatory activity cannot reach the output of the BG network.

In addition, the value of EI remains at a relatively low level when the frequency range is 40 Hz$-100$ Hz. In this paper, we prefer to use lower frequency DBS to achieve therapeutic effect comparable to standard high frequency DBS. It is concluded that the optimum frequency range of the PDP-type DBS on STN-GPe site can be chosen as 40 Hz$-100$ Hz considering EI and $\beta$ power together.

2) STN-GPi DBS

Figure 8 manifests the influence of the stimulation frequency on $\beta$ power of GPi with the PDP-type DBS working on STN-GPi site. It can be found that low frequency also produces greater oscillation power of GPi. In this situation, DBS eliminates $\beta$ oscillation of GPi at a very high level of frequency.

STN and GPi were stimulated concurrently makes GPi neurons extremely active through the STN-GPi pathway. And there is no doubt that in this case, the $\beta$ oscillation in GPi nuclei will be enhanced.

As mentioned in Section 3.1, STN-GPi DBS is not an applicable technique. In view of EI and $\beta$ oscillation outcomes, we no longer take STN-GPi DBS into consideration.

3) GPi-GPe DBS

Figure 9 reveals the effect of the stimulation frequency on $\beta$ power of GPi with the PDP-type DBS on GPi-GPe site. When the frequency is 20 Hz and 30 Hz, the resonance phenomenon mentioned above also occurs in GPi nuclei. Moreover, we can see intuitively that when the frequency is equal to or greater than 40 Hz, the $\beta$ oscillation in GPi always maintains at a low level.

The mechanism of GPi-GPe DBS is similar to STN-GPe DBS, through inhibitory efferent from GPe to GPi, the firing of GPi neurons gradually recover to the regular state. Thus, we can adopt the optimal frequency range of the PDP-type DBS as 40 Hz-100 Hz.

3.3.   Firing rates and patterns

Under diverse circumstances, the firing rates and patterns will be discrepant. Firing rates was used as the third index to measure the effect of the PDP-type DBS. Furthermore, spike rasters for three nuclei in the BG network are also displayed, so as to gain more insights into the mechanism underlying the DBS. In what follows, only the raster diagrams with 40 Hz stimulation and 100 Hz stimulation are shown.

1) STN-GPe DBS

See Fig. 10, as the stimulation frequency increases, the firing rates of STN and GPe neurons also increase, which implies that PDP-type DBS serves as a pacemaker to directly drive each STN and GPe neuron to fire in accordance with the ratio 1:1. While GPi neurons always maintains at a normal level, which is consistent with the phenomenon observed in the experiments.

As illustrated in Fig. 11, we find that while the stimulation frequency change between 40 Hz and 100 Hz, STN neurons will be driven to fire. As the frequency increases, the synchronous firing of neurons in the two nuclei, GPi and GPe, will be enhanced, which in turn promotes the synchronization of the whole network. 40 Hz PDP-type DBS destroys the phasic bursting synchronization shown in the Parkinsonian state, followed by the restoration of the TH relay reliability. In contrast, the neurons stimulated by traditional DBS or 100 Hz PDP-type DBS exhibit a high-frequency firing state.

2) GPi-GPe DBS

From Fig. 12, it is concluded that after simultaneous stimulation on GPi and GPe, GPe neurons will also be driven to fire, and the firing rate will continue to increase. Besides, there is no substantial increase in the firing rate of GPi neurons. However, the firing rate of STN decreases gradually and keeps below the normal state.

Furthermore, Fig. 13 shows that under the influence of the PDP-type DBS on GPi-GPe site, the number of spikes of STN neurons is significantly lower than the healthy level, and the firing activities among most neurons are extremely synchronized. We conclude that it is because that the stimulation on GPe excessively inhibits STN, which makes STN in a state of inactivation. Consequently, when the frequency is increased to 100 Hz, the number of spikes of STN neurons become less, and the abnormal firing behaviour of STN appears. In this case, GPi and GPe exhibit a high-frequency firing state again.

4.   Conclusion

In this study, a new DBS strategy which is called pulse-delay-pulse (PDP) type DBS on multiple-site, is designed to adopt on a computational model of the BG-TH network. On one side, through the application of DBS with lower frequency on two different nuclei, that is, multiple-site stimulation, the therapeutic effect can be compared with that under single-site stimulation, eliminating the risk of using traditional high frequency DBS and reducing side effects. On the other side, as has been studied by many researchers, interfering the DBS waveform with a delay between the cathodic and anodic parts also increases the threshold for the activation of neurons and will make the network less influenced by the anodic phase. In order to obtain the optimal stimulation frequency used in this type of DBS, error index, $\beta$ power of GPi and firing rate of the BG network were introduced. The numerical results showed that the PDP-type DBS on STN-GPe site resulted in better performance at lower frequencies. It can be concluded that the optimal interval of the stimulation frequency is $40-50$ Hz. Stimulation on GPi-GPe causes the spike number of STN to decline and deviate from the healthy status, and it is not a good choice. Furthermore, stimulation on STN-GPi site can only produce fairly good therapeutic effects at higher frequencies (greater than 120 Hz).

In the practical application, PDP-type DBS discussed in this paper is easier to implement in hardware, and it is verified through numerical analysis that multiple-site stimulation can also attain promising recovery outcome. Furthermore, closed-loop control of PDP-type DBS at STN-GPe is also an issue worthy of further consideration. It will be developed to adjust the stimulation frequency with respect to different PD stages. Moreover, we will discuss more waveforms of DBS on multiple-site strategy and test the robustness of the TC relay responses in our future work.