Computational Analysis of Subthalamic Nucleus and Lenticular Fasciculus Activation During Therapeutic Deep Brain Stimulation

Svjetlana Miocinovic, Martin Parent, Christopher R. Butson, Philip J. Hahn, Gary S. Russo, Jerrold L. Vitek, Cameron C. McIntyre


The subthalamic nucleus (STN) is the most common target for the treatment of Parkinson’s disease (PD) with deep brain stimulation (DBS). DBS of the globus pallidus internus (GPi) is also effective in the treatment of PD. The output fibers of the GPi that form the lenticular fasciculus pass in close proximity to STN DBS electrodes. In turn, both STN projection neurons and GPi fibers of passage represent possible therapeutic targets of DBS in the STN region. We built a comprehensive computational model of STN DBS in parkinsonian macaques to study the effects of stimulation in a controlled environment. The model consisted of three fundamental components: 1) a three-dimensional (3D) anatomical model of the macaque basal ganglia, 2) a finite element model of the DBS electrode and electric field transmitted to the tissue medium, and 3) multicompartment biophysical models of STN projection neurons, GPi fibers of passage, and internal capsule fibers of passage. Populations of neurons were positioned within the 3D anatomical model. Neurons were stimulated with electrode positions and stimulation parameters defined as clinically effective in two parkinsonian monkeys. The model predicted axonal activation of STN neurons and GPi fibers during STN DBS. Model predictions regarding the degree of GPi fiber activation matched well with experimental recordings in both monkeys. Only axonal activation of the STN neurons showed a statistically significant increase in both monkeys when comparing clinically effective and ineffective stimulation. Nonetheless, both neural targets may play important roles in the therapeutic mechanisms of STN DBS.


Deep brain stimulation (DBS) has become an established clinical therapy for advanced Parkinson’s disease (PD) (Obeso et al. 2001). Chronic high-frequency electrical stimulation of subcortical structures can provide >50% improvement in clinical ratings of motor symptoms (Walter and Vitek 2004). However, the impressive clinical efficacy achieved by DBS has occurred without a clear understanding of the therapeutic mechanisms of action. In addition, a number of adverse effects can be generated by DBS including sensorimotor impairments, involuntary movements (stimulation-induced dyskinesias), as well as speech, mood, and cognitive disturbances (Okun et al. 2005; Volkmann et al. 2002). Often these side effects can be avoided or alleviated with proper adjustment of the stimulation settings (Krack et al. 2003). However, further improvements in the engineering design and clinical implementation of DBS technology will rely on addressing a number of questions on the effects of DBS on the nervous system. The fundamental goal of this project was to enhance our understanding of the target neural elements of the stimulation.

The subthalamic nucleus (STN) represents the most common anatomical target for DBS treatment of PD (Limousin et al. 1998). Electrodes placed in the STN are surrounded by several neural types (local projection neurons and their axons, fibers of passage, afferent inputs, etc.), but knowledge of the response properties of these different neural types to DBS is limited. Therapeutic DBS electrode contacts are typically located in the region of the dorsal STN, lenticular fasciculus (LF or Forel’s field H2), and zona incerta (Hamel et al. 2003; Nowinski et al. 2005; Saint-Cyr et al. 2002; Starr et al. 2002; Voges et al. 2002; Yelnik et al. 2003; Zonenshayn et al. 2004). STN projection neurons send their highly collateralized axons to the globus pallidus, striatum, and substantia nigra (Sato et al. 2000). The LF courses just dorsal to the STN and is composed of fibers from the internal segment of the globus pallidus (GPi), which carry output from the basal ganglia to the thalamus (Parent and Parent 2004; Parent et al. 2001). Given that GPi DBS provides similar therapeutic benefits as STN DBS (Burchiel et al. 1999; Obeso et al. 2001; Rodriguez-Oroz et al. 2005), both STN projection neurons and pallidothalamic (GPi) fibers represent viable candidates as the therapeutic target of the stimulation (Parent and Parent 2004). However, the axons of local projection neurons and fibers of passage respond at similar extracellular stimulation thresholds (McIntyre and Grill 1999; Nowak and Bullier 1998; Ranck 1975), making it difficult to determine which neuron types are activated during STN DBS.

The extent of neural activation generated by extracellular stimulation depends on the stimulation parameters, electrode and tissue electrical properties, and the position and orientation of neural elements with respect to the electrode (McIntyre et al. 2004a; Ranck 1975; Tehovnik 1996). To address these issues in the context of STN DBS, we developed a comprehensive computer model with detailed representation of the three-dimensional (3D) neuroanatomy, the time-dependent electric field generated by DBS electrodes, and the underlying biophysics that regulate the neural response to stimulation. We tested our hypothesis that both STN projection neurons and GPi fibers of passage are activated during clinically effective STN DBS. This would imply that both neural types could play a role in the therapeutic mechanisms of STN DBS and prompt further study into electrode localization and stimulation parameter selection techniques for optimizing the stimulation to individual subjects.

We customized our modeling framework to analyze neural activation in two parkinsonian macaques implanted with chronic scaled clinical DBS systems (Hashimoto et al. 2003). Our goal was to theoretically reproduce the experimental effects of STN DBS that improved parkinsonian symptoms (bradykinesia and rigidity) in the two monkeys. Our first aim was to quantify the proportion of STN projection neurons and GPi fibers of passage that were activated during clinically effective and ineffective stimulation. We found a significant increase in axonal activation of STN projection neurons during clinically effective compared with ineffective stimulation. Considerable GPi fiber activation was observed in only one of the two monkeys. Single unit extracellular recordings of short latency, presumably antidromic, GPi activation during STN DBS in both monkeys support our model predictions. The second aim was to analyze how changes in electrode location affected neural activation. We found that submillimeter shifts can alter neural activation, highlighting the importance of precise electrode positioning. The final aim was to study the effects of stimulation-induced trans-synaptic inhibition on somatic and axonal firing in STN projection neurons during STN DBS. We found that somatic firing was reduced during DBS compared with the spontaneous prestimulation activity, as seen experimentally (Filali et al. 2004; Meissner et al. 2005; Welter et al. 2004). However, axonal activation was largely unaffected by the somatic inhibition and the axon easily fired at the stimulation frequency (McIntyre et al. 2004a). Preliminary portions of this work have been presented in abstract form (Miocinovic et al. 2004).


