Control analysis of optogenetics and deep brain stimulation targeting basal ganglia for Parkinson's disease

1.   Introduction

Parkinson's disease (PD) is a common neurodegenerative disease that affects many people, especially the elderly. Its clinical manifestations include motor symptoms and non-motor symptoms, such as tremor, abnormal gait, stiffness and so on ^[1]. Rich evidence indicates that the occurrence of PD is closely related to the neurophysiological behaviors of the extrapyramidal system cortex-basal ganglia-thalamus circuit (CBGT) network. Typically, in pathological conditions, neurons in the basal ganglia (BG) networks exhibit abnormal firing patterns and rates, high synchronization between neurons and significant $\beta$ (13-30 $Hz$) band oscillations ^[2]. And the function of thalamus will also be abnormal, which is often manifested as the greatly reduced relay capacity. Notably, these abnormal behaviors are not carried out alone. For example, an increase in abnormal cluster discharge patterns and synchrony often coincide with the abnormal firing rates, just as the intermittent synchronization and periodic cluster firing behaviors can be seen in the recorded GPi neurons of MPTP-treated monkeys ^[2].

The clinical treatments for PD mainly focus on drug therapy, nerve nucleus destruction surgery, DBS, optogenetics and magnetic stimulation ^[3,4]. Among them, DBS, which allows to adjust electrode contact points, pulse width, frequency, and stimulation intensity according to the patients' conditions, is considered to be the most effective treatments for nerve disease. Some researchers explored the effects of the parameters of DBS on PD, in terms of thalamic relay fidelity and the suppression of pathological oscillation respectively ^[5,6,7]. To avoid the specific physical damage caused by continuous stimulation of a single nucleus and reduce energy consumption, there are some researches exploring the effects of combined DBS with the phase deviation on PD about desynchronization and thalamic relay fidelity ^[8]. However, DBS can't precisely target specific neurons because the neuron cell bodies around the electrodes and the axons passing through the electrodes are both activated. Optogenetics, as an advanced neuromodulation method, can control neuronal activity in a cell-specific manner with a high degree of time accuracy, which has significant advantages over traditional technologies such as electrical stimulation and pharmacological methods. Yu et al. performed high frequency optogenetic stimulation in STN neurons and proved it was effective in treating the symptoms of parkinsonism in the 6-hydroxydopamine (6-OHDA) lesion rat ^[9]. After isolating two groups of neurons in GPe before optogenetic intervention, Mastro et al. found that excitatory light stimulation of GPe restore the motor ability of mice in Parkinson's state ^[10]. And there are studies have explored the effects of optogenetic stimulation on the pathological parkinsonian rhythmic neural activity ^[11,12]. Although many attempts to alleviate and control PD have been done, there is currently no complete cure. It is still an urgent problem to explore efficient neural regulation strategies.

In addition, on the one hand, previous researches focusing on DBS only considered the influence of one or two parameters on PD. On the other hand, the current studies on optogenetic stimulation are mainly concentrated in clinical trials, but insufficiently researched on computational models. In order to further explore the regulatory mechanism of optogenetics and DBS on PD, a systematic model study based on basal ganglia-thalamic network model (BGTC) is carried out in this paper ^[6]. And according to experimental research, we select three types of neurons of basal ganglia nuclei (i.e., GPe, STN and GPi) as stimulation targets, so as to study comprehensively the combined effects of three typical parameters(i.e., intensity, frequency and stimulus duration) of DBS and optogenetic stimulation on the thalamic relay ability. Our results further enrich the existing conclusions from another perspective.

2.   Models

2.1.   Description of the BGTC network

The model used in our paper is the modified basal ganglia-thalamic network model (BGTC) developed by So et al. ^[6], which is comprised of thalamic (TC), STN, GPe and GPi and each nucleus contains 100 neurons. The network is connected in a sparse but structured manner. Each GPe projects inhibitory effects on two STN, two GPi, and two GPe neurons. Each STN simultaneously excites two GPe and two GPi neurons. These projections are chosen in accordance with the topographical organization within the basal ganglia. Each TC receives only a single inhibitory projection from one adjacent GPi neuron, and also receives excitatory input from the sensorimotor cortex (SMc), which is modeled as a series of pulses. The connection diagram is shown in Figure 1.

