Analysis of neural clusters due to deep brain stimulation pulses

Abstract

Deep brain stimulation (DBS) is an established method for treating pathological conditions such as Parkinson’s disease, dystonia, Tourette syndrome, and essential tremor. While the precise mechanisms which underly the effectiveness of DBS are not fully understood, several theoretical studies of populations of neural oscillators stimulated by periodic pulses have suggested that this may be related to clustering, in which subpopulations of the neurons are synchronized, but the subpopulations are desynchronized with respect to each other. The details of the clustering behavior depend on the frequency and amplitude of the stimulation in a complicated way. In the present study, we investigate how the number of clusters and their stability properties, bifurcations, and basins of attraction can be understood in terms of one-dimensional maps defined on the circle. Moreover, we generalize this analysis to stimuli that consist of pulses with alternating properties, which provide additional degrees of freedom in the design of DBS stimuli. Our results illustrate how the complicated properties of clustering behavior for periodically forced neural oscillator populations can be understood in terms of a much simpler dynamical system.

Appendices

Appendix A: neuron models

In this appendix, we give details of the neural models used in this paper, specifically the Hodgkin–Huxley model considered in the main text, and the thalamic neuron model considered in Appendix B.

Hodgkin–Huxley neuron model

The full Hodgkin–Huxley model is given by:

$$\begin{aligned} \dot{V}= & {} (I_b -\bar{g}_{Na}h(V-V_{Na})m^3-\bar{g}_K(V-V_K)n^4\\&-\bar{g}_L(V-V_L))/c + u(t) \; , \\ \dot{m}= & {} a_m(V)(1-m)-b_m(V)m \; , \\ \dot{h}= & {} a_h(V)(1-h)-b_h(V)h \; , \\ \dot{n}= & {} a_n(V)(1-n)-b_n(V)n \; , \end{aligned}$$

where

$$\begin{aligned} a_m(V)= & {} 0.1(V+40)/(1-\exp (-(V+40)/10)) \; , \\ b_m(V)= & {} 4\exp (-(V+65)/18) \; , \\ a_h(V)= & {} 0.07\exp (-(V+65)/20) \; , \\ b_h(V)= & {} 1/(1+\exp (-(V+35)/10)) \; , \\ a_n(V)= & {} 0.01(V+55)/(1-\exp (-(V+55)/10)) \; , \\ b_n(V)= & {} 0.125\exp (-(V+65)/80) \; , \end{aligned}$$

The parameters for this model are

$$\begin{aligned} V_{Na}= & {} 50 \; mV \;,\; V_K=-77 \; mV \;,\; V_L=-54.4 \; mV \;, \\ \bar{g}_{Na}= & {} 120 \; mS/cm^2 \; , \; \bar{g}_K=36 \; mS/cm^2 \;,\\ \bar{g}_L= & {} 0.3 \; mS/cm^2 \;,\; I_b=10 \; \mu A/cm^2 \;,\\ c= & {} 1 \; \mu F/cm^2 . \end{aligned}$$

Thalamic neuron model

The full thalamic neuron model is given by:

$$\begin{aligned} \dot{V}= & {} \frac{-I_L-I_{Na}-I_K-I_T+I_b}{C_m}+u(t),\\ \dot{h}= & {} \frac{h_{\infty }-h}{\tau _h},\\ \dot{r}= & {} \frac{r_{\infty }-r}{\tau _r}, \end{aligned}$$

where