The implementation of chronic scaled clinical DBS systems in parkinsonian macaques provides the foundation for detailed study of the therapeutic mechanisms of the stimulation (Elder et al. 2005; Hashimoto et al. 2003). We coupled our experimental results with detailed computer modeling to provide new insight into the cellular effects of STN DBS. We developed computational models customized to two parkinsonian macaques implanted with STN DBS systems. Each model consisted of three fundamental components: 1) a 3D anatomical model of the monkey basal ganglia, 2) a finite element model of the DBS electrode and electric field transmitted to the tissue medium, and 3) multicompartment biophysical models of reconstructed STN projection neurons, GPi fibers of passage, and internal capsule fibers of passage. Populations of the neuron models were placed within context of the 3D anatomical model and the histologically defined DBS electrode positions. The DBS electric field model was applied to the neuron models, thereby allowing for theoretical prediction of the neural response to the stimulation.

Brain atlas

We built a 3D reconstruction of the basal ganglia for monkey R7160 (Macaca mulatta) (Hashimoto et al. 2003). Digital atlas templates of the macaque brain (Martin and Bowden 2000) were warped in Edgewarp v3.28 (Bookstein 1990) to histological brain slices to identify borders of nuclei not visible directly in the Nissl-stained sections (Fig. 1). Nuclei of interest were outlined in the 2D warped atlas slices, spaced in 1-mm increments. 3D volumes were created by interpolating between these contour lines using the graphical modeling program Rhinoceros v3.0 (McNeal and Associates, Seattle, WA). The resulting 3D brain atlas provided an anatomically realistic virtual space to position the DBS electrode and the neuron models (Fig. 1). The DBS electrode trajectory was reconstructed from the histological slices, and it was verified that it matched the electrode position drawing in the top row of Fig. 1 from Hashimoto et al. (2003). The brain of the second monkey (R370) was not available; therefore we used the same 3D atlas, but manually adjusted the DBS electrode position until the sagittal cross-section from the 3D atlas matched the rendering in the bottom row of Fig. 1 from Hashimoto et al. (2003).

FIG. 1.

Three-dimensional (3D) reconstruction of the basal ganglia for monkey R7160. Atlas templates (top right, Martin and Bowden 2000) were warped to Nissl-stained brain slices (top left) using Edgewarp (bottom, Bookstein 1990). Bottom right image shows a warped atlas template that now matches the histological slice, allowing definition of the nuclear borders. Outlines of nuclei from a series of warped atlas templates were interpolated into 3D volumes (right). A 4-contact deep brain stimulation (DBS) electrode was positioned in the posterior subthalamic nucleus (STN) of the customized brain atlas based on the histologically defined electrode location. Each leg of the 3D scale bar is 5 mm (D, dorsal; P, posterior; M, medial). Put, putamen; Caud, caudate nucleus; GPe, globus pallidus pars externa; GPi, globus pallidus pars interna; OT, optic tract; STN, subthalamic nucleus; Th, thalamus including thalamic reticular nucleus.

Neuron geometries

The 3D geometry of a longtailed macaque’s STN projection neuron (Macaca fascicularis) was reconstructed using biotin dextran amine labeling and axonal tracing as described by Sato et al. (2000) (Fig. 2). The neuron geometry was defined using Neurolucida (MicroBrightField, Williston, VT) and converted for display in Rhinoceros and our 3D brain atlas. Axonal tracing experiments (Sato et al. 2000) have revealed that STN projection neurons course either dorsally along the ventral border of the thalamus or ventrally along the lateral border of the STN on their way to the globus pallidus. To account for this anatomical variability, the original neuron reconstruction was used to create two additional STN neuron geometries with alternative axonal paths (Fig. 3A). The GPi axon geometry was based on the description of lenticular fasciculus trajectory from Parent and Parent (2004). The LF fibers emerged dorsomedially from the GPi, crossed the IC at a level dorsal to the STN, turned caudally to run along the dorsal border of the STN in Forel’s field H2, joined the ansa lenticularis in Forel’s field H, continued through Forel’s field H1, and terminated in the ventral thalamus (Fig. 3B).

FIG. 2.

Multicompartment cable model of an STN projection neuron. A: 3D reconstruction of a macaque STN neuron generated in Neurolucida. B: neuron geometry was imported into the NEURON simulation environment. Soma and neuronal processes were divided into compartments and modeled as a series of resistors and capacitors. C: membrane lipid bilayer was represented as a capacitor and ion channels as variable resistors. Na, fast sodium channel; NaP, persistent sodium channel; KDR, potassium delayed rectifier; Kv31, potassium fast rectifier; sKCa, small conductance calcium-activated potassium channel; Ih, hyperpolarization-activated cation channel; CaT, low voltage–activated calcium channel; CaN and CaL, high-voltage activated calcium channel; L, leak channel; Cm, membrane capacitance (see Gillies and Willshaw 2006 for details on the neuron membrane dynamics).

FIG. 3.

Neural populations and DBS electrode in the context of 3D neuroanatomy. A: 3 types of STN projection neurons share the same soma and dendritic tree but have different axonal trajectories. B: GPi fibers form the lenticular fasciculus on their way to the ventral thalamus. C: population of internal capsule fibers. Diameters of neuronal processes were increased 10 times for figure rendering. Each leg of 3D scale bar represents 1 mm (D, dorsal; P, posterior; M, medial).

The internal capsule defines the lateral border of the STN. Consequently, motor evoked responses from activation of the corticospinal tract (CST) can be elicited with relatively low thresholds during STN stimulation. To provide a gross level of model validation and connection to behavioral measurements, we incorporated CST axon trajectories into our anatomical framework. From the level of dorsal thalamus, the CST fibers coursed ventrally at an ∼20° anterior-to-posterior angle (Fig. 3C).