Then, the membrane potential equation of TC, STN and GP are described as below:

Here, $C_{m} = 1 pF/ \mu m^{2}$. $I_{L}$, $I_{N a}$ and $I_{K}$ are leak, sodium, and potassium ion currents, respectively. The current $I_{T}$ is a low-threshold T-type calcium current. $I_{Ca}$ is a high-threshold calcium current. $I_{ahp}$ is a $Ca^{2+}$-activated, voltage-independent afterhyperpolarization $K^{+}$ current. $I_{S M c}$ represents sensorimotor input to the TC neurons and is modeled as a periodic step function of the form:

Here, $H$ is the Heaviside step function, such that $H(x) = 0$ if $x\leq0$ and $H(x) = 1$ if $x > 0$. Note that $\rho_{\mathrm{SMc}}$ is the period of $I_{\mathrm{SMc}}$, $i_{\mathrm{SMc}}$ is the amplititude, and $\delta_{\mathrm{SMc}}$ is the duration of positive input. In addition, such as the form of $I_{G P i \rightarrow T C}$ denotes synaptic current. $I_{app}$ is a constant bias current. For more details, please refer to the references ^[5,6]. And in this paper, $I_{c}$ represents the current applied to STN, GPi or GPe neurons, which can be the DBS current ($I_{DBS}$) or the current generated by light stimulation ($I_{ChR2}$).

Based on the above model, we can simulate the healthy and PD state by changing the value of $I_{app}$. The time series and corresponding power spectral densities are shown in the Figure 2. In the healthy state, there is no obvious beta oscillation (13–30 $Hz$) in the BG nuclei, and TC can transmit SMc signals input faithfully. Conversely, in the PD state, GPe and GPi are bursting oscillating obviously, and the BG nucleus exhibit beta oscillation behavior. Meanwhile, TC neurons could not accurately respond to SMc signals. In this process, we consider three error types: TC neurons fail to respond SMc signals, TC neurons overrespond to SMc signals, and TC neurons spontaneously discharge in the absence of SMc signal input ^[5,6]. So the error index (EI) is proposed to measure the fidelity of TC relay to SMc signals, and EI is calculated as the total number of errors divided by the total number of input stimuli ^[5,6].

2.2.   Description of the four-state model

ChR2 is a sodium ion channel, which can depolarize the cell membrane by blue light stimulation to make cell excitement. And many experiments have reported that the response of neuron expressing ChR2 has three main characteristics: an initial peak with rapid decay, a steady-state plateau, and a decay back to baseline after irradiation ^[13]. In addition, there are currently two main computational models: the three-state model ^[14] and the four-state model ^[15]. The three-state model is not accurate enough for quantitative calculations, moreover, it is too simple to explain the biexponential decay of the light-off current. The four-state model has improved this. Hence, here we use four-state model to simulate ChR2 current dynamics ^[15].

The four-state model simulates the photocurrent kinetics for the light activated ChR2 channel through two sets of intratransitional states: more dark adapted ($C1 \leftrightarrow O1$) and light adapted ($C2 \leftrightarrow O2$). Without light, the ChR2 molecules are assumed to be in the closed state $C1$. However, in the light-adapted state, ChR2 molecules are distributed in all four states. The longer the light duration, the more molecules in the $O2$ state. After the light signal is terminated, the state $C2$ slowly transitions back to the state $C1$, and has a time scale of the recovery period. The percentage of ChR2 molecules in each of the four states of $C1$, $C2$, $O1$, and $O2$ are given by $c1$, $c2$, $o1$ and $o2$, respectively. So the four-state model of ChR2 photocurrent dynamics can be represented by the following differential equations:

where $c1+ c2+ o1+ o2 = 1$. $P_{1}$ and $P_{2}$ are the maximum excitation rates of the first and second closed states, $G_{d1}$ and $G_{d2}$ are the rate constants for the $O1 \rightarrow C1$ and the $O2 \rightarrow C2$ transitions respectively. $e_{12}$ and $e_{21}$ are the rate constants for $O1 \rightarrow O2$ and $O2 \rightarrow O1$ transitions. $G_{r}$ is the recovery rate of the first closed state after the light pulse is turned off. We note that, while the photon absorption and the isomerization of the retinal compound in the rhodopsin is almost instantaneous, the conformational change leading to the open state configuration is a slower process. The function $s$ aims to capture the temporal dynamics of this conformational change in the protein, which is provided in an explicit form ^[6]. In Eq (2.8), $S_{0}(\theta) = 0.5(1+\tanh (120(\theta-0.1)))$ is a sigmoidal function and $\theta(t) = \sum_{i} \Theta\left(t-t_{i_{\mathrm{on}}}\right) \Theta\left(t-t_{i_{\rm {off }}}\right)$ describes the optostimulation protocol. $\Theta(x)$ is the Heaviside function, while $t_{i_{on}}$ and $t_{i_{off}}$ are the onset and offset times of the $i_{th}$ optostimulation pulse, respectively. $\tau_{ChR2}$ is a constant representing the activation time of the ChR2 ion channel with typical values on the order of a few milliseconds. And the ChR2 photocurrent is given by:

Where $\gamma = \frac{g_{2}}{g_{1}}$. $V$ represents the real-time voltage applied to a specific neuron by photocurrent.

All simulations are written and performed under MATLAB environment (The Mathworks, Natick, MA, USA). Ordinary differential equations in our system are integrated using the standard forward Euler method with a fixed temporal resolution of 0.01 ms. Unless otherwise specified, values of model parameters in our simulations follow the baseline values listed in references ^[5,6,15].

3.   Results

DBS and optogenetics, as alternative therapeutic methods for PD, have been proved by some physiological experiments and clinical cases that they can alleviate the disease. STN and GPi are usually selected as DBS targets, which are able to reduce tremor, bradykinesia and stiffness, and then improve quality of PD patients' life ^[16]. Although the surgical experience using GPe as a target is not very rich so far, in view of GPe is a key contributor to motor suppressing pathways in the basal ganglia ^[10], here we apply DBS and ChR2 stimulation on STN, GPe and GPi respectively to explore the inner mechanism of PD.

3.1.   Effects of optogenetic excitatory stimulation targeting basal ganglia

Based on the BGTC model, here we first mainly consider the optogenetic excitatory stimulation of GPe. Intensity $\left(E \right)$, frequency $\left(f \right)$ and pulse width $\left(ws \right)$ of light stimulation are common parameters that can be varied. Ref ^[17] reported that damage to superficial cortical tissue are observed as a result of sustained light stimulation $E > 100 mW/mm^{2}$. So we limit the range of $E$ from $10 mW/mm^{2}$ to $100 mW/mm^{2}$ at $10 mW/mm^{2}$, and take into account the influence of $f$ and $ws$ on EI. As shown in Figure 3, we plot the two-dimensional plane of $\left(f, ws \right)$ with different values of $E$. From Figure 3(a) to Figure 3(c), the values of $E$ are chosen as $30 mW/mm^{2}$, $50 mW/mm^{2}$, and $70 mW/mm^{2}$, respectively. The blue region represents the smaller value of $EI$. It can be seen that EI will be affected only when $f$ and $ws$ reach certain values, such as $E = 50 mW/mm^{2}$ as shown in Figure 3(b). Horizontally, if the pulse width is too small, light stimulation hardly make effects for any frequency, standing still in a higer value of $EI$. But when the width is higher than $5 ms$, the TC neurons will be more sensitive to the frequency of light, emerging a much better improvement of TC relay performance with some different frequencies. From the vertical direction, if light frequency is lower than $30 Hz$, varing the width of each pulse does not yield to a good situation. In contrast, as the gradually increase of $f$ and $ws$, the relay capacity of the TC has been improved to some extent. Moreover, the areas of blue corresponding to the better TC relay performance in Figure 3(a-c) have a significant increasement with $E$ raises. Particularly, when $E = 70mW/mm^{2}$, larger $f$ and $ws$ could make the TC performs perfectly well $\left(EI\approx0 \right)$. Our numerical simulation results confirm the effectiveness of excitatory optogenetic stimulation of GPe, which is consistent with the experimental results ^[10].

