If you don't remember your password, you can reset it by entering your email address and clicking the Reset Password button. You will then receive an email that contains a secure link for resetting your password
If the address matches a valid account an email will be sent to __email__ with instructions for resetting your password
Krembil Brain Institute, University Health Network, Toronto, CanadaInstitute of Biomedical Engineering, University of Toronto, Toronto, CanadaKITE, Toronto Rehabilitation Institute, University Health Network, Toronto, CanadaCRANIA, University Health Network and University of Toronto, Toronto, Canada
Krembil Brain Institute, University Health Network, Toronto, CanadaKITE, Toronto Rehabilitation Institute, University Health Network, Toronto, CanadaCRANIA, University Health Network and University of Toronto, Toronto, CanadaDivision of Neurosurgery, Toronto Western Hospital, University Health Network, Toronto, CanadaDepartment of Surgery, University of Toronto, Toronto, Canada
Krembil Brain Institute, University Health Network, Toronto, CanadaCRANIA, University Health Network and University of Toronto, Toronto, CanadaDivision of Neurosurgery, Toronto Western Hospital, University Health Network, Toronto, CanadaDepartment of Surgery, University of Toronto, Toronto, Canada
Krembil Brain Institute, University Health Network, Toronto, CanadaCRANIA, University Health Network and University of Toronto, Toronto, CanadaDivision of Neurosurgery, Toronto Western Hospital, University Health Network, Toronto, CanadaDepartment of Surgery, University of Toronto, Toronto, Canada
Institute of Biomedical Engineering, University of Toronto, Toronto, CanadaKITE, Toronto Rehabilitation Institute, University Health Network, Toronto, CanadaCRANIA, University Health Network and University of Toronto, Toronto, Canada
CRANIA, University Health Network and University of Toronto, Toronto, CanadaDepartment of Surgery, University of Toronto, Toronto, CanadaDepartment of Physiology, University of Toronto, Toronto, Canada
Krembil Brain Institute, University Health Network, Toronto, CanadaInstitute of Biomedical Engineering, University of Toronto, Toronto, CanadaKITE, Toronto Rehabilitation Institute, University Health Network, Toronto, Canada
Extracellular stimulation is brain-region-specific and frequency-dependent.
•
Neuronal stimulus responses were excitatory in thalamus & inhibitory in basal ganglia.
•
High-frequency stimulation led to neuronal suppression/a loss of site-specificity.
•
Responses to single stimuli could be predicted based on anatomy/local microcircuitry.
•
Frequency-dependant neuronal suppression could be modelled by synaptic depression.
Abstract
Background
Deep brain stimulation is an established therapy for several neurological disorders; however, its effects on neuronal activity vary across brain regions and depend on stimulation settings. Understanding these variable responses can aid in the development of physiologically-informed stimulation paradigms in existing or prospective indications.
Objective
Provide experimental and computational insights into the brain-region-specific and frequency-dependent effects of extracellular stimulation on neuronal activity.
Methods
In patients with movement disorders, single-neuron recordings were acquired from the subthalamic nucleus, substantia nigra pars reticulata, ventral intermediate nucleus, or reticular thalamus during microstimulation across various frequencies (1–100 Hz) to assess single-pulse and frequency-response functions. Moreover, a biophysically-realistic computational framework was developed which generated postsynaptic responses under the assumption that electrical stimuli simultaneously activated all convergent presynaptic inputs to stimulation target neurons. The framework took into consideration the relative distributions of excitatory/inhibitory afferent inputs to model site-specific responses, which were in turn embedded within a model of short-term synaptic plasticity to account for stimulation frequency-dependence.
Results
We demonstrated microstimulation-evoked excitatory neuronal responses in thalamic structures (which have predominantly excitatory inputs) and inhibitory responses in basal ganglia structures (predominantly inhibitory inputs); however, higher stimulation frequencies led to a loss of site-specificity and convergence towards neuronal suppression. The model confirmed that site-specific responses could be simulated by accounting for local neuroanatomical/microcircuit properties, while suppression of neuronal activity during high-frequency stimulation was mediated by short-term synaptic depression.
Conclusions
Brain-region-specific and frequency-dependant neuronal responses could be simulated by considering neuroanatomical (local microcircuitry) and neurophysiological (short-term plasticity) properties.
], the effects of DBS on neuronal activity are not fully understood and neural responses evoked by electrical stimulation have been shown to differ across stimulation targets [
]. As such, the objective of this study was to use microelectrode recordings and stimulation to demonstrate how the neuronal effects of pulsatile stimulation vary depending on the stimulation target region and the frequency of stimulation pulses being delivered. Knowledge of the site-specific and frequency-dependent ability to selectively modulate (e.g. upregulate or downregulate) neuronal output is of importance for stimulation programming and the development of physiologically-informed stimulation paradigms in existing or prospective DBS indications; and can allow the user to leverage DBS in a functionally-specific manner.
It was previously demonstrated that single pulses of electrical stimulation to the substantia nigra pars reticulata (SNr) or globus pallidus internus (GPi) were associated with stimulation-evoked inhibitory responses, likely mediated by local GABA release [
]. High-frequency stimulation (HFS) of the thalamic ventral intermediate nucleus (Vim) on the other hand elicited brief short-latency excitatory responses, likely the result of unsustained local glutamate release [
]. In this study, microelectrode recordings of single-neuron activity across four brain regions (Vim, thalamic reticular nucleus (Rt), subthalamic nucleus (STN), and SNr) were assessed during microstimulation across a range of frequencies. We hypothesized that (i) the effects of individual electrical stimulation pulses would vary with respect to the distribution of afferent inputs converging on target neurons (whether predominantly inhibitory or excitatory), and that based on these local neuroanatomical properties, stimulation pulses would elicit either net inhibitory or excitatory responses. Moreover, based on previous findings of HFS-induced depression of stimulus-evoked field potentials [
], we hypothesized that (ii) suppression of neuronal activity during HFS is mediated by changes in short-term synaptic dynamics (i.e. depression of inhibitory and excitatory synaptic transmission).
In addition to experimental procedures, we developed a computational framework for modelling the site-specific and frequency-dependent neuronal responses to electrical stimulation based on the above hypotheses. Previous theoretical works suggest that individual pulses of extracellular stimulation (i.e. DBS) initiate action potentials which are propagated along axons and/or their terminals [
Cellular, molecular, and clinical mechanisms of action of deep brain stimulation—a systematic review on established indications and outlook on future developments.
]. These axonal activations can in turn mediate synaptic transmission. Based on our first hypothesis, our computational framework considers that the postsynaptic neuronal responses to individual DBS pulses are the result of a simultaneous activation of presynaptic inputs and takes into consideration the site-specific proportions of inhibitory and excitatory inputs converging on target neurons (derived from anatomical literature). However, HFS may reduce synaptic transmission fidelity by way of synaptic depression [
] of short-term synaptic plasticity was embedded within our computational framework in order to account for changes to synaptic transmission fidelity based on the frequency of successive stimuli.
Methods
Experimental: Patients and neurons
115 neurons from patients with Parkinson's disease (n = 47) or essential tremor (n = 11) were included in this study. All experiments conformed to the guidelines set by the Tri-Council Policy on Ethical Conduct for Research Involving Humans and were approved by the University Health Network Research Ethics Board. Moreover, each patient provided written informed consent prior to taking part in the studies.
Experimental: Protocols
Neurophysiological mapping procedures were performed during awake DBS surgeries (OFF-medication) using two closely spaced microelectrodes (600 μm apart, 0.1–0.4 MΩ impedances) [
] neurons have been previously reported. One microelectrode was used for recording single-neuron activity while a second immediately adjacent microelectrode was used to deliver stimulation at different frequencies. Recordings were obtained using two Guideline System GS3000 amplifiers (Axon Instruments, Union City, CA) and signals were digitized at ≥12.5 kHz with a CED 1401 data acquisition system (Cambridge Electronic Design, Cambridge, UK). Microstimulation was delivered using an isolated constant-current stimulator (Neuro-Amp1A, Axon Instruments, Union City, CA) with 0.3 ms biphasic pulses (cathodal followed by anodal).
To generate stimulation frequency response functions, stimulation trains were delivered at 1 Hz (10 pulses), 2 Hz (20 pulses), 3 Hz (60 pulses), 5 Hz (50 pulses), 10 Hz (50 pulses), 20 Hz (60 pulses), 30 Hz (60 pulses), 50 Hz (50 pulses), and 100 Hz (50 pulses) using 100 μA and a 0.3 ms biphasic pulse width. This frequency response protocol was executed at 9 Vim (npatients = 5), 11 Rt (npatients = 11), 27 STN (npatients = 16), and 14 SNr (npatients = 9) recording sites. Data for STN and SNr were previously collected [
], whereas Vim and Rt data for this study were unique. Longer trains (>2 s) of 100Hz stimulation were also delivered to the aforementioned Vim, Rt, and SNr neurons. 100 Hz long train data for STN (44 neurons, npatients = 20) were previously collected [
]. Please refer to Supplementary Table 6 for data summary.
Experimental: Offline analyses and statistics
For artifact removal, data from the start of each stimulation artifact to just after the anodic peak (i.e. from the anodic peak or last saturated value to about 25% of the baseline amplitude) were replaced by a straight line; corresponding to a time window of ∼0.8 ms. Data were then high pass filtered (≥250 Hz) and template matching was done using a principal component analysis method in Spike2 (Cambridge Electronic Design, UK). Artifact subtraction allowed for data to be high-pass filtered without distortion in the time domain as would otherwise occur when filtering a signal containing saturated high-amplitude stimulation artifacts [
Complex locking rather than complete cessation of neuronal activity in the globus pallidus of a 1-methyl-4-phenyl-1,2,3,6-tetrahydropyridine-treated primate in response to pallidal microstimulation.
]. As a single action potential is ≥ 1 ms, then at most one action potential might be lost in the <1 ms artifact subtraction process. With a 0.8 ms artifact removal window, the percentage of data lost during each stimulation train corresponds to: 0.08% (1Hz), 0.16% (2Hz); 0.24% (3Hz); 0.4% (5Hz); 0.8% (10Hz); 1.6% (20Hz); 2.4% (30Hz); 4% (50Hz); 8% (100Hz). To investigate single-pulses responses, peristimulus histograms (120 ms total width, 20 ms offset, 2 ms bins) were created to encompass responses to all 50 stimuli delivered during the 5Hz train, across all neurons. The 20 ms pre-stimulus periods were compared to the 20 ms and 40 ms post-stimulus periods using Bonferroni-corrected (two comparisons) two-tailed paired t-tests, and effect sizes (Cohen's dz) were calculated. For the frequency response protocol (≤60 stimulation pulses delivered at each frequency), firing rates were measured before and during each of the stimulation trains. Kolmogorov-Smirnov tests were used to assess the null hypothesis that the data are normally distributed. One-way repeated measures ANOVA tests (stimulation frequency as a within-subject factor) were carried out, and if significant main effects were found, Bonferroni-corrected (nine comparisons) post-hoc t-tests were used to compare firing rates during the various stimulation trains to pre-stimulation baseline firing. ANOVA effect sizes (η2) and t-test effect sizes (Cohen's dz) were also determined. One neuron from the Vim group and one neuron from the Rt group were excluded from statistical analyses due to incomplete stimulation protocol (i.e. missing data points). Of note, the solid gray lines in Fig. 1B consider that each stimulation pulse generated one action potential on the efferent axon [
], representing a situation in which the overall “neuronal output” is the summation of the somatic firing rate and a stimulus-locked efferent axon activation. However, readers should note that the statistical analyses only consider the action potential firing during periods of time that were not populated by artifacts (i.e. the activity generated at the somatic level). ANOVA analyses were carried out in the same way for both experimental (Fig. 1B) and computational (Fig. 8) results. To investigate possible time-varying responses throughout the stimulation trains, time-series histograms (2–3 s total width, no offset, 50 ms bins) were created for 5 Hz, 10 Hz, 20 Hz, 30 Hz, and long trains of 100 Hz (and 200 Hz long trains for Vim). Of note, the long train (i.e. 3 s) 100 Hz (and 200 Hz) data come from various sources since long trains of high-frequency stimulation were not initially delivered (please refer to Supplementary Table 6 for data summary). The attenuations of excitation over time in Vim and Rt during stimulation trains of ≥20 Hz were fit with double exponential functions. Histograms were also created for the shorter trains (≤1 s) of stimulation at 50 Hz and 100 Hz (0.5–1 s total width, no offset, 20 ms bins; Supplementary Fig. 1). Moreover, to investigate the prominent time-vary effects in Vim and Rt during 3 s, 100 Hz and 200 Hz stimulation trains, baseline firing was compared to the first second of stimulation and the subsequent 2 s of stimulation using Bonferroni-corrected (two comparisons) two-tailed paired t-tests.
Fig. 1Experimental (A) peristimulus responses and (B) frequency response functions. (A) Top panels show an exemplary response to a single stimulation pulse in each structure, whereas bottom panels show groupwise firing rate (mean + standard error) peristimulus time histograms of stimulus-evoked excitatory responses for Vim (n = 9) and Rt (n = 11) and stimulus-evoked inhibitory responses for STN (n = 27) and SNr (n = 14). The average firing rates of the immediate 20 ms and 40 ms periods following stimulations pulses were significantly greater than the 20 ms pre-stimulus periods for Vim and Rt, and significantly lesser for STN and SNr (p-values of Bonferonni-corrected 2-tailed paired t-test are displayed with Cohen's dz effect sizes in parentheses). (B) Stimulation (≤60 pulses) frequency response functions show that average firing rates progressively increased in Vim and Rt as the stimulation frequency became greater, while they progressively decreased in STN and SNr. The average ± standard error baseline firing rates for Vim, Rt, STN, and SNr neurons were 32.0 ± 11 Hz, 8.2 ± 1 Hz, 39.9 ± 3 Hz, and 102.3 ± 16 Hz, respectively (dashed gray lines). Firing rates during the various stimulation trains were compared to the baseline firing rates and the p-values of Bonferroni-corrected post-hoc t-tests (2-tailed, paired) are displayed with Cohen's dz effect sizes in parentheses. ANOVA main effects for stimulation were all significant and are reported in the Results subsection “Experimental: Stimulation frequency response functions”. If one considers that each DBS pulse generates an action potential on the efferent axon, then the overall neuronal output would be the summation of the somatic firing rate and stimulation frequency; this is represented by the solid gray lines in each plot (the values on this line for 100 Hz in Vim, Rt, and STN are 100 (Hz) plus the value on the corresponding coloured line). The right anatomical panels are 12.0 mm and 14.5 mm sagittal sections (Supplementary Fig. 3 shows the locations of the highlighted structures relative to other neuroanatomical landmarks).
To model the effect of DBS pulses on the afferents of the stimulated nuclei, we used a leaky integrate and fire (LIF) single neuron model, together with a TM model of short-term synaptic plasticity [
]. Each model neuron received 500 presynaptic inputs and the proportion of excitatory and inhibitory inputs were obtained using morphological data (detailed below in “Computational: Presynaptic inputs”). In addition to these inputs, the background synaptic activity [
]. In accordance with our first hypothesis, each DBS single pulse simultaneously activated all presynaptic inputs (Fig. 3A). This simultaneous activation was modelled by artificially generating precise spike times which correspond to the arrival of each DBS pulse in the presynaptic inputs. We utilized our modeling framework to recreate the neuronal firing in Vim, STN, and SNr in response to stimulation trains with frequencies of 1, 2, 5, 10, 20, 30, 50, and 100 Hz. Model generation for Rt neurons was omitted to avoid redundancy since the model parameters are identical to Vim except for the parameters which underly the baseline firing rates (this is elaborated upon in detail within the “Computational: Parameter settings” subsection below).
Computational: Presynaptic inputs
The vast majority of inputs to the Vim are glutamatergic projections from the dentate nucleus of the cerebellum [
Distribution and properties of GABAB antagonist [3H]CGP 62349 binding in the rhesus monkey thalamus and basal ganglia and the influence of lesions in the reticular thalamic nucleus.
Organization of projections from the anterior pole of the nucleus reticularis thalami (NRT) to subdivisions of the motor thalamus: light and electron microscopic studies in the Rhesus monkey.
]; like Vim, the majority of afferent inputs to Rt are glutamatergic. The vast majority (∼90%) of inputs to the SNr are GABAergic, projecting from the striatum and globus pallidus externus (GPe) [
]. While the mixed inputs are more homogenous in STN, electron-microscopy work suggests that GABAergic terminals nevertheless outnumber glutamatergic terminals [
]. Based on the cited literature, estimates of the proportions of inhibitory and excitatory inputs were generated (Supplementary Table 1) to be used for the model.
In the model, an ensemble of 500 LIF model neurons produced inputs to the stimulated nuclei. Each neuron received a random input (modelled by OU process of time constant 5 ms) and fired at the rate of about 5 Hz (the total average firing rate across neurons was equal to 5 ± 0.7 Hz). Each of the 500 neurons was labeled either as excitatory or inhibitory based upon estimates of the proportions of excitatory and inhibitory inputs received by Vim, STN, and SNr (Supplementary Table 1); and their spikes were fed to the stimulated nuclei through the TM model. We used an LIF neuron model (see Supplementary Tables 2 and 5 for the LIF parameters) to generate membrane potentials of the stimulated nuclei. The total synaptic current was obtained as a linear combination of presynaptic excitatory (Iexc) and inhibitory currents (Iinh):
Isyn(t) = wexc Iexc(t) + winh Iinh(t)
(1)
where wexc and winh denote the weights of excitatory and inhibitory currents, respectively. These weights, together with the mean and standard deviation of the background synaptic current, were tuned to reproduce the neuronal firing rates at the baseline (DBS-OFF) as well as in response to DBS with different frequencies (Supplementary Table 2).
Computational: Synapse model
We utilized the TM equations to model the function of short-term synaptic plasticity:
(2)
(3)
(4)
where u indicates utilization probability, i.e., the probability of releasing neurotransmitters in synaptic cleft due to presynaptic influx of calcium ions. Upon the arrival of each presynaptic spike, tsp, u increases by and then decays to zero by the facilitation time constant, . U denotes the increment of u produced by each presynaptic spike. A denotes the absolute synaptic efficacy of the synaptic connections. The vesicle depletion process – due to the release of neurotransmitters – was modelled by (2) where r denotes the fraction of available resources after neurotransmitter depletion. In contrast to the increase of u upon the arrival of each presynaptic spike, r drops and then recovers by depression time constant to its steady state value of 1. The competition between the depression ( and facilitation ( time constants determines the dynamics of the synapse. In the TM model, , and are three parameters that determine the type and dynamics of the synapse. In (4), I and indicate the presynaptic current and its time constant, respectively. The time constants of the excitatory and inhibitory inputs were 3 ms and 10 ms, respectively.
], we used OU process of the time constant of 5 ms to represent the effect of synaptic noise. The OU process can be written as:
(5)
where ξ is a random number drawn from a Gaussian distribution with 0 average and unit variance. is the time constant, μ and α indicate the mean and standard deviation of variable x, respectively.
Computational: Neuron model
The membrane potential dynamics in an LIF model can be written as:
(6)
where EL = −70 mV, R = 1 MΩ, and τV = 10 ms. Iinj indicates the total injected current to the model neuron (i.e., Isyn plus background synaptic noise (5)). A spike occurspike occurs when V ≥ Vth, where Vth = - 40 mV and the reset voltage is −90 mV with an absolute refractory period of 1 ms.
Computational: Parameter setting
The proportions of excitatory and inhibitory neurons (Supplementary Table 1), total synaptic current (Supplementary Table 2), parameters of excitatory (Supplementary Table 3) and inhibitory (Supplementary Table 4) synapses, and time constants of membrane dynamics and synaptic currents (Supplementary Table 5) are available in the Supplementary Material. Of note, values were derived from previous experimental work for Supplementary Tables 3 and 4 [
]; parameter setting for Supplementary Tables 1 and 2 are described above. Also, as previously mentioned, modelling of Rt was omitted due to redundancy as all parameters are identical to Vim except for the parameters which mediate the baseline firing rates (i.e. wexc and winh and parameters of background synaptic noise; Supplementary Table 2). Parameters for Rt modelling are nevertheless included within the Supplementary Tables.
Experimental: Responses to single stimulation pulses & stimulation frequency response functions
The average ± standard error baseline firing rates for Vim, Rt, STN, and SNr neurons were 32.0 ± 11Hz, 8.2 ± 1Hz, 39.9 ± 3Hz, and 102.3 ± 16Hz, respectively. The responses to single stimulation pulses (Fig. 1A) showed stimulus-evoked excitatory responses for Vim and Rt, and inhibitory responses for STN and SNr. For Vim, the average firing rates of the immediate 20 ms (181.0Hz ± 33Hz; p = .002) and 40 ms (125.2 ± 25 Hz; p = .003) periods following stimulation pulses were significantly greater than the 20 ms pre-stimulus period. This was also the case for the 20 ms (186.2 ± 29Hz; p < .001) and 40 ms (120.8 ± 20Hz; p <. 001) post-stimulus periods for Rt. For STN, the average firing rates of the 20 ms (22.4 ± 3Hz; p <. 001) and 40 ms (32.6 ± 4Hz; p = .041) post-stimulus periods were significantly less than the 20 ms pre-stimulus period. This was also the case for the 20 ms (7.8 ± 3Hz; p < .001) and 40 ms (21.9 ± 7Hz; p = .003) post-stimulus periods for SNr. All statistics were corrected for multiple comparisons. Cohen's dz effect sizes are depicted in Fig. 1A.
The stimulation frequency response functions (Fig. 1B) show excitatory responses for Vim and Rt, and inhibitory responses for STN and SNr. For Vim, neuronal firing rates progressively increased as the stimulation frequency became greater and a significant main effect of stimulation was found [F = 43.074 (9, 234), p < .001, η2 = 0.624]. Bonferroni-corrected t-tests revealed differences in neuronal firing compared to baseline at stimulation frequencies of 10 Hz (p = .038), 30 Hz (p = .041), and greater (p < .05). For Rt, neuronal firing rates also progressively increased as the stimulation frequency became greater and a significant main effect of stimulation was found [F = 31.170 (9, 117), p < .001, η2 = 0.706]. Statistically significant differences in neuronal firing compared to baseline were found at stimulation frequencies of 30Hz (p = .029) and greater (p < .05). For STN, neuronal firing rates were progressively attenuated as the stimulation frequency became greater and a significant main effect of stimulation was found [F = 26.420 (9, 91), p < .001, η2 = 0.746]. Statistically significant differences in neuronal firing compared to baseline were found at stimulation frequencies of 20 Hz (p = .029) and greater (p < .001). For SNr, neuronal firing rates also progressively attenuated as the stimulation frequency became greater and a significant main effect of stimulation was found [F = 25.890 (9, 63), p < .001, η2 = 0.787]. Statistically significant differences in neuronal firing compared to baseline were found at stimulation frequencies of 3 Hz (p < .05) and greater (p ≤ .01). Detailed post-hoc t-test statistics (all corrected for multiple comparisons within the text and figures) and Cohen's dz effect sizes are depicted in Fig. 1B.
Experimental: Time-domain responses to stimulation
In Vim and Rt, periodic excitatory responses were evident at 5 Hz and 10 Hz (Fig. 2). The strength of the excitatory responses attenuated over time during stimulation trains of ≥20 Hz and were modelled by double exponential decay functions (R2 values within Fig. 2). The time-series histograms for long train ≥100 Hz data (3 s) in Vim and Rt show particularly prominent time-varying responses. These stimulations elicited excitatory responses that were transient in nature and limited to start of stimulation. For Vim at 100Hz (3 s), the firing rate at baseline (41.1 ± 7 Hz) was different from the firing rate during the first 1 s of stimulation (94.8 ± 7 Hz; p = .004), but not for the subsequent 2 s of stimulation (45.4 ± 7 Hz). For Vim at 200 Hz (3 s), the firing rate at baseline (53.1 ± 8 Hz) was not different from the first 1 s of stimulation (35.9 ± 9 Hz; Fig. 2 depicts a very transient initial excitation followed by suppression) but was for the subsequent 2 s (14.0 ± 5 Hz; p = .002). For Rt at 100Hz (3 s), the firing rate at baseline (7.7 ± 1 Hz) was different from the first 1 s of stimulation (90.7 ± 14 Hz; p = .002), but not for the subsequent 2 s (4.3 ± 2 Hz). In SNr, periodic inhibitory responses were evident at 5 Hz and 10 Hz. In STN and SNr, there was an overall stationary neuronal suppressive effect with increasing frequency (rather than an effect which changed dynamically over time as was the case in Vim and Rt). Of note, the data in Fig. 2 for 5–30 Hz stimulation is the same as that presented in Fig. 1B, while the 100 Hz (and 200 Hz) data in Fig. 2 come from various sources since long trains of high-frequency stimulation were not initially delivered (please refer to Supplementary Table 6 for data summary).
Fig. 2Experimental time-domain responses to stimulation trains. For Vim and Rt, firing rates (mean +standard error) progressively increased with increasing stimulation frequencies. Periodic excitatory responses are shown at 5 Hz and 10 Hz, however neuronal excitation declined over time with ≥20 Hz. Excitatory responses with 100 Hz long trains (≥2s) were transient, and a subsequent reduction of neuronal firing is evident after the initial excitation. In Vim, the initial excitatory response at 200 Hz is of shorter duration than at 100 Hz, and the subsequent neuronal suppressive response is stronger at 200 Hz compared to 100 Hz. In STN and SNr, firing rates progressively decreased with increasing stimulation frequencies. In SNr, periodic inhibitory responses are visible at 5 and 10 Hz. Exemplary firing rate raster data from each structure during the various stimulation trains are displayed above each of the panels. This figure is intended to demonstrate the dynamics of the firing rate as a function of time. Of note, the 100 Hz (and 200 Hz) data herein are different than the data presented in Fig. 1B (please refer to Supplementary Table 6 for data summary).
Fig. 3Computational (A) model framework and (B) simulated peristimulus responses. (A) To model the response to single pulses of electrical stimulation, each model neuron was assigned a certain proportion of excitatory and inhibitory presynaptic inputs/weights with proportions derived from anatomical literature. The effect of each DBS single pulse was modelled by simultaneously activating all presynaptic inputs. (B) The corresponding changes to (i) synaptic currents and (ii) somatic firing induced by the simultaneous activations are displayed (i.e. the single-pulse responses). This framework closely replicated the robust stimulus-evoked neuronal excitation in Vim and neuronal inhibition in SNr. In STN, there was a short-latency neuronal excitation which was not observed in the experimental data (though may have been occluded by the stimulation artifact) due to the high speed of excitatory synaptic transmission, followed by an inhibitory period congruent with the experimental data. F: facilitatory; D: depressive; P: pseudolinear; indicating the different types of synapses.
Computational: Responses to single stimulation pulses
The net changes to postsynaptic currents in response to single pulses of stimulation were modelled by simultaneous activations of all presynaptic inputs (Fig. 3Bi). These responses differed across brain regions due to differences in the proportions of excitatory and inhibitory inputs (summarized in the Methods subsection “Computational: Presynaptic inputs” and Supplementary Table 1). The simulated peristimulus firing rate histograms (i.e. the neuronal responses to the aforementioned changes to presynaptic currents) revealed stimulus-evoked excitatory responses for Vim (peak firing rate of 405.9 Hz vs. 245.1 Hz in the experimental data), inhibitory responses for SNr (minimum firing rate of 0 Hz vs. 0.7 Hz in the experimental data), and a short-latency excitatory responses (78.8 Hz peak vs. no peak in the experimental data) followed by a longer latency inhibitory response (8.4 Hz trough vs. 4.8 Hz in the experimental data) for STN (Fig. 3Bii). The lack of short-latency excitation in the experimental data for STN might be explained by discrepancies in the temporal dynamics of excitatory transmission and/or occlusion of the excitatory response by the stimulus artifact.
Computational: Time-domain synaptic currents
Excitatory and inhibitory synaptic currents were generated separately, along with the net (i.e. sum of excitatory and inhibitory) synaptic currents in responses to DBS pulses across a range of frequencies for each of Vim, STN, and SNr (Fig. 4). The TM model accounts for frequency-dependent changes to short-term synaptic dynamics. In all structures, the model suggests frequency-dependent depression of both excitatory and inhibitory synaptic currents. For Vim, sustained periodic excitations were seen with 5 Hz and 10 Hz, while frequency-dependent weakening of the excitatory responses with successive stimuli were observed with frequencies ≥20 Hz. Predominant inhibitory synaptic currents corroborate the strong inhibitions of somatic firing in SNr with low stimulation frequencies; whereas neuronal suppression with higher frequencies was likely the result of frequency-dependent synaptic depression. For STN, the mixed excitatory-inhibitory stimulus-evoked responses likely explain the more net-neutral somatic firing responses in experimental data with lower stimulation frequencies; whereas synaptic depression can explain the frequency-dependent suppression of somatic firing with higher stimulation frequencies.
Fig. 4Computational time-domain synaptic currents. The three figures show excitatory, inhibitory, and total (i.e. sum of excitatory and inhibitory) synaptic currents in responses to DBS pulses across a range of frequencies with an embedded TM model to account for short-term synaptic dynamics. In all cases, the model suggests frequency-dependent depression of both excitatory and inhibitory synaptic currents. For Vim, rather stable periodic excitations are seen with 5 Hz and 10 Hz. Also corroborating experimental data, frequency-dependent weakening of the excitatory responses is observed with frequencies ≥20 Hz. Predominant inhibitory synaptic currents corroborate the strong inhibitions of somatic firing in SNr, together with frequency-dependent synaptic depression. For STN, the mixed excitatory-inhibitory stimulus-evoked responses likely explain the more net-neutral somatic firing responses in experimental data with lower stimulation frequencies, while synaptic depression can explain frequency-dependent suppression of somatic firing.
The membrane potentials of modelled neurons in response to DBS across a range of frequencies were generated for each of Vim (Fig. 5), STN (Fig. 6), and SNr (Fig. 7). The proportions of excitatory and inhibitory inputs (Supplementary Table 1) together with the parameters of the model neurons (Supplementary Table 2) generated baseline (DBS-OFF) firing rates which corresponded to in vivo recordings. The left side of Fig. 5 shows the simulated membrane potential (accounting also for action potential generation) before and during stimulation across a range of frequencies for Vim, whereas the right side shows an exemplary in vivo Vim neuron. Synchronous/periodic neuronal firing due to stimulus entrainment was reproduced by the model neuron for DBS at 20 Hz. The model neuron can moreover partially reproduce the transient excitatory responses at DBS onset with 50 Hz and 100 Hz stimulation and 30 Hz to some degree; however, the transient excitatory responses within the model were of shorter latency. For STN (Fig. 6), the simulated (left) neuronal firing compared to baseline decreases for DBS at ≥30 Hz, corroborating experimental data (exemplary in vivo STN neuron portrayed on the right side). Neuronal firing rates were substantially attenuated with DBS at 50 Hz and 100 Hz (as was the case experimentally) due to synaptic depression. For SNr (Fig. 7), the simulated (left) neuronal firing rate decreases dramatically beginning at 20 Hz due to the dominant inhibitory presynaptic currents, corroborating experimental data (exemplary in vivo SNr neuron portrayed on the right side). The model neuron fails to generate action potentials for DBS ≥50 Hz (as was the case experimentally) due to synaptic depression. Time-domain histograms are also presented in each figure (Fig. 5, Fig. 6, Fig. 7) which were generated by averaging the neuronal firing rates of 10 modelled neurons for each respective structure across 2 s of stimulation at each frequency.
Fig. 5Computational time-domain membrane potential for Vim. The left panels show the membrane potential of a model Vim neuron immediately before (non-shaded) and during (shaded) DBS across a range of frequencies. The right panels are exemplary recordings from an in vivo human Vim neuron (stimulation for 50 Hz was limited to 1 s). The bottom-most panels are time-domain firing rate histograms generated by averaging across 10 model Vim neurons. Synchronous/periodic neuronal firing due to stimulus entrainment was reproduced by the model neuron for DBS at 20 Hz. The model neuron can partially generate the transient excitatory responses at DBS onset with 50 Hz and 100 Hz stimulation and 30 Hz to some degree; however, the transient excitatory responses within the model were of shorter latency.
Fig. 6Computational time-domain membrane potential for STN. The left panels show the membrane potential of a model STN neuron immediately before (non-shaded) and during (shaded) DBS across a range of frequencies. The right panels are exemplary recordings from an in vivo human STN neuron (stimulation for 50 Hz was limited to 1 s). The bottom-most panels are time-domain firing rate histograms generated by averaging across 10 model STN neurons. The neuronal firing rate compared to baseline decreases for DBS at ≥30 Hz, corroborating experimental data. The modelled neuronal firing rates were substantially attenuated with 50 Hz and 100 Hz (as was the case experimentally) due to synaptic depression.
Fig. 7Computational time-domain membrane potential for SNr. The left panels show the membrane potential of a model SNr neuron immediately before (non-shaded) and during (shaded) DBS across a range of frequencies. The right panels are exemplary recordings from an in vivo human SNr neuron (stimulation for 50 Hz was limited to 1 s). The bottom-most panels are time-domain firing rate histograms generated by averaging across 10 model SNr neurons. Due to the dominant inhibitory presynaptic currents, the neuronal firing rate decreases dramatically beginning at 20 Hz, corroborating experimental data. The model neuron fails to generate action potentials for DBS ≥50 Hz (as was the case experimentally) due to synaptic depression.
Fig. 8Computational frequency response functions. The neuronal dynamics for stimulation (≤60 pulses) frequency response functions match experimental results (solid gray lines), though further tuning is required to optimize excitatory responses of Vim more precisely. The average ± standard error baseline firing rates for computational Vim, STN, and SNr neurons were 28.0 ± 0.1 Hz, 30.1 ± 0.2 Hz, and 61.7 ± 0.3 Hz, respectively (dashed gray lines). Firing rates during the various stimulation trains were compared to the baseline firing rates and the p-values of Bonferroni-corrected post-hoc t-tests (2-tailed, paired) are displayed with Cohen's dz effect sizes in parentheses. ANOVA main effects for stimulation were all significant and are reported in the Results subsection “Computational: Stimulation frequency response functions”.
Computational: Stimulation frequency response functions
Similar to experimental results, significant main effects of stimulation were found for Vim [F = 2400.280 (6, 54), p < .001, η2 = 0.996], STN [F = 227.963 (6, 54), p < .001, η2 = 0.962], and SNr [F = 7093.439 (6, 54), p < .001, η2 = 0.999]. 10 modelled neurons were used for each brain structure and the stimulation duration at each frequency was constrained to match experimental data within Fig. 1B. Detailed post-hoc t-test statistics (corrected for multiple comparisons) and Cohen's dz effect sizes are depicted in Fig. 8. Overall, the neuronal dynamics matched experimental results, though further tuning is required to optimize initial excitatory responses of Vim more precisely; a topic for future work.
Discussion
Site-specific and frequency-dependent stimulation effects
At the somatic level, electrical stimulation is both site-specific and frequency-dependent. In Vim and Rt, neuronal activity could be upregulated, whereas in STN and SNr it was downregulated. These mechanistic disparities across brain regions are most likely explained by anatomical differences in local microcircuitries, in that the effects appeared dependent upon the relative distributions of excitatory and inhibitory inputs converging at target neurons [
]. The experimental findings demonstrated that neuronal activity in any brain region could be suppressed either selectively in regions with a high predominance of inhibitory inputs or non-selectively if high enough stimulation frequencies were used. Neuronal excitation, however, could only be achieved when electrical stimulation was delivered to brain regions with a high predominance of glutamatergic inputs. While these bimodal effects (excitatory vs. inhibitory) with low stimulation frequencies were likely attributable to presynaptic activation, the loss of site-specificity and convergence towards neuronal suppression with sustained HFS (≥100 Hz) was most likely attributable to synaptic depression [
]. This phenomenon of short-term synaptic plasticity can be defined as a reversible decrease in synaptic efficacy, caused by the depletion of readily releasable neurotransmitter vesicle pools when successive stimuli are delivered at a fast rate; a reduction of presynaptic calcium conductance; and/or the inactivation of neurotransmitter release sites due to delayed recovery from vesicle fusion events [
Our computational model was designed to test our two main hypotheses (i) that the post-synaptic responses (i.e. neuronal output), to single pulses of electrical stimulation were mediated by the proportions of inhibitory vs. excitatory inputs to the stimulated neuron, and (ii) that weakened synaptic transmission fidelity over time with higher stimulation frequencies was mediated by short-term synaptic plasticity. As such, the biophysical modelling approach takes into consideration both anatomical (local microcircuitry) and physiological (short-term synaptic dynamics) properties. At stimulation frequencies below the threshold for synaptic depression (i.e. <20–30 Hz) [
], our model showed that neuronal responses were the result of a temporal summation of stimulus-evoked responses. In structures with predominantly excitatory inputs, this led to increases in neuronal output, whereas the opposite occurred in structures with predominantly inhibitory inputs. Beyond the threshold for synaptic depression, the strengths of successive stimulus-evoked responses were progressively reduced (i.e. a loss of synaptic transmission fidelity). In the Vim and Rt, with high frequencies, we observed an initial excitatory response which weakened over time. In the SNr, stimulus-evoked inhibitory responses were of sufficient magnitude to induce a substantial amount of neuronal inhibition; however, the SNr is also affected by synaptic depression, as evidenced by our previous work showing progressive, frequency-dependent decreases to the amplitudes of extracellular evoked field potentials in SNr with stimulation frequencies ≥20 Hz [
]. One may then assume that since synaptic depression would weaken the strength of inhibitory synaptic transmission, neuronal firing should increase via disinhibition. However, our model shows non-selective synaptic depression of both inhibitory and excitatory synaptic currents, which is supported by experimental work in rodent STN slices which demonstrated that pharmacologically-isolated excitatory and inhibitory postsynaptic potentials were both depressed during HFS [
]. This would also explain the suppression of somatic firing in STN with higher stimulation frequencies, whereas the stimulus-evoked responses with lower frequencies produced rather weak net inhibitory responses due to the more homogenous distribution of excitatory and inhibitory inputs to STN.
A recent theoretical study incorporated the TM and LIF models to characterize excitatory post-synaptic currents (EPSCs) and action potential signaling of depressive, facilitatory, and pseudolinear synapses being directly activated by DBS [
]. Herein, we have built upon this work by modelling both excitatory and inhibitory postsynaptic currents to generate a site-specific (i.e. dependent upon the proportion of convergent inhibitory/excitatory inputs) and frequency-dependent DBS-mediated net current elicited by each pulse. Thus, our model can capture stimulation-mediated neuronal dynamics across various brain targets and applied stimulation frequencies. Notably, each of our studies suggest that short-term synaptic depression may be a putative mechanism of high-frequency DBS. In line with these findings, previous computational work has suggested that high-frequency DBS may lead to axonal and synaptic failures which suppress the synaptic transfer of firing rate oscillations, synchrony, and rate-coded information from the STN to its synaptic targets [
The primate subthalamic nucleus. III. Changes in motor behavior and neuronal activity in the internal pallidum induced by subthalamic inactivation in the MPTP model of parkinsonism.
], then the overall “neuronal output” should be considered as the summation of the somatic firing (that is influenced by afferent axon/axon terminal activations [
]) plus the direct efferent axonal activations; we have incorporated this summation within Fig. 1B. Thus, HFS applications which completely suppresses somatic firing would replace neuronal output with a more regular pattern of output corresponding to the stimulation frequency; in line with the theory of “decoupling of the axon and soma” [
]. However, we suggest that in cases where somatic firing is not completely suppressed, such as in the STN at lower stimulation frequencies or when stimulating structures such as the Vim and Rt, the effect is a “summation” of axonal and somatic firing, rather than an explicit decoupling.
Translational implications
The selectively bimodal and frequency-dependent somatic responses described here should be taken into consideration in the development of novel stimulation paradigms and DBS indications. In applications of DBS which utilize a high stimulation frequency, suppression of somatic output is likely achieved (though as mentioned above, axonal activations should be considered). Stimulation paradigms which utilize low stimulation frequencies and are applied to areas of the brain with predominantly glutamatergic inputs may depend upon periodic facilitation of somatic firing, with one possible example being low-frequency pedunculopontine-DBS [
]. Low-frequency stimulation in an area of the brain with predominantly inhibitory inputs may on the other hand cause periodic inhibitions. In either case, low-frequency stimulation can induce oscillatory neuronal behaviour (as seen in Supplementary Fig. 2). Knowledge of the site-specific and frequency-dependent properties of DBS can inform the development of novel stimulation paradigms such as closed-loop stimulation for on demand upregulation or downregulation of neuronal firing, or for induction or disruption of neuronal oscillations. Indeed, stimulation parameters are often decided upon empirically. Based on the findings presented here, knowledge of the local microcircuitry (distribution of afferent inputs) inherent to the stimulated brain region (i.e. therapeutic targets of interest for DBS application) may allow us to infer/predict the stimulation frequency response properties. As such, our comprehensive computational model may represent a valuable tool for physiologically-informed stimulation programming and paradigm development in prospective DBS targets and indications, particularly as our model was developed based on in vivo experimental data from the human brain. Further clinical and physiological implications of basal ganglia, Vim, and Rt stimulation are discussed in greater detail within the Supplementary Material.
Considerations and limitations
Although we did not record from any structures downstream of the stimulation site, it is perhaps possible to infer the downstream effects based on the results presented here. For example, in STN stimulation, activation of the glutamatergic subthalamo-pallidal/nigral efferents may cause excitatory responses downstream [
Biochemical and electrophysiological changes of substantia nigra pars reticulata driven by subthalamic stimulation in patients with Parkinson's disease.
], especially if lower stimulation frequencies are used; whereas with higher stimulation frequencies the downstream glutamatergic drive may in fact be weakened [
Electrophysiological and metabolic evidence that high-frequency stimulation of the subthalamic nucleus bridles neuronal activity in the subthalamic nucleus and the substantia nigra reticulata.
] due to synaptic depression. Further studies are warranted in order to better understand the possible orthodromic (and antidromic) network phenomena of DBS [
]. Moreover, studies relating to the downstream and upstream DBS effects would allow us to better understand the mechanisms of DBS applied to white matter tracts (such as forniceal-DBS). Another notable limitation of this study is that the applied stimulation trains were limited to short durations compared to that of hours, days, or longer in clinical applications. Stimulation effects over longer durations are yet to be validated and are to be considered in future work. Moreover, while this study aimed to elucidate differential mechanisms involved in stimulation of various brain structures, behavioural and clinical correlates were not assessed here directly. However, the high-frequency microstimulation applied to Vim was shown to be effective at suppressing tremor [
] confirming that the stimulation parameters used were clinically relevant. Moreover, the stimulation parameters used here (in particular, the 100Hz microstimulation trains) are comparable in terms of total electrical energy delivered during clinically-applied DBS macrostimulation [
] though are of greater current density due to the smaller stimulating surface. Another important limitation of this work is that the explanation of site-specific mechanistic disparities based on the proportions of inhibitory/excitatory afferent inputs does not account for the possible contributions of glia [
] which should be considered in future work. It is also important to note that the experimental data within this study was acquired in context of pathological circuits and may not reflect the typical responses to stimulation in a healthy individual. Moreover, Vim data was acquired from two patients with Parkinson's disease and nine with essential tremor; however, our experimental and computational analyses did not account for possible differences in baseline firing characteristics nor responses to stimulation across disease conditions. Finally, the Vim computational frequency responses require further tuning to more precisely capture stronger initial excitatory responses during HFS; a topic for future work.
Conclusion
The presented results demonstrate the site-specific and frequency-dependent neuronal effects of extracellular stimulation. Neuronal suppression could be achieved either by stimulus-evoked inhibitory events in structures with predominantly GABAergic inputs (STN and SNr) or non-selectively when sustained HFS was delivered. Stimulus-evoked neuronal excitatory responses were exclusive to structures with predominantly glutamatergic inputs (Vim and Rt), particularly with lower stimulation frequencies. Our computational model showed that the bimodal site-specific stimulus-evoked responses could be explained by differences in the distributions of inhibitory and excitatory inputs to the stimulated target structures, whereas convergence towards neuronal suppression with sustained HFS could be explained by synaptic depression.
CRediT authorship contribution statement
Luka Milosevic: Conceptualization, Methodology, Formal analysis, Investigation, Data curation, Writing – original draft. Suneil K. Kalia: Resources, Project administration, Writing – review & editing. Mojgan Hodaie: Resources, Project administration, Writing – review & editing. Andres M. Lozano: Resources, Project administration, Writing – review & editing. Milos R. Popovic: Resources, Project administration, Writing – review & editing. William D. Hutchison: Conceptualization, Investigation, Resources, Project administration, Writing – review & editing. Milad Lankarany: Methodology, Software, Formal analysis, Writing – original draft.
Declaration of competing interest
S.K.K., M.H., W.D.H. have received honoraria, travel funds, and/or grant support from Medtronic (not related to this work). A.M.L. has received honoraria, travel funds, and/or grant support from Medtronic, Boston Scientific, St. Jude-Abbott, and Insightec (not related to this work). M.R.P. is a shareholder in MyndTec Inc. A.M.L. is a co-founder of Functional Neuromodulation Ltd. L.M. and M.L. have no financial disclosures.
Acknowledgements
The authors thank Tameem Al-Ozzi for assistance in data collection and the patients for their participation.
Appendix A. Supplementary data
The following is the Supplementary data to this article:
This work was supported in part by the Natural Sciences and Engineering Research Council: Discovery Grant RGPIN-2016-06358 (M.R.P.), Dean Connor and Maris Uffelmann Donation (M.R.P.), Walter & Maria Schroeder Donation (M.R.P.), and the Dystonia Medical Research Foundation (W.D.H.).
References
Limousin P.
Krack P.
Pollak P.
Benazzouz A.
Ardouin C.
Hoffmann D.
et al.
Electrical stimulation of the subthalamic nucleus in advanced Parkinson's disease.
Cellular, molecular, and clinical mechanisms of action of deep brain stimulation—a systematic review on established indications and outlook on future developments.
Complex locking rather than complete cessation of neuronal activity in the globus pallidus of a 1-methyl-4-phenyl-1,2,3,6-tetrahydropyridine-treated primate in response to pallidal microstimulation.
Distribution and properties of GABAB antagonist [3H]CGP 62349 binding in the rhesus monkey thalamus and basal ganglia and the influence of lesions in the reticular thalamic nucleus.
Organization of projections from the anterior pole of the nucleus reticularis thalami (NRT) to subdivisions of the motor thalamus: light and electron microscopic studies in the Rhesus monkey.
The primate subthalamic nucleus. III. Changes in motor behavior and neuronal activity in the internal pallidum induced by subthalamic inactivation in the MPTP model of parkinsonism.
Biochemical and electrophysiological changes of substantia nigra pars reticulata driven by subthalamic stimulation in patients with Parkinson's disease.
Electrophysiological and metabolic evidence that high-frequency stimulation of the subthalamic nucleus bridles neuronal activity in the subthalamic nucleus and the substantia nigra reticulata.
00092-9&id=gr1.jpg){kind=link}
00092-9&id=gr2.jpg){kind=link}
00092-9&id=gr3.jpg){kind=link}
00092-9&id=gr4.jpg){kind=link}
00092-9&id=gr5.jpg){kind=link}
00092-9&id=gr6.jpg){kind=link}
00092-9&id=gr7.jpg){kind=link}
00092-9&id=gr8.jpg){kind=link}