The three STN neuron geometries, along with the GPi and CST axon trajectories, were placed within the 3D atlas in their respective anatomically realistic positions and orientations (Fig. 3). Neural populations of the STN neurons, GPi fibers, and CST fibers were created by copying the five basic geometries and distributing them randomly within their respective anatomical boundaries. In total, three such populations were created by randomly shifting neurons by ±250 μm in any direction and manually repositioning those that ended up outside their anatomical boundaries. After the populations were established, a DBS electrode was placed within the STN. The cells were defined as damaged if any part of the axon, soma, or substantial portion of the dendritic tree intersected with the electrode. We started with 100 STN neurons, 80 GPi fibers, and 70 CST fibers, removed all the damaged cells, and randomly removed additional cells so that the final count for each of the three populations was 50 cells.

Neuron biophysics

We built multicompartment cable models of STN projection neurons, GPi fibers of passage, and CST fibers of passage using NEURON v5.8 (Hines and Carnevale 1997). The STN soma, initial segment, and dendrites contained channel dynamics of the rat subthalamic projection neuron originally developed by Gillies and Willshaw (2006), where dendritic channel conductance densities scale linearly with distance from the soma (Fig. 2C). To better fit the membrane dynamics to our neuron geometry and the firing pattern of STN neurons in a parkinsonian monkey, we modified the original conductances in the following way: the calcium-activated potassium channel was increased by 80%, the fast potassium rectifier was increased by 20% (soma and initial segment only), and the fast acting sodium channel was decreased by 25% in the soma and initial segment and by 35% in the dendrites. The model temperature was set to 36°C to simulate in vivo conditions. The STN neuron model had a resting potential of approximately −54 mV and spontaneous tonic firing of 32 Hz, consistent with the rate of 36.5 ± 10.8 Hz recorded by Wichmann et al. (2002) in the parkinsonian macaque. The neuron firing rate increased with increasing amplitude of intracellular depolarizing current. The slope of the model frequency-intensity (F-I) curve was 0.26 Hz/pA, lower than slopes recorded in rat brain slices (∼0.54 Hz/pA in Bevan and Wilson 1999). Despite this discrepancy, the model neuron was able to fire at >200 Hz, well above the range of clinically used DBS frequencies (100–185 Hz). The input resistance measured as the slope of the I-V curve while the neuron was in a hyperpolarized state (−59.6 mV resting potential) and injected with small currents (−0.2 to 0.05 nA) was 35 MΩ, twice the value of resistances recorded in rat STN neurons in vivo (mean, 18 MΩ; range, 9–28 MΩ in Kita et al. 1983). This discrepancy could be caused by a lack of synaptic inputs into the dendritic tree, differences between monkey and rat STN neuron morphology, or imperfect measurement of dendritic diameters during histological 3D reconstruction resulting in underestimate of the surface area. Average rat in vitro measurements range from 146 to 200 MΩ (Beurrier et al. 1999; Nakanishi et al. 1987), placing the model neuron within the range of in vivo and in vitro recordings. The membrane time constant estimated from the 1/e point of the membrane potential change induced by a low-intensity hyperpolarizing current pulse was 6 ms, consistent with in vivo rat measurements (6 ± 2 ms in Kita et al. 1983).

The axon of the STN neuron and the GPi and CST axons were based on the myelinated axon model originally described by McIntyre et al. (2002). The fiber diameter was set to 2 μm, and the individual segment dimensions and ion channel conductances were adjusted as previously described (McIntyre et al. 2004a). The axonal resting potential was set to −65 mV. The GPi axon was induced to fire tonically at 80 Hz by injecting short current pulses at the GPi terminal end of the axon to mimic parkinsonian macaque GPi firing rate of 80.1 ± 20.2 Hz (Wichmann et al. 2002). There were no synaptic connections or interactions between any of the different neurons in the model system.

In some simulations, trans-synaptic GABAa input to the STN neurons was added to evaluate the influence of stimulation-induced trans-synaptic conductances on somatic activity during high-frequency stimulation. This was simplistically modeled as an inhibitory postsynaptic current applied only to the central compartment of the cell body. The time-course and amplitude of the GABAa synaptic conductance was modeled with experimentally defined first-order kinetics of the transmitter binding to postsynaptic receptors (Destexhe et al. 1994a,b). Because afferent inputs (i.e., axonal terminals) are typically more excitable than the passing axons (Anderson et al. 2006; McIntyre et al. 2004a), in all applicable simulations, we assumed that synaptic inhibition was always activated by each extracellular stimulus pulse, regardless of the cell position with respect to the electrode.

Electric field model

An axisymmetric finite element model (FEM) of the DBS electrode was created in FEMLAB 3.1 (COMSOL, Burlington, MA) to calculate potentials generated in the tissue medium by the electrode. The stimulating lead was a scaled-down version of the chronic DBS electrode used in humans (Model 3387, Medtronic, Minneapolis, MN) and consisted of four metal contacts each with a diameter of 0.75 mm, height of 0.50 mm, and separation between contacts of 0.50 mm. The volume conductor was 5 × 5 cm in size and grounded at the boundaries. The bulk conductivity of the tissue medium was 0.2 S/m (Ranck 1963), and a 0.25-mm encapsulation sheath (0.1 S/m) surrounded the electrode shaft (Butson et al. 2006; Grill and Mortimer 1994). The electrode shaft was modeled as an insulator (1e−6 S/m) and each electrode contact as a conductor (1e6 S/m). Voltage sources were specified at two electrode contacts for bipolar stimulation. The model had variable resolution mesh with a total of 32,640 mesh elements that increased in size with increasing distance from the electrode. Voltage values within the volume were determined from the Poisson equation, which was solved using direct matrix inversion (UMFPACK solver) (Fig. 4).

FIG. 4.

Field-neuron model of STN DBS. A: finite element model (FEM) voltage solution for 1-V bipolar stimulus overlaid with a sagittal brain cross-section from monkey R7160. A 250-μm-thick encapsulation layer surrounds the DBS electrode. Str, striatum; GPe, globus pallidus pars externa; GPi, globus pallidus pars interna; OT, optic tract; Th, thalamus including reticular thalamic nucleus; STN, subthalamic nucleus; Sn, substantia nigra. B: extracellular potentials generated by the electrode create transmembrane polarization along the STN projection neuron. Neural compartments are colored according to their transmembrane potential at onset of a subthreshold stimulus pulse. Arrows indicate depolarized nodes of Ranvier. Diameters of neuronal processes were thickened for figure rendering.