Secondly, we apply ChR2 on STN, which aims to study the role of STN in mediating the PD symptom-relieving. Figure 3(d–f) displays the distribution of EI within the region $\left(f, ws \right)$ as $E$ is changed in a similar way, and the specific analysis will not be described here. But it's clear that the areas of dark red corresponding to the worse TC relay performance in Figure 3(d–f) have a slightly decrease as $E$ increases. However, it still cannot reach to a healthy level. This is consistent with the experimental result of Gradinaru et al. that direct optogenetic activation of STN neurons fails to relieve PD behavioral symptoms ^[18].

Finally, we discuss the impact of ChR2 stimulation of GPi on EI from the perspective of the model. As shown in Figure 3(g–i), it can be seen that within the range of our parameters, the value of EI is a little lower compared to the pathological state, which may means GPi-ChR2 is effective in improving TC relay capacity and alleviating the symptoms of PD. Based on the above analysis, for each light intensity, we always could find a set of appropriate parameter values in the $\left(f, ws \right)$ plane to obtain the minimum EI, as shown in the Figure 4(a). In summary, the effect of GPe-ChR2 is better than the others, which suggests that GPe can be regarded as a potential stimulus target in subsequent studies to explore the internal mechanism of PD.

Previous studies have shown that in some cases, the effects of the low-frequency with high-pulse-width stimulation mode is different from the high-frequency with low-pulse-width stimulation mode. Here, we fix the pulse duty cycle of light stimulation at $40\%$ to discuss which form of light stimulation is more effective in reducing $EI$. Figure 4(b–d) shows the relationship between $EI$ and $E$ under different pulse width stimulation for GPe, GPi and STN. Specifically, in the case of GPe - ChR2 as illustated in Figure 4(b), the value of $EI$ decreases overall with the increase of $E$ when the stimulus frequency is high and pulse width is low (i.e., $f$ = $100 Hz$, $ws$ = $4 ms$ or $f$ = $80 Hz$, $ws$ = $5 ms$). At the same time, under the situation of $f$ = 50 or $40 Hz$, the effect is relatively good. However, when the stimulus frequency is low enough, such as $f$ = $10Hz$, $ws$ = $40ms$ or $f$ = $20Hz$, $ws$ = $20ms$, the value of $EI$ is large, which means a less effective regulation. The situations of GPi-ChR2 and STN-ChR2 are shown in Figure 4(c, d), we will not repeat the detailed analysis to avoid duplication. One point should be particularly put forward is that, high-frequency with low-pulse-width stimulation mode is much better to control the value of $EI$ in our simulation.

Moreover, for GPe-ChR2, two sets of parameters are selected from the plane region of Figure 3(a-c), corresponding to the controllable and uncontrolled situations of PD symptoms respectively. At the same time, the time series, oscillating grating diagrams and power spectral densities of neurons are presented in Figures 5 and 6. It is important to note that our light stimulation time is from $200ms$ to $1800ms$. Under controllable conditions, ChR2 activates GPe, resulting in high frequency discharge of GPe, which reduces the firing rates of STN and GPi. Meanwhile, TC can transmit sensorimotor cortex (SMc) signal input faithfully. Furthermore, the synchronization between 100 neurons is destroyed, and there is no obvious beta oscillation in the BG nuclei. Conversely, under uncontrolled conditions, GPe and GPi cells are bursting oscillating, while TC neurons could not accurately respond to SMc signals. In addition, the 100 neurons still fire synchronously and the beta oscillations in the BG nuclei are obvious.

3.2.   Effects of single target DBS targeting basal ganglia