$$\begin{aligned} h_\infty= & {} 1/(1+\exp ((V+41)/4)),\\ r_\infty= & {} 1/(1+\exp ((V+84)/4)),\\ \alpha _h= & {} 0.128\exp (-(V+46)/18),\\ \beta _h= & {} 4/(1+\exp (-(V+23)/5)),\\ \tau _h= & {} 1/(\alpha _h+\beta _h),\\ \tau _r= & {} (28+\exp (-(V+25)/10.5)),\\ m_\infty= & {} 1/(1+\exp (-(V+37)/7)),\\ p_\infty= & {} 1/(1+\exp (-(V+60)/6.2)),\\ I_L= & {} g_L(V-e_L),\\ I_{Na}= & {} g_{Na}({m_\infty }^3)h(V-e_{Na}),\\ I_K= & {} g_K((0.75(1-h))^4)(V-e_K),\\ I_T= & {} g_T(p_\infty ^2)r(V-e_T). \end{aligned}$$

The parameters for this model are

$$\begin{aligned} C_m= & {} 1 \; \mu F/cm^2 \;,\; g_L = 0.05 \; mS/cm^2 \;,\; e_L = -70 \; mV \;,\\ g_{Na}= & {} 3 \; mS/cm^2 \;,\; e_{Na} = 50 \; mV, g_K = 5 \; mS/cm^2 \;,\\ e_K= & {} -90 \; mV \;,\; g_T = 5 \; mS/cm^2 \;,\; e_T = 0 \; mV \;,\\ I_b= & {} 5 \; \mu A/cm^2. \end{aligned}$$

Appendix B: results for thalamic neurons

In this appendix, we show simulations and analysis for an (approximately) Type I neuron model, the thalamic neurons from (Rubin and Terman 2004). The full equations are given in Appendix A; for our simulations, we use the corresponding phase model. For reference, for these parameters the thalamic neurons have $\omega = 0.748$ rad/s.

We consider populations of thalamic neurons with the same stimuli (7) with $u_{max}$ corresponding to a current density of 20 $\mu A/cm^2$, $p = 0.5$ ms, and $\lambda = 3$. We simulated 500 thalamic neurons with initial phases evenly spaced between 0 and $2 \pi$, corresponding to a uniform distribution. The stimulation frequency was varied from 70 Hz to 300 Hz in increments of 5 Hz. Figure 20 shows the final phases after 40 periods of stimulation, after transients have decayed. Figure 21 shows the time series of the phases of a population of such neurons for selected frequencies. Here, we again see clustering for some frequencies (such as 250 Hz, where $r_2$ is large and $r_1$ and $r_3$ are small, indicating a 2-cluster solution), and non-clustering behavior for other frequencies (such as 200 Hz, where the Lyapunov exponent $\varLambda$ is positive, corresponding to chaotic dynamics).

Fig. 20

a The final phases $\theta$ of thalamic neurons drawn from an initial uniform distribution as a function of stimulation frequency, after 40 periods of stimulation. Colors correspond to the neurons’ initial phases. b Order parameters $r_1$ (black), $r_2$ (blue), and $r_3$ (red) for the final state as a function of frequency. For the initial uniform distribution, $r_1=r_2=r_3=0$. c Lyapunov exponent $\varLambda$ as a function of stimulation frequency

Fig. 21

Time series showing the phases of thalamic neurons drawn from an initial uniform distribution for frequencies a 200 Hz, and b 250 Hz. For (a), clusters do not form; for (b), there are two clusters after transients decay away

The same analysis techniques can also be used to understand the dynamics of thalamic neurons subjected to periodic pulses. Figure 22a shows the response function $f(\theta )$ for thalamic neurons with the stimulus given by (7) with $u_{max}$ corresponding to a current density of $20 \mu A/cm^2$, $p = 0.5$ ms, and $\lambda = 3$; Fig. 22b shows that there is a stable 2-cluster state for a stimulation frequency of 250 Hz, as expected from Fig. 20.

Fig. 22

a Response function $f(\theta )$ which characterizes the phase response of thalamic neurons to a pulse with $u_{max}$ corresponding to a current of $20 \mu A/cm^2$, $p = 0.5$ ms, and $\lambda =3$. b Map $g^{(2)}(\theta )$ for the thalamic neuron with stimulation frequency 250 Hz, showing two stable fixed points which correspond to a 2-cluster state