The stimulus waveform produced by the Itrell II (Medtronic) implantable pulse generator (IPG), as used in human DBS and the monkey experiments, is a voltage-controlled biphasic asymmetric square pulse. However, the actual stimulus delivered to the brain tissue is modified by the electrode capacitance. To account for the effects of electrode capacitance, a Fourier FEM was used as described in Butson and McIntyre (2005). The four general steps of this method are briefly described here. First, the IPG stimulus waveform was constructed in the time domain. It was converted to frequency domain using discrete Fourier transform (DFT) in Matlab (Mathworks, Natick, MA). Third, the FEM model was solved at each component frequency of the DFT (1,024 frequencies between 0 and 50 kHz). The result at each frequency was scaled and phase shifted using the DFT magnitudes and phases. Finally, an inverse Fourier transform was performed to obtain the stimulus waveform in the time domain. The electrode was modeled as purely capacitive (0.65 μF) and adjusted to account for the smaller surface area of monkey electrode contacts compared with the human DBS electrode (Butson and McIntyre 2005).

Stimulation prediction

We coupled the finite element electric field model with the multicompartment neuron models to enable quantitative stimulation predictions in the context of the 3D neuroanatomy. This coupling was accomplished by applying extracellular voltages from the electric field model to each compartment of each neuron model and simulating the biophysical response (action potential signaling) of each neuron over time (McIntyre et al. 2004a) (Figs. 4 and 5). The magnitude of the applied extracellular voltage was dependent on the stimulus amplitude, stimulus waveform, and the compartment’s distance from the electrode. At each time step of the simulation, the extracellular voltage at each neural compartment was updated to a value determined by the time-dependent stimulus train delivered to the tissue medium (Figs. 4 and 5).

FIG. 5.

STN neuron firing in response to extracellular stimulation. Lowercase letters indicate location in the STN neuron where the transmembrane voltage was recorded. a, soma; b, 1st node of Ranvier; c, 30th node of Ranvier; d, 50th node of Ranvier. A: action potential initiated in the axon and propagated toward the cell body and axonal terminals in the globus pallidus. Traces in the top row represent stimulus voltage waveforms applied to the neuron. The 4 traces below show response to a subthreshold DBS pulse (1.4 V), a suprathreshold DBS pulse (1.8 V), and a suprathreshold DBS pulse (1.8 V) in a model with inhibitory somatic synapses. B: STN neuron spontaneous firing and response to stimulation are stable over time. The 4 traces show neuronal firing before, during, and after 1-s, 136-Hz DBS train in the soma and the distal axon of the same neuron. Results are displayed for models with and without GABAa stimulation induced synaptic input (0.013 μS). Somatic firing rate was lower than the stimulation frequency because action potentials initiated in the axon did not invade the soma for every stimulus pulse. GABAergic synaptic inputs reduce somatic firing, but axonal output was largely unaffected (see also Fig. 9).

Clinical efficacy for various DBS parameter settings was established in two parkinsonian macaques using behavioral tests (for details, see Hashimoto et al. 2003). In both monkeys, R7160 and R370, the electrode was positioned in the posterior STN at a 20° angle in the sagittal plane. In monkey R7160, bipolar stimulation (contact 0 cathode, contact 2 anode) at 136 Hz and 210-μs pulse width produced consistent improvement in rigidity and bradykinesia at 1.8-V amplitude (clinically effective stimulation), but not at 1.4 V (clinically ineffective stimulation). In monkey R370, bipolar stimulation (contact 2 cathode, contact 0 anode) at 136 Hz and 90-μs pulse width was clinically effective at 3-V amplitude and clinically ineffective at 2 V. Thus these various stimulation parameters were applied to our model system. The fundamental differences between the model simulations for the two monkeys were the electrode position and stimulation parameters. In addition, the two electrode positions resulted in somewhat different neural populations because different neurons were destroyed by the electrodes. The model neurons were stimulated with a train of 25 pulses. Longer train durations (1 s; 136 pulses) did not impact neural response to stimulation (Fig. 5B).

We did not observe any model neurons that exhibited a blocking of axonal firing from the direct application of the DBS electric field; therefore we concentrated our analysis on the excitatory response of the stimulation (Fig. 5A). Those that produced orthodromically propagating action potentials in response to >80% of the stimulus pulses were considered to be activated. This percentage was chosen because most neurons responded either to none or to 20 or more of the 25 stimulus pulses. Some activated neurons did not respond to all 25 pulses because, in certain cases, the stimulus pulse was delivered immediately after a spontaneous action potential (originating in the soma), while the axon was still in the refractory state. Percentages of activated neurons were averaged over three randomized populations. Student’s t-tests (1-tailed; P < 0.05) were performed to compare the averages for clinically ineffective and effective stimulation conditions.

Experimental recordings of GPi activity during STN DBS

The experimental recording procedure and data collection from the two monkeys used in this study has been described in detail elsewhere (Elder et al. 2005; Hashimoto et al. 2002, 2003). Briefly, single unit neural activity was recorded extracellularly from GPi-identified neurons using glass-coated platinum–iridium microelectrodes (impedances of 0.4–0.8 MΩ at 2 kHz). Recording penetrations were made in parasagittal planes moving rostral to caudal at an angle of 70° to the orbitomeatal line. Neurons (n = 12 for monkey R7160 and n = 27 for monkey R370) were recoded for 25–35 s during clinically effective STN DBS at 136 Hz. Stimulus artifact template subtraction methods (Hashimoto et al. 2002) and in-house software developed in MATLAB v7.0 (Mathworks) were used for the neural signal analysis. Peristimulus time histograms were constructed with a 0.2-ms bin size. The histograms used 336,373 spikes resulting from 644,929 stimuli in 39 GPi cells. Stimuli for which no spike was recorded in the interstimulus interval were not considered. Probability distributions for R7160 and R370 were compared using a χ2 goodness of fit test (P < 0.05). Cells were further classified as having an early response if >15% of the stimulus responses occurred at a latency of <1.5 ms. Significance of differences in the proportion of early response cells was found by χ2 test (P < 0.05).