As well, we discuss the impact of stimulation frequency ($f$) and pulse width ($\delta$) on EI when the amplitude ($iD$) of DBS starts from $10 \mu A$ to $100 \mu A$ with steps of $10 \mu A$. In this section, we also select three representative amplitudes (i.e., $40 \mu A$, $70 \mu A$ and $100 \mu A$), and then plot two-dimensional planar graphs of $\left(f, \delta\right)$ for each fixed amplitude. The relations between the values of $EI$ and GPe-DBS parameters as seen in Figure 7(a–c). Significantly, the dark blue area represents $EI = 0$. With the increasement of $iD$, the area of the dark blue region gradually increases too, indicating that DBS has a better effect on the recovery of thalamic relay capacity. To be specific, for example, when $iD = 70 \mu A$, either $f$ or $\delta$ is too small, the value of $EI$ can not be reduced to a lower level. Only when $f > 20 Hz$ and $\delta > 0.2 ms$, GPe-DBS can greatly improve the relay performance of TC. Similar results also appear when DBS acts on GPi alone, as shown in the Figure 7(d–f), no further analyses will be given here. For STN-DBS, a larger $iD$ is required to produce a larger range of blue areas in the $\left(f, \delta\right)$ plane, which means a better control effect. However, as can be seen from the Figure 7(g–i), when $iD$ is small enough, such as $iD = 40 \mu A$, STN-DBS can also reduce appropriately the value of $EI$ as $f$ and $\delta$ increase to a certain range. Moreover, from the perspective of reducing the value of $EI$, we can draw a conclusion that STN-DBS largely depends on the stimulation amplitude, frequency and pulse width.

In addition, based on the above discussion, a new indicator is used to further explore the control effect of DBS on three different stimulus targets. As shown in Figure 7, $f \in \left[10,130\right], \delta \in \left[0, 1\right]$, we divide its parameter range into 13 and 11 units on average, so a total of 143 points in a two-dimensional plane are considered. Due to the error of numerical calculation, we believe that when $EI \leq 0.05$, corresponding to a fact that the relay capacity of the TC neurons is restored over $90 \%$, and the control effect of DBS is better. Therefore, we calculate the ratio of the number of points with $EI \leq 0.05$ to the total number of points. Obviously, a larger ratio reflects a better DBS control effect. As shown in Figure 8, with the increase of DBS stimulus amplitude, the ratio gradually increases, and GPe-DBS has the best effect. Similarly, the time series, oscillating grating diagrams and power spectral densities of neurons under two conditions, controllable and uncontrollable, are shown in Figure 9 and Figure 10. We believe that GPe is expected to become a clinically alternative and potential stimulation target in the future.

3.3.   Effects of combined DBS targeting basal ganglia

Considering the specific physical damage, side effects, and energy consumption caused by persistent stimulation of a single neuron causes, in this section, we discuss the influence of dual targets combined DBS on the thalamic reliability. We combine the stimulation targets (i.e., $STN$, $GPe$ and $GPi$) in pairs to obtain three different combined DBS strategies: STN $+$ GPe, STN $+$GPi and GPe $+$ GPi. In addition, the two different DBS protocols are adopted, as shown in Figure 11(a, e). The blue and red rectangles just represent the duration of the DBS pulse received by two different types of neurons. Specifically, the first stimulation mode is to periodically stimulate two nuclei in turn (see Figure 11(a)). The second stimulation mode is displayed in Figure 11(e), in which a complete $1:1$ stimulation is followed by a one period interval.

In order to explore the regulation mechanism of different stimulus strategies, we fix $iD = 100 \mu A$ and use two-dimensional parameter space $\left(f, \delta\right)$ to demonstrate the combined effect on EI. As shown in Figure 11(b–d), in the first stimulus mode, the blue area in the parameter plane is the largest when combined stimulation is GPe $+$ GPi, and $EI \approx 0$ as $f > 50 Hz, \delta > 0.1 ms$, indicating that the effect of stimulus parameters on improving the fidelity of TC relay tends to be stable. Combined stimulation of STN and GPe produce similar results, but the effect is relatively worse. However, when STN and GPi are stimulated at the same time, the area of the blue region is decreased compared with the first two conditions. Furthermore, when $f < 80 Hz$, the value of EI is generally large no matter how to change the value of $\delta$ in the selected range of parameters. To our pleasure, this result is basically consistent with the conclusion of ^[16]. In the second stimulus mode, the simulation results of the three combined stimulation are similar to the above. On the whole, the effect of the second stimulus mode is slightly worse than the first one, which is reasonable because the second stimulus strategy has a periodic intermittent period.