Appendix C: shift properties of the maps

Proposition 1

(Shift properties of $g^{(n)}$): Iterates of the map

$$\begin{aligned} g_\tau (\theta ) = \theta + \omega \tau + f(\theta + \omega \tau ) \end{aligned}$$

(28)

satisfy the property

$$\begin{aligned} g_\sigma ^{(n)}(\theta ) = g_\tau ^{(n)}(\theta + \omega (\sigma - \tau )) + \mathcal{O}(\sigma - \tau ). \end{aligned}$$

(29)

Proof

We will prove this by induction. First, (29) holds for $n=1$ from (19) in the main text (in fact, in this case the $\mathcal{O}(\sigma - \tau )$ correction term vanishes). Next, let us assume that (29) holds for n; we will show that this implies that it also holds for $n+1$. For reference,

$$\begin{aligned} g_\tau ^{(n+1)}(\theta )= & {} g_\tau (g_\tau ^{(n)}(\theta ))\nonumber \\= & {} g_\tau ^{(n)}(\theta ) + \omega \tau + f(g_\tau ^{(n)}(\theta ) + \omega \tau ). \end{aligned}$$

(30)

Now,

$$\begin{aligned} g_\sigma ^{(n+1)}(\theta )= & {} g_\sigma (g_\sigma ^{(n)}(\theta )) \\= & {} g_\sigma ^{(n)}(\theta ) + \omega \sigma + f(g_\sigma ^{(n)}(\theta ) + \omega \sigma ) \\= & {} g_\tau ^{(n)}(\theta + \omega (\sigma - \tau )) + \mathcal{O}(\sigma - \tau ) + \omega \sigma \\&+ f(g_\tau ^{(n)} (\theta + \omega (\sigma - \tau )) + \mathcal{O}(\sigma - \tau ) + \omega \sigma ), \end{aligned}$$

where we have used (29). We now use $\omega \sigma = \omega \tau + \omega (\sigma - \tau ) = \omega \tau + \mathcal{O}(\sigma - \tau )$ to give

$$\begin{aligned} g_\sigma ^{(n+1)}(\theta )= & {} g_\tau ^{(n)}(\theta + \omega (\sigma - \tau )) + \omega \tau + \mathcal{O}(\sigma - \tau ) \\&+ f(g_\tau ^{(n)} (\theta + \omega (\sigma - \tau )) + \omega \tau + \mathcal{O}(\sigma - \tau )). \end{aligned}$$

Next, we Taylor expand the last term about $g_\tau ^{(n)} (\theta + \omega (\sigma - \tau )) + \omega \tau$, treating $(\sigma - \tau )$ as small:

$$\begin{aligned} f(g_\tau ^{(n)}(\theta+ & {} \omega (\sigma - \tau )) + \omega \tau + \mathcal{O}(\sigma - \tau )) \\= & {} f(g_\tau ^{(n)}(\theta + \omega (\sigma - \tau )) + \omega \tau ) + \mathcal{O}(\sigma - \tau ) \\ \end{aligned}$$

Thus,

$$\begin{aligned} g_\sigma ^{(n+1)}(\theta )= & {} g_\tau ^{(n)}(\theta + \omega (\sigma - \tau )) + \omega \tau \\&+ f(g_\tau ^{(n)} (\theta + \omega (\sigma - \tau )) + \omega \tau ) + \mathcal{O}(\sigma - \tau ) \\= & {} g_\tau ^{(n+1)} (\theta + \omega (\sigma - \tau )) + \mathcal{O}(\sigma - \tau ), \end{aligned}$$

where the last equality follows from (30). Therefore, (29) holds for all $n \ge 1$, with no $\mathcal{O}(\sigma - \tau )$ term necessary for $n=1$ from (19).

Proposition 2

(Shift properties of $G^{(n)}$): Iterates of the map