We calculated levels of axonal activation for populations of STN projection neuron, GPi fibers of passage (lenticular fasciculus), and CST fibers of passage during clinically effective and ineffective STN DBS in two parkinsonian macaques. We evaluated four general aspects of our model system. First, activation of CST fibers at muscle contraction thresholds were analyzed to provide a gross degree of model validation. Second, we evaluated the neural response to therapeutic stimulation using three separate randomized populations of STN neurons and GPi fibers, and we correlated the model predictions to single-unit microelectrode recordings from the two monkeys during STN DBS. Third, the sensitivity of neural activation to electrode position was assessed by moving the electrode by 0.25 mm in four directions in the horizontal plane. Finally, the effects of DBS on STN somatic firing were evaluated in a model with stimulation-induced inhibitory trans-synaptic conductances.

Activation of the internal capsule with DBS

To address the experimental predictability of our model system, we evaluated the activation of CST fibers. Visually determined muscle contraction thresholds in the two monkeys, were 3 V for R7160 and 3.5 V for R370. These stimulation parameters resulted in activation of 11 ± 1 and 9 ± 3% (SD) of CST fibers in our R7160 and R370 models, respectively, averaged over three random populations. In addition, clinically effective stimulation parameters resulted in minimal activation of CST fibers: 0% and 5 ± 2% in the R7160 and R370 models, respectively. These results indicate that the model exhibited appreciable increases in CST fiber activation when comparing stimuli that experimentally resulted in no visible muscle contraction (clinically effective) and stimuli where muscle contractions were observed (CST threshold). This provides some evidence that the voltage spread in the tissue around the electrode was accurately predicted by the finite element model and that the neuron models fire at realistic activation thresholds.

Activation of the STN and LF with STN DBS

In the context of STN DBS, both STN projection neurons and GPi fibers of passage in the LF represent viable candidates as the therapeutic target of the stimulation. The GPi is an output nucleus of the basal ganglia and the STN modulates basal ganglia output through excitatory projections into the GPi. We defined DBS induced activation of STN projection neurons as the generation of an action potential anywhere in the neuron, resulting from the applied electric field, which propagated to the axonal terminal. When stimulating STN projection neurons with DBS, action potential initiation always took place in the myelinated axon, and if an action potential was induced by the stimulation, it always propagated to the axon terminal (Fig. 5A).

Our simulations predict activation of both GPi fibers and STN neurons during STN DBS (Fig. 6). These results were consistent over three randomized populations of neurons. Stimulation parameters that failed to improve parkinsonian symptoms in monkey R7160 (1.4 V, 210 μs, 136 Hz) activated 29 ± 2% of STN projection neurons, 9 ± 4% of GPi fibers of passage, and no CST fibers. Clinically effective stimulation (1.8 V) activated 37 ± 4% of STN projection neurons (30% average increase), 18 ± 6% of GPi fibers of passage (107% increase), and no CST fibers (Fig. 6A). In the second monkey, R370, clinically ineffective (2 V, 90 μs, 136 Hz) and clinically effective stimulation (3 V, 90 μs, 136 Hz) activated 31 ± 3% and 49 ± 5% of STN neurons (54% increase), 66 ± 2% and 82 ± 6% of GPi fibers (24% increase), and 1 ± 1% and 5 ± 2% of CST fibers, respectively (Fig. 6B). The increase in activation between clinically ineffective and effective stimulation was statistically significant for STN neurons in both monkeys and for GPi and CST fibers in monkey R370 (P < 0.05; Fig. 6).

FIG. 6.

Neural activation during clinically effective and ineffective DBS. A: STN neuron axons (top) and GPi fibers (bottom) activated by clinically effective DBS in monkey R7160 (1.8 V, 136 Hz, 210 μs, contact 0 cathode, contact 2 anode) are shown in red. Axons that did not respond to ≥80% of the stimulus pulses are shown in gray. Each leg of the 3D scale bar represents 1 mm. B: percent of neurons activated for monkeys R7160 (top) and R370 (bottom) during clinically ineffective, clinically effective, and corticospinal tract threshold DBS, averaged over 3 randomized populations. Only corticospinal tract (CST) fiber activation was evaluated at CST threshold amplitudes. Asterisks indicate significant difference between clinically ineffective and effective stimulation (P < 0.05; t-test).

The therapeutic effects of DBS depend critically on the stimulation frequency, with frequencies >100 Hz being generally beneficial and those <50 Hz sometimes worsening the symptoms (Rizzone et al. 2001). However, the output of our model system was not substantially affected by the stimulation frequency (i.e., approximately the same numbers of neurons were activated by the stimulus train). We did observe small decreases in the number of neurons activated as the frequency decreased from 136 to 2 Hz. For example, monkey R7160 exhibited an 8% reduction in STN neurons and a 2% reduction in GPi fibers activated when comparing 136- and 2-Hz stimulation. More importantly, all neurons activated by a given stimulus pulse fired in synchrony, and as the stimulation frequency changed, the activated neurons were entrained to fire at the given stimulation frequency. These results suggest that changes in stimulation frequency do not dramatically affect the volume of tissue activated, but rather that stimulation frequency exerts its influence on the basal ganglia network functioning (Grill et al. 2004; Montgomery and Baker 2000; Rubin and Terman 2004), which was not addressed in this model.

Given the limited quantitative characterization of STN neuron membrane dynamics, we performed a range of sensitivity analyses to addresses the robustness of our DBS model predictions. Before finalizing the variant of the Gillies and Willshaw (2006) membrane dynamics used in this study, we evaluated a wide range of models of STN neuron membrane dynamics (Otsuka et al. 2004; Wilson et al. 2004). The model used in this study best matched the available in vitro experimental characterization of STN neuron firing. However, nearly every STN neuron model variant that we evaluated generated nearly identical neural activation results in response to STN DBS. This consistency in the model output was primarily linked to the fact that action potential initiation in response to extracellular stimulation takes place in the myelinated axon (McIntyre and Grill 1999; McIntyre et al. 2004a; Nowak and Bullier 1998). In turn, the description of the somatic and dendritic membrane dynamics had little effect on the axonal output generated by DBS. For example, we modified the m (activation) gate of the spike triggering fast sodium channel in the dendrites, soma, and initial segment. Large changes in the time constant (±50%) and small shifts in the steady-state activation voltage (±1 mV) had no effect on the number of STN neurons activated by DBS (data not shown). However, it should be noted that modifications of the steady-state activation voltage strongly affected the spontaneous firing frequency of the STN neurons.

