Modeling shifts in the rate and pattern of subthalamopallidal network activity during deep brain stimulation

Abstract

Deep brain stimulation (DBS) of the subthlamic nucleus (STN) represents an effective treatment for medically refractory Parkinson's disease; however, understanding of its effects on basal ganglia network activity remains limited. We constructed a computational model of the subthalamopallidal network, trained it to fit in vivo recordings from parkinsonian monkeys, and evaluated its response to STN DBS. The network model was created with synaptically connected single compartment biophysical models of STN and pallidal neurons, and stochastically defined inputs driven by cortical beta rhythms. A least mean square error training algorithm was developed to parameterize network connections and minimize error when compared to experimental spike and burst rates in the parkinsonian condition. The output of the trained network was then compared to experimental data not used in the training process. We found that reducing the influence of the cortical beta input on the model generated activity that agreed well with recordings from normal monkeys. Further, during STN DBS in the parkinsonian condition the simulations reproduced the reduction in GPi bursting found in existing experimental data. The model also provided the opportunity to greatly expand analysis of GPi bursting activity, generating three major predictions. First, its reduction was proportional to the volume of STN activated by DBS. Second, GPi bursting decreased in a stimulation frequency dependent manner, saturating at values consistent with clinically therapeutic DBS. And third, ablating STN neurons, reported to generate similar therapeutic outcomes as STN DBS, also reduced GPi bursting. Our theoretical analysis of stimulation induced network activity suggests that regularization of GPi firing is dependent on the volume of STN tissue activated and a threshold level of burst reduction may be necessary for therapeutic effect.

Fig. 1

Subthalamopallidal network architecture. The excitatory connections are shown on the left (functional column 2), and inhibitory connections are shown on the right. Inputs to the model are shown as excitatory cortical inputs to STN and inhibitory striatal inputs to GPe and GPi

Fig. 2

Network input properties. (a) Each cell receives seven external synaptic inputs that are stochastically determined by a common excitatory cortical rhythm. p(i,j) is the probability that input j to cell i responds on any given cycle of the cortical rhythm with an strong increase in activity level. (b) Cumulative probability of response magnitude. The probability of a weak response during an excitatory cortical phase decreases as p increases. Each curve corresponds to a value of p (the probability of a single synapse responding strongly) ranging from 0.05 to 0.55 in steps of 0.1. (c) Examples of total excitatory cortical synaptic input over time in an individual STN cell for three levels of p

Fig. 3

Training algorithm. The training algorithm takes an initial randomized parameter vector (x∈R^m) through a series of steps to minimize the error between the model output values (y∈Rⁿ) and target values (y_T∈Rⁿ). Each step in parameter space is determined by the numerical approximation of partial derivatives of the system

Fig. 4

Network training examples. Columns represent the training of two of the 25 model instances used in the Results. (a, c) Error (||y − y_T||) is reduced in a non-smooth manner through the two phases of training in each trial (open circle—step taken, filled circle—no improvement, star—random bump. The first phase of training considered only spike rates, the second phase of training (identified by black overline) considered both spike rates and burst rates. (b, d) As the algorithm takes steps through parameter space, the distance between the current parameter vector and the final parameter vector for that instance of the model (||x − x_∞||) increases and decreases, eventually reaching zero

Fig. 5

Single GPi neuron model. Neural response to varying levels of cortico-striatal input. (a, c) Inputs derived from the MPTP (p=.25) state of the network. (b, d) Inputs derived from the Normal (p=.05) state of the network. (a, b) Upper trace shows the transmembrane voltage and lower trace shows total synaptic current to the GPi neuron from a single trial. (c, d) One second traces were aligned on the phase of the cortical rhythm and averaged over many trials (smoothed with 20 ms sliding window). The upper trace shows the instantaneous spiking rate and the lower trace shows total synaptic current to the GPi neuron

Fig. 6

Modeling the MPTP state. Raster diagrams show spike times as small vertical hashes. In each row, detected bursts are indicated by horizontal lines through the corresponding burst spikes for that cell. (a) Population raster approximations for experimental data. Each row (1 s duration) corresponds to spike times taken from segments of separate experimental recordings (GPe or GPi) from MPTP treated monkeys (Hashimoto et al. 2003). A₂ expands 500 ms of the four cells indicated. (b) Population spike time rasters for one second of activity from the network model in the MPTP state. B₂ illustrates the transmembrane voltage for 500 ms of the four cells indicated

Fig. 7

Comparison of MPTP and Normal states. (a) Spike rate and (b) burst rate shifts in response to changes in input parameters from the MPTP to Normal state. Each point represents population means for an individual instance of the network model. The target value used for training the network to the MPTP state is indicated by a triangle on each left axis. The mean and standard deviation of the spike and burst rates from normal and MPTP monkeys are displayed (gray markers and error bars—Wichmann and Soares (2006); black markers and error bars—Hahn et al. (2008)). (c) The mean and standard deviation for spike (C₁) and burst (C₂) rates across all instances of the BG network in MPTP (open box) and Normal (grey box) states

Fig. 8

Modeling the STN DBS state. Raster diagrams show spike times as small vertical hashes. In each row, detected bursts are indicated by horizontal lines through the corresponding burst spikes for that cell. (a) Population raster approximations from experimental data. Each row (1 s duration) corresponds to spike times taken from a separate experimental recording (GPe or GPi) from MPTP treated monkeys during STN DBS (Hashimoto et al. 2003). (b) Population rasters for one second of network model activity during STN DBS condition. STN cells with directly stimulated efferent output are indicated by a vertical bar

Fig. 9

Comparison of MPTP and STN DBS states. (a) Spike rate and (b) burst rate shifts for the model in MPTP state with simulated STN DBS. Each point represents population means for an individual instance of the network model. The target value used for training the network to the MPTP state is indicated by a triangle on each left axis. The mean and standard deviation of the spike and burst rates from MPTP monkeys are displayed (gray markers and error bars—Wichmann and Soares (2006); black markers and error bars—Hahn et al. (2008)). (c) The mean and standard deviation for spike (C₁) and burst (C₂) rates across all instances of the network model before (open box) and during (grey box) simulated STN DBS

Fig. 10

GPi burst rate changes. (a) Population mean burst rate (averaged from all 25 model instances) in GPi as a function of STN DBS frequency. (b) Average GPi burst rate from each model instance during STN DBS. (c) Population mean burst rate (averaged from all 25 model instances) in GPi as a function of the percentage of the STN neuron population directly stimulated (grey bars) or lesioned (hatched bars). (d) Average GPi burst rate from each model instance during STN DBS. Color coding for the model instances is consistent for (b) and (d)