$$\begin{aligned} G_{\tau , \tau _2}(\theta )= & {} \theta + \omega \tau + f_2 (\theta + \omega \tau _2) \nonumber \\&+ f(\theta + \omega \tau + f_2(\theta + \omega \tau _2)) \end{aligned}$$

(31)

satisfy the property

$$\begin{aligned} G_{\sigma , \sigma _2}^{(n)}= & {} G_{\tau , \tau _2}^{(n)} (\theta + \omega (\sigma - \tau ) + \omega (\sigma _2 - \tau _2)) \nonumber \\&+ \mathcal{O}(\sigma - \tau ) + \mathcal{O}(\sigma _2 - \tau _2). \end{aligned}$$

(32)

Proof

We will prove this using induction. Let us first show that (32) holds for $n=1$. By definition,

$$\begin{aligned} G_{\sigma , \sigma _2}(\theta )= & {} \theta + \omega \sigma + f_2 (\theta + \omega \sigma _2) \nonumber \\&+ f(\theta + \omega \sigma + f_2(\theta + \omega \sigma _2)). \end{aligned}$$

Letting

$$\begin{aligned} \theta = \hat{\theta } + \omega (\tau - \sigma ) + \omega (\tau _2 - \sigma _2), \end{aligned}$$

and simplifying, we obtain

$$\begin{aligned} G_{\sigma , \sigma _2} (\theta )= & {} \hat{\theta } + \omega \tau - \omega (\sigma _2 - \tau _2) \\&+ f_2(\hat{\theta } + \omega \tau _2 - \omega (\sigma - \tau )) \\&+ f(\hat{\theta } + \omega \tau - \omega (\sigma _2 - \tau _2) \\&\;\;\;\;\; + f_2(\hat{\theta } + \omega \tau _2 - \omega (\sigma - \tau ))) \\= & {} \hat{\theta } + \omega \tau + f_2(\hat{\theta } + \omega \tau _2) \\&+ f(\hat{\theta } + \omega \tau + f_2 (\hat{\theta } + \omega \tau _2)) \\&+ \mathcal{O}(\sigma - \tau ) + \mathcal{O}(\sigma _2 - \tau _2) \\= & {} G_{\tau , \tau _2}(\hat{\theta }) + \mathcal{O}(\sigma - \tau ) + \mathcal{O}(\sigma _2 - \tau _2) \\= & {} G_{\tau , \tau _2} (\theta + \omega (\sigma - \tau ) + \omega (\sigma _2 - \tau _2)) \\&+ \mathcal{O}(\sigma - \tau ) + \mathcal{O}(\sigma _2 - \tau _2). \end{aligned}$$

Here, the second equality follows from Taylor expansion, treating $(\sigma - \tau )$ and $(\sigma _2 - \tau _2)$ as small. Thus, (32) holds for $n=1$.

Now, suppose (32) holds for n; we will show this also implies that it holds for $n+1$. For reference,

$$\begin{aligned} G_{\tau , \tau _2}^{(n+1)}(\theta )= & {} G_{\tau , \tau _2} (G_{\tau , \tau _2}^{(n)}(\theta )) \nonumber \\= & {} G_{\tau , \tau _2}^{(n)}(\theta ) + \omega \tau + f_2(G_{\tau , \tau _2}^{(n)} (\theta ) +\omega \tau _2) \nonumber \\&+ f(G_{\tau , \tau _2}^{(n)}(\theta ) + \omega \tau + f_2 (G_{\tau , \tau _2}^{(n)}(\theta ) +\omega \tau _2)). \nonumber \\ \end{aligned}$$

(33)

Now,