Experimentally recorded short-latency GPi activity during STN DBS

The response of GPi neurons to STN DBS was monitored experimentally by single-unit microelectrode recordings. Using stimulus artifact template subtraction we were able to construct peristimulus time histograms and analyze activity of GPi neurons immediately after each DBS pulse (Fig. 7) (Hashimoto et al. 2002, 2003). Experimentally, the activation of GPi fibers resulted in short latency excitation of GPi cell bodies, which could be interpreted as antidromic activation (Bar-Gad et al. 2004). Using the 3D anatomical brain atlas and GPi axon model, we estimated the propagation latency from DBS activation of the LF to antidromic activation of the GPi to be <1.5 ms. We calculated the number of experimentally recorded GPi cells that exhibited short-latency excitation for 15% or more of the DBS pulses, indicative of reliable antidromic activation (see discussion). In monkey R370, 41% of GPi cells fired within 1.5 ms after a DBS pulse, whereas only 8% of GPi cells fired within the same time period in monkey R7160 (Fig. 7). Consistent with these experimental results, our model predicted a large degree of GPi fiber activation for monkey R370 (82%), but a much smaller activation for monkey R7160 (18%).

FIG. 7.

Experimentally recorded GPi firing during STN DBS. Histograms of monkeys R7160 (A) and R370 (B) show the percentage of total GPi spikes (all spikes recorded from all cells in each animal) during clinically effective stimulation at 136 Hz separated into 0.2-ms bins over the interstimulus interval. Monkey R370 exhibited an increased number of spikes in the 1st 7 bins (<1.5 ms), indicating a higher degree of short latency excitation compared with monkey R7160.

Sensitivity of neural activation to electrode position

Electrode location can have a substantial impact on the therapeutic efficacy and side effects of DBS; however, conflicting opinions exist on the optimal position of STN DBS electrodes (Hamel et al. 2003; Nowinski et al. 2005; Saint-Cyr et al. 2002; Starr et al. 2002; Voges et al. 2002; Yelnik et al. 2003; Zonenshayn et al. 2004). Therefore we examined how small modifications in electrode position affected our simulation results. We compared neural activation for the original electrode position in monkey R7160 to models with the electrode moved either 0.25 mm medial, lateral, anterior, or posterior to the original position (Fig. 8). For clinically ineffective stimulation, the range of activation was 6–16% for GPi fibers, 0–2% for CST fibers, and 12–30% for STN neurons. For clinically effective stimulation, activation extended between 6 and 28% for GPi fibers, between 0 and 10% for CST fibers, and between 18 and 42% for STN neurons. As a result, electrode location within the STN region can substantially affect axonal activation of both STN neurons and GPi fibers.

FIG. 8.

Sensitivity of neural activation to electrode position. Percentage of STN neurons (top) and GPi fibers (bottom) activated during DBS for monkey R7160 with the electrode moved from the original position (O) to a location 0.25 mm medial (M), lateral (L), anterior (A), or posterior (P). The original electrode position was in the sagittal plane, 6.2 mm from the midline, at a 20° anterior-to-posterior angle. The center of the bottom contact (cathode) was 3.3 mm ventral to the AC-PC plane (horizontal plane defined by the anterior and posterior commissures) placing it in the posterior-medial-ventral border of the STN. The 2nd from the top contact (anode) was 1.4 mm below the AC-PC plane, corresponding to the area of LF just dorsal to posterior-lateral border of the STN.

STN neuron somatic firing during STN DBS

In the simulations described above, STN neuron activation was assessed at the distal end of the axon. Because experimental extracellular STN microelectrode recordings register primarily somatic firing rather than the axonal output, we also evaluated firing frequency in the STN soma during STN DBS (Figs. 5 and 9). During extracellular stimulation, action potential initiation occurs in the axon, and as a result, the soma and axon can fire independently (McIntyre and Grill 1999; McIntyre et al. 2004a; Nowak and Bullier 1998). In addition, extracellular stimulation can activate synaptic inputs impinging on local projection neurons (Anderson et al. 2006; Baldissera et al. 1972; Dostrovsky et al. 2000; Filali et al. 2004; Gustafsson and Jankowska 1976; Meissner et al. 2005; Welter et al. 2004). Synaptic inputs to the STN come primarily from the external part of globus pallidus (GPe) and cerebral cortex. Dominant GPe inhibitory afferents converge primarily on the cell body and proximal dendrites (Smith et al. 1990), which we modeled as GABAa trans-synaptic conductance in the central compartment of the soma. The presence of stimulation induced GABAa synaptic input inhibited somatic firing and thus reduced the frequency of spontaneous spikes (Figs. 9A and 5B). The magnitude of somatic firing suppression induced by stimulated GABAa input depended on the neuron position relative to the electrode and the strength of the inhibitory synaptic conductance. In the absence of extracellular stimulation or inhibitory synaptic input, the output of the STN neuron was dictated by the rate of spontaneous somatic spiking (32 Hz). If the STN neuron was too far from the electrode to be directly activated by the extracellular stimulus, but high-frequency inhibitory synaptic inputs were applied, the spontaneous somatic firing could be entirely suppressed, resulting in no axonal output. However, if the neuron was close enough to the electrode to be directly activated by the extracellular stimulus, axonal firing (i.e., neural output) was almost completely dictated by the stimulation frequency and largely unaffected by inhibitory synaptic currents applied to the soma (Fig. 9B). Interestingly, the inhibition of spontaneous somatic firing also increased the reliability of stimulus-evoked action potential generation in the axon from extracellular stimulation because of the elimination of the possible interaction between spontaneous spikes and stimulation induced spikes.

FIG. 9.

STN firing frequency under the influence of stimulation-induced trans-synaptic GABAa inhibitory inputs. Firing frequency was measured in the soma (A) and distal axon (B) of the same cell and averaged for a population of STN neurons (monkey R7160; clinically effective stimulation). Somatic firing, with respect to the no stimulation condition (No DBS), decreased with increasing strength of the GABAa synapse. Distal axon firing, with respect to the 0 conductance condition, was relatively unaffected by somatic inhibition.


