Analysis of Neural Clusters due to Deep Brain Stimulation Pulses
AAnalysis of Neural Clusters due toDeep Brain Stimulation Pulses
Daniel Kuelbs · Jacob Dunefsky · Bharat Monga · Jeff MoehlisAbstract
Deep brain stimulation (DBS) is an estab-lished method for treating pathological conditions suchas Parkinson’s disease, dystonia, Tourette syndrome,and essential tremor. While the precise mechanisms whichunderly the effectiveness of DBS are not fully under-stood, theoretical studies of populations of neural os-cillators stimulated by periodic pulses suggest that thismay be related to clustering, in which subpopulationsof the neurons are synchronized, but the subpopula-tions are desynchronized with respect to each other.The details of the clustering behavior depend on thefrequency and amplitude of the stimulation in a com-plicated way. In the present study, we investigate howthe number of clusters, their stability properties, andtheir basins of attraction can be understood in terms ofone-dimensional maps defined on the circle. Moreover,we generalize this analysis to stimuli that consist ofpulses with alternating properties, which provide addi-tional degrees of freedom in the design of DBS stimuli.Our results illustrate how the complicated properties ofclustering behavior for periodically forced neural oscil-
Daniel KuelbsStanford University, Palo Alto, CA, 94305E-mail: [email protected], Co-first AuthorJacob DunefskyYale University, New Haven, CT, 06520E-mail: [email protected], Co-first AuthorBharat MongaDepartment of Mechanical Engineering, University of Cali-fornia, Santa Barbara, CA 93106E-mail: [email protected]ff MoehlisDepartment of Mechanical Engineering, Program in Dynami-cal Neuroscience, University of California, Santa Barbara, CA93106E-mail: [email protected], Corresponding Author lator populations can be understood in terms of a muchsimpler dynamical system.
Keywords
Neural oscillators, Clustering, Phasemodels, Deep Brain Stimulation
A primary motivation for this study is Parkinson’s dis-ease, which can cause an involuntary shaking that typ-ically affects the distal portion of the upper limbs, anddifficulty initiating motion. For patients with advancedParkinson’s disease who do not respond to drug ther-apy, electrical deep brain stimulation (DBS), an FDA-approved therapeutic procedure, may offer relief [3].Here, a neurosurgeon guides a small electrode into thesub-thalamic nucleus or globus pallidus interna (GPi);the electrode is connected to a pacemaker implanted inthe chest which sends periodic electrical pulses directlyinto the brain tissue. The efficacy of DBS for the treat-ment of Parkinson’s disease has been found to dependon the frequency of stimulation, with high-frequencystimulation (70 to 1000 Hz and beyond) being ther-apeutically effective [3,29,25]. The generally acceptedtherapeutic range is 130-180 Hz [36,14]. Experimentalevidence has suggested that motor symptoms of Parkin-son’s disease are associated with pathological synchro-nization of neurons in the basal ganglia, and that DBSdesynchronizes the neural activity [35,7,12,16,32]. DBShas also shown promising results in treating other neu-rological conditions, for which the stimulation electrodeis implanted in the GPi (for dystonia) or the thalamus(for Tourette syndrome and essential tremor) [31,2].While the precise mechanisms which underly the ef-fectiveness of DBS are not fully understood, theoreti-cal studies have shown that DBS-like stimulation con-sisting of a periodic pulses applied to neural oscillator a r X i v : . [ m a t h . D S ] J u l Daniel Kuelbs et al. populations can lead to chaotic desynchronization [37]or clustering behavior [41], in which subpopulationsof the neurons are synchronized, but the subpopula-tions are desynchronized with respect to each other.Clustering has also been found in theoretical studiesof coordinated reset, in which multiple electrodes de-liver inputs which are separated by a time delay [17,18,19,34]. These studies, along with clinical successeswith coordinated reset [1], point to clustering as anattractive objective for designing stimulation proper-ties; this has motivated the design of single control in-puts which promote clustering [20,21,38], in contrast tomethods which seek to fully desynchronize the neuralactivity [33,26,40,22]. Notably, clustering has at leasttwo important differences from chaotic desynchroniza-tion: clustered states often exist over a much larger pa-rameter range than chaotic desynchronization, a possi-ble explanation why effective DBS parameters are easierto find than chaotic desynchronization would suggest;and clustered states may induce plasticity changes moreeffectively than chaotic desynchronization, which mayexplain why benefits are more persistent for some kindsof stimulation mechanisms than others (cf. [1,21]). Inthis paper, we will focus on clustering which arises froma single stimulation electrode, unlike coordinated resetwhich uses multiple electrodes.Despite substantial data backing the general efficacyof DBS, it can have side effects including disorienta-tion, memory deficits, spatial delayed recall, responseinhibition, episodes of mania, hallucinations, or moodswings, as well as impairment of social functions such asthe ability to recognize the emotional tone of a face [8,6]. Our study develops tools which can help to iden-tify different stimuli that result in the same clusteringbehavior; our hope is that the identification of thesealternatives will allow neurologists to consider differ-ent stimuli in order to find those which are effectiveat treating neurological disorders while minimizing theseverity of side effects.In this paper, we investigate how the details of clus-tering due to periodic pulses of the type used in DBScan be understood in terms of one-dimensional mapsdefined on the circle. As a first step, Section 2 de-scribes phase reduction, a powerful classical techniquefor the analysis of oscillators in which a single variabledescribes the phase of the oscillation with respect tosome reference state. Section 3 shows results from sim-ulations of populations of neural ocsillators stimulatedby periodic pulses of the type used for DBS; this illus-trates the different types of clustering which can occur,and motivates the theoretical analysis. Section 4 de-rives and investigates the one-dimensional maps whichcan be used to understand the types of clusters which occur, their stability properties, and their basins of at-traction. Section 5 then demonstrates how this analysisin terms of maps can be generalized to consider stimulithat consist of pulses with alternating properties, whichprovide additional degrees of freedom for DBS stimulusdesign. Section 6 summarizes the results. The modelsfor the neurons considered in this paper are given inthe Appendix.
A common way to describe the dynamics of neurons isto use conductance-based models such as the Hodgkin-Huxley equations [13]. Such models are typically high-dimensional and contain a large number of parameters,which can make them unwieldy for simulations of largeneural populations. A powerful technique for the analy-sis of oscillatory neurons, whose dynamics are describedby a stable periodic orbit, is the rigorous reduction ofconductance-based models to phase models, with a sin-gle variable θ describing the phase of the oscillationwith respect to some reference state [43,15,23].Suppose that our conductance-based model is de-scribed by the n -dimensional dynamical system d x dt = F ( x ) , x ∈ R n ( n ≥ , (1)with a stable periodic orbit γ ( t ) with period T . For eachpoint x ∗ in the basin of attraction of γ ( t ) there existsa corresponding phase θ ( x ∗ ) such that [11,43]lim t →∞ (cid:12)(cid:12)(cid:12)(cid:12) x ( t ) − γ (cid:18) t + T π θ ( x ∗ ) (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) = 0 , (2)where, under the given vector field, x ( t ) is the trajec-tory of the initial point x ∗ . The asymptotic phase of x , θ ( x ), ranges in value from [0 , π ). In this paper, θ = 0will represent the phase at which the neuron fires anaction potential. Isochrons are level sets of θ ( x ), andwe define isochrons such that the phase of a trajectoryevolves linearly in time both on and off of the periodicorbit [42,43]. As a result, for the entire basin of attrac-tion of the periodic orbit, dθdt = 2 πT ≡ ω. (3)If we now consider the dynamical system d x dt = F ( x ) + U ( t ) , x ∈ R n , (4)where U ( t ) ∈ R n is an infinitesimal control input, phasereduction gives the one-dimensional system [15,5,23] dθdt = ω + U ( t ) T Z ( θ ) . (5) nalysis of Neural Clusters due to Deep Brain Stimulation Pulses 3 In this equation, Z ( θ ) is the gradient of θ evaluatedon the periodic orbit, and is known as the phase re-sponse curve (PRC) [43,10,27]; it represents the changein phase that the control input will cause when appliedat a given phase. In this paper, we consider electricalcurrent inputs which only act in the voltage directiondefined by the unit vector ˆ V , i.e., U ( t ) = u ( t ) ˆ V , withthe corresponding phase reduction dθ i dt = ω + Z ( θ i ) u ( t ) . (6)Here, θ i represents the phase of the i th neuron, ω is thenatural frequency of the neuron in radians per second, Z ( θ ) = ∂θ∂V is the component of the PRC in the voltagedirection, and u ( t ) is the input. For the populationsof neurons considered in this paper, we assume thatthe neurons are identical and they all receive the sameinput, and we will consider uncoupled neurons withoutnoise; these assumptions allow a more detailed analysisto be performed.In the next section, we show simulation results forpopulations of neurons described by such phase modelswith periodic pulses of the type used for DBS. In this section, we show simulation results for popula-tions of neurons stimulated by periodic pulses of thetype used for DBS; these results will inspire the analy-sis in Section 4. To illustrate a range of clustering be-haviors, we show simulations for prototypical systemswhich represent two common types of neurons [28]: as aType I neuron model, we consider the model for thala-mic neurons from [30], and as a Type II model we con-sider the Hodgkin-Huxley equations [13]. These modelsare not meant to correspond to the neurons directly rel-evant to Parkinson’s disease in human patients; rather,they are used to illustrate typical clustering behaviorsfor populations of neural oscillators under DBS-like stim-uli. The full equations and parameters for these modelsare given in the Appendix. For our simulations, we usethe corresponding phase models. For reference, for theseparameters the thalamic neurons have ω = 0 .
748 rad/s,and the Hodgkin-Huxley neurons have ω = 0 .
429 rad/s.The PRC functions for these neurons are shown in Fig-ure 1(a) and (b), respectively. Each PRC was calculatednumerically using XPP [9], and is approximated by aFourier series.The input u ( t ) that we consider, shown in Figure 2and inspired by DBS stimuli [24], is a periodic sequenceof identical charge-balanced pulses parameterized byamplitude u max , period τ (with corresponding frequency (a) (b) Z ( θ ) θ θZ ( θ ) Fig. 1
Panels (a) and (b) show the phase response curves Z ( θ ) of the thalamic (Type I) and Hodgkin-Huxley (Type II)neurons considered in this paper, respectively. t = 0 t = τ t = 2 τp u max u min tu ( t ) λp Fig. 2
Periodic sequence of identical pulses. /τ ), pulse width p , and multiplier λ (the ratio of timethat the pulse is negative to the time that the pulse ispositive). Mathematically, u ( t ) is given by: u ( t ) = u max mod( t, τ ) ≤ pu min ≡ − u max λ p < mod( t, τ ) ≤ ( λ + 1) p . (7)Unless otherwise stated, we will use u max correspond-ing to a current density of 20 µA/cm , p = 0 . λ = 3 in our simulations. We consider different frequen-cies of stimulation between 70-300 Hz, which includesthe typical therapeutic range of 130-180 Hz for DBStreatment of Parkinson’s disease.We simulated 500 Hodgkin-Huxley neurons with ini-tial phases evenly spaced between 0 and 2 π , correspond-ing to an initial uniform phase distribution. The stim-ulation frequency was varied from 70 Hz to 300 Hz inincrements of 5 Hz. Figure 3 shows the final phases af-ter 40 periods of stimulation, after transients have de-cayed away. The colors indicate the initial phases of theneurons. Not all colors are visible for most stimulationfrequencies because the final phases of entire subpop-ulations of neurons are nearly identical, and only onerepresentative initial phase can be seen. All of the neu-rons which have nearly the same final phase are part ofthe same cluster.Figure 4 shows the times series of the phases of apopulation of Hodgkin-Huxley neurons for selected fre-quencies, and helps us to interpret the results shown in Daniel Kuelbs et al. frequency (Hz) θ Fig. 3
The final phases θ of Hodgkin-Huxley neurons drawnfrom an initial uniform distribution as a function of stimula-tion frequency, after 40 periods of stimulation. Colors corre-spond to the neurons’ initial phases. Figure 3. For example, Figure 4(a) shows that for a 100Hz stimulus the neurons separate into three clusters, asis also the case for 250 Hz as shown in Figure 4(e). (No-tice that Figure 3 shows three possible final phases foreach of these frequencies, corresponding to these threeclusters.) Figure 4(b) shows that for a 150 Hz stimulusthey separate into two clusters. For a 180 Hz stimu-lus, there is no clustering; see Figure 4(c). By carefullylooking at Figure 4(d), one sees that for a 185 Hz stim-ulus there are five clusters, and from Figure 4(f) thatfor a 295 Hz stimulus there are seven clusters, as ex-pected from final states shown at these frequencies inFigure 3. Such clustering behavior and non-clustering(chaotic) behavior has been seen in other studies, suchas [37] and [41].Inspired by neural synchrony in Parkison’s patients,we also considered an initial partially synchronized neu-ral population, with phases distributed according to avon Mises distribution [4] centered at θ = 0: ρ ( θ ) = e κ cos θ πI ( κ ) , (8)where I ( κ ) is the modified Bessel function of order 0.This distribution is similar to a Gaussian distribution,but on a circle. We simulated 500 Hodgkin-Huxley neu-rons with initial phases distributed according to the vonMises distribution with κ = 50. As for Figure 3, thestimulation frequency was varied from 70 Hz to 300 Hzin increments of 5 Hz. Figure 5 shows the final phasesafter 40 periods of stimulation, after transients have de-cayed away. We see that the final phases of the neuronsfrom the initial von Mises distribution lie on a subsetof the final phases of the neurons from the intial uni-form distribution. For example, when the stimulationfrequency is 100 Hz, the neurons from the initial vonMises distribution are concentrated in two of the threeclusters which exist for the initial uniform distribution. (d)(a)(b)(c)(e)(f)
185 Hz: 5 clusters250 Hz: 3 clusters290 Hz: 7 clusters θθ tθθθθ
150 Hz: 2 clusters tttt
100 Hz: 3 clusters180 Hz: no clustering t Fig. 4
Time series showing the phases of Hodgkin-Huxleyneurons drawn from an initial uniform distribution for fre-quencies (a) 100 Hz, (b) 150 Hz, (c) 180 Hz, (d) 185 Hz, (e)250 Hz, and (f) 290 Hz. The titles of these panels indicate thenumber of clusters found after transients have decayed away.For (c), clusters do not form. For this and subsequent timeseries figures, t is measured in ms, and the colors indicate theinitial phases of the neurons, with colorbar as in Figure 3.nalysis of Neural Clusters due to Deep Brain Stimulation Pulses 5 frequency (Hz) θ Fig. 5
As a function of stimulation frequency, the final phases θ of Hodgkin-Huxley neurons drawn from an initial von Misesdistribution after 40 periods of stimulation are shown as black ∗ ’s, overlaid on the final phases of Hodgkin-Huxley neuronsdrawn from an initial uniform distribution (as was shown inFigure 3). We also designed an algorithm to detect the sizeof clusters in a population. The algorithm groups thephases of neurons in a population at each timestep intoclusters by sorting the phases in ascending order andchecking if the i th phase is within (cid:15) of the ( i +1)-th phasefor an appropriate small value of (cid:15) . If so, the size of thecurrent cluster is increased by one. If not, the algorithmcreates a new cluster. The process is repeated until allneurons have been grouped into clusters. Figure 6 showsthe number of neurons in the different clusters over arange of frequencies for the intial uniform distribution(for which three clusters are populated) and von Misesdistribution (for which only two clusters are populated).As we will see in Section 4, this figure can be explainedin terms of the basins of attraction of fixed points ofiterates of a one-dimensional map defined on the circle.The initial phase of a given neuron will determine whichcluster it ends up in.We also considered populations of thalamic (TypeI) neurons with the same stimuli (7) with u max corre-sponding to a current density of 20 µA/cm , p = 0 . λ = 3. We simulated 500 thalamic neurons withinitial phases evenly spaced between 0 and 2 π , corre-sponding to a uniform distribution. The stimulation fre-quency was varied from 70 Hz to 300 Hz in incrementsof 5 Hz. Figure 7 shows the final phases after 40 periodsof stimulation, after transients have decayed. Figure 8shows the time series of the phases of a population ofsuch neurons for selected frequencies. Here we again seeclustering for some frequencies (such as 250 Hz), andnon-clustering behavior for other frequencies (such as200 Hz).In the next section, we derive and investigate one-dimensional maps which can be used to understand thetypes of clusters which occur in these simulations, along
190 200 210 220 230 240 250
Frequency (Hz) C l u s t e r S i z e Cluster ICluster II (b)
200 220 240 260
Frequency (Hz) C l u s t e r S i z e Cluster ICluster IICluster III (a)
Fig. 6
The number of Hodgkin-Huxley neurons in differentclusters for a population size of 500, with initial (a) uniformand (b) von Mises distributions. frequency (Hz) θ Fig. 7
The final phases θ of thalamic neurons drawn froman initial uniform distribution as a function of stimulationfrequency, after 40 periods of stimulation. Colors correspondto the neurons’ initial phases. with their stability properties and their basins of attrac-tion. In this section, we show how the clustering behaviorfound in the simulations from Section 3 can be un-derstood in terms of appropriate compositions of one-dimensional maps on the circle.We consider a system of neural oscillators subjectedto a τ -periodic sequence of pulses as shown in Figure 2, Daniel Kuelbs et al. (b)(a)
250 Hz: 2 clusters200 Hz: no clustering θ tθ t
Fig. 8
Time series showing the phases of thalamic neuronsdrawn 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. and described by the dynamics [41]˙ θ i = ω + f ( θ i ) δ (mod( t, τ )) , i = 1 , · · · , N. (9)Here the response function f ( θ ) describes the changein phase due to a single pulse (including the positivecurrent for time p , and the negative current for time λp ). If the pulse was a delta function with unit area, f ( θ ) would be equal to the infinitesimal PRC Z ( θ ); formore general pulses, it can be calculated using a directmethod in which a pulse is applied at a known phase,and the change in phase is deduced from the change intiming of the next action potential [27]. We will thinkof the change in phase due to the pulse as occurringinstantaneously, even though the pulse will typicallyhave a finite duration; this will be a good approximationfor pulses of short duration. Figure 9 shows f ( θ ) for theHodgkin-Huxley neurons considered in this paper forpulses as shown in Figure 2 with u max correspondingto a current density of 20 µA/cm , p = 0 . λ = 3.To understand the clustering behavior, it will beuseful to consider the map which takes the phase ofa neuron to the phase exactly one forcing cycle later,cf. [41]. To find this map, suppose that we start with θ (0 + ) = 0, immediately after the start of a pulse, wherewe assume that we have already accounted for the ef-fect of the pulse according to the function f ( θ ). Thenext pulse comes at time τ . Up until time τ , the phaseevolves according to ˙ θ = ω ; therefore, θ ( τ − ) = θ + ωτ. (10) f () f ( θ ) θ Fig. 9
Response function f ( θ ) which characterizes the phaseresponse of Hodgkin-Huxley neurons to the stimulus, for u max corresponding to a current density of 20 µA/cm , p = 0 . λ = 3. Treating the change in phase due to the next pulse asoccurring instantaneously, we have θ ( τ + ) = θ + ωτ + f ( θ + ωτ ) . (11)The system then evolves for a time τ without stimulus,giving θ (2 τ − ) = θ + 2 ωτ + f ( θ + ωτ ); (12)the next pulse at time 2 τ gives θ (2 τ + ) = θ +2 ωτ + f ( θ + ωτ )+ f ( θ +2 ωτ + f ( θ + ωτ )) , (13)and so on. It is useful to let [41] g ( s ) = s + ωτ + f ( s + ωτ ) , (14)which gives θ ( nτ + ) = g ( n ) ( θ ) , (15)where g ( n ) denotes the composition of g with itself n times, and θ is the initial state of the neuron.We look for fixed points of g ( n ) , that is, solutionsto θ ∗ = g ( n ) ( θ ∗ ); for such solutions, the phase has thesame value after n pulses as where it started. We areparticularly interested in fixed points of g ( n ) which arenot fixed points of g ( m ) for any positive integer m sat-isfying m < n ; then there will be n fixed points of g ( n ) that correspond to points on a period- n orbit of g . If (cid:12)(cid:12)(cid:12)(cid:12) ddθ (cid:12)(cid:12)(cid:12)(cid:12) θ = θ ∗ ( g ( n ) ( θ )) (cid:12)(cid:12)(cid:12)(cid:12) < , (16)then the fixed point θ ∗ of g ( n ) is stable, as is the corre-sponding period- n orbit of g . Neurons which start withinitial phases within the basin of attraction of a givenfixed point of g ( n ) will asymptotically approach thatfixed point under iterations of g ( n ) . The n different fixedpoints will each have a basin of attraction, so a uniform nalysis of Neural Clusters due to Deep Brain Stimulation Pulses 7 intial distribution of neurons will form n clusters, onefor each of these fixed points of g ( n ) , cf. [41].We now illustrate how these maps can be used tounderstand the specific clustering behavior shown inSection 3. As a first example, suppose that a popula-tion of Hodgkin-Huxley neurons is stimulated with fre-quency 150 Hz, corresponding to τ = 6 .
67 ms. Figure 10shows g ( θ ) and g (2) ( θ ). Fixed points of these maps cor-respond to intersections with the diagonal. We see thatthere are two stable fixed points for g (2) ( θ ), at θ = 2 . θ = 5 .
86 (these fixed points are stable because theslope at the intersection is between − g (2) at θ = 1 . θ = 4 . g ( θ ), buta cobweb analysis verifies that there is a period-2 orbit θ = 2 . → . → . → · · · . These fixed points of g (2) correspond a stable 2-clusterstate for a population of oscillators, as shown in Fig-ure 4(b). We note that we can deduce the basin ofattraction for the different stable fixed points of g (2) ;for example, the basin of attraction for the stable fixedpoint at θ = 2 .
86 is the range 1 . < θ < . τ = 10 ms. Figure 11 shows g ( θ ) and g (3) ( θ ). We see that there are three stable fixedpoints for g (3) ( θ ), at θ = 1 . θ = 3 .
37, and θ = 5 . − g ( θ ), but a cobweb analysis verifiesthat there is a period-3 orbit θ = 1 . → . → . → . → · · · . These fixed points of g (3) correspond a stable 3-clusterstate for a population of oscillators, as shown in Fig-ure 4(a).We see that this map can also capture n -clusters forlarger values of n . For example, for frequency 185 Hz,the g and g (5) maps shown in Figure 12 confirm thatthere is a stable period-5 orbit θ = 1 . → . → . → . → . → . → · · · , corresponding to the stable 5-cluster state shown in Fig-ure 4(d).Finally, if we apply these identical stimuli at 300Hz, we obtain a total of four stable fixed points for g (4) , corresponding to a stable period-4 orbit for g ( θ ) θ = 0 . → . → . → . → . → · · · g () g () g (2) ( θ ) θg ( θ ) θ Fig. 10
Maps g ( θ ) and g (2) ( θ ) for Hodgkin-Huxley neuronwith stimulation frequency 150 Hz. Intersections with the di-agonal dashed line indicate fixed points of the respective map.The dotted lines show θ values for the stable fixed points ofthe g (2) map. which is equivalent to two stable period-2 orbits for g (2) θ = 0 . → . → . → · · · ,θ = 2 . → . → . → · · · , and four stable fixed points for g (4) θ = 0 . , θ = 2 . , θ = 3 . , θ = 5 . Daniel Kuelbs et al.
200 Hz stimluus, and 209, 133, and 159 neurons in Clus-ters I, II, and III, respectively, for a 260 Hz stimulus.This is consistent with the results shown in Figure 6(a).The number of neurons in each cluster for Figure 6(b)would be determined by the number of neurons whichare initially in the respective basin of attraction, as de-termined by the initial phase distribution; here, therewere no neurons with initial phases that end up in Clus-ter III.The same analysis techniques can also be used to un-derstand the dynamics of thalamic neurons subjected toperiodic pulses. Figure 15(a) shows the response func-tion f ( θ ) for thalamic neurons with the stimulus givenby (7) with u max corresponding to a current density of20 µA/cm , p = 0 . λ = 3; Figure 15(b) showsthat there is a stable 2-cluster state for a stimulationfrequency of 250 Hz, as expected from Figure 7. In this section, we consider more general stimuli, specif-ically pulses with alternating properties, as shown inFigure 16. Here, the pulses from before, that is with u max corresponding to a current density of 20 µA/cm , g () g () g (3) ( θ ) θg ( θ ) θ Fig. 11
Maps g ( θ ) and g (3) ( θ ) for Hodgkin-Huxley neuronwith stimulation frequency 100 Hz. g () g () g ( θ ) θθg (5) ( θ ) Fig. 12
Maps g ( θ ) and g (5) ( θ ) for Hodgkin-Huxley neuronwith stimulation frequency 185 Hz. p = 0 . λ = 3, will be assumed to occur at times0 , τ, τ, · · · . But now additional pulses with u max cor-responding to a current density of 10 µA/cm , λ = 3, u min = − u max /λ and p = 0 . τ , τ + τ , 2 τ + τ , · · · . Figure 17 showsthat the clustering behavior for such alternating pulseswith τ = τ / τ later. To formulate this map, we need the responsecurves for each type of pulse: the response curve f ( θ )for the pulse with u max corresponding to 20 µA/cm was already shown in Figure 9; the response curve f ( θ )for the pulse with u max corresponding to 10 µA/cm isshown in Figure 18. To find this map, suppose that westart with θ (0 + ) = 0, immediately after the start of apulse, where we assume that we have already accountedfor the effect of the pulse according to the function f ( θ ).The next pulse, of different type, comes at time τ . Upuntil time τ , the phase evolves according to ˙ θ = ω ; nalysis of Neural Clusters due to Deep Brain Stimulation Pulses 9 g () g () g () g (2) ( θ ) θg ( θ ) θθg (4) ( θ ) Fig. 13
For the Hodgkin-Huxley neurons with stimulationfrequency 300 Hz, there is a stable period-4 orbit for g , whichcorresponds to two stable period-2 orbits for g (2) , which inturn correspond to four stable fixed points for g (4) . therefore, θ ( τ − ) = θ + ωτ . (17)Treating the change in phase due to the next pulse asoccurring instantaneously, we have θ ( τ +2 ) = θ + ωτ + f ( θ + ωτ ) . (18)The system then evolves for a time τ − τ without stim-ulus, giving θ ( τ − ) = θ + ωτ + ω ( τ − τ ) + f ( θ + ωτ )= θ + ωτ + f ( θ + ωτ ) . g () g () Cluster III Cluster I Cluster IICluster III Cluster I Cluster II g (3) ( θ ) θθg (3) ( θ ) Fig. 14
Basins of attraction for the different clusters for (a)200 Hz and (b) 260 Hz stimuli. g () (b) f () (a) g (2) ( θ ) θθf ( θ ) Fig. 15 (a) Response function f ( θ ) which characterizes thephase response of thalamic neurons to a pulse with u max cor-responding to a current of 20 µA/cm , p = 0 . λ = 3.(b) Map g (2) ( θ ) for the thalamic neuron with stimulationfrequency 250 Hz, showing two stable fixed points which cor-respond to a 2-cluster state.0 Daniel Kuelbs et al. t = 0 t = τ λp t = τλp t = τ + τ u max pu min u max u ( t ) u min tp Fig. 16
Sequence of alternating pulses. frequency (Hz) θ Fig. 17
The final phases θ of Hodgkin-Huxley neurons drawnfrom an initial uniform distribution as a function of stimu-lation frequency, after 80 periods of pulses with alternatingproperties (to allow transients to decay), as described in thetext. Colors correspond to the neurons’ initial phases. -0.8-0.6-0.4-0.200.20.40.60.811.2 θf ( θ ) Fig. 18
Response function f ( θ ) which characterizes thephase response of a Hodgkin-Huxley neuron to a pulse with u max corresponding to a current density of 10 µA/cm , u min = − u max /
3, and p = 0 . At time τ , we have another pulse of the type thatstarted at t = 0, so θ ( τ + ) = θ + ωτ + f ( θ + ωτ )+ f ( θ + ωτ + f ( θ + ωτ )) . Continuing in this fashion, we obtain θ ( τ + τ − ) = θ + ω ( τ + τ ) + f ( θ + ωτ )+ f ( θ + ωτ + f ( θ + ωτ )) ,θ ( τ + τ +2 ) = θ + ω ( τ + τ ) + f ( θ + ωτ )+ f ( θ + ωτ + f ( θ + ωτ ))+ f ( θ + ω ( τ + τ ) + f ( θ + ωτ )+ f ( θ + ωτ + f ( θ + ωτ ))) θ (2 τ − ) = θ + 2 ωτ + f ( θ + ωτ )+ f ( θ + ωτ + f ( θ + ωτ ))+ f ( θ + ω ( τ + τ ) + f ( θ + ωτ )+ f ( θ + ωτ + f ( θ + ωτ ))) θ (2 τ + ) = θ (2 τ − ) + f ( θ (2 τ − )) . A useful formulation is to let G ( s ) = s + ωτ + f ( s + ωτ )+ f ( s + ωτ + f ( s + ωτ )) , (19)which gives θ ( nτ + ) = G ( n ) ( θ ) . (20)Alternatively, we can view this as a composition of twomaps: θ (0 + ) = θ ,θ ( τ +2 ) = θ + ωτ + f ( θ + ωτ ) ≡ h ( θ ) ,θ ( τ + ) = θ ( τ +2 ) + ω ( τ − τ )+ f ( θ ( τ +2 ) + ω ( τ − τ )) ≡ h ( θ ( τ +2 )) = h ( h ( θ )) = G ( θ ) . Note that we have written G , which is a map over thetime interval τ , as the composition of two maps h and h , that is, G = h ◦ h . Similar to before, we will look for fixed points of G ( n ) ,that is, solutions to θ ∗ = G ( n ) ( θ ∗ ). If (cid:12)(cid:12)(cid:12)(cid:12) ddθ (cid:12)(cid:12)(cid:12)(cid:12) θ = θ ∗ ( G ( n ) ( θ )) (cid:12)(cid:12)(cid:12)(cid:12) < , (21)then the fixed point of G ( n ) is stable. Note that therelationship between fixed points of G ( n ) and clustersis more subtle for pulses with alternating propertiesthan the relationship between fixed points of g ( n ) and n -clusters for identical pulses, because each τ -intervalfor the alternating case contains two pulses. This willbe illustrated in the following examples.Figure 19 shows h ( θ ) for u max corresponding toa current density of 20 µA/cm and h ( θ ) for u max corresponding to a current density of 10 µA/cm , for τ = 10 ms and τ = τ /
2. We notice that these func-tions are quite similar to each other. Next, we show G ( θ ) = h ( h ( θ )) and G (3) ( θ ) in Figure 20. We see thatthere is a stable period-3 orbit for G , corresponding tothree stable fixed points for G (3) . This corresponds toa 3-cluster state, as expected from Figure 17 evaluatedat 100 Hz. Here, the stable fixed points of G (3) cor-respond to a stable period-3 orbit of G , which in turncorresponds to a 3-cluster state. We note that Figure 20looks very similar to Figure 11 for identical pulses with nalysis of Neural Clusters due to Deep Brain Stimulation Pulses 11 θh ( θ ) θh ( θ ) Fig. 19
Functions h ( θ ) for u max corresponding to a currentdensity of 20 µA/cm and h ( θ ) for u max corresponding to acurrent density of 10 µA/cm , for τ = 10 ms and τ = τ/ frequency 100 Hz; however, the sequence of pulses is dif-ferent. For Figure 20, there is a “large” pulse at t = 0,a “small” pulse at t = 5 ms, another large pulse at t = 10 ms, another small pulse at t = 15 ms, anotherlarge pulse at t = 20 ms, etc. For Figure 11, there is alarge pulse at t = 0, another large pulse at t = 10 ms,another large pulse at t = 20 ms, etc, with no smallpulses.Figure 21 shows h ( θ ) for u max corresponding toa current density of 20 µA/cm and h ( θ ) for u max corresponding to a current density of 10 µA/cm , for τ = 6 .
67 ms and τ = τ /
2. As before, h and h arequite similar to each other. Figure 22 shows G ( θ ) = h ( h ( θ )) and G (2) ( θ ). We see that there are four sta-ble fixed points for G (2) ; these actually correspond to a4-cluster state, as shown in Figure 17 evaluated at 150Hz. While at first it might seem surprising that stablefixed points for G (2) correspond to a 4-cluster state, wenote that these results are similar to what we foundfor identical stimuli for a 300 Hz stimulus (or, equiv-alently, for alternating pulses with τ = 6 .
67 ms and u max = u max corresponding to a current density of20 µA/cm , τ = τ / G for alternating pulses with a stimulation frequency g () g () G (3) ( θ ) θθG ( θ ) Fig. 20
Functions G ( θ ) and G (3) ( θ ) for pulses with alternat-ing properties with u max corresponding to a current densityof 20 µA/cm and u max corresponding to a current densityof 10 µA/cm , and τ = 10 ms, τ = τ/ of 150 Hz is similar to g (2) for identical pulses with astimulation frquency of 300 Hz, and G (2) for alternatingpulses for a stimulation frequncy of 150 Hz is similar to g (4) for identical pulses with a stimulation frequency of300 Hz. These results show that we can obtain 4-clustersolutions for a population of oscillators with these al-ternating pulses; see Figure 23(a).Our analytical formalism also allows one to con-sider alternating pulses for which τ (cid:54) = τ /
2. For ex-ample, Figure 24 shows results for u max correspondingto 20 µA/cm , u max corresponding to 10 µA/cm , and τ = 0 . τ and τ = 0 . τ . Interestingly, for τ = 0 . τ there are four fixed points of the G (2) map, correspond-ing to a 4-cluster solution, but for τ = 0 . τ there areonly two fixed points of the G (2) map, correspondingto a 2-cluster solution. Figures 23(b) and (c) show thecorresponding time series for these cases. ComparingFigure 24 with the right panel of Figure 22, we deducethat if τ is treated as a bifurcation parameter, thereis a saddlenode bifurcation (this could also be calleda tangent bifurcation for the G (2) map) for τ slightlylarger than 0 . h () h () h ( θ ) θθh ( θ ) Fig. 21
Functions h ( θ ) for u max corresponding to a currentdensity of 20 µA/cm and h ( θ ) for u max corresponding to acurrent density of 10 µA/cm , for τ = 6 .
67 ms and τ = τ/ Populations of neural oscillators subjected to periodicpulsatile stimuli can display interesting clustering be-havior, in which subpopulations of the neurons are syn-chronized but the subpopulations are desynchronizedwith respect to each other. The details of the clusteringbehavior depend on the frequency and amplitude of thestimuli in a complicated way. Such clustering may be animportant mechanism by which deep brain stimulationcan lead to the alleviation of symptoms of Parkinson’sdisease and other disorders.In this paper, we illustrated how the details of clus-tering for phase models of neurons subjected to peri-odic pulsatile inputs can be understood in terms of one-dimensional maps defined on the circle. In particular,the analysis allows one to predict the number of clus-ters, their stability properties, and their basins of at-traction. Moreover, we generalized our analysis to con-sider stimuli with alternating properties, which provideadditional degrees of freedom in the design of DBS stim-uli. As part of our study, we found multiple ways to getthe same type of clustering behavior, for example by g () g () θG (2) ( θ ) G ( θ ) θ Fig. 22
Functions G ( θ ) and G (2) ( θ ) for pulses with alternat-ing properties with u max corresponding to a current densityof 20 µA/cm and u max corresponding to a current densityof 10 µA/cm , and τ = 6 .
67 ms, τ = τ/ using identical pulses or pulses with alternating prop-erties, or from stimuli with different parameters suchas stimulation frequency or the time spacing betweenpulses with alternating properties. Such clustering oc-curs through the use of a single stimulation electrode,unlike coordinated reset which requires multiple elec-trodes. We expect that the same clustering behavior canalso be obtained for different amplitudes of the pulses,cf. [41]. We believe that the analysis techniques usedin this paper can be useful for identifying a collectionof stimuli which give the same desireable clustering dy-namics for a population of neurons, which will make iteasier to find stimuli which are effective while minimiz-ing the severity of side effects for DBS treatments.Our analysis assumed certain properties of a neuralpopulation: all neurons are identical, they all receivethe same input, they are uncoupled, and there is nonoise. For real neural populations, none of these as-sumptions would be valid. We also assumed that thephase models accurately capture the dynamics of theneurons, which is only true for sufficiently small inputs;see, for example, [39]. However, we believe that the re-sults presented here form an important baseline for the nalysis of Neural Clusters due to Deep Brain Stimulation Pulses 13 (b)(c)(a) τ = 0 . τ : 4 clusters tθθ τ = 0 . τ : 2 clusters ttθ τ = 0 . τ : 4 clusters Fig. 23
Time series showing the phases of Hodgkin-Huxleyneurons drawn from an initial uniform distribution with al-ternating pulses with u max corresponding to 20 µA/cm and u max corresponding to 10 µA/cm , for τ = 6 .
67 ms and (a) τ = 0 . τ , (b) τ = 0 . τ , (c) τ = 0 . τ . Four clusters form for(a) and (b), while only two clusters form for (c). analysis of more realistic neural populations stimulatedby periodic pulses. We note that the effect of noise onperiodically forced neural populations has been consid-ered in [41], which shows that for weak noise and longtimes, the number of neurons in each cluster is roughlythe same. We expect that similar results will hold forneurons in the presence of weak noise subjected to al-ternating stimuli.Our hope is that the techniques in this paper willhelp to guide the design of stimuli for the treatment ofParkinson’s disease and other disorders. We believe thatthe use of pulses with alternating properties is particu-larly worthy of further investigation, since it representsa larger class of stimuli than has been considered inprevious studies. g () g () θG (2) ( θ ) θG (2) ( θ ) Fig. 24
Function G (2) ( θ ) for alternating pulses with u max corresponding to 20 µA/cm and u max corresponding to10 µA/cm , with τ = 6 .
67 ms, and (left) τ = 0 . τ and (right) τ = 0 . τ . This research grew out of the Research Mentorship Pro-gram at the University of California, Santa Barbaraduring summer 2018. We thank Dr. Lina Kim for pro-viding the opportunity for Daniel and Jacob to con-duct this research as high school students, and for TimMatchen for guidance on the project.
Appendix: Neuron Models
In this Appendix, we give details of the neural modelsused in the main text.
Thalamic neuron model
The full thalamic neuron model is given by:˙ V = − I L − I Na − I K − I T + I b C m + u ( t ) , ˙ h = h ∞ − hτ h , ˙ r = r ∞ − rτ r , where h ∞ = 1 / (1 + exp(( V + 41) / ,r ∞ = 1 / (1 + exp(( V + 84) / ,α h = 0 .
128 exp( − ( V + 46) / ,β h = 4 / (1 + exp( − ( V + 23) / ,τ h = 1 / ( α h + β h ) ,τ r = (28 + exp( − ( V + 25) / . ,m ∞ = 1 / (1 + exp( − ( V + 37) / ,p ∞ = 1 / (1 + exp( − ( V + 60) / . ,I L = g L ( v − e L ) ,I Na = g Na ( m ∞ ) h ( v − e Na ) ,I K = g K ((0 . − h )) )( v − e K ) ,I T = g T ( p ∞ ) r ( v − e T ) . The parameters for this model are C m = 1 µF/cm , g L = 0 . mS/cm , e L = − mV ,g Na = 3 mS/cm , e Na = 50 mV, g K = 5 mS/cm ,e K = − mV , g T = 5 mS/cm , e T = 0 mV ,I b = 5 µA/cm . Hodgkin-Huxley neuron model
The full Hodgkin-Huxley model is given by:˙ V = ( I b − ¯ g Na h ( V − V Na ) m − ¯ g K ( V − V K ) n − ¯ g L ( V − V L )) /c + u ( t ) , ˙ m = a m ( V )(1 − m ) − b m ( V ) m , ˙ h = a h ( V )(1 − h ) − b h ( V ) h , ˙ n = a n ( V )(1 − n ) − b n ( V ) n , where a m ( V ) = 0 . V + 40) / (1 − exp( − ( V + 40) / ,b m ( V ) = 4 exp( − ( V + 65) / ,a h ( V ) = 0 .
07 exp( − ( V + 65) / ,b h ( V ) = 1 / (1 + exp( − ( V + 35) / ,a n ( V ) = 0 . V + 55) / (1 − exp( − ( V + 55) / ,b n ( V ) = 0 .
125 exp( − ( V + 65) / , The parameters for this model are V Na = 50 mV , V K = − mV , V L = − . mV , ¯ g Na = 120 mS/cm , ¯ g K = 36 mS/cm ,barg L = 0 . mS/cm , I b = 10 µA/cm ,c = 1 µF/cm . References
1. Adamchic, I., Hauptmann, C., Barnikol, U.B., Pawel-czyk, N., Popovych, O., Barnikol, T.T., Silchenko, A.,Volkmann, J., Deuschl, G., Meissner, W.G., Maarouf,M., Sturm, V., Freund, H.J., Tass, P.A.: Coordinatedreset neuromodulation for Parkinson’s disease: proof-of-concept study. Mov. Disord. (13), 1679–1684 (2014)2. Benabid, A., Benazzous, A., Pollak, P.: Mechanisms ofdeep brain stimulation. Movement Disorders (SUPPL.3), 19–38 (2002)3. Benabid, A.L., Pollak, P., Gervason, C., Hoffmann, D.,Gao, D.M., Hommel, M., Perret, J.E., Rougemont, J.D.:Long-term suppression of tremor by chronic stimulationof the ventral intermediate thalamic nucleus. The Lancet , 403–406 (1991)4. Best, D., Fisher, N.: Efficient simulation of the von Misesdistribution. J. R. Stat. Soc. Ser. C. Appl. Stat. , 152–157 (1979)5. Brown, E., Moehlis, J., Holmes, P.: On the phase reduc-tion and response dynamics of neural oscillator popula-tions. Neural Comp. , 673–715 (2004)6. Buhmann, C., Huckhagel, T., Engel, K., Gulberti, A.,Hidding, U., Poetter-Nerger, M., Goerendt, I., Ludewig,P., Braass, H., Choe, C., et al: Adverse events in deepbrain stimulation: a retrospective long-term analysis ofneurological psychiatric and other occurrences. PLoSOne , e0178984 (2017)7. Chen, C., Litvak, V., Gilbertson, T., Kuhn, A., Lu, C.,Lee, S., Tsai, C., Tisch, S., Limousin, P., Hariz, M.,Brown, P.: Excessive synchronization of basal ganglianeurons at 20 hz slows movement in Parkinson’s disease.Experimental Neurology , 214–221 (2007)8. Cyron, D.: Mental side effects of deep brain stimula-tion (DBS) for movement disorders: the futility of denial.Frontiers in Integrative Neuroscience , 17 (2016)9. Ermentrout, G.: Simulating, Analyzing, and Animat-ing Dynamical Systems: A Guide to XPPAUT for Re-searchers and Students. SIAM, Philadelphia (2002)10. Ermentrout, G.B., Terman, D.H.: Mathematical Founda-tions of Neuroscience. Springer, Berlin (2010)11. Guckenheimer, J.: Isochrons and phaseless sets. J. Math.Biol. , 259–273 (1975)12. Hammond, C., Bergman, H., Brown, P.: Pathological syn-chronization in parkinsons disease: networks, models andtreatments. Trends in Neurosciences , 357–364 (2007)13. Hodgkin, A.L., Huxley, A.F.: A quantitative descriptionof membrane current and its application to conductionand excitation in nerve. J. Physiol. , 500–544 (1952)14. Kuncel, A.M., Grill, W.M.: Selection of stimulus param-eters for deep brain stimulation. Clin. Neurophysiol. (11), 2431–2441 (2004)15. Kuramoto, Y.: Chemical Oscillations, Waves, and Tur-bulence. Springer, Berlin (1984)16. Levy, R., Hutchison, W., Lozano, A., Dostrovsky, J.:High-frequency synchronization of neuronal activity inthe subthalamic nucleus of Parkinsonian patients withlimb tremor. The Journal of Neuroscience , 7766–7775(2000)17. L¨ucken, L., Yanchuk, S., Popovych, O., Tass, P.: Desyn-chronization boost by non-uniform coordinated resetstimulation in ensembles of pulse-coupled neurons. Fron-tiers in Computational Neuroscience (2013)18. Lysyansky, B., Popovych, O., Tass, P.: Desynchronizinganti-resonance effect of m: n ON-OFF coordinated resetstimulation. Journal of Neural Engineering , 036019(2011)nalysis of Neural Clusters due to Deep Brain Stimulation Pulses 1519. Lysyansky, B., Popovych, O., Tass, P.: Optimal numberof stimulation contacts for coordinated reset neuromod-ulation. Frontiers in Neuroengineering (2013). Article520. Matchen, T., Moehlis, J.: Phase model-based neuron sta-bilization into arbitrary clusters. Journal of Computa-tional Neuroscience , 363–378 (2018)21. Monga, B., Moehlis, J.: Phase distribution control of apopulation of oscillators. Physica D , 115–129 (2019)22. Monga, B., Moehlis, J.: Supervised learning algorithmsfor control of underactuated dynamical systems. PhysicaD , 132621 (2020)23. Monga, B., Wilson, D., Matchen, T., Moehlis, J.: Phasereduction and phase-based optimal control for biologicalsystems: a tutorial. Biological Cybernetics , 11–46(2019)24. Montgomery, E.: Deep Brain Stimulation Programming:Principles and Practice. Oxford University Press, Oxford(2010)25. Moro, E., Esselink, R.J., Xie, J., Hommel, M., Benabid,A.L., Pollak, P.: The impact on Parkinson’s disease ofelectrical parameter settings in STN stimulation. Neu-rology (5), 706–713 (2002)26. Nabi, A., Mirzadeh, M., Gibou, F., Moehlis, J.: Minimumenergy desynchronizing control for coupled neurons. J.Comp. Neuro. , 259–271 (2013)27. Netoff, T., Schwemmer, M., Lewis, T.: Experimentallyestimating phase response curves of neurons: theoreticaland practical issues. In: Phase Response Curves in Neu-roscience, pp. 95–129. Springer (2012)28. Rinzel, J., Ermentrout, G.B.: Analysis of neural excitabil-ity and oscillations. In: C. Koch, I. Segev (eds.) Methodsin Neuronal Modeling, pp. 251–291. MIT Press (1998)29. Rizzone, M., Lanotte, M., Bergamasco, B., Tavella, A.,Torre, E., Faccani, G., Melcarne, A., Lopiano, L.: Deepbrain stimulation of the subthalamic nucleus in Parkin-son’s disease: effects of variation in stimulation param-eters. J. Neurol. Neurosurg. Psychiatr. (2), 215–219(2001)30. Rubin, J., Terman, D.: High frequency stimulation ofthe subthalamic nucleus eliminates pathological thala-mic rhythmicity in a computational model. Journal ofComputational Neuroscience (3), 211–235 (2004)31. Savica, R., Stead, M., Mack, K., Lee, K., Klassen, B.:Deep brain stimulation in Tourette syndrome: A descrip-tion of 3 patients with excellent outcome. Mayo ClinicProceedings , 59–62 (2012)32. Schnitzler, A., Gross, J.: Normal and pathological oscil-latory communication in the brain. Nature Review Neu-roscience , 285–296 (2005)33. Tass, P.A.: A model of desynchronizing deep brain stimu-lation with a demand-controlled coordinated reset of neu-ral subpopulations. Biol. Cybern. (2), 81–88 (2003)34. Tass, P.A.: Desynchronization by means of a coordinatedreset of neural sub-populations - A novel technique fordemand-controlled deep brain stimulation. Progress ofTheoretical Physics Supplement , 281–296 (2003)35. Uhlhaas, P., Singer, W.: Neural synchrony in brain disor-ders: relevance for cognitive dysfunctions and pathophys-iology. Neuron , 155–168 (2006)36. Volkmann, J., Herzog, J., Kopper, F., Deuschl, G.: Intro-duction to the programming of deep brain stimulators.Mov. Disord.
17 Suppl 3 , S181–187 (2002)37. Wilson, C.J., Beverlin, B., Netoff, T.: Chaotic desyn-chronization as the therapeutic mechanism of deep brainstimulation. Front. Syst. Neurosci. , 50 (2011) 38. Wilson, D.: Optimal open-loop desynchronization of neu-ral oscillator populations. J. Math. Biol. (2020)39. Wilson, D., Ermentrout, B.: Greater accuracy and broad-ened applicability of phase reduction using isostable coor-dinates. Journal of Mathematical Biology (1-2), 37–66(2018)40. Wilson, D., Moehlis, J.: Optimal chaotic desynchroniza-tion for neural populations. SIAM J. Appl. Dyn. Syst. , 276–305 (2014)41. Wilson, D., Moehlis, J.: Clustered desynchronizationfrom high-frequency deep brain stimulation. PLoS Com-put. Biol. (12), e1004673 (2015)42. Winfree, A.: Biological rhythms and the behavior of pop-ulations of coupled oscillators. J. Theor. Biol.16