$$\begin{aligned} G_{\sigma , \sigma _2}^{(n+1)}(\theta )= & {} G_{\sigma , \sigma _2} (G_{\sigma , \sigma _2}^{(n)}(\theta )) \\= & {} G_{\sigma , \sigma _2}^{(n)}(\theta ) + \omega \sigma + f_2(G_{\sigma , \sigma _2}^{(n)} (\theta ) +\omega \sigma _2) \\&+ f(G_{\sigma , \sigma _2}^{(n)}(\theta ) + \omega \sigma + f_2 (G_{\sigma , \sigma _2}^{(n)}(\theta ) +\omega \sigma _2)) \\= & {} G_{\tau , \tau _2}^{(n)}(\theta + \omega (\sigma - \tau ) +\omega (\sigma _2 - \tau _2)) \\&+ \mathcal{O}(\sigma - \tau ) + \mathcal{O}(\sigma _2 - \tau _2) + \omega \sigma \\&+ f_2(G_{\tau , \tau _2}^{(n)} (\theta + \omega (\sigma - \tau ) + \omega (\sigma _2 - \tau _2)) \\&\;\;\;\;\; + \mathcal{O}(\sigma - \tau ) + \mathcal{O}(\sigma _2 - \tau _2) + \omega \sigma _2) \\&+ f(G_{\tau , \tau _2}^{(n)} (\theta + \omega (\sigma - \tau ) + \omega (\sigma _2 - \tau _2)) \\&\;\;\;\;\;+ \mathcal{O}(\sigma - \tau ) + \mathcal{O}(\sigma _2 - \tau _2) + \omega \sigma \\&\;\;\; + f_2 (G_{\tau , \tau _2}^{(n)} (\theta + \omega (\sigma - \tau ) + \omega (\sigma _2 - \tau _2)) \\&\;\;\;\;\; + \mathcal{O}(\sigma - \tau ) + \mathcal{O}(\sigma _2 - \tau _2) + \omega \sigma _2)), \end{aligned}$$

where the last equality follows from (32). Letting

$$\begin{aligned} \omega \sigma = \omega \tau + \omega (\sigma - \tau ) = \omega \tau + \mathcal{O}(\sigma - \tau ) \end{aligned}$$

and

$$\begin{aligned} \omega \sigma _2= & {} \omega \tau _2 + \omega (\sigma _2 - \tau _2) = \omega \tau _2 + \mathcal{O}(\sigma _2 - \tau _2),\\ G_{\sigma , \sigma _2}^{(n+1)}(\theta )= & {} G_{\tau , \tau _2}^{(n)}(\theta + \omega (\sigma - \tau ) +\omega (\sigma _2 - \tau _2)) \\&+ \omega \tau + \mathcal{O}(\sigma - \tau ) + \mathcal{O}(\sigma _2 - \tau _2) \\&+ f_2(G_{\tau , \tau _2}^{(n)} (\theta + \omega (\sigma - \tau ) + \omega (\sigma _2 - \tau _2)) \\&\;\;\;\;\; + \omega \tau _2 + \mathcal{O}(\sigma - \tau ) + \mathcal{O}(\sigma _2 - \tau _2)) \\&+ f(G_{\tau , \tau _2}^{(n)} (\theta + \omega (\sigma - \tau ) + \omega (\sigma _2 - \tau _2)) \\&\;\;\;\;\; + \omega \tau + \mathcal{O}(\sigma - \tau ) + \mathcal{O}(\sigma _2 - \tau _2) \\&\;\;\; + f_2 (G_{\tau , \tau _2}^{(n)} (\theta + \omega (\sigma - \tau ) + \omega (\sigma _2 - \tau _2)) \\&\;\;\;\;\; + \omega \tau _2 + \mathcal{O}(\sigma - \tau ) + \mathcal{O}(\sigma _2 - \tau _2) )). \end{aligned}$$

Finally, treating $(\sigma - \tau )$ and $(\sigma _2 - \tau _2)$ as small and Taylor expanding f and $f_2$,