The neural response to STN DBS that is responsible for the therapeutic effects of the stimulation has been unclear. This limitation in scientific knowledge has hindered our ability to understand and optimize this medical technology for current and future applications. Using an anatomically and electrically detailed computer model, we evaluated the neural activation generated by clinically effective and ineffective DBS in two parkinsonian macaques. Our simulation results showed axonal activation of both STN projection neurons and GPi fibers of passage during STN DBS. The relative proportion of neural types activated depended strongly on the position of the cathodic contact in the STN region. In monkey R7160, the cathode was in the ventral portion of the STN, which resulted in limited GPi fiber activation (∼10–20%). In monkey R370, the cathode was in the dorsal STN, on the border with LF, resulting in greater GPi fiber recruitment (∼65%), which also significantly increased during therapeutic stimulation conditions (∼80%). These theoretical predictions were supported by the large number of GPi neurons experimentally recorded with short latency excitation in monkey R370, indicative of LF antidromic activation, not seen in monkey R7160. Both monkeys had similar levels of STN neuron activation which showed significant increases with clinically effective stimulation (from ∼30% to 40–50%). These results indicate that activation of approximately one half of the STN is sufficient for the behavioral manifestation of the therapeutic effects of STN DBS. The additional recruitment of GPi fibers may also play an important role in therapeutic outcome, but large-scale activation of GPi fibers may not be necessary.

Model development, analysis, and limitations

The computer model developed in this study used a number of anatomical, physiological, and electrical improvements over previous attempts to theoretically address the cellular effects of DBS (McIntyre et al. 2004a). However, when constructing such a comprehensive model, it is necessary to make a number of assumptions and simplifications. The following section attempts to document these limitations and provide insight into their possible impact on our results.

We dedicated substantial effort toward the development of accurate neuron models for our simulations; however, they were unable to capture the full spectrum of experimentally defined characteristics. We used neuron biophysical properties based on rat STN neurons (Gillies and Willshaw 2006). Rat STN neurons have been extensively characterized with in vitro preparations, making it possible to construct faithful model representations with Hodgkin-Huxley type channel dynamics (Gillies and Willshaw 2006; Otsuka et al. 2004; Terman et al. 2002; Wilson et al. 2004). In addition, we attempted to approximate in vivo conditions by adjusting both the STN projection neurons and GPi fibers of passage to match the spontaneous activity observed in parkinsonian macaques (Wichmann et al. 2002). Nonetheless, our STN DBS model exhibited a simple on-off activation outcome (Fig. 5B), whereas the therapeutic effects of DBS typically evolve over seconds, minutes, and even hours of stimulation (Temperli et al. 2003). This discrepancy in the model predictions and clinical observations may be related to our inability to fully characterize the STN neuron membrane dynamics. Recent in vitro experimental studies have shown slow inactivation of sodium channels may play an important role in STN somatic activity during high-frequency stimulation (Beurrier et al. 2001; Do and Bean 2003). In addition, there exists a long list of factors influencing neural plasticity that could be affected by DBS (e.g., gene expression changes, channel density changes, synaptic strength changes) that were not included in our model system. However, our previous experience suggests that the specific details of the somatic ion channel membrane dynamics are of limited importance in the global effects predicted by computer models of extracellular stimulation (see sensitivity analyses in McIntyre and Grill 1999, 2000; McIntyre et al. 2004a). Possibly of greater importance is the use of realistic neural morphologies that directly interact with the electric field and determine polarization profiles of each neuron (Figs. 4B and 5A). Although our model included a full 3D reconstruction of a macaque STN projection neuron (Sato et al. 2000), as well as GPi and CST axonal trajectories based on documented morphologies (Parent et al. 2001), all neurons in each population were copies of a generic geometry. To compensate for this limitation, we added diversity in our STN projection neurons by creating two additional axonal trajectories based on histological axonal tracing experiments (Sato et al. 2000).

The results suggest that electrode position with respect to the surrounding anatomy is an important factor in the stimulation outcome (Figs. 6 and 8). Using histological brain slices, we developed a 3D basal ganglia reconstruction for one of the monkeys (R7160) and were able to precisely recreate the electrode position in the tissue. Unfortunately, we did not have access to R370’s brain, so the same anatomical model was used for both monkeys. However, the neuroanatomical differences are likely to be small because the two monkeys were both rhesus macaques of similar size and the same sex. Furthermore, the electrode for R370 was positioned using the descriptions and drawings from Hashimoto et al. (2003). It should also be noted that the 3D complex tissue anisotropy around the STN can affect the neural response to DBS (McIntyre et al. 2004b). However, in this study, we were unable to obtain diffusion tensor imaging data to estimate the 3D tissue conductivity properties specific to these animals, so an isotropic conductivity was used for the bulk tissue medium. To address some of these issues, a gross level model validation was attempted by comparing the experimentally defined CST thresholds in the two monkeys with the corresponding model predictions (Fig. 6). The models showed minimal CST activation during clinically effective stimulation (∼0–5% fibers activated), but more substantial activation at experimentally determined CST thresholds (∼10%). The exact percentage of CST fibers needed to be activated to generate a noticeable muscle contraction is unknown. However, we believe that the increase in CST activation shown by the model is substantial enough to justify the assumption that it would correspond to a perceptible effect experimentally. This is particularly true for monkey R7160 whose CST activation went from 0 to 11 ± 1%, and whose brain we explicitly reconstructed in 3D and were able to determine electrode location with a high degree of certainty.

We evaluated somatic firing during high-frequency stimulation (HFS) by including a highly simplistic representation of stimulation-induced trans-synaptic inhibition. STN microelectrode recording studies in monkey (Filali et al. 2004; Meissner et al. 2005) and human (Welter et al. 2004) have shown that STN HFS reduces the somatic firing rate of STN neurons, which are supported by our results (Fig. 9A). By selecting the appropriate synaptic conductance value, our model is in close agreement with Meissner et al. (2005), who documented an ∼50% decrease in somatic firing frequency with STN HFS. Our simulations also show that axonal output is largely unaffected by the inhibition of somatic firing (Fig. 9B). This disconnect between axonal and somatic firing has previously been noted in a model of DBS of a thalamocortical relay neuron (McIntyre et al. 2004a) and we now show it for a population of STN projection neurons.