Finally, the indicator proposed in Sec. 3.2 is also used here to compare the stimulation effects of these two stimulus modes. As shown in Figure 12, the combined stimulus effect of GPe and GPi is significantly better than that of STN and GPe no matter which stimulation mode is applied, while the combined stimulus effect of STN and GPi is the worst. Moreover, it is obvious that the first stimulus strategy has a better performance.

4.   Discussion

Based on experimental research, we select three types of neurons of basal ganglia nuclei (i.e., GPe, STN and GPi) as stimulation targets. And we mainly focus on the thalamic relay ability to define the control effect, which is closely related to PD's motor symptoms ^[5]. Results show that the excitatory optogenetic stimulation of these three targets can improve the thalamic relay ability to a certain extent. Among them, the control effect of excitatory light stimulation targeting GPe is the best, which far exceeds that of GPi and STN. Most of these results are consistent with existing clinical trials and model conclusions ^[9,10,18]. For example, Mastro et al. showed that excitatory light stimulation of GPe could restore movement in dopamine-depleted mice ^[10]. Yu et al. found that frequency-specific optogenetic deep brain stimulation of STN improves PD's motor behaviors ^[9]. It is worth noting that although experiments have shown that optogenetic inactivation of GPi improves forelimb akinesia ^[19], the effect of excitatory light stimulation targeting GPi is still unclear. It is hoped that our modeling results can provide some feasible guidance.

GPi and STN are the most common targets for DBS in PD. In this paper, we also demonstrate that high-frequency DBS of GPi and STN can significantly reduce the error rate of thalamic relay and thus alleviate symptoms. To our surprise, DBS targeting GPe can also achieve perfect control effects, which is even better than targeting GPi. This finding supports current modeling and physiological findings that DBS targeting GPe is much better than GPi at eliminating Parkinson's motor symptoms ^[20], although surgical experience with GPe is still limited. To avoid the specific physical damage caused by continuous stimulation of a single nucleus and reduce energy consumption, we also explore the combined DBS of two nuclei with two stimulation protocols. Three different stimulation target combination are STN and GPe, STN and GPi, as well as GPi and GPe. Results show that the combined DBS of GPi and GPe is superior to the other two combined stimulations, which is consistent with the conclusion obtained by Yu et al. ^[8].

Furthermore, it should be pointed out the evaluation of the effect of PD treatment needs to further combine the pathological $\beta$ oscillation behavior and the synchronization of the network. And we only attempts to use the most frequently used activated channel protein ChR2 to excite neurons. More excitatory or even inhibitory channel proteins are worthy of further exploration for the regulatory effect of PD, such as ChETA, NpHR and so on ^[15].

5.   Conclusions

Neural stimulation has been the research focus in the field of neuroscience over the past century. In this paper, we systematically study the therapeutic effect of DBS and photogenetic stimulation. Simulation results demonstrate that the excitatory optogenetic stimulation and high-frequency DBS of these three targets can decrease the thalamic errors in relaying cortical information to a certain extent. Among them, the control effect of both stimulation targeting GPe is the best. Moreover, it is shown that the combined DBS of GPi and GPe is superior to the other two combined stimulations. Thus, Optostimulation with its optimal parameters can be considered as an alternative stimulation method to improve the Parkinsonian state. In addtion, our work shows that GPe can be a potential target for clinical stimulation.

Although these numerical results remains to be clinically studied, they also provide a theoretical basis for brain stimulation methods to further explore the pathogenesis and assist clinical treatment of PD and other neuropsychiatric disorders of motor and cognitive abnormalities.

Acknowledgments

This research was supported by the Natural Science Foundation of China (Nos. 11872304, 12072265), the Double First-class Construction Special Fund of Northwestern Polytechnical University (No. 0650021GH0201171) and the National Training Program of Innovation and Entrepreneurship for Undergraduates (No. S202010699537).

Conflict of interest

The authors declare there is no conflicts of interest.