$$\begin{aligned} G_{\sigma , \sigma _2}^{(n+1)}(\theta )= & {} G_{\tau , \tau _2}^{(n)}(\theta + \omega (\sigma - \tau ) +\omega (\sigma _2 - \tau _2)) + \omega \tau \\&+ f_2(G_{\tau , \tau _2}^{(n)} (\theta + \omega (\sigma - \tau ) \nonumber \\&+ \omega (\sigma _2 - \tau _2)) + \omega \tau _2) \\&+ f(G_{\tau , \tau _2}^{(n)} (\theta + \omega (\sigma - \tau ) + \omega (\sigma _2 - \tau _2)) + \omega \tau \\&\;\;\; + f_2 (G_{\tau , \tau _2}^{(n)} (\theta + \omega (\sigma - \tau ) \nonumber \\&+ \omega (\sigma _2 - \tau _2))+ \omega \tau _2)) \\&+ \mathcal{O}(\sigma - \tau ) + \mathcal{O}(\sigma _2 - \tau _2) \\= & {} G_{\tau , \tau _2}^{(n+1)} (\theta + \omega (\sigma - \tau ) + \omega (\sigma _2 - \tau _2)) \\&+ \mathcal{O}(\sigma - \tau ) + \mathcal{O}(\sigma _2 - \tau _2), \end{aligned}$$

as desired, where the last equality follows from (33). Thus, (32) holds for all $n \ge 1$.

Proposition 3

(Shift properties of $h_1$ and $h_2$):

The maps

$$\begin{aligned} h_{1 \tau , \tau _2}(\theta )= & {} \theta + \omega (\tau - \tau _2) + f(\theta + \omega (\tau - \tau _2)) \end{aligned}$$

(34)

$$\begin{aligned} h_{2 \tau _2}(\theta )= & {} \theta + \omega \tau _2 + f_2 (\theta + \omega \tau _2) \end{aligned}$$

(35)

satisfy the properties

$$\begin{aligned} h_{1 \sigma , \sigma _2}(\theta )= & {} h_{1 \tau , \tau _2}(\theta + \omega (\sigma - \tau ) + \omega (\tau _2 - \sigma _2)), \end{aligned}$$

(36)

$$\begin{aligned} h_{2 \sigma _2}(\theta )= & {} h_{2 \tau _2} (\theta + \omega (\sigma _2 - \tau _2)). \end{aligned}$$

(37)

Proof

First, consider

$$\begin{aligned} h_{1 \sigma , \sigma _2}(\theta ) = \theta + \omega (\sigma - \sigma _2) + f(\theta + \omega (\sigma - \sigma _2)). \end{aligned}$$

Letting

$$\begin{aligned} \theta = \hat{\theta } + \omega (\tau - \sigma ) + \omega (\sigma _2 - \tau _2) \end{aligned}$$

and simplifying,

$$\begin{aligned} h_{1 \sigma , \sigma _2}(\theta )= & {} \hat{\theta } + \omega (\tau - \tau _2) + f(\hat{\theta } + \omega (\tau - \tau _2)) \nonumber \\= & {} h_{1 \tau , \tau _2}(\hat{\theta }) \\= & {} h_{1 \tau , \tau _2}(\theta + \omega (\sigma - \tau ) + \omega (\tau _2 - \sigma _2)). \end{aligned}$$

Now, consider

$$\begin{aligned} h_{2 \sigma _2}(\theta ) = \theta + \omega \sigma _2 + f_2(\theta + \omega \sigma _2). \end{aligned}$$

Letting

$$\begin{aligned} \theta = \hat{\theta } + \omega (\tau _2 - \sigma _2) \end{aligned}$$

and simplifying,

$$\begin{aligned} h_{2 \sigma _2}(\theta )= & {} \hat{\theta } + \omega \tau _2 + f(\hat{\theta } + \omega \tau _2) = h_{2 \tau _2}(\hat{\theta }) \\= & {} h_{2 \tau _2} (\theta + \omega (\sigma _2 - \tau _2)). \end{aligned}$$