It should also be noted that the STN is surrounded by many other fiber tracts that we ignored in this study. The reciprocal STN-GPe projections, STN-substantia nigra pars reticulata connections, and nigrostriatal fibers may also be directly affected by STN DBS (Lee et al. 2006). While it is possible that these tracts play a role in therapeutic mechanisms of DBS, we limited this study to evaluate the two most probable candidates (based on location of optimal therapeutic contacts in human patients), namely STN projection neurons and pallidothalamic (GPi) fibers. We also ignored possible physiologic effects that could result from antidromic activation of afferent inputs projecting to the STN (mostly from GPe and cerebral cortex), as well as the potential role of glial cells in the regulation of the extracellular environment (e.g., extracellular ionic concentrations and neurotransmitter levels). In turn, our model had a number of limitations, but we do not believe they impact our fundamental conclusions. In addition, the model system used in this study represents one of the most anatomically and electrically accurate computer models of DBS ever created and never before has such explicit connection between modeling and experimental DBS results been attempted.

Neural target of STN DBS

In current clinical practice, the STN is the target of choice for DBS treatment of PD, even though GPi DBS is similarly effective (Burchiel et al. 1999; Obeso et al. 2001; Rodriguez-Oroz et al. 2005). Several studies have shown that clinically effective STN DBS electrode contacts are located in the dorsal STN or the white matter above the STN formed by the LF and zona incerta (Hamel et al. 2003; Nowinski et al. 2005; Saint-Cyr et al. 2002; Starr et al. 2002; Voges et al. 2002; Yelnik et al. 2003; Zonenshayn et al. 2004). This has introduced the hypothesis that activation of either STN projection neurons and/or pallidothalamic fibers could be responsible for the therapeutic effects of STN DBS.

The results of this study suggest that both STN neurons and GPi fibers are activated during clinically effective STN DBS. Recent histological tracing studies have shown that the LF is composed of axons from the entire GPi (Parent and Parent 2004; Parent et al. 2001) and not just its medial (nonmotor) part as previously believed (Kim et al. 1976; Kuo and Carpenter 1973). Consequently, stimulating LF through contacts located in the dorsal STN would mechanistically be equivalent to direct GPi stimulation, at least when considering GPi output to thalamus. Because the GPi fibers of the LF are running in a tightly packed bundle, they can be effectively stimulated with lower charge injection than with electrodes in the GPi where the neurons are widely dispersed. However, our simulations show that, for monkey R7160, where the cathode was located in the ventral STN, relatively few GPi fibers were activated, yet the stimulation was clinically effective. Overall, our results suggest that large-scale activation of GPi fibers may not be necessary for therapeutic benefit, and a recent clinical study found that stimulating in the white matter above the STN was less effective than stimulation of the dorsolateral STN border (Herzog et al. 2004).

The large difference in model predicted GPi fiber activation in the two monkeys was supported by experimental single unit recordings performed during STN DBS (Hashimoto et al. 2003). The proportion of GPi cells exhibiting short latency, presumably antidromic activation, during STN DBS was similar to our theoretical predictions in both monkeys (R7160-18% vs. R370-82% in the model; R7160-8% vs. R370-41% in the experiments). The discrepancies observed between the model and experimental results could be related to a number of factors. Our analysis of the experimental data required that ≥15% of the DBS pulses must exhibit short latency activation in the GPi cell to consider it antidromically activated. This somewhat arbitrary threshold was necessary to perform our analysis for a number of reasons. First, no cell will respond in exactly the same manner for all stimulus pulses because of stochastic biological variability in the neural response and inherent problems with extracellular microelectrode recordings during DBS (artifact subtraction, signal fluctuations, unit isolation, spike sorting, etc.). Second, GPi neurons fire at relatively high frequencies with a mix of antidromic and orthodromic action potentials during STN DBS (Hashimoto et al. 2003). Orthodromic activity in the GPi cells will collide with antidromic signals generated by stimulation of the LF. Therefore antidromic excitation will not be recorded for every action potential evoked by the stimulation. In addition, our experimental sample of GPi cells was limited, and approximately one half of the GPi exit fibers would be expected to travel via the ansa lenticularis, a fiber tract rostral to the LF and outside the reach of the DBS voltage spread. Nonetheless, the experimental results support our model conclusions that therapeutic STN DBS generated significantly different activation of GPi fibers in the two monkeys.

The importance of activation of either STN projection neurons or GPi fibers of passage on the therapeutic outcome of STN DBS may be fundamentally dictated by the precise electrode location within the STN region. This concept is supported by our simulations where small changes (±0.25 mm) in electrode position could strongly affect GPi fiber and STN neuron activation (Fig. 8). Clinical studies that found LF/ZI stimulation to be effective could have placed DBS electrodes in locations more suitable for maximal GPi fiber activation. However, postoperative electrode localization techniques used in human studies are not accurate to a submillimeter level, making it difficult to determine electrode position with high certainty. Furthermore, it is possible that more than one optimal contact and set of stimulation parameters exists, but because of time-consuming programming methods currently used, it is not possible to definitively address every clinically beneficial setting. The results of this study lay the foundation for a prospective and coupled, theoretical and experimental, analysis of STN DBS where electrode location and stimulation parameter selection can be systematically manipulated to address the impact of STN neuron and/or GPi fiber activation on therapeutic outcome. Such future developments in DBS research will enable a more complete understanding of the optimal electrode location and give better insight into the neural targets of STN DBS.


This project was financially supported by National Institutes of Health Grants T32 GM-07250, R01 NS-47388, and R01 NS-37019, IRSC Grant MOP-5781, and the Ohio BRTT/WCI Partnership.


The authors thank Drs. Taka Hashimoto for experimental contributions, Andre Parent for providing the stained STN neuron histology, Douglas Bowden for help with the warping of atlas templates in Edgewarp, and Andrew Gilles for assistance with the STN neuron membrane dynamics.


  • The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “advertisement” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.


View Abstract